text:
mutable struct Iterator{K,V}
    # We used to have path::Vector{Tuple{INode{K,V},Int}} here, but to reduce allocations in
    # older versions of julia it was changed to two vectors, one of inodes, and one of indexes
    # in the inodes. In the future we can revert this change.
    path_inode::Vector{INode{K,V}}
    path_i::Vector{Int}
    leaf::Leaf{K,V}
    leaf_index::Int
    pathmask::UInt64
    pathhash::UInt64
end

function Base.iterate(dict::MuttsDict{K,V}) where {K,V}
    if dict.n == 0
        return nothing
    end

    iter = _make_iterator(dict)
    return _iterate_next(iter)
end

function Base.iterate(::MuttsDict{K,V}, iter::Iterator{K,V}) where {K,V}
    return _iterate_next(iter)
end

# Avoid allocating when iterating tiny dictionaries that contain just a single leaf.
@generated function _empty_vector(::Type{T}) where {T}
    vec = T[]
    return vec
end

# Returns an Iterator{K,V} positioned before the first record. Must call _next(iter)
# before attempting to get the key/value pair.
function _make_iterator(dict::MuttsDict{K,V}) where {K,V}
    root = dict.root
    pathmask = UInt64(0)
    pathhash = UInt64(0)
    if isa(root, Leaf{K,V})
        leaf = root::Leaf{K,V}
        # Note: it is safe to pass the shared `_empty_path(K,V)` instance here, because
        # we know it will not be mutated in _next(), since root isa Leaf.
        path_inode = _empty_vector(INode{K,V})
        path_i = _empty_vector(Int)
        return Iterator{K,V}(path_inode, path_i, leaf, 0, pathmask, pathhash)
    end

    # Descend to leftmost leaf
    path_inode = INode{K,V}[]
    path_i = Int[]
    node = root::INode{K,V}
    while isa(node, INode{K,V})
        push!(path_inode, node)
        push!(path_i, 1)
        pathmask = (pathmask << UInt64(8)) | (length(node.entries)-1)
        pathhash = (pathhash << UInt64(8))
        if isa(node.entries, Vector{Leaf{K,V}})
            # Specialize this path to avoid dynamic typing
            node = @inbounds node.entries[1]
            break
        end
        node = @inbounds node.entries[1]
    end
    pathmask = pathmask << UInt64(16)
    pathhash = pathhash << UInt64(16)
    return Iterator{K,V}(path_inode, path_i, node::Leaf{K,V}, 0, pathmask, pathhash)
end

function _iterate_next(iter::Iterator{K,V}) where {K,V}
    if _next(iter)
        i = iter.leaf_index
        (ki,vi) = @inbounds iter.leaf.entries[i]
        return ((ki,vi),iter)
    else
        return nothing
    end
end

function _next(iter::Iterator{K,V}) where {K,V}
    height = length(iter.path_inode)
    while true
        # Scan for the next (k,v) pair in the leaf
        num_slots = length(iter.leaf.entries)
        leaf = iter.leaf
        for i in iter.leaf_index+1:num_slots
            if _slot_occupied(leaf, reinterpret(UInt64,i), num_slots)
                (ki,vi) = @inbounds leaf.entries[i]
                h = hash(ki)
                if (h & iter.pathmask) == iter.pathhash
                    iter.leaf_index = i
                    return true
                end
            end
        end

        # Lowest 16 bits are for leaf
        iter.pathmask >>= UInt64(16)
        iter.pathhash >>= UInt64(16)

        # Ascend until we are in a child with a right sibling
        while length(iter.path_inode) > 0
            node = pop!(iter.path_inode)
            i = pop!(iter.path_i)
            iter.pathmask >>= UInt64(8)
            iter.pathhash >>= UInt64(8)
            if node.entries isa Vector{INode{K,V}}
                # This branch is specialized for Vector{INode{K,V}}
                if i < length(node.entries)
                    push!(iter.path_inode, node)
                    push!(iter.path_i, i+1)
                    iter.pathmask = (iter.pathmask << UInt64(8)) | (length(node.entries)-1)
                    iter.pathhash = (iter.pathhash << UInt64(8)) | i
                    break
                end
            else
                # This branch is specialized for Vector{Leaf{K,V}}
                if i < length(node.entries)
                    push!(iter.path_inode, node)
                    push!(iter.path_i, i+1)
                    iter.pathmask = (iter.pathmask << UInt64(8)) | (length(node.entries)-1)
                    iter.pathhash = (iter.pathhash << UInt64(8)) | i
                    break
                end
            end
        end

        # Have reached root? If so, done.
        if length(iter.path_inode) == 0
            iter.leaf_index = 0
            return false
        end

        # Descend through leftmost children
        node = iter.path_inode[end]
        i = iter.path_i[end]
        while true
            if node.entries isa Vector{Leaf{K,V}}
                (node,i) = (node.entries[i],1)
                break
            end
            (node,i) = (node.entries[i],1)
            push!(iter.path_inode, node)
            push!(iter.path_i, i)

            if node.entries isa Vector{Leaf{K,V}}
                # This branch is specialized for Vector{Leaf{K,V}}
                entries_mask = reinterpret(UInt64,length(node.entries)-1)
                iter.pathmask = (iter.pathmask << UInt64(8)) | entries_mask
                iter.pathhash = (iter.pathhash << UInt64(8))
            else
                # This branch is specialized for Vector{INode{K,V}}
                entries_mask = reinterpret(UInt64,length(node.entries)-1)
                iter.pathmask = (iter.pathmask << UInt64(8)) | entries_mask
                iter.pathhash = (iter.pathhash << UInt64(8))
            end
        end

        iter.pathmask <<= UInt64(16)
        iter.pathhash <<= UInt64(16)
        iter.leaf = node::Leaf{K,V}
        iter.leaf_index = 0
    end
end

function _get_key(iter::Iterator{K,V})::K where {K,V}
    @assert iter.leaf_index != 0
    (k,_) = @inbounds iter.leaf.entries[iter.leaf_index]
    return k
end

function _get_value(iter::Iterator{K,V})::V where {K,V}
    @assert iter.leaf_index != 0
    (_,v) = @inbounds iter.leaf.entries[iter.leaf_index]
    return v
end

# More efficient version that avoids allocations from the
# return value of Base.iterate().
function Base.foreach(f, dict::MuttsDict{K,V}) where {K,V}
    iter = _make_iterator(dict)
    while _next(iter)
        k = _get_key(iter)
        v = _get_value(iter)
        f(k, v)
    end
end

"""
    keys(dict)

Returns an iterator over the keys in the dict.
"""
function Base.keys(dict::MuttsDict)
    return (k for (k, _) in dict)
end


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

name = "LCIO_Julia_Wrapper"
version = v"0.13.2"
julia_versions = [v"1.6", v"1.7", v"1.8"]

# Collection of sources required to build LCIOWrapBuilder
sources = [
	GitSource("https://github.com/jstrube/LCIO_Julia_Wrapper.git", "e28132bfdc0664faf9724a74b2ae33803c26dc5a")
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/LCIO_Julia_Wrapper/
mkdir build && cd build
cmake -DJulia_PREFIX=${prefix} -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_FIND_ROOT_PATH=${libdir}/cmake -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ..
make -j${nproc}
make install
install_license $WORKSPACE/srcdir/LCIO_Julia_Wrapper/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")

platforms = expand_cxxstring_abis(vcat(libjulia_platforms.(julia_versions)...))
filter!(!Sys.isfreebsd, platforms)
filter!(!Sys.iswindows, platforms)
filter!(p -> arch(p) ∉ ("armv6l", "armv7l") , platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("liblciowrap", :lciowrap)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="libcxxwrap_julia_jll")),
    Dependency(PackageSpec(name="LCIO_jll")),
    BuildDependency(PackageSpec(name="libjulia_jll"))
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"8", julia_compat="1.6")

ext: jl

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function ellpic_bulirsch(n::T,k2::T) where {T <: Real}
# Computes the complete elliptical integral of the third kind using
# the algorithm of Bulirsch (1965):
kc = sqrt(1.0-k2)
ca = sqrt(eps(n))
if (kc*(n+1.0)) == 0
  ellpic = NaN
else
  ee = kc
  m0 = 1.0
  if n > -1.0
    c = 1.0
    p = sqrt(n+1.0)
    d = 1.0/p
  else
    g = -n
    f = -k2-n
    p = sqrt(f/g)
    d = -k2/(g*p)
    c = 0.0
  end

  f = copy(c)
  c += d/p
  g = ee/p
  d = (f*g+d)*2.0
  p += g
  g = m0
  m0 += kc
  while abs(g-kc) > (ca*g)
    kc = 2.0*sqrt(ee)
    ee = kc*m0
    f = copy(c)
    c += d/p
    g = ee/p
    d = (f*g+d)*2.0
    p += g
    g = m0
    m0 += kc
  end
  ellpic= 0.5*pi*(c*m0+d)/(m0*(m0+p))
end
return ellpic
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
using ReservoirComputing

train_data = ones(Int, 2, 10)
res_size = 10
generations = 5
permutations = 2

reca = RECA_TwoDim(train_data, res_size, generations, permutations)

@test isequal(reca.train_data, train_data)
@test isequal(reca.res_size, res_size)
@test isequal(reca.generations, generations)
@test isequal(reca.permutations, permutations)


text:
using MyExample
using Test

#import MyExample
#@test my_f(2,1) == 7

@test_throws TestSetException @testset "MyExample.jl" begin
    # 2x + 3y
    @test my_f(2,1) == 7
    @test my_f(2,3) == 13
    @test my_f(2,3) == 14
end

@testset "MyExample.jl" begin
    # 2x + 3y
    @test my_f(2,1) == 7
    @test my_f(2,3) == 13
    # @test my_f(2,3) == 14
end


text:
module KevinsPC

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Kevin's PC" => Ahorn.EntityPlacement(
        Maple.KevinsPC
    )
)

spritePC = "objects/kevinspc/pc"
spriteSpectogram = "objects/kevinspc/spectogram"

# We don't need the back texture, it fits inside
function Ahorn.selection(entity::Maple.KevinsPC)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(spritePC, x, y, jx=0.5, jy=1.0)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.KevinsPC, room::Maple.Room)
    Ahorn.drawSprite(ctx, spritePC, 0, 0, jx=0.5, jy=1.0)
    Ahorn.drawImage(ctx, spriteSpectogram, -16, -39, 0, 0, 32, 18)
end

end


text:
import Dolo: Model

import Dolo
import Dolo: time_iteration
import Dolo: EmptyDomain

model = Model("examples/models/rbc.yaml")

m,s,x,p = model.calibration[:exogenous, :states, :controls, :parameters]

Dolo.transition(model, m,s,x,m,p)
Dolo.perturb(model)

import Dolo: Domain

Dolo.time_iteration(model)

function get_infos(model::AModel)
    if "infos" in keys(model.data)
        return model.data["infos"]
    else
        return Dict()
    end
end

# ####################################################
# ###### TODO   conversion to SymExpr incorrect ######
# ####################################################

# model = Model("examples/models/rbc.yaml")

import Dolang
# ff_a = get_factory(model, "transition")

ff, ff_l, ff_u = Dolo.get_factory(model, "arbitrage")

source_l = Dolang.gen_generated_gufun(ff_l)

l = Base.eval(Dolang, source_l)

# m,s,x,p = model.calibration[:exogenous, :states, :controls, :parameters]

# source_g = Dolang.gen_generated_gufun(ff_a)
# source_f = Dolang.gen_generated_gufun(ff_s)

# g = Base.eval(Dolang, source_g)
# f = Base.eval(Dolang, source_f)

# g(m,s,x,m,p)
# f(m,s,x,m,s,x,p)


text:
module FrostHelperFlutterBird

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/CustomFlutterBird" FHCustomFlutterbird(x::Integer, y::Integer, colors::String="89fbff,f0fc6c,f493ff,93baff")

const placements = Ahorn.PlacementDict(
    "Custom Flutterbird (Frost Helper)" => Ahorn.EntityPlacement(
        FHCustomFlutterbird
    ),
)

const sprite = "scenery/flutterbird/idle00.png"

const colors = [
    (137, 251, 255, 255) ./ 255,
    (240, 252, 108, 255) ./ 255,
    (244, 147, 255, 255) ./ 255,
    (147, 186, 255, 255) ./ 255
]

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

    return Ahorn.getSpriteRectangle(sprite, x, y, jx=0.5, jy=1.0)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FHCustomFlutterbird, room::Maple.Room)
    rng = Ahorn.getSimpleEntityRng(entity)
    color = rand(rng, colors)

    Ahorn.drawSprite(ctx, sprite, 0, 0, jx=0.5, jy=1.0, tint=color)
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
module MMDQueues

# TODO: I don't think the extra slot at the end of the queue is needed anymore.

#=
This module defines a custom priority queue optimised for implementing
the mmd heuristic.

Retrieving the highest priority element and decrementing the priority
of elements by 1 are both order 1 operations in complextity. mmd does
not need an insertion operation so it can be omitted.

The queue is represented by a sorted array of (vertex, degree) pairs.
An array of keys or indices for the vertices is maintained for order
1 look up of a vertex in the queue.

A bag structure is also maintained to help keep the queue sorted
when the degree of a vertex is decremented. ie. the queue is
partitioned into bags, with all vertices in bag d having degree d.
This is maintained by an array containing the positions of the left
most boundary of a bag. The vertices of degree d begin at bags[d+1]
inclusively and end at bags[d+2] exclusively.

            Bag 0 <----- empty Bags ------>  Bag 3
            |___|                            |___|
              |--Bag 1-|-------Bag 2-----------|--Bag 4-|--Bag 5-|
              |        |                       |        |        |
              |        |                       |        |        |
queue: [ ][ ][ ][ ][ ][ ][ ][ ][ ][ ][ ][ ][ ][ ][ ][ ][ ][ ][ ][ ]
        |  |  |
        |  |   \__
removed----       \--min
=#

export MMDQueue, initialise_mmdqueue!, pull!, decrement!

"""
A custom priority queue optimised for implementing the mmd heuristic.

Retrieving the highest priority element and decrementing the priority
of elements by 1 are both order 1 operations in complextity. mmd does
not need an insertion operation so it can be omitted.

See comments above docstring for implementation details.
"""
mutable struct MMDQueue
    queue::Vector{Pair{UInt16, UInt16}}
    keys::Vector{UInt16}
    bags::Vector{UInt16}
    min::UInt16
end

function MMDQueue(vertices::Vector{UInt16}, degrees::Vector{UInt16})
    N = length(vertices)

    queue = Array{Pair{UInt16, UInt16}, 1}(undef, N+1)
    keys = similar(vertices)
    bags = fill(UInt16(N+1), N)

    pq = MMDQueue(queue, keys, bags, 0x0001)
    initialise_mmdqueue!(pq, vertices, degrees)

    pq
end

"""
Initialises the given MMDQueue to hold the given vertices
in order of increasing degree.

I is assumed to hold permutation of indices which puts
the vertices into the correct order.
"""
function initialise_mmdqueue!(pq::MMDQueue,
                            vertices::AbstractArray,
                            degrees::Vector{UInt16})
    queue = pq.queue
    keys = pq.keys
    bags = pq.bags
    N = UInt16(length(vertices))

    # Initialise the queue.
    for i = 0x0001:N
        @inbounds vi = vertices[i]
        @inbounds queue[i] = vi => degrees[vi]
    end

    # Sort the queue (Use quicksort as it's non-allocating).
    q = @view queue[0x0001:N]
    sort!(q; alg=QuickSort, by=p -> p.second)

    # Initialise keys and bag structure
    bag = 0x0000
    @inbounds for i = 0x0001:N
        (v, d) = queue[i]
        keys[v] = i

        if d >= bag
            bags[bag+1:d+1] .= i
            bag = d + 0x0001
        end
    end

    pq.min = 0x0001
    nothing
end

"""
    pull!(q::MMDQueue)

Return the element with highest priority and increment all
bag boundaries to point to the next highest priority element.
"""
function pull!(q::MMDQueue)
    @inbounds p = q.queue[q.min]
    q.min += 0x0001
    @inbounds q.bags[0x0001:p.second+0x0001] .= q.min
    p
end

#=
decrement! is implemented by swapping v, in bag d, with the vertex
at the boundary of bag d. The degree of v is decremented and the
boundary of bag d is incremented so that v is no longer in bag d
and now in bag d - 1.

Example:

Initial state:
               ---to swap --
               |           |
               |           |
degree: [ ][ ][ ][ ][ ][ ][d][ ][ ]
queue : [ ][ ][p][ ][ ][ ][v][ ][ ]
bags  : [ ][b][ ][ ][ ][ ][ ][ ][ ]
keys  : [ ][ ][ ][ ][i][ ][ ][ ][ ]
ind   :     d  b     v     i

Final state: (D = d - 1, B = b + 1)

                  --- new boundary of bag d
                  |
                  |
degree: [ ][ ][D][|][ ][ ][ ][ ][ ]
queue : [ ][ ][v][|][ ][ ][p][ ][ ]
bags  : [ ][B][ ][|][ ][ ][ ][ ][ ]
keys  : [ ][ ][ ][|][b][ ][ ][ ][ ]
ind   :     d  b  |  v     i
=#

"""
    decrement!(q::MMDQueue, v::UInt16)

Decrease the degree of v by 1 and update the queue to maintain
the correct priority for vertices and bag structure.

See comments above docstring for implementation details.
"""
function decrement!(q::MMDQueue, v::UInt16)
    @inbounds i = q.keys[v]           # index of v.
    @inbounds if i >= q.min
        d = q.queue[i].second   # the degree of v.
        b = q.bags[d + 0x0001]  # index of bag d boundary.

        # Move the pair at the boundary to where v was.
        p = q.queue[b]
        q.queue[i] = p
        q.keys[p.first] = i

        # Place new pair for v at the boundary and
        # increment the boundary position.
        q.queue[b] = v => d - 0x0001
        q.keys[v] = b
        q.bags[d + 0x0001] += 0x0001
    end

    nothing
end

###
### Base.show implementation for MMDQueues.
###

function Base.show(io::IO, q::MMDQueue)
    if length(q.queue) <= 32
        show_verbose(io, q)
    else
        show_compact(io, q)
    end
end

"""
Print some of the fields of the given MMDQueue.
"""
function show_compact(io::IO, q::MMDQueue)
    ks = [p.first for p in q.queue]
    ds = [p.second for p in q.queue]
    println(io, "degree: ", Int.(ds))
    println(io, "vertex: ", Int.(ks))
    println(io, "bag   : ", Int.(q.bags))
    println(io, "key   : ", Int.(q.keys))
end

"""
Draw a verbose image of the given MMDQueue.
"""
function show_verbose(io::IO, q::MMDQueue)
    println()
    bar = "  |  "
    space = "     "

    # Draw the bag structure of the queue.
    display_array = -1*ones(Int, length(q.queue))
    for (i, b) in enumerate(q.bags)
        display_array[b] = i - 1
    end

    Bl = ""
    Bu = ""
    for i in display_array
        if i == -1
            Bl *= space
            Bu *= space
        else
            Bl *= bar
            Bu *= "  $i  "
        end
    end

    println(io, Bu)
    println(io, Bl)

    # Draw the queue.
    Q = ""
    D = ""
    for i = 1:length(q.queue)
        v, d = q.queue[i]
        Q *= "[" * lpad(v, 3) * "]"
        D *= "[" * lpad(d, 3) * "]"
    end

    println(io, Q)
    println(io, D)

    # Draw the arrow pointing to the next item of the queue.
    m = q.min
    min_arrow = space^(m - 1) * "  ^  "
    min_bar   = space^(m - 1) * bar
    min_base  = space^(m - 1) * " min "

    println(io, min_arrow)
    println(io, min_bar)
    println(io, min_base)
end

end


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

"""
    CatmullClark()

Catmull-Clark refinement of polygonal meshes.

Strictly speaking, the Catmull-Clark algorithm
is used for subdivision surface modeling, not
just mesh refinement. At each step of refinement,
the vertices are adjusted to approximate a smooth
surface.

## References

* Catmull & Clark. 1978. [Recursively generated
  B-spline surfaces on arbitrary topological meshes]
  (https://www.sciencedirect.com/science/article/abs/pii/0010448578901100)
"""
struct CatmullClark <: RefinementMethod end

function refine(mesh, ::CatmullClark)
  # retrieve geometry and topology
  points = vertices(mesh)
  connec = topology(mesh)

  # convert to half-edge structure
  t = convert(HalfEdgeTopology, connec)

  # add centroids of elements
  ∂₂₀ = Boundary{2,0}(t)
  epts = map(1:nelements(t)) do elem
    ps = view(points, ∂₂₀(elem))
    cₒ = sum(coordinates, ps) / length(ps)
    Point(cₒ)
  end

  # add midpoints of edges
  ∂₁₂ = Coboundary{1,2}(t)
  ∂₁₀ = Boundary{1,0}(t)
  fpts = map(1:nfacets(t)) do edge
    ps = view(epts, ∂₁₂(edge))
    qs = view(points, ∂₁₀(edge))
    ∑p = sum(coordinates, ps)
    ∑q = sum(coordinates, qs)
    M  = length(ps) + length(qs)
    Point((∑p + ∑q) / M)
  end

  # move original vertices
  ∂₀₂ = Coboundary{0,2}(t)
  ∂₀₀ = Adjacency{0}(t)
  vpts = map(1:nvertices(t)) do u
    # original point
    P = coordinates(points[u])

    # average of centroids
    ps = view(epts, ∂₀₂(u))
    F  = sum(coordinates, ps) / length(ps)

    # average of midpoints
    vs = ∂₀₀(u)
    n  = length(vs)
    R  = sum(vs) do v
      uv = view(points, [u,v])
      sum(coordinates, uv) / 2
    end / n

    Point((F + 2R + (n - 3)P) / n)
  end

  # new points in refined mesh
  newpoints = [vpts; epts; fpts]

  offset₁ = length(vpts)
  offset₂ = offset₁ + length(epts)

  # connect vertices into new quadrangles
  ∂₂₁ = Boundary{2,1}(t)
  newconnec = Connectivity{Quadrangle,4}[]
  for elem in 1:nelements(t)
    verts = CircularVector(∂₂₀(elem))
    edges = CircularVector(∂₂₁(elem))
    for i in 1:length(edges)
      u = elem       + offset₁
      v = edges[i]   + offset₂
      w = verts[i+1]
      z = edges[i+1] + offset₂
      quad = connect((u, v, w, z))
      push!(newconnec, quad)
    end
  end

  SimpleMesh(newpoints, newconnec)
end


text:
function renderloop(screen::Screen; framerate = 1/30, prerender = () -> nothing)
    try
        while isopen(screen)
            t = time()
            GLFW.PollEvents() # GLFW poll
            prerender()
            make_context_current(screen)
            render_frame(screen)
            GLFW.SwapBuffers(to_native(screen))
            diff = framerate - (time() - t)
            if diff > 0
                sleep(diff)
            else # if we don't sleep, we need to yield explicitely
                yield()
            end
        end
    catch e
        destroy!(screen)
        rethrow(e)
    end
    destroy!(screen)
    return
end

function setup!(screen)
    glEnable(GL_SCISSOR_TEST)
    if isopen(screen)
        glScissor(0, 0, widths(screen)...)
        glClearColor(1, 1, 1, 1)
        glClear(GL_COLOR_BUFFER_BIT)
        for (id, rect, clear, visible, color) in screen.screens
            if visible[]
                a = rect[]
                rt = (minimum(a)..., widths(a)...)
                glViewport(rt...)
                if clear[]
                    c = color[]
                    glScissor(rt...)
                    glClearColor(red(c), green(c), blue(c), alpha(c))
                    glClear(GL_COLOR_BUFFER_BIT)
                end
            end
        end
    end
    glDisable(GL_SCISSOR_TEST)
    return
end

const selection_queries = Function[]

"""
Renders a single frame of a `window`
"""
function render_frame(screen::Screen)
    nw = to_native(screen)
    GLAbstraction.is_context_active(nw) || return
    fb = screen.framebuffer
    wh = Int.(framebuffer_size(nw))
    resize!(fb, wh)
    w, h = wh
    glDisable(GL_STENCIL_TEST)
    #prepare for geometry in need of anti aliasing
    glBindFramebuffer(GL_FRAMEBUFFER, fb.id[1]) # color framebuffer
    glDrawBuffers(2, [GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1])
    glClearColor(0,0,0,0)
    glClear(GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT | GL_COLOR_BUFFER_BIT)
    setup!(screen)
    glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE)
    GLAbstraction.render(screen, true)
    # transfer color to luma buffer and apply fxaa
    glBindFramebuffer(GL_FRAMEBUFFER, fb.id[2]) # luma framebuffer
    glDrawBuffer(GL_COLOR_ATTACHMENT0)
    glViewport(0, 0, w, h)
    glClearColor(0,0,0,0)
    glClear(GL_COLOR_BUFFER_BIT)
    GLAbstraction.render(fb.postprocess[1]) # add luma and preprocess

    glBindFramebuffer(GL_FRAMEBUFFER, fb.id[1]) # transfer to non fxaa framebuffer
    glViewport(0, 0, w, h)
    glDrawBuffer(GL_COLOR_ATTACHMENT0)
    GLAbstraction.render(fb.postprocess[2]) # copy with fxaa postprocess

    #prepare for non anti aliased pass
    glDrawBuffers(2, [GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1])

    GLAbstraction.render(screen, false)
    #Read all the selection queries
    glReadBuffer(GL_COLOR_ATTACHMENT1)
    for query_func in selection_queries
        query_func(fb.objectid, w, h)
    end
    glBindFramebuffer(GL_FRAMEBUFFER, 0) # transfer back to window
    glViewport(0, 0, w, h)
    glClearColor(0, 0, 0, 0)
    glClear(GL_COLOR_BUFFER_BIT)
    GLAbstraction.render(fb.postprocess[3]) # copy postprocess
    return
end

function id2rect(screen, id1)
    # TODO maybe we should use a different data structure
    for (id2, rect, clear, color) in screen.screens
        id1 == id2 && return true, rect
    end
    false, IRect(0,0,0,0)
end

function GLAbstraction.render(screen::Screen, fxaa::Bool)
    for (zindex, screenid, elem) in screen.renderlist
        found, rect = id2rect(screen, screenid)
        found || continue
        a = rect[]
        glViewport(minimum(a)..., widths(a)...)
        if fxaa && elem[:fxaa][]
            render(elem)
        end
        if !fxaa && !elem[:fxaa][]
            render(elem)
        end
    end
    return
end


text:
# # Fractional Kelvin-Voigt
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/model_kv.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/model_kv.ipynb)

using RHEOS
## include a helper function for plotting
include("assets/plothelper.jl");
#-

Fract_KelvinVoigt

# #### Constitutive Equation
# ```math
# \sigma(t) = c_\alpha \frac{d^{\alpha} \epsilon(t)}{dt^{\alpha}}+ c_{\beta} \frac{d^\beta \epsilon(t)}{dt^\beta}
# ```

# ```math
# \text{for}\; \ 0 \leq \beta \leq \alpha \leq 1
# ```

# #### Relaxation Modulus
# ```math
# G(t) =\frac{c_{\alpha} }{\Gamma(1-\alpha)} t^{-\alpha}+\frac{c_{\beta} }{\Gamma(1-\beta)} t^{-\beta}
# ```

# #### Creep Modulus
# ```math
# J(t) = \frac{t^{\alpha}}{c_\alpha} E_{\alpha-\beta,1+\alpha}\left(-\frac{c_\beta}{c_\alpha} t^{\alpha-\beta}\right)
# ```

# #### Storage Modulus
# ```math
# G^{\prime}(\omega) = c_\alpha \omega^\alpha \cos\left(\alpha \frac{\pi}{2}\right) + c_\beta \omega^\beta \cos \left( \beta \frac{\pi}{2}\right)
# ```

# #### Loss Modulus
# ```math
# G^{\prime\prime}(\omega) = c_\alpha \omega^\alpha \sin\left(\alpha \frac{\pi}{2}\right) + c_\beta \omega^\beta \sin \left( \beta \frac{\pi}{2}\right)
# ```

# # Fractional (Spring) Kelvin-Voigt

FractS_KelvinVoigt
#-

models = Vector{RheoModel}()

## plot moduli for varying α
for alpha in [0.1, 0.3, 0.5, 0.7, 0.9]

    push!(models, RheoModel(FractS_KelvinVoigt, (cₐ = 1.0, a = alpha, k = 1.0)))

end

plotmodel(models, ymaxG = 5.0)

# # Fractional (Dashpot) Kelvin-Voigt

FractD_KelvinVoigt
#-

models = Vector{RheoModel}()

## plot moduli for varying β
for beta in [0.1, 0.3, 0.5, 0.7, 0.9]

    push!(models, RheoModel(FractD_KelvinVoigt, (η = 10, cᵦ= 1.0, β = beta)))

end

plotmodel(models, ymaxG = 3.0)

# # Kelvin-Voigt

KelvinVoigt
#-

models = Vector{RheoModel}()

## plot moduli for varying k
for η in [0.1, 1.0, 5.0]

    push!(models, RheoModel(KelvinVoigt, (η = η, k = 1)))

end

plotmodel(models)


text:
# -*- coding: utf-8 -*-
# dx9adl_demo_0002

VERSION >= v"0.4.0-dev+6521" && __precompile__()
module Demo_0002

include("../src/lib_fnc_defs.jl") # macros/functions with import Relocator: _mf

import WCharUTF8
import DirectX: D3DMatrix, D3DMatrixBits, array_to, Q_D3DMatrix, QQMatrix
import DirectX: VERTEX_GLYPH, FT_PosBits, FT_VectorBits, GLYPH_VTX, GLYPH_TBL
import DirectX: D3DVector, D3DVectorBits, D9F_Vecs, D9Foundation
import DirectX: RenderD3DItemsState, Dx9adl
import DirectX: PSO_SPRITE, PSO_STRING

const TXSRC, TXDST, VtxGlp = 0:2
const txs = ["_col_4.png", "_D_00.png", "_D_01.dds", "_D_02.png", "_D_go.png"]
const face_default = ("mikaP.ttf",)

m_transform = D3DMatrix()
m_rotation = D3DMatrix(
  1.0,         0.,         0.,  0.,
   0.,  1./1.4142,  1./1.4142,  0.,
   0., -1./1.4142,  1./1.4142,  0.,
   0.,         0.,         0., 1.0)
m_scale = D3DMatrix(
  1.5,  0.,  0.,  0.,
   0., 1.5,  0.,  0.,
   0.,  0., 1.5,  0.,
   0.,  0.,  0., 1.0)
m_translate = D3DMatrix(
  1.0,  0.,  0.,  0.,
   0., 1.0,  0.,  0.,
   0.,  0., 1.0,  0.,
  -4., -1., -2., 1.0)
qqm = QQMatrix(C_NULL, C_NULL, C_NULL, C_NULL) # re-set later
vg = VERTEX_GLYPH(C_NULL, C_NULL, C_NULL, 0) # re-set later
gv = GLYPH_VTX(C_NULL, 0, 0, C_NULL, 0, 0, # re-set later
  reinterpret(FT_VectorBits, Int32[0, 0])[],
  reinterpret(FT_PosBits, Int32[0])[], reinterpret(FT_PosBits, Int32[0])[],
  0, 0)
gt = GLYPH_TBL(C_NULL, C_NULL, C_NULL, C_NULL, # re-set later
  Float32(6000.), Float32(25.), 0, 0,
  0, 1024, 256, 256, C_NULL, C_NULL,
  C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL)

function initD3DItems(pIS::Ptr{RenderD3DItemsState})
  ist = @juliaobj pIS
  debugout("initD3DItems: %08X\n", ist.stat)
  d9f = @juliaobj ist.d9fnd # expect Julia structure
  d9 = @juliaobj ist.parent # expect Julia structure
  imp = replace(d9.respath * "/" * txs[1], "/", "\\") # only for Windows
  D3DXCreateTextureFromFileA(d9f.pDev, pointer(imp.data), PtrPtrU(pIS, TXSRC))
  D3DXTXB_CreateTexture(d9f.pDev, 256, 256, PtrPtrU(pIS, TXDST))
  if PtrUO(pIS, TXSRC) == C_NULL
    debugout("not found ? [%s]\n", pointer(imp.data))
    return 0::Cint
  end
  return 1::Cint
end

function cleanupD3DItems(pIS::Ptr{RenderD3DItemsState})
  ist = @juliaobj pIS
  debugout("cleanupD3DItems: %08X\n", ist.stat)
  return 1::Cint
end

function renderD3DItems(pIS::Ptr{RenderD3DItemsState})
  ist = @juliaobj pIS # expect Julia structure
  if ist.stat & 0x00000001 != 0; return 0::Cint; end
  d9f = @juliaobj ist.d9fnd # expect Julia structure
  if ist.stat & 0x00008000 != 0
    pSprite = PtrSO(pIS, PSO_SPRITE)
    D3DXTXB_RewriteTexture(PtrPtrU(pIS, TXDST), PtrPtrU(pIS, TXSRC))
    BltTexture(pIS, 0xFFFFFFFF, PtrPtrU(pIS, TXDST), 0, 0, 256, 256,
      0., 0., 0., 10., 100., 1.)
    BltTexture(pIS, 0xFFFFFFFF, PtrPtrS(pIS, PSO_STRING), 0, 0, 512, 512,
      0., 0., 0., 10., 10., .5)
    BltString(pIS, 0xFF808080, "BLTSTRING", 2, 192, 32, 0.1)
  else
    menv = @juliaobj d9f.pMenv # expect Julia structure
    D3DXMatrixRotationY(menv.world, ist.nowTime * 0.06 * pi / 180)
    SetupMatrices(pIS)
    DrawAxis(pIS)
    if ist.nowTime - ist.prevTime <= 1000
      if ist.nowTime - ist.prevTime <= 5
        debugout("renderD3DItems %02d %08X\n", ist.fps, ist.nowTime)
      end
      t = (75. - 60. * ((ist.nowTime >> 4) % 256) / 256) * pi / 180;
      gt.pIS = pIS
      qqm.transform = @ptr m_transform
      qqm.rotation = @ptr m_rotation
      qqm.scale = @ptr m_scale
      qqm.translate = @ptr m_translate
      vg.pQQM = @ptr qqm
      vg.ppTexture = PtrPtrU(pIS, TXDST)
      vg.ppVtxGlyph = PtrPtrU(pIS, VtxGlp)
      ReleaseNil(vg.ppVtxGlyph)
      vg.szGlyph = 0;
      gt.pVG = @ptr vg
      d9 = @juliaobj ist.parent
      facepath = replace(d9.respath * "/" * face_default[1], "/", "\\")
      # debugout("[%s]\n", pointer(facepath.data))
      gt.facename = pointer(facepath.data)
      nt = ist.nowTime & 0x0000FFFF
      u8 = @sprintf "%02d%c%04X%s" ist.fps 'ぁ' nt "\u3047\u3046" # UTF-8
      gt.utxt = pointer(WCharUTF8.UTF8toWCS(u8, eos=true))
      gt.ratio = Float32(6000.)
      gt.angle = Float32(25.)
      gt.reserved1 = gt.reserved0 = 0
      gt.mode = ((ist.nowTime >> 12) % 2) != 0 ? 4 : 8
      # gt.mode |= 0x40000000
      gt.td = 32 * (((ist.nowTime >> 4) % 256) + 1) # about 0.8sec
      gt.tw = 256
      gt.th = 256
      gt.rct = C_NULL
      gt.glyphBmp = C_NULL
      m_rotation.cc = m_rotation.bb = cos(t)
      m_rotation.bc = - (m_rotation.cb = sin(t))
      @q3m m_transform m_rotation m_scale m_translate
      gt.matrix = qqm.transform
      gv.col = ist.fgc
      gv.bgc = ist.bgc
      gt.vtx = @ptr gv
      gt.funcs = C_NULL
      gt.glyphContours = _mf(:d3dxglyph, :D3DXGLP_GlyphContours)
      gt.glyphAlloc = _mf(:d3dxfreetype2, :D3DXFT2_GlyphAlloc)
      gt.glyphFree = _mf(:d3dxfreetype2, :D3DXFT2_GlyphFree)
      D3DXFT2_GlyphOutline(@ptr gt)
    end
    D3DXGLP_DrawGlyph(@ptr gt)
    DrawString(pIS, ist.fgc, "DRAWSTRING", 2, 0.5, 0.5, 0.1, -3., 1., -2.)
  end
  return 1::Cint
end

function demo_0002(d9)
  debugout("adl_test &d9.istat = %08X\n", @ptr d9.istat)
  hInst = GetModuleHandleA(C_NULL)
  nShow = 1 # 1: SW_SHOWNORMAL or 5: SW_SHOW
  className = "juliaClsDx9ADLtest"
  appName = "juliaAppDx9ADLtest"
  r = InitD3DApp(
    hInst, nShow, className, appName, (@ptr d9.istat),
    cfunction(initD3DItems, Cint, (Ptr{RenderD3DItemsState},)),
    cfunction(cleanupD3DItems, Cint, (Ptr{RenderD3DItemsState},)),
    cfunction(renderD3DItems, Cint, (Ptr{RenderD3DItemsState},)))
  debugout("[%08X]\n", r)
  debugout("in\n")
  try
    r = MsgLoop(@ptr d9.istat)
  catch err
    debugout("err\n")
    println(err)
  finally
    debugout("out\n")
  end
  debugout("[%08X]\n", r)
  return r
end

end


text:
QUERY_FILES = """
  query Files {
    files {
      id
      filename
    }
  }"""

QUERY_CALIBRATION = """
  query Calibration(\$fileId: ID!) {
    calibration(fileId: \$fileId) {
      placeholder
    }
  }"""

MUTATION_ADDVARIABLE = """
  mutation addVariable (\$variable: FactorGraphInput!) {
    addVariable(variable: \$variable)
  }"""

MUTATION_ADDFACTOR = """
  mutation addFactor (\$factor: FactorGraphInput!) {
    addFactor(factor: \$factor)
  }"""

MUTATION_ADDSESSIONDATA = """
  mutation addSessionData (\$sessionData: SessionData!) {
    addSessionData(sessionData: \$sessionData)
  }"""

MUTATION_SOLVESESSION = """
  mutation solveSession (\$client: ClientInput!) {
    solveSession(client: \$client)
  }"""

MUTATION_SOLVEFEDERATED = """
  mutation solveFederated (\$scope: ScopeInput!) {
    solveFederated(scope: \$scope)
  }"""

MUTATION_DEMOCANONICALHEXAGONAL = """
  mutation demoCanonicalHexagonal (\$client: ClientInput!) {
    demoCanonicalHexagonal(client: \$client)
  }"""

MUTATION_CREATE_UPLOAD = """
  mutation CreateUpload (\$file: FileInput!, \$parts: Int) {
    createUpload(file: \$file, parts: \$parts) {
      uploadId
      file {
        id
        filename
        filesize
      }
      parts {
        partNumber
        url
      }
    }
  }"""

MUTATION_ABORT_UPLOAD = """
  mutation AbortUpload (\$fileId: ID!, \$uploadId: ID!) {
    abortUpload(fileId: \$fileId, uploadId: \$uploadId)
  }"""

MUTATION_COMPLETE_UPLOAD = """
  mutation CompleteUpload (\$fileId: ID!, \$completedUpload: CompletedUploadInput!) {
    completeUpload(fileId: \$fileId, completedUpload: \$completedUpload)
  }"""

MUTATION_PROC_CALIBRATION = """
  mutation ProcessCalibration (\$fileId: ID!) {
    procCalibration(fileId: \$fileId)
  }"""

SUBSCRIPTION_UPDATES = """
  subscription TrackEvents(\$client: ClientInput!) {
    mutationUpdate(client: \$client) {
      requestId,
      action,
      state,
      timestamp
    }
  }"""


text:
using Diversity
using Diversity.ShortNames
using Diversity.API
using LinearAlgebra
using DataFrames
using EcoBase

"""
    generalisedrichness(level::DiversityLevel, proportions::AbstractArray,
                        Z::AbstractMatrix)
    generalisedrichness(level::DiversityLevel, proportions::AbstractArray,
                        sim::AbstractTypes)

Calculates species richness (diversity at q = 0) of a series of
columns representing subcommunity counts, allowing a similarity matrix
for the types / species.

# Arguments:

- `level`: DiversityLevel to calculate at (e.g. subcommunityDiversity)

- `proportions`: population proportions

- `Z`: similarity matrix or
- `sim`: instance of AbstractTypes

# Returns:

- diversity (at ecosystem level) or diversities (of subcommunities)
"""
generalisedrichness(level::DiversityLevel,
                    proportions::AbstractArray,
                    Z::AbstractMatrix = Matrix(1.0I, size(proportions, 1), size(proportions, 1))) =
    generalisedrichness(level, proportions, GeneralTypes(Z))
function generalisedrichness(level::DiversityLevel,
                             proportions::AbstractArray,
                             sim::AbstractTypes)
    if (level == subcommunityDiversity)
        dm = ᾱ
    elseif (level == metacommunityDiversity)
        dm = Gamma
    else
        error("Can't calculate richness for $level")
    end
    gr=level(dm(Metacommunity(proportions, sim)), 0)
    gr[!,:measure] .= "Richness"
    return gr
end

"""
    richness(proportions::AbstractMatrix)

Calculates species richness (diversity at q = 0) of a series of
columns representing independent subcommunity counts.

#### Arguments:
- `proportions`: population proportions

#### Returns:
- diversities of subcommunities
"""
function richness(proportions::AbstractVecOrMat)
    gr = generalisedrichness(subcommunityDiversity, proportions,
                             UniqueTypes(size(proportions, 1)))
    gr[!,:diversity] .= Int.(round.(gr[!,:diversity]))
    return gr
end
function richness(asm::EcoBase.AbstractAssemblage)
    hassimilarity(asm) && error("function cannot run with $(typeof(gettypes(asm))) types as contains similarity")
    return richness(occurrences(asm))
end

"""
    generalisedshannon(level::DiversityLevel, proportions::AbstractArray,
                       Z::AbstractMatrix)
    generalisedshannon(level::DiversityLevel, proportions::AbstractArray,
                       sim::AbstractTypes)

Calculates Shannon entropy (log of diversity at q = 1) of a series of
columns representing independent subcommunity counts, allowing a
similarity matrix for the types / species.

# Arguments:
- `level`: DiversityLevel to calculate at (e.g. subcommunityDiversity)

- `proportions`: population proportions

- `Z`: similarity matrix or
- `sim`: instance of AbstractTypes

#### Returns:
- entropy (at metacommunity level) or entropies (of subcommunities)
"""
function generalisedshannon end
generalisedshannon(level::DiversityLevel,
                   proportions::AbstractArray,
                   Z::AbstractMatrix = Matrix(1.0I, size(proportions, 1), size(proportions, 1))) =
    generalisedshannon(level, proportions, GeneralTypes(Z))

function generalisedshannon(level::DiversityLevel,
                            proportions::AbstractArray,
                            sim::AbstractTypes)
    if (level == subcommunityDiversity)
        dm = ᾱ
    elseif (level == metacommunityDiversity)
        dm = Gamma
    else
        error("Can't calculate richness for $level")
    end
    gs = level(dm(Metacommunity(proportions, sim)), 1)
    gs[!,:diversity] .= log.(gs[!,:diversity])
    gs[!,:measure] .= "Shannon"
    select!(gs, Not(:q))
    return gs
end

"""
    shannon(proportions::AbstractVecOrMat)

Calculates shannon entropy (log of diversity at q = 1) of a series of
columns representing independent subcommunity counts.

#### Arguments:
- `proportions`: population proportions

#### Returns:
- entropies of subcommunities
"""
shannon(proportions::AbstractVecOrMat) =
    generalisedshannon(subcommunityDiversity, proportions,
                       UniqueTypes(size(proportions, 1)))
function shannon(asm::EcoBase.AbstractAssemblage)
    hassimilarity(asm) && error("function cannot run with $(typeof(gettypes(asm))) types as contains similarity")
    return shannon(occurrences(asm))
end

"""
    generalisedsimpson(level::DiversityLevel, proportions::AbstractArray,
                       Z::AbstractMatrix)
    generalisedsimpson(level::DiversityLevel, proportions::AbstractArray,
                       sim::AbstractTypes)

Calculates Simpson's index (1 / diversity at q = 2) of a series of
columns representing independent subcommunity counts, allowing a
similarity matrix for the types / species.

#### Arguments:
- `level`: DiversityLevel to calculate at (e.g. subcommunityDiversity)

- `proportions`: population proportions

- `Z`: similarity matrix or
- `sim`: instance of AbstractTypes

#### Returns:
- concentration (at ecosystem level) or concentrations (of subcommunities)
"""
function generalisedsimpson end
generalisedsimpson(level::DiversityLevel,
                   proportions::AbstractArray,
                   Z::AbstractMatrix = Matrix(1.0I, size(proportions, 1), size(proportions, 1))) =
    generalisedsimpson(level, proportions, GeneralTypes(Z))

function generalisedsimpson(level::DiversityLevel,
                            proportions::AbstractArray,
                            sim::AbstractTypes)
    if (level == subcommunityDiversity)
        dm = ᾱ
    elseif (level == metacommunityDiversity)
        dm = Gamma
    else
        error("Can't calculate richness for $level")
    end
    gs = level(dm(Metacommunity(proportions, sim)), 2)
    gs[!,:diversity] .= gs[!,:diversity] .^ -1
    gs[!,:measure] .= "Simpson"
    select!(gs, Not(:q))
    return gs
end

"""
    simpson(proportions::AbstractMatrix)

Calculates Simpson's index (1 / diversity at q = 2) of a series of
columns representing independent subcommunity counts.

# Arguments:

- `proportions`: population proportions

# Returns:

- concentrations of subcommunities
"""
simpson(proportions::AbstractVecOrMat) =
    generalisedsimpson(subcommunityDiversity, proportions,
                       UniqueTypes(size(proportions, 1)))
function simpson(asm::EcoBase.AbstractAssemblage)
    hassimilarity(asm) && error("function cannot run with $(typeof(gettypes(asm))) types as contains similarity")
    return simpson(occurrences(asm))
end

"""
    generalisedjaccard(proportions::AbstractArray, qs, Z::AbstractMatrix)
    generalisedjaccard(proportions::AbstractArray, qs, sim::AbstractTypes)

Calculates a generalisation of the Jaccard index of two columns
representing the counts of two subcommunities. This evaluates to raw
alpha / gamma - 1 for a series of orders, repesented as a vector of qs
(or a single number). It also includes an optional similarity matrix
for the species. This gives a measure of the distinctness of the
subcommunities, though we believe that beta and normalised beta have
better properties.

# Arguments:

- `proportions`: population proportions

- `qs`: single number or vector of values of parameter q

- `Z`: similarity matrix or
- `sim`: instance of AbstractTypes

# Returns:

- Jaccard-related distinctivess measures
"""
function generalisedjaccard end

generalisedjaccard(proportions::AbstractMatrix, qs,
                   Z::AbstractMatrix =
                   Matrix(1.0I, size(proportions, 1),
                                size(proportions, 1))) =
    generalisedjaccard(proportions, qs, GeneralTypes(Z))

function generalisedjaccard(proportions::AbstractMatrix, qs,
                            sim::AbstractTypes)
    meta = Metacommunity(proportions, sim)
    countsubcommunities(meta) == 2 ||
    error("Can only calculate Jaccard index for 2 subcommunities")
    ab = metadiv(α(meta), qs)
    g = metadiv(Γ(meta), qs)
    j = innerjoin(ab, g, on=[:q, :type_level, :type_name,
                        :partition_level, :partition_name, :div_type],
             makeunique=true)
    j[!,:diversity] .= j[!,:diversity] ./ j[!,:diversity_1] .- 1
    j[!,:measure] .= "Jaccard"
    select!(j, Not([:diversity_1, :measure_1]))
    return j
end

"""
    jaccard(proportions::AbstractMatrix)

Calculates Jaccard index (Jaccard similarity coefficient) of two
columns representing independent subcommunity counts, which is
normmetaalpha(proportions, 0) / metagamma(proportions, 0) - 1

# Arguments:

- `proportions`: population proportions

# Returns:

- the Jaccard index
"""
jaccard(proportions::AbstractMatrix) =
    generalisedjaccard(proportions, 0,
                       UniqueTypes(size(proportions, 1)))
function jaccard(asm::EcoBase.AbstractAssemblage)
    hassimilarity(asm) && error("function cannot run with $(typeof(gettypes(asm))) types as contains similarity")
    return jaccard(occurrences(asm))
end


text:
using Dash

app = dash()

app.layout = html_div() do
    html_h1("hello us")

end

# run_server(app, "0.0.0.0", debug = true)
port = parse(Int64, ENV["PORT"])
run_server(app, "0.0.0.0", port)


text:
type MUtilTimeOutEx <: Exception
end

type SyncObjRef{T}
    rr::RemoteRef
    SyncObjRef() = new(RemoteRef())
end

function call_on_owner(f, h::SyncObjRef, args...; kw...)
    if h.rr.where == myid()
        f(h, args...; kw...)
    else
        remotecall_fetch(h.rr.where, () -> f(h, args...; kw...))
    end
end

# Just a wrapper for a RemoteRef - so that we can support our own
# methods similar to that supported for RemoteRefs

type SyncObjData
    cantake::Condition   # waiting for a value
    canput::Condition  # waiting for value to be removed

    SyncObjData() = new (Condition(), Condition())
end

# data stored by the owner of a RemoteRef
abstract AbstractRemoteSyncObj

function syncobj(T::Type, args...; kw...)
    h = SyncObjRef{T}()
    put!(h.rr, T(args...; kw...))
    h
end

function syncobj_create(pid, T, args...; kw...)
    remotecall_fetch(pid, ()->syncobj(T, args...; kw...))
end

function setup_waittimer(condvar, timeout)
    t1 = time()
    if isv4
        t = Timer(x -> notify(condvar), timeout, 0.0)
    else
        t = Timer(x -> notify(condvar))
        start_timer(t, timeout, 0.0)
    end
    (t1, t)
end

function kwextract(kw, fld, default)
    for (k,v) in kw
        if (k == fld)
            return v
        end
    end
    return default
end

function wait_cantake (rv::AbstractRemoteSyncObj, args...; kw...)
    timeout = kwextract(kw, :timeout, 0.0)
    (t1, t) = (timeout > 0.0) ? setup_waittimer(rv.so.cantake, timeout) : (0, nothing)

    try
        while !rv.cantake(rv, args...)
            # Multiple tasks are waiting on the same condition variable.
            # Consequently, we may get woken up but not find any data for ourselves.
            (timeout > 0.0) && ((time() - t1) > timeout) && throw(MUtilTimeOutEx())
            wait( rv.so.cantake )
        end
    finally
        if isv4
            (t != nothing) && close(t)
        else
            (t != nothing) && stop_timer(t)
        end
    end
    return ( rv.fetch(rv, args...) )
end

function wait_canput(rv::AbstractRemoteSyncObj, args...; kw...)
    timeout = kwextract(kw, :timeout, 0.0)
    (t1, t) = (timeout > 0.0) ? setup_waittimer(rv.so.canput, timeout) : (0, nothing)

    try
        while !rv.canput(rv, args...)
            (timeout > 0.0) && ((time() - t1) > timeout) && throw(MUtilTimeOutEx())
            wait(rv.so.canput)
        end
    finally
        if isv4
            (t != nothing) && close(t)
        else
            (t != nothing) && stop_timer(t)
        end
    end
    return nothing
end

wait_ref(h, args...; kw...) = (wait_cantake(fetch(h.rr), args...; kw...); nothing)
wait(h::SyncObjRef, args...; kw...) = (call_on_owner(wait_ref, h, args...; kw...); h)

fetch_ref(h, args...; kw...) = wait_cantake(fetch(h.rr), args...; kw...)
fetch(h::SyncObjRef, args...; kw...) = call_on_owner(fetch_ref, h, args...; kw...)

# storing a value to a Ref
function put!(rv::AbstractRemoteSyncObj, args...; kw...)
    wait_canput(rv, args...; kw...)
    rv.put(rv, args...; kw...)
    notify(rv.so.cantake, rv.fetch(rv))
end

put_ref(h, args...; kw...) = put!(fetch(h.rr), args...; kw...)
put!(h::SyncObjRef, args...; kw...) = (call_on_owner(put_ref, h, args...; kw...); h)

function take!(rv::AbstractRemoteSyncObj, args...; kw...)
    wait_cantake(rv, args...; kw...)
    val = rv.take(rv, args...)
    notify(rv.so.canput)
    val
end

take_ref(h; kw...) = take!(fetch(h.rr); kw...)
take_ref(h, key; kw...) = take!(fetch(h.rr), key; kw...)
take!(h::SyncObjRef; kw...) = call_on_owner(take_ref, h; kw...)
take!(h::SyncObjRef, key; kw...) = call_on_owner(take_ref, h, key; kw...)

query_ref(h, args...) = query(fetch(h.rr), args...)
query(h::SyncObjRef, args...) = call_on_owner(query_ref, h, args...)
query(rv::AbstractRemoteSyncObj, args...) = rv.query(rv, args...)

function isready(h::SyncObjRef, args...)
    remotecall_fetch(h.rr.where, (rr, rargs...) -> begin rv = fetch(rr); rv.cantake(rv, rargs...) end, h.rr, args...)
end


text:
module TestExtensions

import Test: @test, test_expr!, get_test_result, do_test, Returned

"""
    @test test_expr error_expr
    @test f(args...) error_expr key=val...

Exactly like `@test`, but evaluates `error_expr` if `@test` did not `Pass`.

# Examples
```jldoctest
julia> @test true 37+5
Test Passed
julia> @test false 37+5
[ Info: 42
Test Failed at REPL[34]:1
  Expression: false
ERROR: There was an error during testing
```
"""
macro test(ex, msg, kws...)
    test_expr!("@test", ex, kws...)
    orig_ex = Expr(:inert, ex)
    result = get_test_result(ex, __source__)
    quote
        if !(isa($result, Returned) && isa($result.value, Bool) && $result.value)
            @info $msg
        end
        do_test($result, $orig_ex)
    end
end

end # module


text:
struct CompressedBuffer
    data::Vector{UInt8}
    uncompressedlength::Int64
end

struct Compressed{Z, A}
    data::A
    buffers::Vector{CompressedBuffer}
    len::Int64
    nullcount::Int64
    children::Vector{Compressed}
end

Base.length(c::Compressed) = c.len
Base.eltype(c::Compressed{Z, A}) where {Z, A} = eltype(A)
getmetadata(x::Compressed) = getmetadata(x.data)
compressiontype(c::Compressed{Z}) where {Z} = Z

function compress(Z::Meta.CompressionType, x::Array)
    GC.@preserve x begin
        y = unsafe_wrap(Array, convert(Ptr{UInt8}, pointer(x)), sizeof(x))
        if Z == Meta.CompressionType.LZ4_FRAME
            return CompressedBuffer(transcode(LZ4FrameCompressor, y), length(y))
        elseif Z == Meta.CompressionType.ZSTD
            return CompressedBuffer(transcode(ZstdCompressor, y), length(y))
        else
            throw(ArgumentError("unsupported compression type: $Z"))
        end
    end
end

compress(Z::Meta.CompressionType, x) = compress(Z, convert(Array, x))

abstract type ArrowVector{T} <: AbstractVector{T} end

Base.IndexStyle(::Type{A}) where {A <: ArrowVector} = Base.IndexLinear()
Base.similar(::Type{A}, dims::Dims) where {T, A <: ArrowVector{T}} = Vector{T}(undef, dims)
validitybitmap(x::ArrowVector) = x.validity
nullcount(x::ArrowVector) = validitybitmap(x).nc
getmetadata(x::ArrowVector) = x.metadata

function toarrowvector(x, de=DictEncoding[], meta=getmetadata(x); compression::Union{Nothing, Symbol}=nothing, kw...)
    @debug 2 "converting top-level column to arrow format: col = $(typeof(x)), compression = $compression, kw = $kw"
    @debug 3 x
    A = arrowvector(x, de, meta; compression=compression, kw...)
    if compression === :lz4
        A = compress(Meta.CompressionType.LZ4_FRAME, A)
    elseif compression === :zstd
        A = compress(Meta.CompressionType.ZSTD, A)
    end
    @debug 2 "converted top-level column to arrow format: $(typeof(A))"
    @debug 3 A
    return A
end

function arrowvector(x, de, meta; dictencoding::Bool=false, dictencode::Bool=false, kw...)
    if !(x isa DictEncode) && !dictencoding && (dictencode || (x isa AbstractArray && DataAPI.refarray(x) !== x))
        x = DictEncode(x)
    end
    T = eltype(x)
    S = maybemissing(T)
    return arrowvector(S, T, x, de, meta; kw...)
end

# conversions to arrow types
arrowvector(::Type{Dates.Date}, ::Type{S}, x, de, meta; kw...) where {S} =
    arrowvector(converter(DATE, x), de, meta; kw...)
arrowvector(::Type{Dates.Time}, ::Type{S}, x, de, meta; kw...) where {S} =
    arrowvector(converter(TIME, x), de, meta; kw...)
arrowvector(::Type{Dates.DateTime}, ::Type{S}, x, de, meta; kw...) where {S} =
    arrowvector(converter(DATETIME, x), de, meta; kw...)
arrowvector(::Type{P}, ::Type{S}, x, de, meta; kw...) where {P <: Dates.Period, S} =
    arrowvector(converter(Duration{arrowperiodtype(P)}, x), de, meta; kw...)

# fallback that calls ArrowType
function arrowvector(::Type{S}, ::Type{T}, x, de, meta; kw...) where {S, T}
    if ArrowTypes.istyperegistered(S)
        meta = meta === nothing ? Dict{String, String}() : meta
        arrowtype = ArrowTypes.getarrowtype!(meta, S)
        return arrowvector(converter(arrowtype, x), de, meta; kw...)
    end
    return arrowvector(ArrowType(S), S, T, x, de, meta; kw...)
end

arrowvector(::NullType, ::Type{Missing}, ::Type{Missing}, x, de, meta; kw...) = MissingVector(length(x))
compress(Z::Meta.CompressionType, v::MissingVector) =
    Compressed{Z, MissingVector}(v, CompressedBuffer[], length(v), length(v), Compressed[])

struct ValidityBitmap <: ArrowVector{Bool}
    bytes::Vector{UInt8} # arrow memory blob
    pos::Int # starting byte of validity bitmap
    ℓ::Int # # of _elements_ (not bytes!) in bitmap (because bitpacking)
    nc::Int # null count
end

compress(Z::Meta.CompressionType, v::ValidityBitmap) =
    v.nc == 0 ? CompressedBuffer(UInt8[], 0) : compress(Z, view(v.bytes, v.pos:(v.pos + bitpackedbytes(v.ℓ) - 1)))

Base.size(p::ValidityBitmap) = (p.ℓ,)
nullcount(x::ValidityBitmap) = x.nc

function ValidityBitmap(len::Integer)
    return ValidityBitmap(UInt8[], 1, len, 0)
end

function ValidityBitmap(x)
    len = length(x)
    blen = bitpackedbytes(len)
    bytes = Vector{UInt8}(undef, blen)
    st = iterate(x)
    i = 0
    nc = 0
    for k = 1:blen
        b = 0x00
        for j = 1:8
            if (i + j) <= len
                y, state = st
                if !y
                    nc += 1
                    b = setbit(b, false, j)
                else
                    b = setbit(b, true, j)
                end
                st = iterate(x, state)
            end
        end
        i += 8
        @inbounds bytes[k] = b
    end
    return ValidityBitmap(nc == 0 ? UInt8[] : bytes, 1, nc == 0 ? 0 : len, nc)
end

@propagate_inbounds function Base.getindex(p::ValidityBitmap, i::Integer)
    # no boundscheck because parent array should do it
    # if a validity bitmap is empty, it either means:
    #   1) the parent array null_count is 0, so all elements are valid
    #   2) parent array is also empty, so "all" elements are valid
    p.nc == 0 && return true
    # translate element index to bitpacked byte index
    a, b = fldmod1(i, 8)
    @inbounds byte = p.bytes[p.pos + a - 1]
    # check individual bit of byte
    return getbit(byte, b)
end

@propagate_inbounds function Base.setindex!(p::ValidityBitmap, v, i::Integer)
    x = convert(Bool, v)
    p.ℓ == 0 && !x && throw(BoundsError(p, i))
    a, b = fldmod1(i, 8)
    @inbounds byte = p.bytes[p.pos + a - 1]
    @inbounds p.bytes[p.pos + a - 1] = setbit(byte, x, b)
    return v
end

struct Primitive{T, S, A} <: ArrowVector{T}
    arrow::Vector{UInt8} # need to hold a reference to arrow memory blob
    validity::ValidityBitmap
    data::A
    ℓ::Int64
    metadata::Union{Nothing, Dict{String, String}}
end

Primitive(::Type{T}, b::Vector{UInt8}, v::ValidityBitmap, data::A, l::Int, meta) where {T, A} =
    Primitive{T, eltype(A), A}(b, v, data, l, meta)

Base.size(p::Primitive) = (p.ℓ,)

function arrowvector(::PrimitiveType, ::Type{T}, ::Type{S}, x, de, meta; kw...) where {T, S}
    len = length(x)
    if T !== S
        validity = ValidityBitmap((y !== missing for y in x))
    else
        validity = ValidityBitmap(len)
    end
    return Primitive{S, T, typeof(x)}(UInt8[], validity, x, len, meta)
end

function Base.copy(p::Primitive{T, S}) where {T, S}
    if T === S
        return copy(p.data)
    else
        return convert(Array, p)
    end
end

@propagate_inbounds function Base.getindex(p::Primitive{T, S}, i::Integer) where {T, S}
    @boundscheck checkbounds(p, i)
    if T >: Missing
        return @inbounds (p.validity[i] ? ArrowTypes.arrowconvert(T, p.data[i]) : missing)
    else
        return @inbounds ArrowTypes.arrowconvert(T, p.data[i])
    end
end

@propagate_inbounds function Base.setindex!(p::Primitive{T, S}, v, i::Integer) where {T, S}
    @boundscheck checkbounds(p, i)
    if T >: Missing
        if v === missing
            @inbounds p.validity[i] = false
        else
            @inbounds p.data[i] = convert(S, v)
        end
    else
        @inbounds p.data[i] = convert(S, v)
    end
    return v
end

function compress(Z::Meta.CompressionType, p::P) where {P <: Primitive}
    len = length(p)
    nc = nullcount(p)
    validity = compress(Z, p.validity)
    data = compress(Z, p.data)
    return Compressed{Z, P}(p, [validity, data], len, nc, Compressed[])
end

struct Offsets{T <: Union{Int32, Int64}} <: ArrowVector{Tuple{T, T}}
    arrow::Vector{UInt8} # need to hold a reference to arrow memory blob
    offsets::Vector{T}
end

Base.size(o::Offsets) = (length(o.offsets) - 1,)

@propagate_inbounds function Base.getindex(o::Offsets, i::Integer)
    @boundscheck checkbounds(o, i)
    @inbounds lo = o.offsets[i] + 1
    @inbounds hi = o.offsets[i + 1]
    return lo, hi
end

struct List{T, O, A} <: ArrowVector{T}
    arrow::Vector{UInt8} # need to hold a reference to arrow memory blob
    validity::ValidityBitmap
    offsets::Offsets{O}
    data::A
    ℓ::Int
    metadata::Union{Nothing, Dict{String, String}}
end

Base.size(l::List) = (l.ℓ,)

function arrowvector(::ListType, ::Type{T}, ::Type{S}, x, de, meta; largelists::Bool=false, kw...) where {T, S}
    len = length(x)
    if S !== T
        validity = ValidityBitmap((y !== missing for y in x))
    else
        validity = ValidityBitmap(len)
    end
    flat = ToList(x; largelists=largelists)
    offsets = Offsets(UInt8[], flat.inds)
    if eltype(flat) == UInt8
        data = flat
    else
        data = arrowvector(flat, de, nothing; lareglists=largelists, kw...)
    end
    return List{S, eltype(flat.inds), typeof(data)}(UInt8[], validity, offsets, data, len, meta)
end

@propagate_inbounds function Base.getindex(l::List{T}, i::Integer) where {T}
    @boundscheck checkbounds(l, i)
    @inbounds lo, hi = l.offsets[i]
    if ArrowTypes.isstringtype(T)
        if Base.nonmissingtype(T) !== T
            return l.validity[i] ? ArrowTypes.arrowconvert(T, unsafe_string(pointer(l.data, lo), hi - lo + 1)) : missing
        else
            return ArrowTypes.arrowconvert(T, unsafe_string(pointer(l.data, lo), hi - lo + 1))
        end
    elseif Base.nonmissingtype(T) !== T
        return l.validity[i] ? ArrowTypes.arrowconvert(T, view(l.data, lo:hi)) : missing
    else
        return ArrowTypes.arrowconvert(T, view(l.data, lo:hi))
    end
end

# @propagate_inbounds function Base.setindex!(l::List{T}, v, i::Integer) where {T}

# end

function compress(Z::Meta.CompressionType, x::List{T, O, A}) where {T, O, A}
    len = length(x)
    nc = nullcount(x)
    validity = compress(Z, x.validity)
    offsets = compress(Z, x.offsets.offsets)
    buffers = [validity, offsets]
    children = Compressed[]
    if eltype(A) == UInt8
        push!(buffers, compress(Z, x.data))
    else
        push!(children, compress(Z, x.data))
    end
    return Compressed{Z, typeof(x)}(x, buffers, len, nc, children)
end

struct FixedSizeList{T, A <: AbstractVector} <: ArrowVector{T}
    arrow::Vector{UInt8} # need to hold a reference to arrow memory blob
    validity::ValidityBitmap
    data::A
    ℓ::Int
    metadata::Union{Nothing, Dict{String, String}}
end

Base.size(l::FixedSizeList) = (l.ℓ,)

function arrowvector(::FixedSizeListType, ::Type{T}, ::Type{S}, x, de, meta; kw...) where {T, S}
    len = length(x)
    if S !== T
        validity = ValidityBitmap((y !== missing for y in x))
    else
        validity = ValidityBitmap(len)
    end
    flat = ToFixedSizeList(x)
    if eltype(flat) == UInt8
        data = flat
    else
        data = arrowvector(flat, de, nothing; kw...)
    end
    return FixedSizeList{S, typeof(data)}(UInt8[], validity, data, len, meta)
end

@propagate_inbounds function Base.getindex(l::FixedSizeList{T}, i::Integer) where {T}
    @boundscheck checkbounds(l, i)
    N = ArrowTypes.getsize(Base.nonmissingtype(T))
    off = (i - 1) * N
    if Base.nonmissingtype(T) !== T
        return l.validity[i] ? ArrowTypes.arrowconvert(T, ntuple(j->l.data[off + j], N)) : missing
    else
        return ArrowTypes.arrowconvert(T, ntuple(j->l.data[off + j], N))
    end
end

@propagate_inbounds function Base.setindex!(l::FixedSizeList{T}, v::T, i::Integer) where {T}
    @boundscheck checkbounds(l, i)
    if v === missing
        @inbounds l.validity[i] = false
    else
        N = ArrowTypes.getsize(Base.nonmissingtype(T))
        off = (i - 1) * N
        foreach(1:N) do j
            @inbounds l.data[off + j] = v[j]
        end
    end
    return v
end

function compress(Z::Meta.CompressionType, x::FixedSizeList{T, A}) where {T, A}
    len = length(x)
    nc = nullcount(x)
    validity = compress(Z, x.validity)
    buffers = [validity]
    children = Compressed[]
    if eltype(A) == UInt8
        push!(buffers, compress(Z, x.data))
    else
        push!(children, compress(Z, x.data))
    end
    return Compressed{Z, typeof(x)}(x, buffers, len, nc, children)
end

struct Map{T, O, A} <: ArrowVector{T}
    validity::ValidityBitmap
    offsets::Offsets{O}
    data::A
    ℓ::Int
    metadata::Union{Nothing, Dict{String, String}}
end

Base.size(l::Map) = (l.ℓ,)

function arrowvector(::MapType, ::Type{T}, ::Type{S}, x, de, meta; largelists::Bool=false, kw...) where {T, S}
    len = length(x)
    if S !== T
        validity = ValidityBitmap((y !== missing for y in x))
    else
        validity = ValidityBitmap(len)
    end
    flat = ToMap(x; largelists=largelists)
    offsets = Offsets(UInt8[], flat.inds)
    data = arrowvector(flat, de, nothing; lareglists=largelists, kw...)
    return Map{S, eltype(flat.inds), typeof(data)}(validity, offsets, data, len, meta)
end

@propagate_inbounds function Base.getindex(l::Map{T}, i::Integer) where {T}
    @boundscheck checkbounds(l, i)
    @inbounds lo, hi = l.offsets[i]
    if Base.nonmissingtype(T) !== T
        return l.validity[i] ? ArrowTypes.arrowconvert(T, Dict(x.key => x.value for x in view(l.data, lo:hi))) : missing
    else
        return ArrowTypes.arrowconvert(T, Dict(x.key => x.value for x in view(l.data, lo:hi)))
    end
end

function compress(Z::Meta.CompressionType, x::A) where {A <: Map}
    len = length(x)
    nc = nullcount(x)
    validity = compress(Z, x.validity)
    offsets = compress(Z, x.offsets.offsets)
    buffers = [validity, offsets]
    children = Compressed[]
    push!(children, compress(Z, x.data))
    return Compressed{Z, A}(x, buffers, len, nc, children)
end

struct Struct{T, S} <: ArrowVector{T}
    validity::ValidityBitmap
    data::S # Tuple of ArrowVector
    ℓ::Int
    metadata::Union{Nothing, Dict{String, String}}
end

Base.size(s::Struct) = (s.ℓ,)

function arrowvector(::StructType, ::Type{T}, ::Type{S}, x, de, meta; kw...) where {T, S}
    len = length(x)
    if S !== T
        validity = ValidityBitmap((y !== missing for y in x))
    else
        validity = ValidityBitmap(len)
    end
    data = Tuple(arrowvector(ToStruct(x, i), de, nothing; kw...) for i = 1:fieldcount(T))
    return Struct{S, typeof(data)}(validity, data, len, meta)
end

@propagate_inbounds function Base.getindex(s::Struct{T}, i::Integer) where {T}
    @boundscheck checkbounds(s, i)
    NT = Base.nonmissingtype(T)
    if ArrowTypes.structtype(NT) === ArrowTypes.NAMEDTUPLE
        if NT !== T
            return s.validity[i] ? NT(ntuple(j->s.data[j][i], fieldcount(NT))) : missing
        else
            return NT(ntuple(j->s.data[j][i], fieldcount(NT)))
        end
    elseif ArrowTypes.structtype(NT) === ArrowTypes.STRUCT
        if NT !== T
            return s.validity[i] ? NT(ntuple(j->s.data[j][i], fieldcount(NT))...) : missing
        else
            return NT(ntuple(j->s.data[j][i], fieldcount(NT))...)
        end
    end
end

@propagate_inbounds function Base.setindex!(s::Struct{T}, v::T, i::Integer) where {T}
    @boundscheck checkbounds(s, i)
    if v === missing
        @inbounds s.validity[i] = false
    else
        NT = Base.nonmissingtype(T)
        N = getN(NT)
        foreach(1:N) do j
            @inbounds s.data[j][i] = v[j]
        end
    end
    return v
end

function compress(Z::Meta.CompressionType, x::A) where {A <: Struct}
    len = length(x)
    nc = nullcount(x)
    validity = compress(Z, x.validity)
    buffers = [validity]
    children = Compressed[]
    for y in x.data
        push!(children, compress(Z, y))
    end
    return Compressed{Z, A}(x, buffers, len, nc, children)
end

struct DenseUnion{T, S} <: ArrowVector{T}
    arrow::Vector{UInt8} # need to hold a reference to arrow memory blob
    arrow2::Vector{UInt8} # need to hold a reference to arrow memory blob
    typeIds::Vector{UInt8}
    offsets::Vector{Int32}
    data::S # Tuple of ArrowVector
    metadata::Union{Nothing, Dict{String, String}}
end

Base.size(s::DenseUnion) = size(s.typeIds)
nullcount(x::DenseUnion) = 0

arrowvector(U::Union, ::Type{S}, x, de, meta; denseunions::Bool=true, kw...) where {S} =
    arrowvector(denseunions ? DenseUnionVector(x) : SparseUnionVector(x), de, meta; denseunions=denseunions, kw...)

function arrowvector(::UnionType, ::Type{UT}, ::Type{S}, x, de, meta; kw...) where {UT <: UnionT, S}
    if unionmode(UT) == Meta.UnionMode.Dense
        x = x isa DenseUnionVector ? x.itr : x
        typeids, offsets, data = todense(UT, x)
        data2 = map(y -> arrowvector(y, de, nothing; kw...), data)
        return DenseUnion{UT, typeof(data2)}(UInt8[], UInt8[], typeids, offsets, data2, meta)
    else
        x = x isa SparseUnionVector ? x.itr : x
        typeids = sparsetypeids(UT, x)
        data3 = Tuple(arrowvector(ToSparseUnion(fieldtype(eltype(UT), i), x), de, nothing; kw...) for i = 1:fieldcount(eltype(UT)))
        return SparseUnion{UT, typeof(data3)}(UInt8[], typeids, data3, meta)
    end
end

@propagate_inbounds function Base.getindex(s::DenseUnion{T}, i::Integer) where {T}
    @boundscheck checkbounds(s, i)
    @inbounds typeId = s.typeIds[i]
    @inbounds off = s.offsets[i]
    @inbounds x = s.data[typeId + 1][off + 1]
    return x
end

@propagate_inbounds function Base.setindex!(s::DenseUnion{UnionT{T, typeIds, U}}, v, i::Integer) where {T, typeIds, U}
    @boundscheck checkbounds(s, i)
    @inbounds typeId = s.typeIds[i]
    typeids = typeIds === nothing ? (0:(fieldcount(U) - 1)) : typeIds
    vtypeId = Int8(typeids[isatypeid(v, U)])
    if typeId == vtypeId
        @inbounds off = s.offsets[i]
        @inbounds s.data[typeId +1][off + 1] = v
    else
        throw(ArgumentError("type of item to set $(typeof(v)) must match existing item $(fieldtype(U, typeid))"))
    end
    return v
end

function compress(Z::Meta.CompressionType, x::A) where {A <: DenseUnion}
    len = length(x)
    nc = nullcount(x)
    typeIds = compress(Z, x.typeIds)
    offsets = compress(Z, x.offsets)
    buffers = [typeIds, offsets]
    children = Compressed[]
    for y in x.data
        push!(children, compress(Z, y))
    end
    return Compressed{Z, A}(x, buffers, len, nc, children)
end

struct SparseUnion{T, S} <: ArrowVector{T}
    arrow::Vector{UInt8} # need to hold a reference to arrow memory blob
    typeIds::Vector{UInt8}
    data::S # Tuple of ArrowVector
    metadata::Union{Nothing, Dict{String, String}}
end

Base.size(s::SparseUnion) = size(s.typeIds)
nullcount(x::SparseUnion) = 0

@propagate_inbounds function Base.getindex(s::SparseUnion{T}, i::Integer) where {T}
    @boundscheck checkbounds(s, i)
    @inbounds typeId = s.typeIds[i]
    @inbounds x = s.data[typeId + 1][i]
    return x
end

@propagate_inbounds function Base.setindex!(s::SparseUnion{UnionT{T, typeIds, U}}, v, i::Integer) where {T, typeIds, U}
    @boundscheck checkbounds(s, i)
    typeids = typeIds === nothing ? (0:(fieldcount(U) - 1)) : typeIds
    vtypeId = Int8(typeids[isatypeid(v, U)])
    @inbounds s.typeIds[i] = vtypeId
    @inbounds s.data[vtypeId + 1][i] = v
    return v
end

function compress(Z::Meta.CompressionType, x::A) where {A <: SparseUnion}
    len = length(x)
    nc = nullcount(x)
    typeIds = compress(Z, x.typeIds)
    buffers = [typeIds]
    children = Compressed[]
    for y in x.data
        push!(children, compress(Z, y))
    end
    return Compressed{Z, A}(x, buffers, len, nc, children)
end

mutable struct DictEncoding{T, A} <: ArrowVector{T}
    id::Int64
    data::A
    isOrdered::Bool
end

Base.size(d::DictEncoding) = size(d.data)

@propagate_inbounds function Base.getindex(d::DictEncoding{T}, i::Integer) where {T}
    @boundscheck checkbounds(d, i)
    return @inbounds ArrowTypes.arrowconvert(T, d.data[i])
end

# convenience wrapper to signal that an input column should be
# dict encoded when written to the arrow format
struct DictEncodeType{T} end
getT(::Type{DictEncodeType{T}}) where {T} = T

struct DictEncode{T, A} <: AbstractVector{DictEncodeType{T}}
    data::A
end

DictEncode(x::A) where {A} = DictEncode{eltype(A), A}(x)
Base.IndexStyle(::Type{<:DictEncode}) = Base.IndexLinear()
Base.size(x::DictEncode) = (length(x.data),)
Base.iterate(x::DictEncode, st...) = iterate(x.data, st...)
Base.getindex(x::DictEncode, i::Int) = getindex(x.data, i)
ArrowTypes.ArrowType(::Type{<:DictEncodeType}) = DictEncodedType()

struct DictEncoded{T, S, A} <: ArrowVector{T}
    arrow::Vector{UInt8} # need to hold a reference to arrow memory blob
    validity::ValidityBitmap
    indices::Vector{S}
    encoding::DictEncoding{T, A}
    metadata::Union{Nothing, Dict{String, String}}
end

DictEncoded(b::Vector{UInt8}, v::ValidityBitmap, inds::Vector{S}, encoding::DictEncoding{T, A}, meta) where {S, T, A} =
    DictEncoded{T, S, A}(b, v, inds, encoding, meta)

Base.size(d::DictEncoded) = size(d.indices)

isdictencoded(d::DictEncoded) = true
isdictencoded(x) = false
isdictencoded(c::Compressed{Z, A}) where {Z, A <: DictEncoded} = true

signedtype(::Type{UInt8}) = Int8
signedtype(::Type{UInt16}) = Int16
signedtype(::Type{UInt32}) = Int32
signedtype(::Type{UInt64}) = Int64

indtype(d::D) where {D <: DictEncoded} = indtype(D)
indtype(::Type{DictEncoded{T, S, A}}) where {T, S, A} = signedtype(S)
indtype(c::Compressed{Z, A}) where {Z, A <: DictEncoded} = indtype(A)

getid(d::DictEncoded) = d.encoding.id
getid(c::Compressed{Z, A}) where {Z, A <: DictEncoded} = c.data.encoding.id

arrowvector(::DictEncodedType, ::Type{DictEncodeType{T}}, ::Type{DictEncodeType{S}}, x, de, meta; kw...) where {T, S} =
    arrowvector(DictEncodedType(), maybemissing(T), S, x, de, meta; kw...)

function arrowvector(::DictEncodedType, ::Type{T}, ::Type{S}, x, de, meta; dictencode::Bool=false, dictencodenested::Bool=false, kw...) where {T, S}
    @assert x isa DictEncode
    x = x.data
    len = length(x)
    if S !== T
        validity = ValidityBitmap((y !== missing for y in x))
    else
        validity = ValidityBitmap(len)
    end
    if x isa AbstractArray && DataAPI.refarray(x) !== x
        inds = copy(DataAPI.refarray(x))
        # adjust to "offset" instead of index
        for i = 1:length(inds)
            @inbounds inds[i] -= 1
        end
        pool = DataAPI.refpool(x)
        data = arrowvector(pool, de, nothing; dictencode=dictencodenested, dictencodenested=dictencodenested, dictencoding=true, kw...)
        encoding = DictEncoding{S, typeof(data)}(0, data, false)
    else
        # need to encode ourselves
        y = PooledArray(x)
        inds = DataAPI.refarray(y)
        # adjust to "offset" instead of index
        for i = 1:length(inds)
            @inbounds inds[i] = inds[i] - 1
        end
        data = arrowvector(DataAPI.refpool(y), de, nothing; dictencode=dictencodenested, dictencodenested=dictencodenested, dictencoding=true, kw...)
        encoding = DictEncoding{S, typeof(data)}(0, data, false)
    end
    push!(de, encoding)
    return DictEncoded(UInt8[], validity, inds, encoding, data.metadata)
end

@propagate_inbounds function Base.getindex(d::DictEncoded, i::Integer)
    @boundscheck checkbounds(d, i)
    @inbounds valid = d.validity[i]
    !valid && return missing
    @inbounds idx = d.indices[i]
    return @inbounds d.encoding[idx + 1]
end

@propagate_inbounds function Base.setindex!(d::DictEncoded{T}, v, i::Integer) where {T}
    @boundscheck checkbounds(d, i)
    if v === missing
        @inbounds d.validity[i] = false
    else
        ix = findfirst(d.encoding.data, v)
        if ix === nothing
            push!(d.encoding.data, v)
            @inbounds d.indices[i] = length(d.encoding.data) - 1
        else
            @inbounds d.indices[i] = ix - 1
        end
    end
    return v
end

function Base.copy(x::DictEncoded{T, S}) where {T, S}
    pool = copy(x.encoding.data)
    valid = x.validity
    inds = x.indices
    if T !== S
        refs = Vector{S}(undef, length(inds))
        @inbounds for i = 1:length(inds)
            refs[i] = ifelse(valid[i], inds[i] + one(S), missing)
        end
    else
        refs = copy(inds)
        @inbounds for i = 1:length(inds)
            refs[i] = refs[i] + one(S)
        end
    end
    return PooledArray(PooledArrays.RefArray(refs), Dict{T, S}(val => i for (i, val) in enumerate(pool)), pool)
end

function compress(Z::Meta.CompressionType, x::A) where {A <: DictEncoded}
    len = length(x)
    nc = nullcount(x)
    validity = compress(Z, x.validity)
    inds = compress(Z, x.indices)
    return Compressed{Z, A}(x, [validity, inds], len, nc, Compressed[])
end

max_stars_repo_name: StefanKarpinski/arrow

max_issues_repo_name: StefanKarpinski/arrow

max_forks_repo_name: StefanKarpinski/arrow


text:
module SchwarzSolver

using Printf
using DataStructures
using LinearAlgebra
import LinearAlgebra: BLAS
using MathOptInterface
using Plasmo
using Ipopt

export schwarz_solve

include("schwarz_solver.jl")

end # module


text:
using Documenter, PlasmaDispersionRelations

makedocs(sitename="PlasmaDispersionRelations.jl")

deploydocs(
    repo = "github.com/jwscook/PlasmaDispersionFunctions.jl.git",
)


text:
#!/usr/bin/env julia

using SpectralElements, Flux, LinearAlgebra, Plots, UnPack
import Zygote
#----------------------------------#
function setIC!(u,x,y,t)
    u .= 0.0
    return
end

function setBC!(ub,x,y,t)
    ub .= @. 0+0*x
    return
end

function setForcing!(f,x,y,t)
    f .= @. 1+0*x
    return
end

function setVisc!(ν,x,y,t)
    ν .= @. 1+0*x
    return
end

function callback!(dfn::Diffusion)
    @unpack fld,msh = dfn
    @unpack time, istep = dfn.tstep

    u = fld.u
    println("iter= $(istep[1]), time=$(time[1])")
    plt = meshplt(u,msh)
    display(plt)
    return
end

#----------------------------------#
Ex = 5; nr1 = 8;
Ey = 5; ns1 = 8;

ifperiodic = [false,true]
m1 = Mesh(nr1,ns1,Ex,Ey,ifperiodic,SpectralElements.annulus)
bc = ['D','D','N','N']
diffuseU = Diffusion(bc,m1,Tf=0.0,dt=0.00)

ps = 0

simulate!(diffuseU,callback!,setIC!,setBC!,setForcing!,setVisc!)

#----------------------------------#
#utrue = diffuseU.fld.u
#
#nu = [0.8]
#ps = Flux.params(nu)
#
#function model()
#    dU = Zygote.ignore() do
#        Diffusion(bc,m1,Tf=0.0,dt=0.0)
#    end
#    function varVisc!(ν,x,y,t)
#        ν .= @. nu + 0*x
#        return ν
#    end
#    simulate!(dU,callback!,setIC!,setBC!,setForcing!,varVisc!)
#    return dU.fld.u
#end
#function loss()
#    upred = model()
#    return Flux.Losses.mse(upred,utrue)
#end
#
#opt = Flux.Optimise.ADAM(1e-2)
#gs = Flux.gradient(()->loss(),ps)

#Flux.train!(loss,ps,Iterators.repeated((),100),opt)
#----------------------------------#


text:
datasets() = subtypes(Dataset)

function init_dataset(dataset::Type{<:Dataset})
    register(DataDep(
        dataset_name(dataset),
        dataset_message(dataset),
        dataset_remote_path(dataset),
        dataset_checksum(dataset);
        post_fetch_method=dataset_preprocess(dataset),
    ))
end

function init_dataset(::Type{Dataset})
    for dataset in datasets()
        init_dataset(dataset)
    end
end


text:
@testset "LRNs" begin
    @testset "constructor" begin
        # default arguments
        l = LRN()
        @test l == LRN(1.0, 1.0, 0.5, 5)
        @test l.b == 1.0
        @test l.α == 1.0
        @test l.β == 0.5
        @test l.k == 5
        # custom arguments
        l = LRN(b=2.0, α=1.5, β=0.7, k=3)
        @test l.b == 2.0
        @test l.α == 1.5
        @test l.β == 0.7
        @test l.k == 3
        # custom arguments, standard constructor
        l = LRN(2.0, 1.5, 0.7, 3)
        @test l.b == 2.0
        @test l.α == 1.5
        @test l.β == 0.7
        @test l.k == 3
    end # @testset "constructor"

    @testset "apply" begin
        # the first two dimensions are treated as image dims, the third is the
        # feature dimension and the last is batch.
        #
        # For the first few tests, we keep the third dimension at 1, so we can
        # ignore k.
        # We tests on zeros and ones with different biases and alphas:
        arr0 = zeros(5, 5, 1, 4)
        arr1 = ones(5, 5, 1, 4)
        l = LRN()
        @test l(arr0) == arr0 # multiplicative scaling does not affect zeros
        @test l(arr1) ≈ arr1 ./ sqrt(2.0)
        scaling = rand()
        b = rand()
        l = LRN(b=b)
        @test l(arr0) == arr0
        @test l(arr1) ≈ arr1 ./ sqrt(1.0 + b)
        @test l(scaling.* arr1) ≈ scaling .* arr1 ./ sqrt(scaling^2 + b)
        α = rand()
        l = LRN(α=α)
        @test l(arr0) == arr0
        @test l(arr1) ≈ arr1 ./ sqrt(1.0 + α)
        @test l(scaling.* arr1) ≈ scaling .* arr1 ./ sqrt(1.0 + α * scaling^2)
        # Next, test that the summation over the feature dimension happens
        # correctly:
        arr1 = ones(1, 1, 10, 2)
        l = LRN(k=5)
        summed = arr1 .* 5.0
        summed[:, :, 1, :] .= 3.0
        summed[:, :, 2, :] .= 4.0
        summed[:, :, 9, :] .= 4.0
        summed[:, :, 10, :] .= 3.0
        @test l(arr1) ≈ arr1 ./ sqrt.(1.0 .+ summed)
    end # @testset "apply"
end # @testset "LRNs"


text:
@doc doc"""
Résolution des problèmes de minimisation sous contraintes d'égalités

# Syntaxe
```julia
Lagrangien_Augmente(algo,fonc,contrainte,gradfonc,hessfonc,grad_contrainte,
            hess_contrainte,x0,options)
```

# Entrées
      * **algo**    : (String) l'algorithme sans contraintes à utiliser:
    - **"newton"**  : pour l'algorithme de Newton
    - **"cauchy"**  : pour le pas de Cauchy
    - **"gct"**     : pour le gradient conjugué tronqué
      * **fonc**    : (Function) la fonction à minimiser
  * **contrainte**   : (Function) la contrainte [x est dans le domaine des contraintes ssi ``c(x)=0``]
  * **gradfonc**       : (Function) le gradient de la fonction
  * **hessfonc**    : (Function) la hessienne de la fonction
  * **grad_contrainte** : (Function) le gradient de la contrainte
  * **hess_contrainte** : (Function) la hessienne de la contrainte
      * **x0**    : (Array{Float,1}) la première composante du point de départ du Lagrangien
  * **options**   : (Array{Float,1})
    1. **epsilon**    : utilisé dans les critères d'arrêt
    2. **tol**         : la tolérance utilisée dans les critères d'arrêt
    3. **itermax**    : nombre maximal d'itération dans la boucle principale
    4. **lambda0**   : la deuxième composante du point de départ du Lagrangien
    5. **mu0,tho**   : valeurs initiales des variables de l'algorithme

# Sorties
* **xmin**   : (Array{Float,1}) une approximation de la solution du problème avec contraintes
* **fxmin**    : (Float) ``f(x_{min})``
* **flag**   : (Integer) indicateur du déroulement de l'algorithme
   - **0**    : convergence
   - **1**    : nombre maximal d'itération atteint
   - **(-1)** : une erreur s'est produite
* **niters**    : (Integer) nombre d'itérations réalisées

# Exemple d'appel
```julia
using LinearAlgebra
f(x)=100*(x[2]-x[1]^2)^2+(1-x[1])^2
gradf(x)=[-400*x[1]*(x[2]-x[1]^2)-2*(1-x[1]) ; 200*(x[2]-x[1]^2)]
hessf(x)=[-400*(x[2]-3*x[1]^2)+2  -400*x[1];-400*x[1]  200]
algo = "gct" # ou newton|gct
x0 = [1; 0]
options = []
contrainte(x) =  (x[1]^2) + (x[2]^2) -1.5
grad_contrainte(x) = [2*x[1] ;2*x[2]]
hess_contrainte(x) = [2 0;0 2]
output = Lagrangien_Augmente(algo,f,contrainte,gradf,hessf,grad_contrainte,hess_contrainte,x0,options)
```
"""
function Lagrangien_Augmente(algo, fonc::Function, contrainte::Function, gradfonc::Function,
    hessfonc::Function, grad_contrainte::Function, hess_contrainte::Function, x0, options)

    if options == []
        epsilon = 1e-8
        tol = 1e-15
        itermax = 000
        lambda0 = 2
        mu0 = 100
        tho = 2
    else
        epsilon = options[1]
        tol = options[2]
        itermax = options[3]
        lambda0 = options[4]
        mu0 = options[5]
        tho = options[6]
    end

    xk = x0
    xk_1 = x0
    mu = mu0
    lambda = lambda0
    betha = 0.9
    alpha = 0.1
    etha_c = 0.1258925
    epsilon0 = 1
    etha = etha_c / (mu ^ alpha)
    flag = 0
    i = 0

    for i in 1 : itermax
        Lx(x) = fonc(x) + lambda' * contrainte(x) + (mu / 2) * (norm(contrainte(x))) ^ 2
        grad_x_Lx(x) = gradfonc(x) + lambda' * grad_contrainte(x) + mu * grad_contrainte(x) * contrainte(x)
        Hess_x_Lx(x) = hessfonc(x) + lambda' * hess_contrainte(x) + mu * hess_contrainte(x) * contrainte(x) + mu * grad_contrainte(x) * grad_contrainte(x)'

        while (norm(grad_x_Lx(xk)) > epsilon)
            a = norm(grad_x_Lx(xk))
            if algo == "newton"
                xk, ~ = Algorithme_De_Newton(Lx, grad_x_Lx, Hess_x_Lx, xk, [])
            elseif algo == "cauchy" || algo == "gct"
                xk, ~ = Regions_De_Confiance(algo, Lx, grad_x_Lx, Hess_x_Lx, xk, [])
            end
            if a == norm(grad_x_Lx(xk))
                break
            end
        end

        gradL(x,l) = gradfonc(x) + l' * grad_contrainte(xk)
        if norm(gradL(xk, lambda)) <= max(tol, tol * norm(gradL(x0, lambda0)))
            flag = 0
            break
        elseif norm(contrainte(xk)) <= max(tol, tol * norm(contrainte(x0)))
            flag = 1
            break
        elseif i == itermax
            flag = 3
            break
        else
            if norm(contrainte(xk)) <= etha
                lambda = lambda + mu * contrainte(xk)
                epsilon = epsilon / mu
                etha = etha / (mu ^ betha)
            else
                mu = tho * mu
                epsilon = epsilon0 / mu
                etha = etha_c / (mu ^ alpha)
            end # if
        end # if

    end

    xmin = xk
    fxmin = fonc(xmin)

    return xmin, fxmin, flag, i
end

max_stars_repo_name: MOUDDENEHamza/ENSEEIHT

max_issues_repo_name: MOUDDENEHamza/ENSEEIHT

max_forks_repo_name: MOUDDENEHamza/ENSEEIHT


text:
"""
Atomistic solution for an PDAE.

### Fields

* `t`: time of current time step
* `t̅`: time of previous time step
* `q`: current solution of q
* `q̅`: previous solution of q
* `q̃`: compensated summation error of q
* `p`: current solution of p
* `p̅`: previous solution of p
* `p̃`: compensated summation error of p
* `λ`: current solution of λ
* `λ̅`: previous solution of λ
* `v`: vector field of q
* `v̅`: vector field of q̅
* `f`: vector field of p
* `f̅`: vector field of p̅
* `u`: projective vector field of q
* `u̅`: projective vector field of q̅
* `g`: projective vector field of p
* `g̅`: projective vector field of p̅
"""
mutable struct AtomicSolutionPDAE{DT,TT} <: AtomicSolution{DT,TT}
    t::TT
    t̅::TT

    q::Vector{DT}
    q̅::Vector{DT}
    q̃::Vector{DT}

    p::Vector{DT}
    p̅::Vector{DT}
    p̃::Vector{DT}

    λ::Vector{DT}
    λ̅::Vector{DT}

    v::Vector{DT}
    v̅::Vector{DT}
    f::Vector{DT}
    f̅::Vector{DT}

    u::Vector{DT}
    u̅::Vector{DT}
    g::Vector{DT}
    g̅::Vector{DT}

    function AtomicSolutionPDAE{DT, TT}(nd) where {DT <: Number, TT <: Real}
        new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd),
                                zeros(DT, nd), zeros(DT, nd), zeros(DT, nd),
                                zeros(DT, nd), zeros(DT, nd),
                                zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd),
                                zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd))
    end
end

AtomicSolutionPDAE(DT, TT, nd) = AtomicSolutionPDAE{DT, TT}(nd)

function CommonFunctions.set_solution!(asol::AtomicSolutionPDAE, sol)
    t, q, p, λ = sol
    asol.t  = t
    asol.q .= q
    asol.p .= p
    asol.λ .= λ
    asol.v .= 0
    asol.f .= 0
end

function CommonFunctions.get_solution(asol::AtomicSolutionPDAE)
    (asol.t, asol.q, asol.p, asol.λ)
end

function CommonFunctions.reset!(asol::AtomicSolutionPDAE, Δt)
    asol.t̅  = asol.t
    asol.q̅ .= asol.q
    asol.p̅ .= asol.p
    asol.λ̅ .= asol.λ
    asol.v̅ .= asol.v
    asol.f̅ .= asol.f
    asol.u̅ .= asol.u
    asol.g̅ .= asol.g
    asol.t += Δt
end

function update!(asol::AtomicSolutionPODE{DT}, v::Vector{DT}, f::Vector{DT}, λ::Vector{DT}) where {DT}
    for k in eachindex(v,f)
        update!(asol, v[k], f[k], λ[k])
    end
end

function update!(asol::AtomicSolutionPODE{DT}, v::DT, f::DT, λ::DT, k::Int) where {DT}
    asol.q[k], asol.q̃[k] = compensated_summation(v, asol.q[k], asol.q̃[k])
    asol.p[k], asol.p̃[k] = compensated_summation(f, asol.p[k], asol.p̃[k])
    asol.λ[k] = λ
end


text:
export MtlEvent, MtlSharedEvent

abstract type MtlAbstractEvent end

const MTLEvent = Ptr{MtEvent}
const MTLSharedEvent = Ptr{MtSharedEvent}
const MTLSharedEventHandle = Ptr{MtSharedEventHandle}

mutable struct MtlEvent <: MtlAbstractEvent
	handle::MTLEvent
	device::MtlDevice
end

mutable struct MtlSharedEvent <: MtlAbstractEvent
	handle::MTLEvent
	device::MtlDevice
end

Base.convert(::Type{MTLEvent}, ev::MtlAbstractEvent) = ev.handle
Base.unsafe_convert(::Type{MTLEvent}, ev::MtlAbstractEvent) = convert(MTLEvent, ev.handle)

Base.:(==)(a::MtlAbstractEvent, b::MtlAbstractEvent) = a.handle == b.handle
Base.hash(ev::MtlAbstractEvent, h::UInt) = hash(ev.handle, h)

function unsafe_destroy!(fun::MtlAbstractEvent)
	fun.handle !== C_NULL && mtEventRelease(fun)
end

function MtlEvent(dev::MtlDevice)
	handle = mtDeviceNewEvent(dev)
	obj = MtlEvent(handle, dev)
	finalizer(unsafe_destroy!, obj)
	return obj
end

function MtlSharedEvent(dev::MtlDevice)
	handle = mtDeviceNewSharedEvent(dev)
	obj = MtlSharedEvent(handle, dev)
	finalizer(unsafe_destroy!, obj)
	return obj
end

device(ev::MtlAbstractEvent) = ev.device
function label(l::MtlAbstractEvent)
	ptr = mtEventLabel(l)
	return ptr == C_NULL ? "" : unsafe_string(ptr)
end

# shared event
value(ev::MtlSharedEvent) = mtSharedEventSignaledValue(ev)

mutable struct MtlSharedEventHandle
	handle::MTLSharedEventHandle
	event::MtlSharedEvent
end
function MtlSharedEventHandle(event::MtlSharedEvent)
	handle = mtSharedEventNewHandle(event)
	obj = MtlSharedEventHandle(handle, event)
	finalizer(unsafe_destroy!, obj)
	return obj
end
function unsafe_destroy!(fun::MtlSharedEventHandle)
	fun.handle !== C_NULL && mtSharedEventHandleRelease(fun)
end
Base.convert(::Type{MTLSharedEventHandle}, lib::MtlSharedEventHandle) = lib.handle
Base.unsafe_convert(::Type{MTLSharedEventHandle}, lib::MtlSharedEventHandle) = convert(MTLSharedEventHandle, lib.handle)

Base.:(==)(a::MtlSharedEventHandle, b::MtlSharedEventHandle) = a.handle == b.handle
Base.hash(lib::MtlSharedEventHandle, h::UInt) = hash(lib.handle, h)

## FENCES
const MTLFence = Ptr{MtFence}

mutable struct MtlFence
	handle::MTLFence
	device::MtlDevice
end

Base.convert(::Type{MTLFence}, fen::MtlFence) = fen.handle
Base.unsafe_convert(::Type{MTLFence}, fen::MtlFence) = convert(MTLFence, fen.handle)

Base.:(==)(a::MtlFence, b::MtlFence) = a.handle == b.handle
Base.hash(ev::MtlFence, h::UInt) = hash(ev.handle, h)

function MtlFence(dev::MtlDevice)
	handle = mtDeviceNewFence(dev)
	obj = MtlFence(handle, dev)
	finalizer(unsafe_destroy!, obj)
	return obj
end

device(fen::MtlFence) = fen.device


text:
############################################################################################
# Constants
"RNG for sampling based solutions"
const _rng = Random.MersenneTwister(1)   # shorthand
Random.seed!(_rng, 1)

"Tolerance for stochastic solutions"
const _TOL = 1.0e-6

"Number of samples"
N = 10^3
N_SMC2 = 10^2

############################################################################################
# Kernel and AD backends
include("mcmc.jl")
include("smc.jl")


text:
# Fig_GS_Roessler_Lorenz_v2.jl
# =============================

using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, OrdinaryDiffEq
using StaticArrays
using Random
using Statistics
using Neighborhood
using LinearAlgebra

function ccm(x, y, d, τ)
    ψ = copy(y); Mx = embed(x, d, τ); tw = Theiler(τ);
    # bulk search of d+1 neighbors for all i points:
    idxs=bulkisearch(KDTree(Mx), Mx, NeighborNumber(d+1), tw)
    for i in 1:length(Mx)
      nni = idxs[i]; xi = Mx[i]; n1 = norm(xi - Mx[nni[1]])
      # implement equation (7.5):
      u = [exp(-norm(xi - Mx[j])/n1) for j ∈ nni]
      w = u ./ sum(u)
      ψ[i] = sum(w[k]*y[j] for (k, j) in enumerate(nni))
    end
    return cor(y, ψ) # Pearson correlation coef.
  end

function Roessler_Lorenz_aux_ODE(du,u,p,t)
    a, b, c, alpha,  cc = p

    du[1] = -u[2] - u[3]
    du[2] = u[1] + a*u[2]
    du[3] = b + u[3]*(u[1]-c)

    du[1:3] = alpha .* du[1:3]   # adjust time scales

    du[4] = 10. * (-u[4] + u[5])
    du[5] = 28*u[4] - u[5] - u[4]*u[6]  + cc*u[2]
    du[6] = u[4]*u[5] - 2.666666 * u[6]

    du[7] = 10. * (-u[7] + u[8])
    du[8] = 28*u[7] - u[8] - u[7]*u[9]  + cc*u[2]
    du[9] = u[7]*u[8] - 2.666666 * u[9]

end

function Roessler_Lorenz_CLE_ODE(du,u,p,t)
    a, b, c, alpha,  cc = p

    du[1] = -u[2] - u[3]
    du[2] = u[1] + a*u[2]
    du[3] = b + u[3]*(u[1]-c)

    du[1:3] = alpha .* du[1:3]   # adjust time scales

    du[4] = 10. * (-u[4] + u[5])
    du[5] = 28*u[4] - u[5] - u[4]*u[6] + cc*u[2]
    du[6] = u[4]*u[5] - 2.666666 * u[6]

    # linearized equation
    du[7] =       -10 * u[7]  +   10 * u[8]
    du[8] = (28-u[6]) * u[7]         - u[8]     -  u[4] * u[9]
    du[9] =      u[5] * u[7]  + u[4] * u[8]  - 2.666666 * u[9]
    return nothing
end

#  ------------------ main program --------------------------

@time begin

# parameters and initial values
# -----------------------------
a = 0.2
b = 0.2
c = 5.7
alpha = 6

ccmin = 0.
ccmax = 12.
ncc = 100

ccvec = LinRange(ccmin,ccmax,ncc)
GSvec = zeros(ncc)
cxyvec = zeros(ncc)
cyxvec = zeros(ncc)

dxyvec = zeros(ncc)
dyxvec = zeros(ncc)

CLEvec = zeros(ncc)

ttrans = 2000  # transient time
tattr = 50000   # time for averaging on attractor

d = 6
τ = 1  #  3
tsamp = 0.1 # 0.1

# coupling parameter loop
# -----------------------
for icc = 1:ncc
    @show icc
    cc = ccvec[icc]
    p = a, b, c, alpha,  cc
    uinit = [-1.589656  3.0478997  0.823033  0.11  0.19  0.1 -0.1 -0.2 0.1]

   # transient to attractor
   # ----------------------
     tspan = (0., ttrans)
     prob = ODEProblem(Roessler_Lorenz_aux_ODE, uinit, tspan, p)
     sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, maxiters = 1e7)

   # attractor
   # ---------

    tspan = (0., tattr)
    y_init = sol.u[end]
    prob = ODEProblem(Roessler_Lorenz_aux_ODE, y_init, tspan, p)
    sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, maxiters = 1e7, saveat = tsamp)

    tvec = sol.t
    ures  = sol.u
    nts = length(tvec)

    xpts = zeros(3,nts)
    ypts = zeros(3,nts)
    zpts = zeros(3,nts)

    for its = 1:nts     # .... this can be done more effciently
        uu = ures[its]
        xpts[1:3,its] = uu[1:3]
        ypts[1:3,its] = uu[4:6]
        zpts[1:3,its] = uu[7:9]
    end

   # Synchronization error (auxiliary system method)
   # -----------------------------------------------
    GSerror = 0.
    for its = 1:nts # .... can this be done without a loop?
        GSerror = GSerror +  norm(ypts[:,its] - zpts[:,its] )
    end
    GSerror = GSerror / nts
    GSvec[icc] = GSerror

    # convergent cross mapping
    # ------------------------
    xx = xpts[1,:]
    yy = ypts[1,:]

    cxyvec[icc] = ccm(xx, yy, d, τ)
    cyxvec[icc] = ccm(yy, xx, d, τ)

    if icc == 1
    fig = plt.figure("Fig Roessler Lorenz d=3 τ=1",figsize=(figx,figy))

    subplot(1,2,1)
    plt.scatter(xx[1:end-1],xx[2:end],s=1)

    subplot(1,2,2)
    plt.scatter(yy[1:end-1],yy[2:end],s=1)
    end

   # CLEs
   # ----
   tsamp = 1. # 0.1
   tspan = (0., tattr)
   y_init = sol.u[end]
   y_init[7:9] = [1 0 0]
   ncle = 50000 # Lyapunov steps
   csum = 0.

   for icle = 1:ncle
       tspan = (0., tsamp)
       prob = ODEProblem(Roessler_Lorenz_CLE_ODE, y_init, tspan, p)
       sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, maxiters = 1e7, saveat = tsamp)

       tvec = sol.t
       y_init = sol.u[end]

       ynorm = sqrt( norm( y_init[7:9] ) )
       y_init[7:9] = y_init[7:9] ./ ynorm   # normalization of tangent vector

       if icle > 300
          csum = csum + log(ynorm)
       end

   nts = length(tvec)

   end

   CLE = csum / (ncle-300)
   CLEvec[icc] = CLE

# println("c=",cc,"  GS=",GSerror, "  ccm_xy=", cxyvec[icc],"  ccm_yx=",cyxvec[icc],"  dxy=",dxyvec[icc],"  dyx=",dyxvec[icc],"  CLE=", CLE)
println("c=",cc,"  GS=",GSerror, "  ccm_xy=", cxyvec[icc],"  ccm_yx=",cyxvec[icc],"  CLE=", CLE)

end

end  # of time

# %% plot

fig = plt.figure(figsize=(0.5*figx,figy))

plt.plot(ccvec,GSvec ./maximum(GSvec),color=COLORS[1])  # ,linewidth=1)
plt.plot(ccvec,CLEvec ./ maximum(CLEvec),color=COLORS[2], ls = "--")

plt.plot(ccvec,cxyvec,color=COLORS[3])
plt.plot(ccvec,cyxvec,color=COLORS[4], ls = "-.")

plt.xlabel(L"k"; labelpad = -20)

ax = gca()
ax.set_xlim([ccmin,ccmax])
ax.set_xticks(1:2:11)
# ax.set_ylim([-2.1,2.6])
ax.legend([L"E"; L"\lambda_1^{\rm{CLE}}";
           "CCM: \$x_1 \\to y_1\$"; "CCM: \$y_1 \\to x_1\$"],
         bbox_to_anchor=(0., 1.02, 1., .102), loc="lower left",
        ncol=2, mode="expand", borderaxespad=0, handlelength=2,
)

fig.tight_layout(pad=0.3)
wsave(plotsdir("9", "generalized_sync"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using DiffEqProblemLibrary, Sundials, Test
using DiffEqProblemLibrary.DAEProblemLibrary: importdaeproblems;
importdaeproblems();
using DiffEqProblemLibrary.DAEProblemLibrary: prob_dae_resrob

# Test DAE
mutable struct precflags
    prec_used::Bool
    psetup_used::Bool
end

p = precflags(false, false)
prob = DAEProblem(prob_dae_resrob.f, prob_dae_resrob.du0, prob_dae_resrob.u0, prob_dae_resrob.tspan, p)

dt = 1000
saveat = float(collect(0:dt:100000))
sol = solve(prob, IDA())
@info "Multiple abstol"
sol = solve(prob, IDA(), abstol = [1e-9, 1e-8, 1e-7])
@info "Band solver"
sol2 = solve(prob, IDA(linear_solver = :Band, jac_upper = 2, jac_lower = 2))
@info "GMRES solver"
sol3 = solve(prob, IDA(linear_solver = :GMRES))
#sol4 = solve(prob,IDA(linear_solver=:BCG)) # Fails but doesn't throw an error?
@info "TFQMR solver"
sol5 = solve(prob, IDA(linear_solver = :TFQMR))
@info "FGMRES solver"
sol6 = solve(prob, IDA(linear_solver = :FGMRES))
@info "PCG solver"
sol7 = solve(prob, IDA(linear_solver = :PCG)) # Requires symmetric linear
#@info "KLU solver"
#sol8 = solve(prob,IDA(linear_solver=:KLU)) # Requires Jacobian

sol9 = solve(prob, IDA(linear_solver = :LapackBand, jac_upper = 2, jac_lower = 2))
sol10 = solve(prob, IDA(linear_solver = :LapackDense))
sol11 = solve(prob, IDA(linear_solver = :Dense))

# Test identity preconditioner
prec = (z, r, p, t, y, fy, resid, gamma, delta) -> (p.prec_used = true; z .= r)
psetup = (p, t, resid, u, du, gamma) -> (p.psetup_used = true)
@info "GMRES for identity preconditioner"
sol4 = solve(prob, IDA(linear_solver = :GMRES, prec = prec))
@test p.prec_used
@info "GMRES with pset"
sol4 = solve(prob, IDA(linear_solver = :GMRES, prec = prec, psetup = psetup))
@test p.psetup_used

@info "IDA with saveat"
sol = solve(prob, IDA(), saveat = saveat)

@test sol.t == saveat
@info "IDA with saveat everystep"
sol = solve(prob, IDA(), saveat = saveat, save_everystep = true)

@test sol.t != saveat
@test intersect(sol.t, saveat) == saveat
@info "IDA with tstops"
sol = solve(prob, IDA(), tstops = [0.9])

@test 0.9 ∈ sol.t

prob = deepcopy(prob_dae_resrob)
prob2 = DAEProblem(prob.f, prob.du0, prob.u0, (1.0, 0.0))
sol = solve(prob2, IDA())
@test maximum(diff(sol.t)) < 0 # Make sure all go negative

function f!(res, du, u, p, t)
    res[1] = du[1] - 1.01
    return
end

u0 = [0.0]
du0 = [1.01]
tspan = (0.0, 10.0)

println("With consistent initial conditions:")

dae_prob = DAEProblem(f!, du0, u0, tspan, differential_vars = [true])
sol = solve(dae_prob, IDA())

println("With inconsistent initial conditions:")

du0 = [0.0]

dae_prob = DAEProblem(f!, du0, u0, tspan, differential_vars = [true])
sol = solve(dae_prob, IDA())


text:
module TestInteractive
using Debug

println(Debug.UI.helptext)
println()
println("Type an expression to evaluate it in the current scope.")

@debug begin
    function sumpow(n::Int, p::Int = 1)
        @bp
        if n <= 0; return 0; end
        n^p + sumpow(n-1, p)
    end
    sumpow(10)
end

end # module