text:
## 2D Panels
#==========================================================================================#

abstract type AbstractPanel2D <: AbstractPanel end

struct Panel2D{T <: Real} <: AbstractPanel2D
    p1 :: SVector{2,T}
    p2 :: SVector{2,T}
end

struct WakePanel2D{T <: Real} <: AbstractPanel2D
    p1 :: SVector{2,T}
    p2 :: SVector{2,T}
end

Base.length(:: Panel2D) = 1

Panel2D(p1 :: FieldVector{2,T}, p2 :: FieldVector{2,T}) where T <: Real = Panel2D{T}(p1, p2)
WakePanel2D(p1 :: FieldVector{2,T}, p2 :: FieldVector{2,T}) where T <: Real = WakePanel2D{T}(p1, p2)

a :: AbstractPanel2D + b :: AbstractPanel2D = Panel2D(p1(a) + p1(b), p2(a) + p2(b))
a :: AbstractPanel2D - b :: AbstractPanel2D = Panel2D(p1(a) - p1(b), p2(a) - p2(b))

collocation_point(panel :: AbstractPanel2D, a = 0.5) = a * (panel.p1 + panel.p2)
panel_vector(panel :: AbstractPanel2D) = panel.p2 - panel.p1
panel_length(panel :: AbstractPanel2D) = (norm ∘ panel_vector)(panel)

function transform_panel_points(panel_1 :: AbstractPanel2D, panel_2 :: AbstractPanel2D)
    x1, y1, x2, y2 = panel_2.p1, panel_2.p2
    xs, ys = panel_1.p1

    xp1, yp1 = affine_2D(x1, y1, xs, ys, panel_angle(panel_1))
    xp2, yp2 = affine_2D(x2, y2, xs, ys, panel_angle(panel_1))

    xp1, yp1, xp2, yp2
end

function transform_panel(panel :: AbstractPanel2D, x, y)
    xs, ys = panel.p1
    affine_2D(x, y, xs, ys, panel_angle(panel))
end

transform_panel(panel :: AbstractPanel2D, point :: SVector{2,<: Real}) = transform_panel(panel, point.x, point.y)

function panel_angle(panel :: AbstractPanel2D)
    (x, y) = panel_vector(panel)
    θ = atan(y, x)
    # ifelse(θ < 0., θ + 2π, θ) # Branch cut
end

panel_normal(panel :: AbstractPanel2D) = let (x, y) = panel_vector(panel); normalize([-y; x]) end

tangent_vector(panel :: AbstractPanel2D) = rotation(1., 0., -panel_angle(panel))

normal_vector(panel :: AbstractPanel2D) = inverse_rotation(0., 1., panel_angle(panel))

panel_location(panel :: AbstractPanel2D) = let angle = panel_angle(panel); ifelse((π/2 <= angle <= π) || (-π <= angle <= -π/2), "lower", "upper") end

panel_points(panels) = [ p1.(panels); [(p2 ∘ last)(panels)] ]

reverse_panel(panel :: AbstractPanel2D) = Panel2D(panel.p2, panel.p1)

trailing_edge_panel(panels) = Panel2D((p2 ∘ last)(panels), (p1 ∘ first)(panels))

bisector(f_pan, l_pan) = (normalize(panel_vector(l_pan)) + normalize(-panel_vector(f_pan))) / 2

function trailing_edge_info(panels)
    # Make trailing edge
    te_panel  = trailing_edge_panel(panels)
    s    = panel_vector(reverse_panel(te_panel))
    p    = normalize(s)

    t1 = -normalize(panel_vector(panels[1]))
    t2 = normalize(panel_vector(panels[end]))
    t  = normalize((t1 + t2) / 2)

    h_TE = -s[1] * t[2] + s[2] * t[1]
    tcp  = abs(t[1] * p[2] - t[2] * p[1])
    tdp  = dot(p, t)
    dtdx = t1[1] * t2[2] - t2[1] * t1[2]

    te_panel, h_TE, tcp, tdp, dtdx
end

function wake_panel(panels, bound, α)
    _, firsty        = (p1 ∘ first)(panels)
    lastx, lasty     = (p2 ∘ last)(panels)
    y_mid            = (firsty + lasty) / 2
    y_bound, x_bound = bound .* sincos(α)
    WakePanel2D(SVector(lastx, y_mid), SVector(x_bound * lastx, y_bound * y_mid))
end

function wake_panels(panels, chord, length, num)
    _, firsty = (p1 ∘ first)(panels)
    _, lasty  = (p2 ∘ last)(panels)
    y_mid     = (firsty + lasty) / 2
    bounds    = sine_spacing(chord + length / 2, length, num)
    @. WakePanel2D(SVector(bounds[1:end-1], y_mid), SVector(bounds[2:end], y_mid))
end

function panel_scalar(scalar_func, strength, panel :: AbstractPanel2D, x, y)
    # Transform point to local panel coordinates
    xp, yp = transform_panel(panel, SVector(x, y))
    scalar_func(strength, xp, yp, 0., panel_length(panel))
end

function panel_scalar(scalar_func, strength, panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D, point = 0.5)
    x, y = collocation_point(panel_i, point)
    panel_scalar(scalar_func, strength, panel_j, x, y)
end

function panel_velocity(velocity_func, strength, panel :: AbstractPanel2D, x, y)
    # Transform point to local panel coordinates
    xp, yp = transform_panel(panel, SVector(x, y))

    # Compute velocity in local panel coordinates
    u, w = velocity_func(strength, xp, yp, 0., panel_length(panel))

    # Transform velocity to original coordinates
    inverse_rotation(u, w, panel_angle(panel))
end

function panel_velocity(velocity_func, strength, panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D, point = 0.5)
    x, y = collocation_point(panel_i, point)
    u, w = panel_velocity(velocity_func, strength, panel_j, x, y)
end

panel_velocity(f1, f2, str1, str2, panel :: AbstractPanel2D, x, y) = panel_velocity(f1, str1, panel, x, y) + panel_velocity(f2, str2, panel, x, y)

panel_velocity(f1, f2, str_j1, str_j2, panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D) = panel_velocity(f1, str_j1, panel_j, panel_i) .+ panel_velocity(f2, str_j2, panel_j, panel_i)

get_surface_values(panels, vals, surf = "upper", points = false) = partition(x -> (panel_location ∘ first)(x) == surf, (collect ∘ zip)(panels[2:end], vals), x -> ((first ∘ ifelse(points, p1, collocation_point) ∘ first)(x), last(x)))


text:
### FlatMap and FlatFourier types
struct FlatMap{P<:Flat,T<:Real,M<:AbstractMatrix{T}} <: Field{Map,S0,P,T}
    Ix :: M
end
struct FlatFourier{P<:Flat,T<:Real,M<:AbstractMatrix{Complex{T}}} <: Field{Fourier,S0,P,Complex{T}}
    Il :: M
end
const FlatS0{P,T,M} = Union{FlatMap{P,T,M},FlatFourier{P,T,M}}

### convenience constructors
for (F, X, T) in [
        (:FlatMap,     :Ix, :T),
        (:FlatFourier, :Il, :(Complex{T})),
    ]
    doc = """
        # main constructor:
        $F($X::AbstractMatrix[, θpix={resolution in arcmin}, ∂mode={fourier∂ or map∂})

        # more low-level:
        $F{P}($X::AbstractMatrix) # specify pixelization P explicilty
        $F{P,T}($X::AbstractMatrix) # additionally, convert elements to type $T
        $F{P,T,M<:AbstractMatrix{$T}}($X::M) # specify everything explicilty

    Construct a `$F` object. The top form of the constructor is most convenient
    for interactive work, while the others may be more useful for low-level code.
    """
    @eval begin
        @doc $doc $F
        $F($X; kwargs...) = $F{Flat(Nside=size($X,2);kwargs...)}($X)
        $F{P}($X::M) where {P,T,M<:AbstractMatrix{$T}} = $F{P,T,M}($X)
        $F{P,T}($X::AbstractMatrix) where {P,T} = $F{P}($T.($X))
    end
    T!=:T && @eval $F{P}($X::M) where {P,T,M<:AbstractMatrix{T}} = $F{P,T}($X)
end

### array interface
size(f::FlatS0) = (length(firstfield(f)),)
lastindex(f::FlatS0, i::Int) = lastindex(f.Ix, i)
size_2d(::Type{<:FlatMap{<:Flat{N}}}) where {N} = (N,N)
size_2d(::Type{<:FlatFourier{<:Flat{N}}}) where {N} = (N÷2+1,N)
@propagate_inbounds @inline getindex(f::FlatS0, I...) = getindex(firstfield(f), I...)
@propagate_inbounds @inline setindex!(f::FlatS0, X, I...) = (setindex!(firstfield(f), X, I...); f)
adapt_structure(to, f::F) where {P,F<:FlatS0{P}} = basetype(F){P}(adapt(to,firstfield(f)))
function similar(f::F,::Type{T},dims::Dims) where {P,F<:FlatS0{P},T<:Number}
    @assert size(f)==dims "Tried to make a field similar to $F but dims should have been $(size(f)), not $dims."
    basetype(F){P}(similar(firstfield(f),T))
end

### broadcasting
struct FlatS0Style{F,M} <: AbstractArrayStyle{1} end
(::Type{FS})(::Val{1}) where {FS<:FlatS0Style} = FS()
@generated BroadcastStyle(::Type{F}) where {P,T,M,F<:FlatS0{P,T,M}} = FlatS0Style{basetype(F){P},basetype(M)}()
BroadcastStyle(::FlatS0Style{<:Field{B1}}, ::FlatS0Style{<:Field{B2}}) where {B1,B2} = invalid_broadcast_error(B1,B2)
BroadcastStyle(S::FieldTupleStyle, ::FlatS0Style) = S
BroadcastStyle(S::FieldOrOpArrayStyle, ::FlatS0Style) = S
similar(::Broadcasted{FS}, ::Type{T}) where {T<:Number,FS<:FlatS0Style} = similar(FS,T)
similar(::Type{FlatS0Style{F,M}}, ::Type{T}) where {F<:FlatS0,M,T<:Number} = F(basetype(M){eltype(F{real(T)})}(undef,size_2d(F)...))
@inline preprocess(dest::F, bc::Broadcasted) where {F<:FlatS0} = Broadcasted{DefaultArrayStyle{2}}(bc.f, preprocess_args(dest, bc.args), map(OneTo,size_2d(F)))
preprocess(dest::F, arg) where {F<:FlatS0} = broadcastable(F, arg)
broadcastable(::Type{F}, f::FlatS0{P}) where {P,F<:FlatS0{P}} = firstfield(f)
broadcastable(::Type{F}, f::AbstractVector) where {P,F<:FlatS0{P}} = reshape(f, size_2d(F))
broadcastable(::Any, x) = x

### basis conversion
Fourier(f::FlatMap{P}) where {P} = FlatFourier{P}(fieldinfo(f).FFT * f.Ix)
Map(f::FlatFourier{P}) where {P} = FlatMap{P}(fieldinfo(f).FFT \ f.Il)

### inplace conversion
Fourier(f′::FlatFourier, f::FlatMap) =  (mul!(f′.Il, fieldinfo(f).FFT, f.Ix); f′)
Map(f′::FlatMap, f::FlatFourier)     = (ldiv!(f′.Ix, fieldinfo(f).FFT, f.Il); f′)

### properties
getproperty(f::FlatS0, s::Symbol) = getproperty(f,Val(s))
getproperty(f::FlatS0, ::Val{s}) where {s} = getfield(f,s)
getproperty(f::FlatS0, ::Val{:I}) = f
function getindex(f::FlatS0, k::Symbol)
    k in [:I, :Ix,:Il] || throw(ArgumentError("Invalid FlatS0 index: $k"))
    k == :I ? f : getproperty((k == :Ix ? Map : Fourier)(f),k)
end

### dot products
# do in Map space for simplicity, and use sum_kbn to reduce roundoff error
dot(a::FlatS0{P}, b::FlatS0{P}) where {P} = sum_kbn(Map(a).Ix .* Map(b).Ix) * fieldinfo(a).Δx^2

### simulation and power spectra
function white_noise(::Type{F}) where {N,P<:Flat{N},T,M,F<:FlatS0{P,T,M}}
    FlatMap{P}(randn!(basetype(M){T}(undef,N,N)) / fieldinfo(F).Δx)
end
function Cℓ_to_Cov(::Type{P}, ::Type{T}, ::Type{S0}, Cℓ::InterpolatedCℓs) where {P,T}
    Diagonal(FlatFourier{P}(Cℓ_to_2D(P,T,Cℓ)))
end

function cov_to_Cℓ(L::DiagOp{<:FlatS0})
    ii = sortperm(fieldinfo(L.diag).kmag[:])
    InterpolatedCℓs(fieldinfo(L.diag).kmag[ii], real.(unfold(L.diag.Il))[ii], concrete=false)
end

function get_Cℓ(f::FlatS0{P}, f2::FlatS0{P}=f; Δℓ=50, ℓedges=0:Δℓ:16000, Cℓfid=ℓ->1, err_estimate=false) where {P}
    @unpack Nside,Δx,kmag = fieldinfo(f)
    α = (Nside/Δx)^2

    L = kmag[:]
    CLobs = real.(dot.(unfold(f[:Il]),unfold(f2[:Il])))[:] ./ α
    w = @. nan2zero((2*Cℓfid(L)^2/(2L+1))^-1)

    sum_in_ℓbins(x) = fit(Histogram, L, Weights(x), ℓedges).weights

    A  = sum_in_ℓbins(w)
    Cℓ = sum_in_ℓbins(w .* CLobs) ./ A
    ℓ  = sum_in_ℓbins(w .* L)     ./ A

    if err_estimate
        N   = sum_in_ℓbins(one.(w)) / 2
        σℓ  = sqrt.((sum_in_ℓbins(w .* CLobs.^2) ./ A .- Cℓ.^2) ./ N)
        InterpolatedCℓs(ℓ,  Cℓ .± σℓ)
    else
        InterpolatedCℓs(ℓ,  Cℓ)
    end
end

"""
    ud_grade(f::Field, θnew, mode=:map, deconv_pixwin=true, anti_aliasing=true)

Up- or down-grades field `f` to new resolution `θnew` (only in integer steps).
Two modes are available specified by the `mode` argument:

* `:map`     — Up/downgrade by replicating/averaging pixels in map-space
* `:fourier` — Up/downgrade by extending/truncating the Fourier grid

For `:map` mode, two additional options are possible. If `deconv_pixwin` is
true, deconvolves the pixel window function from the downgraded map so the
spectrum of the new and old maps are the same. If `anti_aliasing` is true,
filters out frequencies above Nyquist prior to down-sampling.

"""
function ud_grade(f::FlatS0{P,T,M}, θnew; mode=:map, deconv_pixwin=(mode==:map), anti_aliasing=(mode==:map)) where {T,M,θ,N,∂mode,P<:Flat{N,θ,∂mode}}
    θnew==θ && return f
    (mode in [:map,:fourier]) || throw(ArgumentError("Available modes: [:map,:fourier]"))

    fac = θnew > θ ? θnew÷θ : θ÷θnew
    (round(Int, fac) ≈ fac) || throw(ArgumentError("Can only ud_grade in integer steps"))
    fac = round(Int, fac)
    Nnew = round(Int, N * θ ÷ θnew)
    Pnew = Flat{Nnew,θnew,∂mode}

    if deconv_pixwin
        @unpack Δx,k = fieldinfo(Pnew,T,M)
        Wk =  @. T(pixwin(θnew, k) / pixwin(θ, k))
    end

    if θnew>θ
        # downgrade
        if anti_aliasing
            kmask = ifelse.(abs.(fieldinfo(P,T,M).k) .> fieldinfo(Pnew,T,M).nyq, 0, 1)
            AA = Diagonal(FlatFourier{P}(kmask[1:N÷2+1] .* kmask'))
        else
            AA = 1
        end
        if mode==:map
            fnew = FlatMap{Pnew}(mapslices(mean,reshape((AA*f)[:Ix],(fac,Nnew,fac,Nnew)),dims=(1,3))[1,:,1,:])
            deconv_pixwin ? FlatFourier{Pnew}(fnew[:Il] ./ Wk' ./ Wk[1:Nnew÷2+1]) : fnew
        else
            FlatFourier{Pnew}((AA*f)[:Il][1:(Nnew÷2+1), [1:(isodd(Nnew) ? Nnew÷2+1 : Nnew÷2); (end-Nnew÷2+1):end]])
        end
    else
        # upgrade
        if mode==:map
            fnew = FlatMap{Pnew}(permutedims(hvcat(N,(x->fill(x,(fac,fac))).(f[:Ix])...)))
            deconv_pixwin ? FlatFourier{Pnew}(fnew[:Il] .* Wk' .* Wk[1:Nnew÷2+1]) : fnew
        else
            fnew = FlatFourier{P}(zeros(Nnew÷2+1,Nnew))
            broadcast_setindex!(fnew.Il, f[:Il], 1:(N÷2+1), [findfirst(fieldinfo(fnew).k .≈ fieldinfo(f).k[i]) for i=1:N]');
            fnew
        end
    end
end

ext: jl


text:
"""Computes a one-dimensional Morlet wavelet in the Fourier domain.
A Morlet wavelet of center frequency `ξ` and of variance `σ` looks almost like
a Gaussian bell curve. To ensure that the wavelet has a vanishing moment, we
substract a corrective term around the zeroth frequency."""
function AbstractFilter{T<:FFTW.fftwReal,G<:LineGroups}(ψmeta::ΨMeta1D,
        spec::Spec1D{T,FourierDomain{1},G,Morlet})
    """1. **Gaussian denominator `den = 2σ²`**
    The FWHM (full width at half maximum) bw of a Gaussian bell curve of
    variance `σ` is defined by the equation
        `g(±bw/2) = exp(- (bw/2)²/(2σ²)) = 1/sqrt(2)`
    which leads to
        `bw² = 2 log(2) * 2σ²`.
    The denominator `den = 2σ²` of the Gaussian is thus equal to
        `den = 2σ² = bw² / (2 log(2))`."""
    @inbounds log2_length = spec.log2_size
    N = 1 << log2_length
    halfN = N >> 1
    center = N * T(ψmeta.centerfrequency)
    bw = N * T(ψmeta.bandwidth)
    den = @fastmath bw * bw / T(2.0 * log(2.0))
    """2. **Number of periods**"""
    halfsupport = sqrt(den * log(inv(spec.ɛ)))
    firstω = max(center - halfsupport, -5N/2)
    lastω = min(center + halfsupport, 5N/2 - 1)
    nPeriods = 1 + ceil(Int, (lastω-halfN) / N)
    """3. **Call to morlet**"""
    y = morlet(FourierDomain(1), center, den, N, nPeriods)
    """4. **Trimming to true support boundaries**"""
    return AbstractFilter(y, spec)
end

function AbstractFilter{T<:FFTW.fftwReal,G<:LineGroups}(ϕmeta::ΦMeta,
        spec::Spec1D{T,FourierDomain{1},G,Morlet})
    """1. **Gaussian denominator `den = 2σ²`**"""
    @inbounds log2_length = spec.log2_size
    N = 1 << log2_length
    halfN = N >> 1
    bw = N * T(ϕmeta.bandwidth)
    den = @fastmath bw * bw / T(2.0 * log(2.0))
    """2. **Call to `morlet`**"""
    lastω = round(Int, min(bw * sqrt(2.0 / spec.ɛ), halfN - 1))
    leg = T[ gauss(ω, den) for ω in 1:lastω ]
    """3. **Trimming to true support boundaries**"""
    leg = leg[1:findlast(leg .> spec.ɛ)]
    return FourierSymmetric1DFilter(leg, one(T))
end

"""Computes `gauss{T<:Real}(ω, den::T) = exp(- ω*ω/den)`.
The Gaussian bell curve is defined as `gauss(ω) = exp(- ω² / 2σ²)`.
For performance reasons, we memoize the denominator `2σ²`,
which is computed only once in the caller `morlet`.
Also note that the exponentiation `ω²` is replaced by the explicit
product `ω*ω`.
"""
gauss{T<:Real}(ω, den::T) = @fastmath Base.convert(T, exp(- ω*ω/den))::T

function morlet{T<:FFTW.fftwReal}(::FourierDomain{1},
        center::T, den::T, N::Int, nPeriods::Int)
    """1. **Compute range of frequencies with nonneglible magnitude**"""
    halfN = N >> 1
    pstart = - ((nPeriods-1)>>1)
    pstop = (nPeriods-1)>>1 + iseven(nPeriods)
    ωstart = - halfN + pstart * N
    ωstop = halfN + pstop * N - 1
    """2. **Compute Gaussians**"""
    @inbounds begin
        gauss_center = T[ gauss(ω-center, den) for ω in ωstart:ωstop ]
        gauss_0 = T[ gauss(ω, den)
            for ω in (ωstart + pstart*N):(ωstop + pstop*N) ]
        corrective_gaussians = T[ gauss_0[1 + ω + p*N]
            for ω in 0:(N*nPeriods-1), p in 0:(nPeriods-1) ]
    end
    """3. **Solve linear system to find out the corrective factors**"""
    b = T[ gauss(p*N - center, den) for p in pstart:pstop ]
    A = T[ gauss((q-p)*N, den)
        for p in 0:(nPeriods-1), q in 0:(nPeriods-1) ]
    corrective_factors = A \ b
    """4. **Periodize in Fourier domain**"""
    y = gauss_center - corrective_gaussians * corrective_factors
    y = reshape(y, N, nPeriods)
    y = squeeze(sum(y, 2), 2)
end


text:
mutable struct LinkedStack{T} <: AbstractStack{T}
    top::LinkNode{T}
    size::Int
    LinkedStack{T}() where T = new(LinkNode{T}(), 0)
end

isempty(l::LinkedStack) = l.size == 0
length(l::LinkedStack) = l.size
topvalue(l::LinkedStack) = isempty(l) ? throw(ArgumentError("stack must be non-empty")) : l.top.data

function push!(l::LinkedStack, item)
    l.top = LinkNode(l.top, item)
    l.size += 1
end

function pop!(l::LinkedStack)
    (l.size == 0) && throw(ArgumentError("stack must be non-empty"))
    item = l.top.data
    l.top = l.top.next
    l.size -= 1
    return item
end


text:
struct DescriptorPoolState
    allocated::Dictionary{DescriptorType,Int}
    limits::Dictionary{DescriptorType,Int}
    max::Int
end

struct PoolAllocations
    set::DescriptorSet
    pool::DescriptorPool
    allocated::Dictionary{DescriptorType,Int}
end

"""
Entity which manages descriptor pools and descriptor sets.
"""
struct DescriptorAllocator
    set_layout::DescriptorSetLayout
    device::Device
    pools::Dictionary{DescriptorPool,DescriptorPoolState}
    allocations::Dictionary{DescriptorSet,PoolAllocations}
end

Base.broadcastable(da::DescriptorAllocator) = Ref(da)

DescriptorAllocator(device) = DescriptorAllocator(device, Dictionary{DescriptorPool,DescriptorPoolState}(), Dictionary{DescriptorSet,PoolAllocations}())

function new_pool!(da::DescriptorAllocator, allocations::Dictionary{DescriptorType,Int}; next = C_NULL, flags = DescriptorPoolCreateFlag(0))
    limits = [DescriptorPoolSize(t, max(allocations[t], 20)) for t in keys(allocations)]
    maxsize = max(sum(allocations), 50)
    create_info = DescriptorPoolCreateInfo(maxsize, limits; next, flags = flags | DESCRIPTOR_POOL_CREATE_FREE_DESCRIPTOR_SET_BIT)
    new = Created(unwrap(create_descriptor_pool(da.device, create_info)), create_info)
    limits_dict = Dictionary(getproperty.(limits, :type), Int.(getproperty.(limits, :descriptor_count)))
    insert!(da.pools, handle(new), DescriptorPoolState(Dictionary{DescriptorType,Int64}(), limits_dict, maxsize))
    new
end

function compute_allocations(da::DescriptorAllocator, layouts::AbstractVector{DescriptorSetLayoutCreateInfo})
    allocations = Dictionary{DescriptorType,Int}()
    for layout in layouts
        for binding in layout.bindings
            set!(allocations, binding.descriptor_type, get(allocations, binding.descriptor_type, 0) + binding.descriptor_count)
        end
    end
    allocations
end

function compute_allocations(da::DescriptorAllocator, layout::DescriptorSetLayoutCreateInfo)
    dictionary(map(layout.bindings) do binding
        binding.descriptor_type => Int(binding.descriptor_count)
    end)
end

function compute_allocations(da::DescriptorAllocator, layouts::AbstractVector{Created{DescriptorSetLayout, DescriptorSetLayoutCreateInfo}})
    compute_allocations(da, info.(layouts))
end

function compute_allocations(da::DescriptorAllocator, layouts::AbstractVector{DescriptorSetLayout})
    compute_allocations(da, [da.layouts[layout] for layout in layouts])
end

function find_pool!(da::DescriptorAllocator, layouts)
    find_pool!(da, compute_allocations(da, layouts))
end

function find_pool!(da::DescriptorAllocator, allocations::Dictionary{DescriptorType,Int})
    for (pool, state) in zip(keys(da.pools), da.pools)
        n = 0
        has_room = sum(allocations) < sum(state.max) && all(keys(state.allocated)) do t
            state.allocated[t] + allocations[t] < state.limits[t]
        end
        if has_room
            return pool
        end
    end

    # no pool could be reused
    new_pool!(da, allocations)
end

function allocate_descriptor_sets!(da::DescriptorAllocator, layouts::Vector{Created{DescriptorSetLayout,DescriptorSetLayoutCreateInfo}}; next=C_NULL)
    pool = find_pool!(da, layouts)
    allocations = compute_allocations(da, info.(layouts))

    for (type, count) in pairs(allocations)
        allocated = da.pools[pool].allocated
        set!(allocated, type, get(allocated, type, 0) + count)
    end

    allocate_info = DescriptorSetAllocateInfo(next, pool, layouts)
    sets = unwrap(allocate_descriptor_sets(da.device, allocate_info))

    # record allocations for ulterior deletion
    foreach(zip(layouts, sets)) do (layout, set)
        insert!(da.allocations, set, PoolAllocations(set, pool, compute_allocations(da, info(layout))))
    end

    Created.(sets, allocate_info)
end

struct DescriptorSetCache
    used::Dictionary{UInt,Vector{DescriptorSet}}
    unused::Dictionary{DescriptorSet,Int8}
    recycled::Dictionary{UInt,DescriptorSet}
    allocator::DescriptorAllocator
end

function DescriptorSetCache(device::Device)
    DescriptorSetCache(
        Dictionary{UInt,Vector{DescriptorSet}}(),
        Dictionary{DescriptorSet,Int8}(),
        Dictionary{UInt,DescriptorSet}(),
        DescriptorAllocator(device),
    )
end

function get_descriptor_set!(cache::DescriptorSetCache, key::UInt)

end

function cycle!(cache::DescriptorSetCache)

end


text:
@diffusion_process FitzHughDiffusionAux begin
    :dimensions
    process --> 2
    wiener --> 1

    :parameters
    (ϵ, s, γ, β, σ, t, T) --> Float64
    (u, v) --> ℝ{2}

    :additional
    constdiff --> true
    linear --> true
end

function B(t, P::FitzHughDiffusionAux)
    @SMatrix [
        1/P.ϵ-3*P.v^2/P.ϵ   -1/P.ϵ;
        P.γ                 -1.0
    ]
end

function β(t, P::FitzHughDiffusionAux)
    ℝ{2}(
        P.s/P.ϵ+2*P.v^3/P.ϵ,
        P.β
    )
end

σ(t, P::FitzHughDiffusionAux) = ℝ{2}(0.0, P.σ)


text:
######################################################################
# Pia - an "internet" of peers.
# -----
# Licensed under Apache License 2.0
module Pia
include("./Errors.jl")

include("./Futures.jl")
include("./Lockables.jl")
include("./Logging.jl")
include("./Packets.jl")
include("./Signals.jl")
include("./Voidables.jl")

include("./utils.jl")

include("./Pia.Interface.jl")
include("./PacketFactories.jl")

include("./Pia.Applications.jl")
include("./Pia.LinkProtocol.jl")
include("./Pia.PortalProtocol.jl")

using Reexport
@reexport using .LinkProtocol
@reexport using .PortalProtocol
@reexport using .Applications
import .Logging

const DEFAULT_PORT = 27472
      DEFAULT_TIMEOUT = 10
const DEFAULT_DESIRED_CONNECTIONS = 32
const DEFAULT_CONNECTION_HARDCAP = 100
const DEFAULT_PEER_LIFETIME = 60

const LogConnectivity = Logging.LogCategory("Pia Connectivity")
const LogDiscovery    = Logging.LogCategory("Pia Discovery")

function set_default_timeout(timeout::Real)
    global DEFAULT_TIMEOUT
    DEFAULT_TIMEOUT = timeout
end

end # module


text:
using SparseArrays
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

V, (VV,EV,FV,CV) = Lar.simplex(3, true)
tetra = V, EV,FV,CV
twotetra = Lar.Struct([ tetra, Lar.t(0.25,0.25,0.25), tetra ])
V,EV,FV,CV = Lar.struct2lar(twotetra)
GL.VIEW([ GL.GLGrid(V,FV, GL.Point4d(1,1,1,0.2)) ]);

cop_EV = Lar.coboundary_0(EV::Lar.Cells);
cop_EW = convert(Lar.ChainOp, cop_EV);
cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
W = convert(Lar.Points, V');

V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement( W::Lar.Points, cop_EW::Lar.ChainOp, cop_FE::Lar.ChainOp)

EV = Lar.cop2lar(copEV)
FE = [findnz(copFE[k,:])[1] for k=1:size(copFE,1)]
FV = [collect(Set(cat(EV[e] for e in FE[f]))) for f=1:length(FE)]
FV = convert(Lar.Cells, FV)
W = convert(Lar.Points, V')
WW = [[k] for k=1:size(W,2)]

GL.VIEW(GL.numbering(.5)((W,[WW,EV]) ));

triangulated_faces = Lar.triangulate(V, [copEV, copFE])
FVs = convert(Array{Lar.Cells}, triangulated_faces)
V = convert(Lar.Points, V')
GL.VIEW(GL.GLExplode(V,FVs,1.5,1.5,1.5,99));

EVs = Lar.FV2EVs(copEV, copFE) # polygonal face fragments
GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));


text:
# # Gaussian profile

struct GaussLaserConstantCache{IC,W,K,T,Z,E}
    inv_c::IC
    ω::W
    k::K
    T₀::T
    z_R::Z
    E₀::E
end

function GaussLaserConstantCache(;c, λ, w₀, a₀, mₑ, q)
    ω = 2π * c / λ
    k = 2π / λ
    T = 2π / ω
    z_R = w₀^2 * k / 2
    E₀ = a₀ * mₑ * c * ω / abs(q)

    return GaussLaserConstantCache(inv(c), ω, k, T, z_R, E₀)
end

@auto_hash_equals mutable struct GaussLaserCache{X,CE}
    x::X
    y::X
    wz::X
    Ex::CE
    Ey::CE
    Ez::CE
end

function GaussLaserCache(λ, E)
    ThreadLocal(GaussLaserCache(
        zero(λ),        # x
        zero(λ),        # y
        zero(λ),        # wz
        zero(E*im),     # Ex
        zero(E*im),     # Ey
        zero(E*im)      # Ez
    ))
end

function Base.fill!(cache::GaussLaserCache, x::AbstractVector)
    cache.x = x[1]
    cache.y = x[2]

    return nothing
end

@doc """
    struct GaussLaser <: AbstractLaser

The `GaussLaser` is defined by the `units` of the laser
(a positional argument that can be `:SI`, `:atomic` or `:SI_unitful` and `:atomic_unitful`)
and the following parameters
- `λ` is the laser wavelangth
- `a₀` is the normalized vector potential (defined as ``a_0=\\frac{eA}{m_e c^2}``)
- `ϕ₀` is the initial phase with the default value 0.0
- `w₀` is the beam radius at the Rayleigh range or [beam waist](https://en.wikipedia.org/wiki/Gaussian_beam#Beam_waist)
- `ξx` and `ξy` give the polarization and have the default value `1.0 + 0im` and `0.0 + 0im`
- `orientation` specifies how the laser is oriented with respect to the
default coordinate system (default `(:x, :z)`). See [`LaserGeometry`](@ref) for more details.
- `propagation_dir` is the propagation direction of the wave in the
intrinsic coordinate system (default `:z`). See [`LaserGeometry`](@ref) for more details.
- `profile` is the temporal profile of the pulse and the default one is constant (infinite pules duration)

During the initialization of the `GaussLaser` type, some useful derived values are
also computed
- `ω` is the angular frequency and is given by ``2π c / λ``
- `k` is the wavenumber and is given by ``2π / λ``
- `z_R` is the [Rayleigh range](https://en.wikipedia.org/wiki/Rayleigh_length) and is given by ``k w_0^2 / 2``
- `T₀` is the laser period and is given by ``2π / ω``
- `E₀` is the amplitude of the electric field and is given by ``a_0\\frac{m_q c \\omega}{q}``
"""
GaussLaser

struct GaussLaser{PhysicsConst,
                  Derived,
                  Cache,
                  Geometry,
                  Polarization,
                  Profile,
                  L,
                  F} <: AbstractLaser
    constants::PhysicsConst
    derived::Derived
    cache::Cache
    geometry::Geometry
    polarization::Polarization
    profile::Profile
    # laser parameters
    λ::L
    a₀::F
    ϕ₀::F
    w₀::L
end

function GaussLaser(units;
        λ,
        a₀,
        ϕ₀ = 0.0,
        w₀,
        ξx = 1.0+0im,
        ξy = 0,
        orientation = (:x, :z),
        profile = ConstantProfile()
    )

    ξx, ξy = promote(ξx, ξy)
    a₀, ϕ₀ = promote(a₀, ϕ₀)
    λ, w₀ = promote(λ, w₀)
    @assert hypot(ξx, ξy) ≈ 1 "Invalid ξx and ξy"
    constants = FundamentalConstants(units)
    @unpack mₑ, c, q = constants

    derived = GaussLaserConstantCache(; c, λ, w₀, a₀, mₑ, q)
    E₀ = derived.E₀

    cache = GaussLaserCache(λ, E₀)

    geometry = LaserGeometry(orientation)

    polarization = LaserPolarization(ξx, ξy)

    return GaussLaser(
        constants,
        derived,
        cache,
        geometry,
        polarization,
        profile,
        λ,
        a₀,
        ϕ₀,
        w₀
    )
end

update_cache!(::GaussLaser{P,D,Cache}, ::Any, ::Any) where {P,D,Cache<:Nothing} = nothing

function required_coords(::GaussLaser, r)
    CylindricalFromCartesian()(r)
end

function Ex(laser::GaussLaser, coords)
    @unpack E₀, k, z_R = immutable_cache(laser)
    ξx = polarization(laser, :ξx)
    @unpack w₀, ϕ₀ = laser
    @unpack r, z = coords

    wz = w(z, laser)
    update_cache!(laser, :wz, wz)

    ξx * E₀ * w₀/wz * exp(-(r/wz)^2 + im*(-(r^2*z)/(z_R*wz^2) + atan(z, z_R) - k*z + ϕ₀))
end

function Ez(laser::GaussLaser, coords)
    @unpack k, z_R = immutable_cache(laser)
    @unpack Ex, Ey, x, y, wz = mutable_cache(laser)
    z = coords.z

    2im/(k*wz^2)*(1 + im*(z/z_R))*(x*Ex + y*Ey)
end

function Bz(laser::GaussLaser, coords)
    @unpack k, z_R = immutable_cache(laser)
    @unpack Ex, Ey, x, y, wz = mutable_cache(laser)
    c = fundamental_constants(laser, :c)
    z = coords.z

    2im/(k*c*wz^2)*(1 + im*(z/z_R))*(y*Ex - x*Ey)
end


text:
# Neighbour finders

export
    NoNeighbourFinder,
    find_neighbours!,
    DistanceNeighbourFinder

"""
    NoNeighbourFinder()

Placeholder neighbour finder that returns no neighbours.
When using this neighbour finder, ensure that `nl_only` for the interactions is
set to `false`.
"""
struct NoNeighbourFinder <: NeighbourFinder end

"""
    find_neighbours!(simulation, neighbour_finder, step_n; parallel=true)

Obtain a list of close atoms in a system.
Custom neighbour finders should implement this function.
"""
function find_neighbours!(s::Simulation,
                            ::NoNeighbourFinder,
                            ::Integer;
                            kwargs...)
    return
end

"""
    DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff)
    DistanceNeighbourFinder(nb_matrix, n_steps)

Find close atoms by distance.
"""
struct DistanceNeighbourFinder{T} <: NeighbourFinder
    nb_matrix::BitArray{2}
    n_steps::Int
    dist_cutoff::T
end

function DistanceNeighbourFinder(nb_matrix::BitArray{2},
                                 n_steps::Integer)
    return DistanceNeighbourFinder(nb_matrix, n_steps, 1.2)
end

function find_neighbours!(s::Simulation,
                            nf::DistanceNeighbourFinder,
                            step_n::Integer;
                            parallel::Bool=true)
    !iszero(step_n % nf.n_steps) && return

    neighbours = s.neighbours
    empty!(neighbours)
    sqdist_cutoff = nf.dist_cutoff ^ 2

    if parallel && nthreads() > 1
        nl_threads = [Tuple{Int, Int}[] for i in 1:nthreads()]

        @threads for i in 1:length(s.coords)
            nl = nl_threads[threadid()]
            ci = s.coords[i]
            nbi = nf.nb_matrix[:, i]
            for j in 1:(i - 1)
                r2 = sum(abs2, vector(ci, s.coords[j], s.box_size))
                if r2 <= sqdist_cutoff && nbi[j]
                    push!(nl, (i, j))
                end
            end
        end

        for nl in nl_threads
            append!(neighbours, nl)
        end
    else
        for i in 1:length(s.coords)
            ci = s.coords[i]
            nbi = nf.nb_matrix[:, i]
            for j in 1:(i - 1)
                r2 = sum(abs2, vector(ci, s.coords[j], s.box_size))
                if r2 <= sqdist_cutoff && nbi[j]
                    push!(neighbours, (i, j))
                end
            end
        end
    end
end


text:
const nucleic = ["GUN","ADE","CYT","THY","URA","DA","DC","DG","DT","DU","A","C","G","T","U",
                "AMP","ADP","ATP","CDP","CTP","GMP","GDP","GTP","TMP","TTP","UMP","UDP","UTP"]
const water = ["HOH","H2O"]
const aas = ["ALA","ARG","ASN","ASP","CYS","GLN","GLU","GLY","HIS","ILE","LEU","LYS","MET","PHE","PRO","SER","THR","TRP","TYR","VAL"]
const lipids = []
#and more

function isheader(linestr::String)
    @inbounds return linestr[1] == '_' || linestr[1] == '#'
end

veclength2(ax,bx,ay,by,az,bz) = (ax-bx)^2 + (ay-by)^2 + (az-bz)^2

const _default_delims = ' '
function splitmmcif(str::AbstractString,atomid,resn,chain,resid,inscode,xlist,ylist,zlist,strs,coords)
    #strs = fill("",5)
    #coords = zeros(Float64,3)
    splitter = _default_delims

    i = start(str)
    n = endof(str)
    r = search(str,splitter,i)
    j, k = first(r), nextind(str,last(r))
    colno = 1

    while 0 < j <= n
        if i < k
            if i < j
                if colno == atomid
                    strs[1] = str[i:j-1]
                elseif colno == resn
                    strs[2] = str[i:j-1]
                elseif colno == chain
                    strs[3] = str[i:j-1]
                elseif colno == resid
                    strs[4] = str[i:j-1]
                elseif colno == inscode
                    strs[5] = str[i:j-1]
                elseif colno == xlist
                    coords[1] = coordsfloat(str[i:j-1])
                elseif colno == ylist
                    coords[2] = coordsfloat(str[i:j-1])
                elseif colno == zlist
                    coords[3] = coordsfloat(str[i:j-1])
                end
                colno += 1
            end
            i = k
        end
        (k <= j) && (k = nextind(str,j))
        r = search(str,splitter,k)
        j, k = first(r), nextind(str,last(r))
    end
    if !done(str,i)
        if colno == atomid
            strs[1] = str[i:j-1]
        elseif colno == resn
            strs[2] = str[i:j-1]
        elseif colno == chain
            strs[3] = str[i:j-1]
        elseif colno == resid
            strs[4] = str[i:j-1]
        elseif colno == inscode
            strs[5] = str[i:j-1]
        elseif colno == xlist
            coords[1] = coordsfloat(str[i:j-1])
        elseif colno == ylist
            coords[2] = coordsfloat(str[i:j-1])
        elseif colno == zlist
            coords[3] = coordsfloat(str[i:j-1])
        end
    end
    #return strs,coords
end

function coordsfloat(str)
    ptr = pointer(str)
    #parse str(xxxxx.xxxxx) -> float
    negaflag = unsafe_load(ptr,1) == 0x2d ? -1 : 1
    n = endof(str)
    i = (3 - negaflag) >> 1
    decspace = 0
    isdot = false
    val = 0.0

    for idx = i:n
        digit = unsafe_load(ptr,idx)
        if digit == 0x2e
            decspace = 0
            isdot = true
            continue
        end
        val *= 10
        val += digit - 0x30
        decspace += 1
    end
    if isdot
        val /= 10^decspace
    end

    val *= negaflag
    return val
end

#--------------------------------------MAIN
function read_mmcif(datalines,thresholdlength=5.0::Float64,targetmoleculelists=nucleic) #change here or specify the target explicitly
    const thresholdlengthpow2 = thresholdlength^2

    #find the _atom_site-loop and memo column no of coords
    const atomcolumn = findfirst(x -> contains(x,"_atom_site."),datalines) -1 #line of "loop_"
    const datastart = findnext(!isheader,datalines,atomcolumn+1)
    const datalast = findnext(isheader,datalines,datastart) -1 #previous line of "#"

    const atomid = findnext(x -> contains(x,"_atom_site.auth_atom_id"),datalines,atomcolumn) - atomcolumn
    const resn = findnext(x -> contains(x,"_atom_site.auth_comp_id"),datalines,atomcolumn) - atomcolumn
    const chain = findnext(x -> contains(x,"_atom_site.auth_asym_id"),datalines,atomcolumn) - atomcolumn
    const resid = findnext(x -> contains(x,"_atom_site.auth_seq_id"),datalines,atomcolumn) - atomcolumn
    const inscode = findnext(x -> contains(x,"_atom_site.pdbx_PDB_ins_code"),datalines,atomcolumn) - atomcolumn
    const xlist = findnext(x -> contains(x,"_atom_site.Cartn_x"),datalines,atomcolumn) - atomcolumn
    const ylist = findnext(x -> contains(x,"_atom_site.Cartn_y"),datalines,atomcolumn) - atomcolumn
    const zlist = findnext(x -> contains(x,"_atom_site.Cartn_z"),datalines,atomcolumn) - atomcolumn

    #prepare array to store (textbase/coordinary) data of (amino acids/ligands)
    aalist = fill("",(5,(datalast-datastart)))
    aacoords = Array{Float64,2}(3,(datalast-datastart))
    liglist = fill("",(5,(datalast-datastart)))
    ligcoords = Array{Float64,2}(3,(datalast-datastart))

    aacount = 1
    ligcount = 1
    aap = pointer(aacoords)
    ligp = pointer(ligcoords)
    strpart = fill("",5)
    coordspart = zeros(Float64,3)

    for lineno = datastart:datalast
        #split one line text to array
        splitmmcif(datalines[lineno],atomid,resn,chain,resid,inscode,xlist,ylist,zlist,strpart,coordspart)

        if strpart[2] in aas #canonical aa.
            #I can't use pointer for aa/lig-list, because something cause error.
            @inbounds aalist[:,aacount] = strpart
            @inbounds aacoords[:,aacount] = coordspart
            aacount += 1
        end
        if strpart[2] in targetmoleculelists
            @inbounds liglist[:,ligcount] = strpart
            @inbounds ligcoords[:,ligcount] = coordspart
            ligcount += 1
        end
    end

    aacount -= 1
    ligcount -= 1

    #cut unnecessary region of coords matrix.
    #text-base information matrixes may not be needed to cut, so I don't.
    @inbounds aacoords = aacoords[:,1:aacount]
    @inbounds ligcoords = ligcoords[:,1:ligcount]

    #-------------------------------------------O(N) search from here
    @inbounds xmin = min(minimum(aacoords[1,:]),minimum(ligcoords[1,:]))
    @inbounds xmax = max(maximum(aacoords[1,:]),maximum(ligcoords[1,:]))
    @inbounds ymin = min(minimum(aacoords[2,:]),minimum(ligcoords[2,:]))
    @inbounds ymax = max(maximum(aacoords[2,:]),maximum(ligcoords[2,:]))
    @inbounds zmin = min(minimum(aacoords[3,:]),minimum(ligcoords[3,:]))
    @inbounds zmax = max(maximum(aacoords[3,:]),maximum(ligcoords[3,:]))

    #make the grid, it's like a hash function of 3d-space of the mmcif
    xgrid = ceil(Int,(xmax-xmin)/thresholdlength)
    ygrid = ceil(Int,(ymax-ymin)/thresholdlength)
    zgrid = ceil(Int,(zmax-zmin)/thresholdlength)
    numgrid = xgrid * ygrid * zgrid

    #using sortedbuffer (, count) and  indexes to find of aacoords
    count = zeros(Int64,numgrid) #the no of aa. atoms in each grid
    indexes = zeros(Int64,numgrid+1) #the 1st index of each grid in sortedbuffer
    particlepos = zeros(Int64,aacount) #the grid of each aa. atom
    sortedbuffer = zeros(Int64,aacount) #the list of aa. atoms sorted by grid no.

    for i = 1:aacount
        @inbounds ix = floor(Int, (aacoords[1,i] - xmin) / thresholdlength)
        @inbounds iy = floor(Int, (aacoords[2,i] - ymin) / thresholdlength)
        @inbounds iz = floor(Int, (aacoords[3,i] - zmin) / thresholdlength)

        #grid no. index  "+1": julia is 1-origin.
        idx = floor(Int, ix + iy * xgrid + iz * xgrid * ygrid) + 1
        @inbounds count[idx] += 1
        @inbounds particlepos[i] = idx
    end

    sum = 0
    for i in 1:numgrid
        @inbounds sum += count[i]
        @inbounds indexes[i+1] = sum
    end

    register = zeros(Int64,numgrid) #temp. values how many aa. atoms are already put into sortedbuffer

    #SORT
    for i in 1:aacount
        @inbounds pos = particlepos[i]
        @inbounds j = indexes[pos] + register[pos] + 1
        @inbounds sortedbuffer[j] = i
        @inbounds register[pos] += 1
    end

    #redefine pointers
    aap = pointer(aacoords)
    ligp = pointer(ligcoords)

    #--------------------------------------------------------
    for tgtidx = 1:ligcount

        @inbounds target = liglist[:,tgtidx]
        #@inbounds tgtcoords = ligcoords[:,tgtidx]
        tx = unsafe_load(ligp,3*tgtidx-2)
        ty = unsafe_load(ligp,3*tgtidx-1)
        tz = unsafe_load(ligp,3*tgtidx)

        #find the grid to where the target atom should be belong
        ix = floor(Int, (tx - xmin) / thresholdlength)
        iy = floor(Int, (ty - ymin) / thresholdlength)
        iz = floor(Int, (tz - zmin) / thresholdlength)

        #Because grid length is equal to thresholdlength,
        #the sought aa. atoms within the threshold length must be in the same or the 26-adjacent grid(s).
        for bufx = -1:1
            for bufy = -1:1
                for bufz = -1:1
                    meshsearch(thresholdlengthpow2,aalist,aap,target,tx,ty,tz,sortedbuffer,count,indexes,xgrid,ygrid,zgrid,ix+bufx,iy+bufy,iz+bufz)
                end
            end
        end
    end

end

function meshsearch(thresholdlengthpow2,aalist,aap,target,tx,ty,tz,sortedbuffer,count,indexes,xgrid,ygrid,zgrid,ix,iy,iz)
    if ix < 0 || xgrid <= ix || iy < 0 || ygrid <= iy || iz < 0 || zgrid <= iz
        return 1
    end

    #index of searching grid
    id2 = ix + iy * xgrid + iz * xgrid * ygrid + 1
    @inbounds for m in indexes[id2]+1:indexes[id2]+count[id2]
        @inbounds aaidx = sortedbuffer[m]

        ax = unsafe_load(aap,3*aaidx-2)
        ay = unsafe_load(aap,3*aaidx-1)
        az = unsafe_load(aap,3*aaidx)

        #slight faster than bother to find root
        #lenpow2 = veclength2(vectors)
        lenpow2 = veclength2(ax,tx,ay,ty,az,tz)
        if lenpow2 <= thresholdlengthpow2
            @inbounds aainfo1 = aalist[1,aaidx]
            @inbounds aainfo2 = aalist[2,aaidx]
            @inbounds aainfo3 = aalist[3,aaidx]
            @inbounds aainfo4 = aalist[4,aaidx]
            @inbounds aainfo5 = aalist[5,aaidx]
            #only print 5 decimal space: you can change .[5]f
            @printf "%s.%s.%s.%s.%s:%s.%s.%s.%s:%.5f \n" aainfo2 aainfo3 aainfo5 aainfo4 aainfo1 target[2] target[3] target[4] target[1] (lenpow2^0.5)
        end
    end
end

if length(ARGS) == 1
    readcif = readlines(ARGS[1])
    read_mmcif(readcif)
elseif length(ARGS) == 2
    readcif = readlines(ARGS[1])
    threlen = parse(Float64, ARGS[2])
    read_mmcif(readcif,threlen)
else
    print("Input MMcif file needed: \n /your/path/to/julia mmcif.jl *.cif (threshold_length[Ang.] : optional, default = 5.0)")
end
#=
hoge = readlines("1vy4.cif")
read_mmcif(hoge)
for i=1:5
    @time read_mmcif(hoge)
end
#Profile.print()
=#

max_stars_repo_name: BILAB/scripts

max_issues_repo_name: BILAB/scripts

max_forks_repo_name: BILAB/scripts


text:
using AcFDTD
using Base.Test

srand(123)

include("schemes_test.jl")
include("geometry_test.jl")
include("fdtd_test.jl")
include("adjoint_test.jl")

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


text:
using NPZ
using Distributions
using LinearAlgebra

#irit

"""
    generate_gaussian_data(N::Int64, D::Int64, K::Int64,MixtureVar::Number)

Generate `N` observations, generated from `K` `D` dimensions Gaussians, with the Gaussian means sampled from a `Normal` distribution with mean `0` and `MixtureVar` variance.

Returns `(Samples, Labels, Clusters_means, Clusters_cov)`

# Example
```julia
julia> x,y,clusters = generate_gaussian_data(10000,2,6,100.0)
[3644, 2880, 119, 154, 33, 3170]
...
```
"""
function generate_gaussian_data(N::Int64, D::Int64, K::Int64,MixtureVar::Number)
	x = randn(Float32,D,N)
	tpi = rand(Dirichlet(ones(K)))
	tzn = rand(Multinomial(N,tpi))
	tz = zeros(Float32,N)

	tmean = zeros(Float32,D,K)
	tcov = zeros(Float32,D,D,K)

	ind = 1
	println(tzn)
	for i=1:length(tzn)
		indices = ind:ind+tzn[i]-1
		tz[indices] .= i
		tmean[:,i] .= rand(MvNormal(zeros(Float32,D), MixtureVar*Matrix{Float32}(I, D, D)))
		tcov[:,:,i] .= rand(InverseWishart(D+2, Matrix{Float32}(I, D, D)))
		#println(minimum(eigvals(tcov[:,:,i]))>0)
		d = MvNormal(tmean[:,i], tcov[:,:,i])
		for j=indices
			x[:,j] = rand(d)
		end
		ind += tzn[i]
	end
	return x, tz, tmean, tcov
end

"""
     generate_mnmm_data(N::Int64, D::Int64, K::Int64, trials::Int64)

Generate `N` observations, generated from `K` `D` features Multinomial vectors, with `trials` draws from each vector.

Returns `(Samples, Labels, Vectors)`

# Example
```julia
julia> generate_mnmm_data(10000, 10, 5, 100)
...
```
"""
function generate_mnmm_data(N::Int64, D::Int64, K::Int64, trials::Int64)
	clusters = zeros(Float64,D,K)
	x = zeros(Float32,D,N)
	labels = rand(1:K,(N,))
	for i=1:K
		alphas = rand(1:20,(D,))
		alphas[i] = rand(30:100)
		clusters[:,i] = rand(Dirichlet(alphas))
	end
	for i=1:N
		x[:,i] = rand(Multinomial(trials,clusters[:,labels[i]]))
	end
	return x, labels, clusters
end


text:
include("utils.jl")

@testset "SymmetricMatrices" begin
    M=SymmetricMatrices(3,ℝ)
    A = [1 2 3; 4 5 6; 7 8 9]
    A_sym = [1 2 3; 2 5 -1; 3 -1 9]
    A_sym2 = [1 2 3; 2 5 -1; 3 -1 9]
    B_sym = [1 2 3; 2 5 1; 3 1 -1]
    M_complex = SymmetricMatrices(3,ℂ)
    C = [1+im 1 im; 1 2 -im; im -im -1-im]
    D = [1 0; 0 1];
    X = zeros(3,3)
    @testset "Real Symmetric Matrices Basics" begin
        @test representation_size(M) == (3,3)
        @test check_manifold_point(M,B_sym)===nothing
        @test_throws DomainError is_manifold_point(M,A,true)
        @test_throws DomainError is_manifold_point(M,C,true)
        @test_throws DomainError is_manifold_point(M,D,true)
        @test_throws DomainError is_manifold_point(M_complex, [:a :b :c; :b :d :e; :c :e :f],true)
        @test check_tangent_vector(M,B_sym,B_sym)===nothing
        @test_throws DomainError is_tangent_vector(M,B_sym,A,true)
        @test_throws DomainError is_tangent_vector(M,A,B_sym,true)
        @test_throws DomainError is_tangent_vector(M,B_sym,D,true)
        @test_throws DomainError is_tangent_vector(M,B_sym, 1*im * zero_tangent_vector(M,B_sym),true)
        @test_throws DomainError is_tangent_vector(M_complex, B_sym, [:a :b :c; :b :d :e; :c :e :f],true)
        @test manifold_dimension(M) == 6
        @test manifold_dimension(M_complex) == 12
        @test A_sym2 == project_point!(M,A_sym)
        @test A_sym2 == project_tangent(M,A_sym,A_sym)
    end
    types = [
        Matrix{Float64},
        MMatrix{3,3,Float64},
        Matrix{Float32},
    ]
    bases = (ArbitraryOrthonormalBasis(), ProjectedOrthonormalBasis(:svd))
    for T in types
        pts = [convert(T,A_sym),convert(T,B_sym),convert(T,X)]
        @testset "Type $T" begin
            test_manifold(
                M,
                pts,
                test_injectivity_radius = false,
                test_reverse_diff = isa(T, Vector),
                test_project_tangent = true,
                test_musical_isomorphisms = true,
                test_vector_transport = true,
                basis_types_vecs = bases,
                basis_types_to_from = bases
            )
            test_manifold(
                M_complex,
                pts,
                test_injectivity_radius = false,
                test_reverse_diff = isa(T, Vector),
                test_project_tangent = true,
                test_musical_isomorphisms = true,
                test_vector_transport = true,
                basis_types_vecs = (ArbitraryOrthonormalBasis(),),
                basis_types_to_from = (ArbitraryOrthonormalBasis(),)
            )
            @test isapprox(-pts[1], exp(M, pts[1], log(M, pts[1], -pts[1])))
        end # testset type $T
    end # for
end # test SymmetricMatrices


text:
using Base.Test

srand(87687634)
@test length(pow2edges(0.1,1:100))==1025
@test isapprox(first(pow2edges(0.1,1:100)),1)
@test isapprox(last(pow2edges(0.1,1:100)),100)

known_dc = LinearDriftCorrect(0,0.0001)
n=10000
indicator = randn(n)*10+known_dc.indicator_median
ideal_corrected = 1000+randn(n)
uncorrected = ideal_corrected./(1+indicator*known_dc.param)
found_dc = drift_correct(indicator, uncorrected, 0.5)
@test abs(found_dc.indicator_median-known_dc.indicator_median)<1.0
@test abs(1-found_dc.param/known_dc.param)<0.05

corrected = applycorrection(found_dc, indicator, uncorrected)
known_dc = LinearDriftCorrect(0,0.0001)
known_corrected = applycorrection(known_dc, indicator, uncorrected)

@test all(known_corrected-ideal_corrected.<1e-8)


text:
# ============================================================
# Structs of Interface Fluxes
# compatible with solution structs
# ============================================================

mutable struct Flux1D{A,B} <: AbstractFlux1D

    w::A
    fw::B

    function Flux1D(w::AbstractArray, fw::AbstractArray)
        new{typeof(w),typeof(fw)}(w, fw)
    end

end

mutable struct Flux1D1F{A,B,C} <: AbstractFlux1D

    w::A
    fw::B
    ff::C

    function Flux1D1F(w::AbstractArray, fw::AbstractArray, ff::AbstractArray)
        new{typeof(w),typeof(fw),typeof(ff)}(w, fw, ff)
    end

end

mutable struct Flux1D2F{A,B,C} <: AbstractFlux1D

    w::A
    fw::B
    fh::C
    fb::C

    function Flux1D2F(
        w::AbstractArray,
        fw::AbstractArray,
        fh::AbstractArray,
        fb::AbstractArray,
    )
        new{typeof(w),typeof(fw),typeof(fh)}(w, fw, fh, fb)
    end

end

mutable struct Flux2D{A,B,C} <: AbstractFlux2D

    n1::A
    w1::B
    fw1::C

    n2::A
    w2::B
    fw2::C

    function Flux2D(
        n1::AbstractArray,
        w1::AbstractArray,
        fw1::AbstractArray,
        n2::AbstractArray,
        w2::AbstractArray,
        fw2::AbstractArray,
    )
        new{typeof(n1),typeof(w1),typeof(fw1)}(n1, w1, fw1, n2, w2, fw2)
    end

end

mutable struct Flux2D1F{A,B,C,D} <: AbstractFlux2D

    n1::A
    w1::B
    fw1::C
    ff1::D

    n2::A
    w2::B
    fw2::C
    ff2::D

    function Flux2D1F(
        n1::AbstractArray,
        w1::AbstractArray,
        fw1::AbstractArray,
        ff1::AbstractArray,
        n2::AbstractArray,
        w2::AbstractArray,
        fw2::AbstractArray,
        ff2::AbstractArray,
    )
        new{typeof(n1),typeof(w1),typeof(fw1),typeof(ff1)}(
            n1,
            w1,
            fw1,
            ff1,
            n2,
            w2,
            fw2,
            ff2,
        )
    end

end

mutable struct Flux2D2F{A,B,C,D} <: AbstractFlux2D

    n1::A
    w1::B
    fw1::C
    fh1::D
    fb1::D

    n2::A
    w2::B
    fw2::C
    fh2::D
    fb2::D

    function Flux2D2F(
        n1::AbstractArray,
        w1::AbstractArray,
        fw1::AbstractArray,
        fh1::AbstractArray,
        fb1::AbstractArray,
        n2::AbstractArray,
        w2::AbstractArray,
        fw2::AbstractArray,
        fh2::AbstractArray,
        fb2::AbstractArray,
    )
        new{typeof(n1),typeof(w1),typeof(fw1),typeof(fh1)}(
            n1,
            w1,
            fw1,
            fh1,
            fb1,
            n2,
            w2,
            fw2,
            fh2,
            fb2,
        )
    end

end


text:
# utils for Graff / Caesar style BigData

function getPointCloudFromKinect(data, dcamjl, imshape)
  ri,ci = imshape[1], imshape[2] # TODO -- hack should be removed since depth is array and should have rows and columns stored in Mongo
  arr = bin2arr(data, dtype=Float32) # should also store dtype for arr in Mongo
  img = reshape(arr, ci, ri)'
  reconstruct(dcamjl, Array{Float64,2}(img))
end

function getPointCloudFromBSON(data)
  buf = IOBuffer(data)
  st = takebuf_string(buf)
  bb = BSONObject(st)
  return map(x -> convert(Array, x), bb["pointcloud"])
end


text:
module TestNaive
	using Base.Test
    using Distance
	using NearestNeighbors

    X = readcsv(Pkg.dir("NearestNeighbors", "test", "iris.csv"))

	t = NaiveNeighborTree(X)

	v = X[:, 21]

	inds, dists = nearest(t, v, 1)
	inds, dists = nearest(t, v, 2)
	inds, dists = nearest(t, v, 3)

	inds, dists = nearest(t, v, 1, 21)
	inds, dists = nearest(t, v, 2, 21)
	inds, dists = nearest(t, v, 3, 21)

	inds, dists = inball(t, v, 0.5)
	inds, dists = inball(t, v, 0.5, 21)
end


text:
using Brownian
using Test

@testset "Brownian" begin

t = LinRange(0, 1, 5)
p = BrownianMotion(t)

rand(p)
rand([p, p])

@test true

end


text:
mutable struct Handle
    fd::IOStream

    # Open the os file - attach a finalizer so the file is appropriately closed.
    function Handle(str::String)
        handle = new(open(str; read = true, write = true))
        finalizer(close, handle)
        return handle
    end
end
Base.close(P::Handle) = close(P.fd)

# Construct the Handle from a group, bus, device, and offset
#
# This is largely reverse engineered from `pci.cpp` from `https://github.com/opcm/pcm`
Handle(bus, device, fn) = Handle(pcipath(bus, device, fn))
function pcipath(bus, device, fn)
    bus_str    = string(bus; base = 16, pad = 2)
    device_str = string(device; base = 16, pad = 2)
    fn_str     = string(fn; base = 16, pad = 1)
    path = joinpath(
        "/proc/bus/pci",
        bus_str,
        "$device_str.$fn_str"
    )

    return path
end

Base.read(P::Handle, ::Type{T}) where {T <: Integer} = read(P.fd, T)
Base.write(P::Handle, v) = write(P.fd, v)
Base.seek(P::Handle, offset::IndexZero) = seek(P.fd, value(offset))

# Wrapping the `pread` calls is MUCH faster because it only requires a single
# transition into kernel code.
function Base.unsafe_read(
        P::Handle,
        ::Type{T},
        offset::IndexZero;
        buffer = Vector{UInt8}(undef, sizeof(T))
    ) where {T}

    # Make a direct `pread` call
    nb = ccall(
        :pread,
        Csize_t,
        (Cint, Ptr{Cvoid}, Csize_t, Cint),
        fd(P.fd),
        buffer,
        sizeof(T),
        value(offset)
    )

    return first(reinterpret(T, buffer))
end


text:
"""
    SubDataset{D, T, <:AbstractDataset{D,T}, <:SubArray{<:SVector{D,T},1}} <: AbstractDataset{D,T}

A view of an `AbstractDataset`, as returned by the `view` function
or the `@view` macro on a dataset instance. A `SubDataset` is an `AbstractDataset`
of the same type as its parent, so indexing, iteration, and most other functions
can be expected to work in the same way for both the parent and the view.
"""
struct SubDataset{D, T, P<:AbstractDataset{D,T}, S<:SubArray{<:SVector{D,T},1}} <: AbstractDataset{D,T}
    parent::P
    data::S
    function SubDataset(par, data)
        @assert parent(data) === par.data
        P = typeof(par)
        S = typeof(data)
        SV = eltype(P)
        T = eltype(SV)
        D = length(SV)
        new{D,T,P,S}(par, data)
    end
end

function Base.summary(sd::SubDataset{D, T}) where {D, T}
    N = length(sd)
    return "$N-element view of $D-dimensional Dataset{$(T)}"
end

Base.parent(sd::SubDataset) = sd.parent
Base.parentindices(sd::SubDataset) = parentindices(sd.data)

"""
    view(d::Dataset, indices)

Return a view into the parent dataset `d`, as a [`SubDataset`](@ref)
that contains the datapoints of `d` referred to by `indices`.
"""
Base.view(d::AbstractDataset, i) = SubDataset(d, view(d.data, i))

function Base.view(::AbstractDataset, ::Any, ::Any, ::Vararg)
    throw(ArgumentError("Dataset views only accept indices on one dimension"))
end


text:
@testset "Conjugacy" begin
  a = fmpz[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
           1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
           0, -9604, -7056, -12076, -2392, -2253, 952, 46, -16, 0, 0, 0, 0, 0,
           0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, -4900, -140]

  b = fmpz[-4645900, -49391, -3848404, -16744, -15771, 6664, 17066, 470484,
           33488, 3779643, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1,
           0, -2, 0, -26, 0, -8, 0, 1, 0, 0, 1, 0, 9, -55750240, -592692,
           -46180848, -200928, -189252, 79969, 204792, 5645808, 396956,
           45355576, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -10, 0, -8, 0, 0, 0, 0, 1,
           0, 8, 101541620, 1079546, 84115116, 365984, 344709, -145656,
           -373022, -10283436, -692768, -82611077, -4, 0, 0, 0, 0, 0, 0, 0, 0,
           1, -18583040, -197564, -15393616, -66976, -63084, 26656, 68264,
           1881936, 129052, 15118432]

  A = matrix(ZZ, 10, 10, a)
  B = matrix(ZZ, 10, 10, b)

  fl, C = Hecke.isGLZ_conjugate(A, B)
  @test fl
  @test C * A == B * C

  A = matrix(QQ, 10, 10, a)
  B = matrix(QQ, 10, 10, b)

  fl, C = Hecke.isGLZ_conjugate(A, B)
  _C = map_entries(QQ, C)
  @test fl
  @test _C * A == B * _C

  A = fmpq(1, 10) * A
  B = fmpq(1, 10) * B
  fl, C = Hecke.isGLZ_conjugate(A, B)
  _C = map_entries(QQ, C)
  @test fl
  @test _C * A == B * _C

  fl, _ = Hecke.isGLZ_conjugate(QQ[1 2; 3 4], QQ[1 2 3; 4 5 6; 7 8 9])
  @test !fl

  R = ResidueRing(FlintZZ, fmpz(7))
  x, y = R(6), R(6)
  q, r = divrem(x, y)
  @test y * q + r == x
  @test iszero(r)
end


text:
typedef struct Bad *Dab;

struct Bad {
  Dab d;
}

int main() {
  Dab a = new Bad;
  Dab b = a->d;
  b->y = new Bad;
  return 0;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
POMDPs.discount(v::AbstractVector) = discount(v[1])

function log_multitask_performances!(𝒮, tasks, logfn=log_undiscounted_return)
    push!(𝒮.log.extras, logfn([Sampler(t, 𝒮.π) for t in tasks]))
end

function continual_learning(tasks, solver_generator)
    solvers = Solver[]
    history = []
    𝒮 = solver_generator(i=1, tasks=tasks[1:1], history=history)
    for i in 1:length(tasks)
        solve(𝒮, tasks[i])
        push!(solvers, deepcopy(𝒮))
        if i < length(tasks)
            𝒮 = solver_generator(tasks=tasks[1:i+1], solvers=solvers, i=i+1, history=history)
        end
    end
    solvers
end


text:
function char2Str(charTuple::Tuple{Vararg{UInt8}})
    cs = collect(charTuple)
    idx = 1
    for i in eachindex(cs)
        if cs[i] == 0
            break
        end
        idx = i
    end
    cs = Char.(cs[1:idx])

    return String(cs)
end

function rotsFromTuple(rotsTuple::Array{NTuple{9,Int32},1}, nop::Integer)
    r = Array{Int64,3}(undef, 3, 3, nop)
    for i in 1:nop
        r[:,:,i] = reshape([Base.convert(Int64, e) for e in rotsTuple[i]], 3, 3)
    end
    return r
end

function transFromTuple(transTuple::Array{NTuple{3,Float64}}, nop::Integer)
    t = Array{Float64,2}(undef, 3, nop)
    for i in 1:nop
        t[:,i] = [e for e in transTuple[i]]
    end
    return t
end


text:
using QMTK
using Compat.Test
using Compat

@testset "Check Type Property" begin

    # check default eltype
    @test eltype(Bit) == Float32
    @test eltype(Spin) == Float32
    @test eltype(Half) == Float32

    # check up down type
    for each in (Bit, Spin, Half)
        @test typeof(up(each)) == eltype(each)
    end

end

@testset "Check Initialization" begin
    # check initial
    SiteBit = Sites(Bit, 3, 3)
    SiteSpin = Sites(Spin, 3, 3)
    SiteHalf = Sites(Half, 3, 3)

    @test SiteBit.data == zeros(eltype(Bit), 3, 3)
    @test SiteSpin.data == -ones(eltype(Spin), 3, 3)
    @test SiteHalf.data == eltype(Half)(-0.5) * ones(eltype(Half), 3, 3)
end

@testset "Check Array Interface" begin

    SHAPE = (3, 4)
    LENGTH = prod(SHAPE)
    SiteA = Sites(Bit, SHAPE)

    @test eltype(SiteA) == eltype(Bit)
    @test length(SiteA) == LENGTH
    @test ndims(SiteA) == length(SHAPE)
    @test size(SiteA) == SHAPE
    @test size(SiteA, 1) == SHAPE[1]
    @test size(SiteA, 2) == SHAPE[2]
    @test Compat.axes(SiteA) == Compat.axes(SiteA.data)
    @test Compat.axes(SiteA, 1) == Compat.axes(SiteA.data, 1)
    @test eachindex(SiteA) == eachindex(SiteA.data)
    @test stride(SiteA, 1) == stride(SiteA.data, 1)
    @test strides(SiteA) == strides(SiteA.data)
    @test getindex(SiteA, 1, 1) == SiteA.data[1, 1]

    setindex!(SiteA, 2, 1, 1)
    @test SiteA[1, 1] == 2
end

@testset "Check Site Interface" begin
    SHAPE = (3, 4)
    SiteA = Sites(Bit, SHAPE)
    SiteB = copy(SiteA)

    # check if this copy is shallow
    rand!(SiteA)
    @test SiteB[1, 1] == 0

    reset!(SiteA)
    @test SiteA.data == zeros(eltype(Bit), SHAPE)

    SiteB[1, 1] = up(Bit)
    @test flip!(SiteA, 1, 1) == SiteB
end

@testset "Check Site Conversions" begin

    for LT in (Bit, Spin, Half)
        A = Sites(LT, 4, 4)
        set!(A, 142)

        @test convert(Int, A) == 142

        # will throw InexactError if Sites is too large
        @test_throws InexactError convert(Int8, A)
    end

end


text:
"""
tools.jl déclare divers utilitaires
"""

# retourne le nom du node correspondant
function getNameNode(indexNode::Int)
    return string("node_",indexNode)
end

max_stars_repo_name: LucasBaussay/ProjetODC

max_issues_repo_name: LucasBaussay/ProjetODC

max_forks_repo_name: LucasBaussay/ProjetODC


text:
# index sets and vector scatters

###########################################################################
export IS, ISBlock # index sets
# Note: we expose a 1-base Julian index interface, but internally
# PETSc's indices are 0-based.
#TODO: support block versions
mutable struct IS{T}
  p::C.IS{T}
  function IS{T}(p::C.IS{T}) where {T}
    o = new{T}(p)
    if (!deactivate_finalizers)
      finalizer(PetscDestroy,o)
    end
    return o
  end
end

comm(a::IS{T}) where{T} = MPI.Comm(C.PetscObjectComm(T, a.p.pobj))

function PetscDestroy(o::IS{T}) where{T}
  PetscFinalized(T) || C.ISDestroy(Ref(o.p))
end

# internal constructor, takes array of zero-based indices:
function IS_(::Type{T}, idx::Array{PetscInt}; comm::MPI.Comm=MPI.COMM_WORLD) where {T<:Scalar}
  is_c = Ref{C.IS{T}}()
  chk(C.ISCreateGeneral(comm, PetscInt(length(idx)), idx, C.PETSC_COPY_VALUES, is_c))
  return IS{T}(is_c[])
end

IS(::Type{T}, idx::AbstractArray{I}; comm::MPI.Comm=MPI.COMM_WORLD) where {I<:Integer, T<:Scalar} =
  IS_(T, PetscInt[i-1 for i in idx]; comm=comm)

function IS(::Type{T}, idx::AbstractRange{I}; comm::MPI.Comm=MPI.COMM_WORLD) where {I<:Integer, T<:Scalar}
  is_c = Ref{C.IS{T}}()
  chk(C.ISCreateStride(comm, PetscInt(length(idx)), PetscInt(first(idx)-1), PetscInt(step(idx)), is_c))
  return IS{T}(is_c[])
end

# there is no Strided block index set, so convert everything to an array
function ISBlock(::Type{T}, bs::Integer,  idx::AbstractArray{I}; comm=MPI.COMM_WORLD) where {I<:Integer, T<:Scalar}
  idx_0 = PetscInt[i-1 for i in idx]
  return ISBlock_(T, bs, idx_0, comm=comm)
end

function ISBlock_(::Type{T}, bs::Integer, idx::AbstractArray{PetscInt};comm=MPI.COMM_WORLD) where {T}
  is_c = Ref{C.IS{T}}()
  chk(C.ISCreateBlock(comm, PetscInt(bs), PetscInt(length(idx)), idx, C.PETSC_COPY_VALUES, is_c))
  return IS{T}(is_c[])
end

function Base.copy(i::IS{T}) where {T}
  is_c = Ref{C.IS{T}}()
  chk(C.ISDuplicate(i.p, is_c))
  return IS{T}(is_c[])
end

function Base.length(i::IS)
  len = Ref{PetscInt}()
  chk(C.ISGetSize(i.p, len))
  return Int(len[])
end

function lengthlocal(i::IS)
  len = Ref{PetscInt}()
  chk(C.ISGetLocalSize(i.p, len))
  return Int(len[])
end

import Base.==
function ==(i::IS{T}, j::IS{T}) where {T}
  b = Ref{PetscBool}()
  chk(C.ISEqual(i.p, j.p, b))
  return b[] != 0
end

function Base.sort!(i::IS)
  chk(C.ISSort(i.p))
  return i
end
Base.sort(i::IS) = sort!(copy(i))

function Base.issorted(i::IS)
  b = Ref{PetscBool}()
  chk(C.ISSorted(i.p, b))
  return b[] != 0
end

function Base.union(i::IS{T}, j::IS{T}) where {T}
  is_c = Ref{C.IS{T}}()
  chk(C.ISExpand(i.p, j.p, is_c))
  return IS{T}(is_c[])
end

function Base.setdiff(i::IS{T}, j::IS{T}) where {T}
  is_c = Ref{C.IS{T}}()
  chk(C.ISDifference(i.p, j.p, is_c))
  return IS{T}(is_c[])
end

function Base.extrema(i::IS) where {T}
  min = Ref{PetscInt}()
  max = Ref{PetscInt}()
  chk(C.ISGetMinMax(i.p, min, max))
  return (Int(min[])+1, Int(max[])+1)
end
Base.minimum(i::IS) = extrema(i)[1]
Base.maximum(i::IS) = extrema(i)[2]

function Base.convert(::Type{Vector{T}}, idx::IS) where {T<:Integer}
  pref = Ref{Ptr{PetscInt}}()
  chk(C.ISGetIndices(idx.p, pref))
  inds = Int[i+1 for i in unsafe_wrap(Array, pref[], lengthlocal(idx))]
  chk(C.ISRestoreIndices(idx.p, pref))
  return inds
end
Base.Set(i::IS) = Set(Vector{Int}(i))
Base.Vector{T}(i::IS) where T<:Integer = convert(Vector{T},i)
export set_blocksize, get_blocksize

function set_blocksize(is::IS, bs::Integer)
  chk(C.ISSetBlockSize(is.p, PetscInt(bs)))
end

function get_blocksize(is::IS)
  bs = Ref{PetscInt}()
  chk(C.ISGetBlockSize(is.p, bs))
  return Int(bs[])
end

function petscview(is::IS{T}) where {T}
  viewer = C.PetscViewer{T}(C_NULL)
  chk(C.ISView(is.p, viewer))
end

###############################################################################
# we expose a 1 based API, but internally ISLoalToGlobalMappings are zero based

export ISLocalToGlobalMapping

mutable struct ISLocalToGlobalMapping{T}
  p::C.ISLocalToGlobalMapping{T}
  function ISLocalToGlobalMapping{T}(p::C.ISLocalToGlobalMapping{T}) where {T}
    o = new{T}(p)
    if (!deactivate_finalizers)
      finalizer(PetscDestroy,o)
    end
    return o
  end
end

# zero based, not exported
function _ISLocalToGlobalMapping(::Type{T}, indices::AbstractArray{PetscInt}, bs=1; comm=MPI.COMM_WORLD, copymode=C.PETSC_COPY_VALUES) where {T}

  isltog = Ref{C.ISLocalToGlobalMapping{T}}()
  chk(C.ISLocalToGlobalMappingCreate(comm, PetscInt(bs), PetscInt(length(indices)), indices, copymode, isltog))

  return ISLocalToGlobalMapping{T}(isltog[])
end

# one based, exported
#TODO: add a data argument to ISLocalToGlobalMapping, to store intermediate
# array for copymode = don't copy
function ISLocalToGlobalMapping(::Type{T}, indices::AbstractArray{I}, bs=1; comm=MPI.COMM_WORLD, copymode=C.PETSC_COPY_VALUES) where {T, I <: Integer}

  indices_0 = PetscInt[ i-1 for i in indices]
  return _ISLocalToGlobalMapping(T, indices_0, bs; comm=comm, copymode=copymode)

end

function ISLocalToGlobalMapping(is::IS{T})  where {T}

  isltog = Ref{C.ISLocalToGlobalMapping{T}}()
  chk(C.ISLocalToGlobalMappingCreateIS(is.p, isltog))
  return ISLocalToGlobalMapping{T}(isltog[])
end

comm(a::ISLocalToGlobalMapping{T}) where {T} = MPI.Comm(C.PetscObjectComm(T, a.p.pobj))

function PetscDestroy(o::ISLocalToGlobalMapping{T}) where {T}
  PetscFinalized(T) || C.ISLocalToGlobalMappingDestroy(Ref(o.p))
end

function petscview(is::ISLocalToGlobalMapping{T}) where {T}
  viewer = C.PetscViewer{T}(C_NULL)
  chk(C.ISLocalToGlobalMappingView(is.p, viewer))
end


text:
function get_unique_yes_questions_in_group(group)
    unique_questions = Set{Char}([])
    for line in split(group, "\n")
        new_questions = Set{Char}(collect(strip(line)))
        union!(unique_questions, new_questions)
    end
    unique_questions
end

function get_shared_yes_questions_in_group(group)
    shared_questions = nothing
    for line in split(group, "\n")
        if shared_questions === nothing
            shared_questions = Set{Char}(collect(strip(line)))
        else
            new_questions = Set{Char}(collect(strip(line)))
            intersect!(shared_questions, new_questions)
        end
    end
    shared_questions
end

function main()
    input_file = "day6_input.txt"
    total_unique_questions = open("day6_input.txt") do file
        file_string = read(file, String)
        groups = split(file_string, "\r\n\r\n")
        total = 0
        for group in groups
            new_questions = get_shared_yes_questions_in_group(group)
            total += length(new_questions)
        end
        total
    end
end

main()


text:
#! /usr/bin/julia

# Rosetta Code, Program name

prog = basename(Base.source_path())
println("This program file is \"", prog, "\".")

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using JointPosteriors

struct BinaryClassification{T} <: parameters{T,1}
  x::Vector{T}
  p::ProbabilityVector{3,T}
end

function log_density(Θ::BinaryClassification, data::Data)

  log_π = data.αm_m1 * log(Θ.p[3]) + data.βm_m1 * log(1 - Θ.p[3]) + data.αp_m1*log(Θ.p[2]) + data.βp_m1 * log(1 - Θ.p[2]) + data.ατ_m1 * log(Θ.p[1]) + data.βτ_m1 * log(1 - Θ.p[1])

  for i ∈ eachindex(data.X)
    log_π += data.freq[i] * log( Θ.p[1] * (1 - Θ.p[3])^data.X[i] * Θ.p[3]^(data.NmX[i]) + (1-Θ.p[1]) * Θ.p[2]^data.X[i] * (1-Θ.p[2])^data.NmX[i] )
  end

  log_π

end

struct BinaryClassificationData <: Data
  X::Array{Int64,1}
  freq::Array{Int64,1}
  NmX::Array{Int64, 1}
  N::Int64
  αm_m1::Float64
  βm_m1::Float64
  αp_m1::Float64
  βp_m1::Float64
  ατ_m1::Float64
  βτ_m1::Float64

end

function BinaryClassificationData(x::Array{Int64, 1}, Freq::Array{Int64,1}; αm::Real = 1, βm::Real = 1, αp::Real = 1, βp::Real = 1, ατ::Real = 1, βτ::Real = 1)
  N = x[end]
  non_zero_freq = sum( Freq .!= 0 )
  X = zeros( Int64, non_zero_freq )
  freq = zeros( Int64, non_zero_freq )
  j = 1
  for i ∈ 1:non_zero_freq
    while Freq[j] == 0
      j += 1
    end
    X[i] = x[j]
    freq[i] = Freq[j]
    j += 1
  end
  BinaryClassificationData(X, freq, N .- X, N, αm - 1, βm - 1, αp - 1, βp - 1, ατ - 1, βτ - 1)
end

X = [i for i ∈ 0:9];
freq = [10, 2, 2, 1, 2, 0, 0, 3, 2, 16];
data = BinaryClassificationData(X, freq, βm = 2, βp = 2);

bc_model = Model(BinaryClassification)

jp = JointPosterior(bc_model, data);

τ(Θ::BinaryClassification) = Θ.p[1]
θ_plus(Θ::BinaryClassification) = Θ.p[2]
θ_minus(Θ::BinaryClassification) = Θ.p[3]

marginal_τ = marginal(jp, τ)

marginal_θ_plus = marginal(jp, θ_plus)
marginal_θ_minus = marginal(jp, θ_minus)

f(x) = x.p[1]
d = Distributions.Normal
q = SparseQuadratureGrids.GenzKeister

values_positive = Array{Float64}(length(jp.weights_positive));
values_negative = Array{Float64}(length(jp.weights_negative));
for i ∈ eachindex(values_positive)
  values_positive[i] = f(JointPosteriors.transform!(jp.M.Θ, jp.M.Grid.nodes_positive[:,i], jp.Θ_hat, jp.L, q))
end
for i ∈ eachindex(values_negative)
  values_negative[i] = f(JointPosteriors.transform!(jp.M.Θ, jp.M.Grid.nodes_negative[:,i], jp.Θ_hat, jp.L, q))
end
μ = dot(jp.weights_positive, values_positive) + dot(jp.weights_negative, values_negative)
σ = sqrt(dot(jp.weights_positive, values_positive .^ 2) + dot(jp.weights_negative, values_negative .^ 2) - μ^2)

poly = JointPosteriors.polynomial_interpolation(jp.weights_positive, jp.weights_negative, values_positive, values_negative, d)

cdf_error(zeros(6), poly)

function prec(n::Int, ρ::Real)
  p = zeros(n,n)
  np1 = n+1
  np1h2 = (np1/2)^2
  β2 = 1 / ( np1h2 - n )
  β1 = -np1*β2
  β0 = 2 + np1h2*β2
  p[1] = p[end] = 3
  det_im1 = 3
  det_im2 = 1
  a0 = 3
  for i ∈ 2:n
    p[i,i] = a = (β0 + β1*i + β2*i^2)
    b = a*a0*ρ/4
    p[i,i-1] = p[i-1,i] = -sqrt(b)
    det_i = a*det_im1-b*det_im2
    det_im2 = det_im1
    det_im1 = det_i
    a0 = a
  end
  p, det_im1
end
function cet(δ::Vector, n::Int, ρ::Real)
  np1 = n+1
  np1h2 = (np1/2)^2
  β2 = 1 / ( np1h2 - n )
  β1 = -np1*β2
  β0 = 2 + np1h2*β2
  det_im1 = 3.0
  det_im2 = 1.0
  out = 3.0 * δ[1]^2
  a0 = 3.0
  for i ∈ 2:n
    a = (β0 + β1*i + β2*i^2)
    b = a*a0*ρ/4
    out += δ[i]^2*a - 2δ[i-1]*δ[i]*√b

    det_i = a*det_im1-b*det_im2
    det_im2 = det_im1
    det_im1 = det_i
    a0 = a
  end
  out, det_im1
end

itp = JointPosteriors.InterpolateIntegral(jp.weights_positive, jp.weights_negative, values_positive, values_negative, d)

using JointPosteriors
struct HiWorld{T} <: parameters{T,1}
  x::Vector{T}
  β::RealVector{3,T}
  σ::PositiveVector{1,T}
end

struct HiWorldData <: Data
  X::Array{Float64,2}
  y::Vector{Float64}
end

function log_density(Θ::HiWorld, data::Data)
  lpdf_normal(Θ.β, 0, 10) + lpdf_normal(Θ.σ[1], 0, 1) + lpdf_normal(data.y, data.X * Θ.β, Θ.σ[1])
end

# True parameter values
sigma = 1;
beta = [1, 1, 2.5];

# Size of dataset
size = 100;

# Predictor variable
X = hcat(ones(size), randn(size,2));
X[:,3] .*= 0.2;

# Simulate outcome variable
y = X*beta .+ randn(size) .* sigma;

HW_data = HiWorldData(X, y);
HW_model = Model(HiWorld);
hw_jp = JointPosterior(HW_model, HW_data);

marginal(hw_jp, x -> x.β[1])
marginal(hw_jp, x -> x.β[2])
marginal(hw_jp, x -> x.β[3])
marginal(hw_jp, x -> x.σ[1])

const hw_stan = "data {
  int N;
  vector[N] y;
  matrix[N, 2] x;
}
parameters {
  real alpha;
  vector[2] beta;
  real<lower=0> sigma;
}
model {
  alpha ~ normal(0, 10);
  beta ~ normal(0, 10);
  sigma ~ normal(0, 1);
  y ~ normal(alpha + x * beta, sigma);
}"
using Stan, Mamba
set_cmdstan_home!("/mnt/ssd/cmdstan-2.14.0")
#set_cmdstan_home!("/home/christel/Downloads/cmdstan-2.14.0")

hw_stan_data = Dict("N" => length(y), "x" => X[:,2:3], "y" => y)
hw_stan_model = Stanmodel(Sample(), name = "HelloWorld", model = hw_stan, monitors = ["alpha", "beta.1", "beta.2", "sigma"]);
hw_stan_res = stan(hw_stan_model, [hw_stan_data])
describe(hw_stan_res[2])

function run_hw()
  hw_jp = JointPosterior(HW_model, HW_data);
  marginal(hw_jp, x -> x.β[1])
  marginal(hw_jp, x -> x.β[2])
  marginal(hw_jp, x -> x.β[3])
  marginal(hw_jp, x -> x.σ[1])
end
using BenchmarkTools
@benchmark run_hw()

using JointPosteriors

include("/mnt/ssd/MSA/Univariate/dataGen.jl")

tv = true_values_uni(15, √99, √0.6, √0.3, √.1);
y, yp, yo = gen_data(tv, 40, 12, 12);

d = LogDensities.TF_RE_ANOVA_Data(y, yp, yo);
anova = Model(LogDensities.TF_RE_ANOVA);

jp = JointPosterior(anova, d);

function rGT(Θ::LogDensities.TF_RE_ANOVA)
  σg2 = sum(Θ.σ[2:end])
  sqrt(σg2 / (σg2 + Θ.σ[1]))
end

marginal_rGT = marginal(jp, rGT)

function run_anova()
  jp = JointPosterior(anova, d)
  marginal_rGT = marginal(jp, rGT)
end
73, 67, 75, 81) seconds, 4.9 minutes total
Sampling took (70, 103, 69, 72)
const tfanova = "
data {
  int Np1;
  int p;
  int o;
  int POp1;
  matrix[p, o] cat_count;
  matrix[p, o] y_bars;
  real ns2h;
  real cauchy_spread;
}
transformed data {
  real n;
  real po;
  matrix[p, o] invRootCounts;
  invRootCounts = inv_sqrt(cat_count);
  n = sum(cat_count);
  po = p * o;
}
parameters {
  real mu_G;
  vector[p] theta_p;
  vector[o] theta_o;
  matrix[p,o] theta_po;
  real<lower = 0> sigma_p;
  real<lower = 0> sigma_o;
  real<lower = 0> sigma_po;
  real<lower = 0> sigma_r;
}
model {
  matrix[p,o] quad;
  for (i in 1:p){
    theta_po[i, ] ~ normal(0, sigma_po);
    for (j in 1:o){
      y_bars[i,j] ~ normal( mu_G + theta_p[i] + theta_o[j] + theta_po[i, j], invRootCounts[i,j] * sigma_r);
    }
  }
  mu_G ~ normal(0, 50);
  theta_p ~ normal(0, sigma_p);
  theta_o ~ normal(0, sigma_o);
  sigma_o ~ cauchy(0, cauchy_spread);
  target += ns2h / sigma_r^2 - (n - po) * log(sigma_r) ;
}
generated quantities {
  real sigma_g2;
  real sigma_g;
  real sigma_t;
  real rGT;
  sigma_g2 = sigma_o^2 + sigma_po^2 + sigma_r^2;
  sigma_g = sqrt(sigma_g2);
  sigma_t = sqrt(sigma_g2 + sigma_p^2);
  rGT = sigma_g / sigma_t;
}
"

using Stan, Mamba#, Gadfly
set_cmdstan_home!("/mnt/ssd/cmdstan-2.14.0")
anova_stan = Stanmodel(Sample(), name = "ANOVA", model = tfanova, monitors = ["sigma_g", "sigma_t", "rGT"]);

function StanDataANOVA(data::TF_RE_ANOVA_Data)
  Dict( "Np1" => data.N + 1, "p" => data.P, "Pp1" => data.P + 1, "o" => data.O, "POp1" => data.PO + 1, "cat_count" => fill(data.R, data.P, data.O), "y_bars" => data.δ .+ data.μ_hat, "ns2h" =>  - data.s2 / 2.0)
end

res = stan(anova_stan, [StanDataANOVA(data)])
describe(res[2])


text:
using PortHamiltonianBenchmarkSystems, MAT

for n_cells in [50, 500, 10_000]
  J, R, Q, B = gugercin_pH_msd_chain(n_cells = n_cells)
  filename = "GugercinN$(2n_cells)"
  writeHDF5(filename*".h5"; J, R, Q, B)
  writeMAT(filename*".mat"; J, R, Q, B)
end


text:
module ParticipationData

using Cascadia
using DataFrames
using Gumbo
using SQLite
using XLSX

#=
Common contribution and comment structure

Contribution:
id                    internal, primary key
number                id from data source
author
content
ref
created_at
modified_at
pos_ratings
neg_ratings
assessment_decision
assessment_content

Comment:
id              internal, primary key
number          id from data source
author
content
reply_to
created_at
pos_ratings
neg_ratings
=#

export liqd,
        mwp

include("liqd.jl")
include("mwp.jl")

end # module ParticipationData


text: "foo" && true

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
mutable struct Concept
    id::Symbol
    defn::String
    predicate::Function
    interest::Float64
    objects_to_try::Array
    example_found::Array
    unused_predicates::Array
    predicates_used_in_spec::Array
    predicates_used_in_desc::Array
    number_found::Int64
    number_tried::Int64
    parent::Union{Concept, Missing}
    subconcepts::Array{Concept}
    function Concept(W::World,Desc,pred,intr,spec::Array,desc::Array,par,subc)
        new(gensym(), # Use Julia's build in name generator
            Desc,
            pred,
            intr,
            copy(W.universe),
            [],
            copy(W.predicates),
            spec,
            desc,
            0,
            0,
            par,
            subc)
    end
end

function Base.show(io::IO,  P::Concept)
    print(io, P.defn)
end

"""
    make_specialization(C::Concept,W::World)

    ; A task of type MAKE-SPECIALIZATION requires that the system
    ;  attempt to create a representation for a new concept.
    ; To create a new concept, the function MAKE-SPECIALIZATION
    ;  creates a structure containing  the following:
    ; - A unique id which is a symbol, e.g., C1, C2, etc.
    ; - A definition of the concept in terms of the parent concept.
    ;   suitable for an explanatory printout.
    ; - A predicate that can be applied to any object
    ;   to determine whether it is an example of this concept.
    ; - An interest value for the concept computed using a rule
    ;   which takes into account the interest of the parent concept
    ;   and the interest of the predicate used to form the restriction.
    ; - A list of objects that have not yet been tried as possible examples,
    ;   initially the whole "UNIVERSE".
    ; - A list of examples found.
    ; - A list of the predicates NOT used in the definition of this concept -
    ;   this simplifies the procedure MAKE-SPECIALIZATION.
    ; - A list of the original (provided) predicates that have been used
    ;   along the path from OBJECT to this concept (used in DISPLAY-CONCEPTS).
    ; - A list of the predicates that have been used to create specializations
    ;   of this concept.
    ; - The number of examples found (so far).
    ; - The number of objects tried (so far).
    ; - A parent concept.
    ; - A list of subconcepts.
    ; These items are put on the property list of the atom under the
    ; types: ID, DEFN, PREDICATE, INTEREST, OBJECTS-TO-TRY, EXAMPLES-FOUND,
    ;	UNUSED-PREDICATES, PREDICATES-USED-IN-SPEC,
    ;	PREDICATES-USED-IN-DESC, NUMBER-FOUND, NUMBER-TRIED,
    ;       PARENT, SUBCONCEPTS.
"""
function make_specialization(C::Concept,W::World)
    unused_pred = setdiff(C.unused_predicates,C.predicates_used_in_spec)
    # Select a predicate not already involved in the parent
    # and not already used for a specialization of C.

    if isempty(unused_pred)
        print("Cannot further specialize concept $C")
        return missing ##TODO maybe false? not sure
    end
    pred = first(unused_pred)

    #Indicate that the selected predicate is no longer available
    push!(C.predicates_used_in_spec,pred)

    #Create a new concept structure
    newc = Concept(
        W,
        #Formulate the definition
        "$(C.defn) with $pred",
        #Create the predicate which tests an object to see if
        #it is an example of the new concept:
        (x)->(C.predicate(x) & pred(x)),
        #interest = interest of parent concept.
        C.interest,
        #inherit specializations used from parent
        #not so in original code
        copy(C.predicates_used_in_spec),
        push!(copy(C.predicates_used_in_desc),pred),
        #parent is C
        C,
        #no childreen yet
        []
        )
        push!(C.subconcepts,newc)

        println("Created $(newc.id): $(newc.defn)")
    if !(haskey(W.agenda,(find_examples_of,newc)))
        enqueue!(W.agenda,(find_examples_of,newc)=>C.interest)
    end
    if !(haskey(W.agenda,(make_specialization,C)))
        enqueue!(W.agenda,(make_specialization,C)=>spec_task_interest(C))
    end

    newc.unused_predicates=filter((x)->(x!=pred),C.unused_predicates)
end

"""
    find_examples_of(C::Concept)

    Tests objects for being examples of concept C, and updates
    the agenda according to the findings

;;; To find examples, the procedure FIND-EXAMPLES-OF takes the list of
;;; objects not yet tried and tries a fixed number of them, (3 of them).
;;; It puts any examples found on the list of examples for the concept,
;;; and it updates the list of objects left to try.
;;; This procedure also does the following:
;;; It updates the interest value for the concept in accordance to
;;; the results of looking for examples.
"""
function find_examples_of(C::Concept,W::World)
    n = min(3,length(C.objects_to_try))
    for i=1:1:n
        C.number_tried += 1;
        x, C.objects_to_try = C.objects_to_try[1],C.objects_to_try[2:end]

        if (C.predicate(x))
            println("$x is an example of $(C.id)")
            C.number_found += 1;
            push!(C.example_found,x)
        else
            println("$x is not an example of $(C.id)")
        end
    end

    # The example-checking part of this task is over.
    # Now update the interest value for C:

    C.interest = compute_concept_interest(C)

    # If there are still objects not yet tried, enter a new
    # task on the agenda to try 3 more objects.

    if !(isempty(C.objects_to_try)) & !(haskey(W.agenda,(find_examples_of,C)))
        enqueue!(W.agenda,(find_examples_of,C)=>examples_task_interest(C))
    else
        if (W.reporting)
            print("All objects are tested for $(C.id)")
        end
    end
    # If there is at least one example of the concept and no
    # specializations for this concept have yet been created,
    # and no tasks for such specialization are already on the
    # agenda, create a new task to make a specialization of C:

    if (C.number_found>0 &
        isempty(C.subconcepts) )

        if (W.reporting)
            println((make_specialization,C))
        end
        if !(haskey(W.agenda,(make_specialization,C)))
            enqueue!(W.agenda,(make_specialization,C)=>spec_task_interest(C))
        end
    end
end

"""
    compute_concept_interest(C)

    ;;; CONCEPT-INTEREST computes the current interest value
    ;;; for concept C using a formula that involves the hit ratio.
"""
function compute_concept_interest(C)
    if C.number_tried == 0
        if (ismissing(C.parent))
            error("Can't compute interest for $C")
        else
            return C.parent.interest
        end
    else
        r = C.number_found / C.number_tried
        return 400.0 * (r-r^2)
    end
end

"""
    examples_task_interest(C::Concept)

    ;;; EXAMPLES-TASK-INTEREST computes the interest of a task to
    ;;; find examples of C as a weighted sum of the interests of
    ;;; C and it parent.
"""
function examples_task_interest(C::Concept)
    if (ismissing(C.parent))
        return 50
    else
        return C.parent.interest*0.8+0.2*C.interest
    end
end

"""
    spec_task_interest(C::Concept)

    ;;; SPEC-TASK-INTEREST computes the interest value
    ;;; for a specialization task, according to the formula:
    ;;;  value = 10 times the parent's hit ratio.
    ;;; We add one to the denominator to avoid the possibility
    ;;; of division by zero.
"""
function spec_task_interest(C::Concept)
    return (10*C.number_found)/(1 + C.number_tried)
end


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

name = "Singular"
version = v"4.1.3"

# Collection of sources required to build normaliz
sources = [
    GitSource("https://github.com/Singular/Sources.git", "4c1fc06f2d81e12a8e75a5419444cb157fbe45e9"),
]

# Bash recipe for building across all platforms
script = raw"""
cd Sources
./autogen.sh
./configure --prefix=$prefix --host=$target --build=${MACHTYPE} \
    --with-libparse \
    --enable-shared \
    --disable-static \
    --enable-p-procs-static \
    --disable-p-procs-dynamic \
    --disable-gfanlib \
    --with-readline=no \
    --with-gmp=$prefix \
    --with-flint=$prefix \
    --without-python

#    --with-ntl=$prefix

make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = filter!(p -> !(p isa Windows), platforms)
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libpolys", :libpolys),
    LibraryProduct("libSingular", :libsingular),
    # LibraryProduct("customstd", :customstd),
    # LibraryProduct("subsets", :subsets),
    ExecutableProduct("Singular", :Singular),
    ExecutableProduct("libparse", :libparse),
    # LibraryProduct("syzextra", :syzextra),
    # LibraryProduct("interval", :interval),
    LibraryProduct("libfactory", :libfactory),
    LibraryProduct("libsingular_resources", :libsingular_resources),
    LibraryProduct("libomalloc", :libomalloc)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("FLINT_jll"),
    Dependency("GMP_jll"),
    Dependency("MPFR_jll"),
    #Dependency("ntl_jll"),
]

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


text:
using PkgBenchmark
using DiffEqBase: DiscreteProblem
using LyapunovExponents: TangentRenormalizer, MLERenormalizer,
    post_evolve!, get_integrator, objname

""" Fake LESolver for benchmarking `post_evolve!`. """
function noop_tangent_integrator(N, dim_lyap)
    tangent_prob = DiscreteProblem((x...) -> nothing,  # irrelevant
                                   zeros(N, 1+ dim_lyap),
                                   (0, 1))  # irrelevant
    return get_integrator(tangent_prob)
end

suite = BenchmarkGroup()
for renormalizer in [TangentRenormalizer, MLERenormalizer]
    bg = suite[objname(renormalizer)] = BenchmarkGroup()
    for N in 2.^(6:7), dim_lyap in (renormalizer === TangentRenormalizer ?
                                    floor.(Int, [0.1, 0.5, 0.9, 1] .* N) :
                                    [1])
        stage = renormalizer(noop_tangent_integrator(N, dim_lyap),
                             1, 1)
        bg["N=$N", "dim_lyap=$dim_lyap"] =
            @benchmarkable(post_evolve!($stage),
                           setup = (randn!($stage.tangent_state)))
    end
end

suite


text:
# set power limit
# SP
export psp_setplimit

"""
    psp_setplimit(io, power_limit)

Set the power limit.
"""
function psp_setplimit(io_psp::IO,plim::Real)
  @assert plim<1000 "power limit must be < 1000W"
  plimstring = @sprintf("%03d",plim)
  set_power_limit_string = "SP "*plimstring*"\r"
  write(io_psp, set_power_limit_string)
  return nothing
end


text:
module YARS

export yars_start, yars_send_reset, yars_send_quit, yars_send_message
export yars_send_actuator_commands, yars_read_sensors

export yars_read_float,  yars_send_float
export yars_read_string, yars_send_string
export yars_read_int,    yars_send_int

export yars_read_float_vector, yars_send_float_vector
export yars_read_int_vector,   yars_send_int_vector

function print_yars(st)
  while true
    print(readline(st))
  end
end

###########################################################################
#                        basic data types - float                         #
###########################################################################

function yars_send_float_vector(hd, mv::Vector{Float64})
  c = bytes_of_size_int(length(mv))
  s = vcat(convert(Vector{Uint8}, ['D']), convert(Vector{Uint8}, c))
  write(hd,s)
  for f in mv
    write(hd, f)
  end
end
send_float_vector(hd, mv) = yars_send_float_vector(hd, mv)

function yars_read_float_vector(hd)
  s = []
  s = readbytes(hd, 1)
  s = convert(Char, s[1])
  if s != 'D'
    println("expected 'D' but received '", s, "'")
  end
  nr = read(hd,Int32)
  data = zeros(nr)
  for i=1:nr
    data[i] = read(hd, Float64)
  end
  data
end
read_float_vector(hd) = yars_read_float_vector(hd)

function yars_send_float(hd, value::Float64)
  write(hd, value)
end
send_float(hd, value) = yars_send_float(hd, value)

function yars_read_float(hd)
  read(hd, Float64)
end
read_float(hd) = yars_read_float(hd)

###########################################################################
#                         basic data types - int                          #
###########################################################################

function yars_send_int_vector(hd, mv::Vector{Int64})
  c = bytes_of_size_int(length(mv))
  s = vcat(convert(Vector{Uint8}, ['D']), convert(Vector{Uint8}, c))
  write(hd,s)
  for f in mv
    write(hd, f)
  end
end
send_int_vector(hd, mv) = yars_send_int_vector(hd, mv)

function yars_read_int_vector(hd)
  s = []
  s = readbytes(hd, 1)
  s = convert(Char, s[1])
  if s != 'D'
    println("expected 'D' but received '", s, "'")
  end
  nr = read(hd,Int32)
  data = zeros(nr)
  for i=1:nr
    data[i] = read(hd, Int64)
  end
  data
end
read_int_vector(hd) = yars_read_int_vector(hd)

function yars_send_int(hd, value::Int64)
  write(hd, value)
end
send_int(hd, value) = yars_send_int(hd, value)

function yars_read_int(hd)
  read(hd, Int64)
end
read_int(hd) = yars_read_int(hd)

###########################################################################
#                         communication functions                         #
###########################################################################

function yars_start(working_dir::String, options::Vector{ASCIIString})
  hd = nothing
  pr = nothing
  opened = false
  while opened == false
    try
      current_dir = pwd()
      cd("$working_dir")
      st = []
      pr = []
      (st, pr) = open(`yars $options`,"r")
      s = "ABC"
      while ismatch(r"port",s) == false && length(s) > 0
        s = readline(st)
      end
      if length(s) == 0
        println("error opening yars")
        exit(-1)
      end
      s    = split(s)
      port = int(s[end])
      cd("$current_dir")
      sleep(1)
      hd = connect(port)
      sleep(1)
      r = @spawn print_yars(st)
      opened = true
    catch
      opened = false
    end
  end
  return hd, pr
end
start(working_dir::String, options::Vector{ASCIIString}) = yars_start(working_dir, options)

function bytes_of_size_int(i::Int64)
  [(i >> (j * 8)) & 0xFF for j=0:3]
end

function yars_send_string(hd, str::String)
  c = bytes_of_size_int(length(str))
  s = vcat(convert(Vector{Uint8}, ['s']),
           convert(Vector{Uint8}, c),
           convert(Vector{Uint8}, str))
  write(hd, s)
end
send_string(hd, str::String) = yars_send_string(hd, str)

function yars_send_reset(hd)
  yars_send_string(hd, "RESET")
end
reset(hd) = yars_send_reset(hd)

function yars_send_quit(hd)
  yars_send_string(hd, "QUIT")
end
quit(hd) = yars_send_quit(hd)

function yars_send_message(hd, str)
  yars_send_string(hd, "MESSAGE")
  yars_send_string(hd, str)
end
message(hd, str) = yars_send_message(hd, str)

function yars_read_sensors(hd)
  yars_send_string(hd, "SENSORS")
  yars_read_float_vector(hd)
end
sensors(hd) = yars_read_sensors(hd)

function yars_send_actuator_commands(hd, mv::Vector{Float64})
  yars_send_string(hd, "ACTUATORS")
  yars_send_float_vector(hd, mv)
end
actuators(hd, mv) = yars_send_actuator_commands(hd, mv)

function yars_read_string(hd)
  s = readbytes(hd, 1)
  s = convert(Char, s[1])
  if s != 's'
    println("expected 's' but received '", s, "'")
  end
  nr = read(hd,Int32)
  r = ""
  for i = 1:nr
    r = string(r, read(hd, Char))
  end
  return r
end
read_string(hd) = yars_read_string(hd)

type Entity
  name::String
  dimension::Int64
  internal::Vector{(Float64, Float64)}
  external::Vector{(Float64, Float64)}
end

type Robot
  sensors::Vector{Entity}
  actuators::Vector{Entity}
  sensor_dim::Int64
  actuator_dim::Int64
end

function parse_configuration_strings(strings::Array{ASCIIString,1})

  sensors   = []
  actuators = []

  current   = nothing

  for str in strings
    if contains(str, "BEGIN SENSOR") || contains(str, "BEGIN ACTUATOR")
      current = Entity("", 0, [(0.0, 0.0)], [(0.0, 0.0)])
    end

    if contains(str, "END SENSOR")
      sensors = [sensors, current]
      current = nothing
    end

    if contains(str, "END ACTUATOR")
      actuators = [actuators, current]
      current   = nothing
    end

    if contains(str,"NAME") && current != nothing
      current.name = strip(str[5:end])
    end

    if contains(str,"DIMENSION") && current != nothing
      values            = split(str)
      current.dimension = int64(values[2])
    end

    if contains(str,"INTERNAL DOMAIN")
      values           = split(str)
      v1               = float(values[end-1])
      v2               = float(values[end])
      current.internal = [current.internal, (v1, v2)]
    end

    if contains(str,"EXTERNAL DOMAIN")
      values           = split(str)
      v1               = float(values[end-1])
      v2               = float(values[end])
      current.external = [current.external, (v1, v2)]
    end
  end

  s = 0
  for sen in sensors
    s = s + sen.dimension
  end

  a = 0
  for act in actuators
    a = a + act.dimension
  end

  return Robot(sensors, actuators, s, a)
end

function yars_get_configuration(hd)
  configuration_strings = []
  yars_send_string(hd, "CONFIGURATION")
  s = ""
  while s != "END CONFIGURATION"
    s = yars_read_string(hd)
    configuration_strings = [configuration_strings, s]
  end
  return parse_configuration_strings(configuration_strings)
end
configuration(hd) = yars_get_configuration(hd)

end # module


text:
module RocketUppercaseOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: uppercase()" begin

    println("Testing: operator uppercase()")

    run_proxyshowcheck("Uppercase", uppercase())

    run_testset([
        (
            source = from("Hello, world") |> uppercase(),
            values = @ts(['H', 'E', 'L', 'L', 'O', ',', ' ', 'W', 'O', 'R', 'L', 'D', c])
        ),
        (
            source = completed() |> uppercase(),
            values = @ts(c)
        ),
        (
            source      = faulted(String, "e") |> uppercase(),
            values      = @ts(e("e")),
            source_type = String
        ),
        (
            source = never() |> uppercase(),
            values = @ts()
        )
    ])

end

end


text:
using Gen
using Gen2DAgentMotion
using Test
import Random

include("motion.jl")