text:
# This file is a part of Julia. Except for the rem_pio2_kernel, and
# cody_waite_* methods (see below) license is MIT: https://julialang.org/license

# rem_pio2_kernel and cody_waite_* methods are heavily based on FDLIBM code:
# __ieee754_rem_pio2, that is made available under the following licence:

## Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
##
## Developed at SunPro, a Sun Microsystems, Inc. business.
## Permission to use, copy, modify, and distribute this
## software is freely granted, provided that this notice
## is preserved.

# Bits of 1/2π
#   1/2π == sum(x / 0x1p64^i for i,x = enumerate(INV_2PI))
# Can be obtained by:
#
#    setprecision(BigFloat, 4096)
#    I = 0.5/big(pi)
#    for i = 1:19
#        I *= 0x1p64
#        k = trunc(UInt64, I)
#        @printf "0x%016x,\n" k
#        I -= k
#    end
const INV_2PI = UInt64[
    0x28be_60db_9391_054a,
    0x7f09_d5f4_7d4d_3770,
    0x36d8_a566_4f10_e410,
    0x7f94_58ea_f7ae_f158,
    0x6dc9_1b8e_9093_74b8,
    0x0192_4bba_8274_6487,
    0x3f87_7ac7_2c4a_69cf,
    0xba20_8d7d_4bae_d121,
    0x3a67_1c09_ad17_df90,
    0x4e64_758e_60d4_ce7d,
    0x2721_17e2_ef7e_4a0e,
    0xc7fe_25ff_f781_6603,
    0xfbcb_c462_d682_9b47,
    0xdb4d_9fb3_c9f2_c26d,
    0xd3d1_8fd9_a797_fa8b,
    0x5d49_eeb1_faf9_7c5e,
    0xcf41_ce7d_e294_a4ba,
    0x9afe_d7ec_47e3_5742,
    0x1580_cc11_bf1e_daea]

@inline function cody_waite_2c_pio2(x::Float64, fn, n)
    pio2_1 = 1.57079632673412561417e+00
    pio2_1t = 6.07710050650619224932e-11

    z = muladd(-fn, pio2_1, x) # x - fn*pio2_1
    y1 = muladd(-fn, pio2_1t, z) # z - fn*pio2_1t
    y2 = muladd(-fn, pio2_1t, (z - y1)) # (z - y1) - fn*pio2_1t
    n, DoubleFloat64(y1, y2)
end

@inline function cody_waite_ext_pio2(x::Float64, xhp)
    pio2_1 = 1.57079632673412561417e+00
    pio2_1t = 6.07710050650619224932e-11
    pio2_2 = 6.07710050630396597660e-11
    pio2_2t = 2.02226624879595063154e-21
    pio2_3 = 2.02226624871116645580e-21
    pio2_3t = 8.47842766036889956997e-32

    fn = round(x*(2/pi)) # round to integer
    # on older systems, the above could be faster with
    # rf = 1.5/eps(Float64)
    # fn = (x*(2/pi)+rf)-rf

    r  = muladd(-fn, pio2_1, x) # x - fn*pio2_1
    w  = fn*pio2_1t # 1st round good to 85 bit
    j  = xhp>>20
    y1 = r-w
    high = highword(y1)
    i = j-((high>>20)&0x7ff)
    if i>16  # 2nd iteration needed, good to 118
        t  = r
        w  = fn*pio2_2
        r  = t-w
        w  = muladd(fn, pio2_2t,-((t-r)-w))
        y1 = r-w
        high = highword(y1)
        i = j-((high>>20)&0x7ff)
        if i>49 # 3rd iteration need, 151 bits acc
            t  = r # will cover all possible cases
            w  = fn*pio2_3
            r  = t-w
            w  = muladd(fn, pio2_3t, -((t-r)-w))
            y1 = r-w
        end
    end
    y2 = (r-y1)-w
    return unsafe_trunc(Int, fn), DoubleFloat64(y1, y2)
end

"""
    fromfraction(f::Int128)

Compute a tuple of values `(z1,z2)` such that
    ``z1 + z2 == f / 2^128``
and the significand of `z1` has 27 trailing zeros.
"""
function fromfraction(f::Int128)
    if f == 0
        return (0.0,0.0)
    end

    # 1. get leading term truncated to 26 bits
    s = ((f < 0) % UInt64) << 63     # sign bit
    x = abs(f) % UInt128             # magnitude
    n1 = 128-leading_zeros(x)         # ndigits0z(x,2)
    m1 = ((x >> (n1-26)) % UInt64) << 27
    d1 = ((n1-128+1021) % UInt64) << 52
    z1 = reinterpret(Float64, s | (d1 + m1))

    # 2. compute remaining term
    x2 = (x - (UInt128(m1) << (n1-53)))
    if x2 == 0
        return (z1, 0.0)
    end
    n2 = 128-leading_zeros(x2)
    m2 = (x2 >> (n2-53)) % UInt64
    d2 = ((n2-128+1021) % UInt64) << 52
    z2 = reinterpret(Float64,  s | (d2 + m2))
    return (z1,z2)
end

function paynehanek(x::Float64)
    # 1. Convert to form
    #
    #    x = X * 2^k,
    #
    # where 2^(n-1) <= X < 2^n  is an n-bit integer (n = 53, k = exponent(x)-52 )

    # Computations are integer based, so reinterpret x as UInt64
    u = reinterpret(UInt64, x)
    # Strip x of exponent bits and replace with ^1
    X = (u & significand_mask(Float64)) | (one(UInt64) << significand_bits(Float64))
    # Get k from formula above
    # k = exponent(x)-52
    k = Int((u & exponent_mask(Float64)) >> significand_bits(Float64)) - exponent_bias(Float64) - significand_bits(Float64)

    # 2. Let α = 1/2π, then:
    #
    #    α*x mod 1 ≡ [(α*2^k mod 1)*X] mod 1
    #
    # so we can ignore the first k bits of α. Extract the next 3 64-bit parts of α.
    #
    # i.e. equivalent to
    #     setprecision(BigFloat,4096)
    #     α  = 1/(2*big(pi))
    #     A  = mod(ldexp(α,k), 1)
    #     z1 = ldexp(A,64)
    #     a1 = trunc(UInt64, z1)
    #     z2 = ldexp(z1-a1, 64)
    #     a2 = trunc(UInt64, z2)
    #     z3 = ldexp(z2-a2, 64)
    #     a3 = trunc(UInt64, z3)

    # This is equivalent to
    #     idx, shift = divrem(k, 64)
    # but divrem is slower.
    idx = k >> 6

    shift = k - (idx << 6)
    if shift == 0
        @inbounds a1 = INV_2PI[idx+1]
        @inbounds a2 = INV_2PI[idx+2]
        @inbounds a3 = INV_2PI[idx+3]
    else
        # use shifts to extract the relevant 64 bit window
        @inbounds a1 = (idx < 0 ? zero(UInt64) : INV_2PI[idx+1] << shift) | (INV_2PI[idx+2] >> (64 - shift))
        @inbounds a2 = (INV_2PI[idx+2] << shift) | (INV_2PI[idx+3] >> (64 - shift))
        @inbounds a3 = (INV_2PI[idx+3] << shift) | (INV_2PI[idx+4] >> (64 - shift))
    end

    # 3. Perform the multiplication:
    #
    #      X.  0  0  0
    #   ×  0. a1 a2 a3
    #   ==============
    #      _.  w  w  _
    #
    # (i.e. ignoring integer and lowest bit parts of result)

    w1 = UInt128(X*a1) << 64 # overflow becomes integer
    w2 = widemul(X,a2)
    w3 = widemul(X,a3) >> 64
    w = w1 + w2 + w3         # quotient fraction after division by 2π

    # adjust for sign of x
    w = flipsign(w,x)

    # 4. convert to quadrant, quotient fraction after division by π/2:
    q = (((w>>125)%Int +1)>>1) # nearest quadrant
    f = (w<<2) % Int128 # fraction part of quotient after division by π/2, taking values on [-0.5,0.5)

    # 5. convert quotient fraction to split precision Float64
    z_hi,z_lo = fromfraction(f)

    # 6. multiply by π/2
    pio2 = 1.5707963267948966
    pio2_hi = 1.5707963407039642
    pio2_lo = -1.3909067614167116e-8
    y_hi = (z_hi+z_lo)*pio2
    y_lo = (((z_hi*pio2_hi - y_hi) + z_hi*pio2_lo) + z_lo*pio2_hi) + z_lo*pio2_lo
    return q, DoubleFloat64(y_hi, y_lo)
end

"""
    rem_pio2_kernel(x)

Return the remainder of `x` modulo π/2 as a double-double pair, along with a `k`
such that ``k \\mod 3 == K \\mod 3`` where ``K*π/2 = x - rem``. Note, that it is
only meant for use when ``|x|>=π/4``, and that ``π/2`` is always subtracted or
added for ``π/4<|x|<=π/2`` instead of simply returning `x`.
"""
@inline function rem_pio2_kernel(x::Float64)
    xhp = poshighword(x)
    #  xhp <= highword(5pi/4) implies |x| ~<= 5pi/4,
    if xhp <= 0x400f6a7a
        #  last five bits of xhp == last five bits of highword(pi/2) or
        #  highword(2pi/2) implies |x| ~= pi/2 or 2pi/2,
        if (xhp & 0xfffff) == 0x921fb # use precise Cody Waite scheme
            return cody_waite_ext_pio2(x, xhp)
        end
        # use Cody Waite with two constants
        #  xhp <= highword(3pi/4) implies |x| ~<= 3pi/4
        if xhp <= 0x4002d97c
            if x > 0.0
                return cody_waite_2c_pio2(x, 1.0, 1)
            else
                return cody_waite_2c_pio2(x, -1.0, -1)
            end
        # 3pi/4 < |x| <= 5pi/4
        else
            if x > 0.0
                return cody_waite_2c_pio2(x, 2.0, 2)
            else
                return cody_waite_2c_pio2(x, -2.0, -2)
            end
        end
    end
    #  xhp <= highword(9pi/4) implies |x| ~<= 9pi/4
    if xhp <= 0x401c463b
        #  xhp <= highword(7pi/4) implies |x| ~<= 7pi/4
        if xhp <= 0x4015fdbc
            #  xhp == highword(3pi/2) implies |x| ~= 3pi/2
            if xhp == 0x4012d97c # use precise Cody Waite scheme
                return cody_waite_ext_pio2(x, xhp)
            end
            # use Cody Waite with two constants
            if x > 0.0
                return cody_waite_2c_pio2(x, 3.0, 3)
            else
                return cody_waite_2c_pio2(x, -3.0, -3)
            end
        # 7pi/4 < |x| =< 9pi/4
        else
            #  xhp == highword(4pi/2) implies |x| ~= 4pi/2
            if xhp == 0x401921fb # use precise Cody Waite scheme
                return cody_waite_ext_pio2(x, xhp)
            end
            # use Cody Waite with two constants
            if x > 0.0
                return cody_waite_2c_pio2(x, 4.0, 4)
            else
                return cody_waite_2c_pio2(x, -4.0, -4)
            end
        end
    end
    #  xhp < highword(2.0^20*pi/2) implies |x| ~< 2^20*pi/2
    if xhp < 0x413921fb # use precise Cody Waite scheme
        return cody_waite_ext_pio2(x, xhp)
    end
    # if |x| >= 2^20*pi/2 switch to Payne Hanek
    return paynehanek(x)
end

## Float32
@inline function rem_pio2_kernel(x::Float32)
    pio2_1 = 1.57079631090164184570e+00
    pio2_1t = 1.58932547735281966916e-08
    inv_pio2 = 6.36619772367581382433e-01
    xd = convert(Float64, x)
    absxd = abs(xd)
    # it is assumed that NaN and Infs have been checked
    if absxd <= pi*5/4
        if absxd <= pi*3/4
            if x > 0
                return 1, DoubleFloat32(xd - pi/2)
            else
                return -1, DoubleFloat32(xd + pi/2)
            end
        end
        if x > 0
            return 2, DoubleFloat32(xd - pi)
        else
            return -2, DoubleFloat32(xd + pi)
        end
    elseif absxd <= pi*9/4
        if absxd <= pi*7/4
            if x > 0
                return 3, DoubleFloat32(xd - pi*3/2)
            else
                return -3, DoubleFloat32(xd + pi*3/2)
            end
        end
        if x > 0
            return 4, DoubleFloat32(xd - pi*4/2)
        else
            return -4, DoubleFloat32(xd + pi*4/2)
        end
    end
    #/* 33+53 bit pi is good enough for medium size */
    if absxd < Float32(pi)/2*2.0f0^28 # medium size */
        # use Cody Waite reduction with two coefficients
        fn = round(xd*inv_pio2)
        r  = xd-fn*pio2_1
        w  = fn*pio2_1t
        y = r-w;
        return unsafe_trunc(Int, fn), DoubleFloat32(y)
    end
    n, y = rem_pio2_kernel(xd)
    return n, DoubleFloat32(y.hi)
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
module LHEF

using LightXML, LorentzVectors

export parse_lhe

struct EventHeader
    nup::UInt8          # Number of particles
    ldprup::UInt8       # Process type?
    xwgtup::Float64     # Event Wight
    scalup::Float64     # Scale
    αem::Float64        # AQEDUP
    αs::Float64         # AQCDUP
end

struct Particle
    particle::Int32
    status::Int8
    mothup::NTuple{2,UInt8}
    color::NTuple{2,UInt16}
    pμ::LorentzVector{Float64}
    m::Float64
    vtimup::Float64
    spinup::Float64
end

struct Event
    header::EventHeader
    data::Vector{Particle}
end

function parse_lhe(filename; format=nothing)
    if format === nothing
        # Format not declared, inferring from extension
        fparts = split(basename(filename), ".")
        format = if fparts[end] == "lhe"
            :lhe
        elseif fparts[end] == "gz" && fparts[end - 1] == "lhe"
            :lhegz
        end
    end

    lhefile = parse_file(filename)
    lhenode = root(lhefile)

    (name(lhenode) == "LesHouchesEvents") || error("Invalid root node")
    events = lhenode["event"]
    return [
        begin
            data = content(first(child_nodes(event)))
            lines = split(data, '\n'; keepempty=false)
            headerdata = split(lines[1], ' '; keepempty=false)
            header = EventHeader(
                parse(UInt8, headerdata[1]),
                parse(UInt8, headerdata[2]),
                parse(Float64, headerdata[3]),
                parse(Float64, headerdata[4]),
                parse(Float64, headerdata[5]),
                parse(Float64, headerdata[6]),
            )
            data = [
                begin
                    fields = split(line, ' '; keepempty=false)
                    p = Particle(
                        parse(Int32, fields[1]),
                        parse(Int8, fields[2]),
                        (parse(UInt8, fields[3]), parse(UInt8, fields[4])),
                        (parse(UInt16, fields[5]), parse(UInt16, fields[6])),
                        LorentzVector(
                            parse(Float64, fields[10]),
                            parse(Float64, fields[7]),
                            parse(Float64, fields[8]),
                            parse(Float64, fields[9]),
                        ),
                        parse(Float64, fields[11]),
                        parse(Float64, fields[12]),
                        parse(Float64, fields[13]),
                    )
                    p
                end for line in lines[2:end]
            ]
            Event(header, data)
        end for event in events
    ]
end

end # module

ext: jl


text:
@doc raw"""
`VDAE`: Variational Differential Algebraic Equation *EXPERIMENTAL*

Defines an implicit initial value problem
```math
\begin{align*}
\dot{q} (t) &= v(t) + λ(t), &
q(t_{0}) &= q_{0} , \\
\dot{p} (t) &= f(t, q(t), v(t)) + g(t, q(t), λ(t)) + \obar{g}(t, q(t), μ(t)) , &
p(t_{0}) &= p_{0} , \\
p(t) &= ϑ(t, q(t), v(t)) , \\
0 &= \phi (t, q(t), p(t)) , \\
0 &= \psi (t, q(t), p(t), \dot{q}(t), \dot{p}(t)) ,
\end{align*}
```
with vector field ``f``, the momentum defined by ``p``, initial conditions ``(q_{0}, p_{0})`` and the solution
``(q,p)`` taking values in ``\mathbb{R}^{d} \times \mathbb{R}^{d}`` and
the algebraic variables ``(v, \lambda, \mu)`` taking values in
``\mathbb{R}^{d} \times \mathbb{R}^{d} \times \mathbb{R}^{m}``.
This is a special case of a differential algebraic equation with dynamical
variables ``(q,p)`` and algebraic variables ``v``, ``\lambda`` and ``\mu``.

### Fields

* `d`: dimension of dynamical variables ``q`` and ``p`` as well as the vector fields ``f`` and ``p``
* `ϑ`: function determining the momentum
* `f`: function computing the vector field
* `g`: function determining the primary projection, usually given by ∇ϑ(q)⋅λ
* `g̅`: function determining the secondary projection, usually given by λ⋅∇ϑ(q)
* `ϕ`: primary constraints, usually given by p-ϑ(q)
* `ψ`: secondary constraints, usually given by ṗ-q̇⋅∇ϑ(q)
* `v`: function computing an initial guess for the velocity field (optional)
* `Ω`: symplectic matrix (optional)
* `∇H`: gradient of the Hamiltonian (optional)
* `t₀`: initial time (optional)
* `q₀`: initial condition for `q`
* `p₀`: initial condition for `p`
* `λ₀`: initial condition for `λ` (optional)
* `μ₀`: initial condition for `μ` (optional)

The functions `ϑ` and `f` must have the interface
```julia
    function ϑ(t, q, v, p)
        p[1] = ...
        p[2] = ...
        ...
    end
```
and
```julia
    function f(t, q, v, f)
        f[1] = ...
        f[2] = ...
        ...
    end
```
where `t` is the current time, `q` is the current solution vector, `v` is the
current velocity and `f` and `p` are the vectors which hold the result of
evaluating the functions ``f`` and ``ϑ`` on `t`, `q` and `v`.
The funtions `g` and `v` are specified by
```julia
    function g(t, q, λ, g)
        g[1] = ...
        g[2] = ...
        ...
    end
```
and
```julia
    function v(t, q, p, v)
        v[1] = ...
        v[2] = ...
        ...
    end
```
"""
struct VDAE{dType <: Number, tType <: Number, ϑType <: Function,
            fType <: Function, gType <: Function, g̅Type <: Function,
            ϕType <: Function, ψType <: Function, vType <: Function,
            ΩType <: Union{Function,Nothing}, ∇HType <: Union{Function,Nothing},
            pType <: Union{Tuple,Nothing}, N} <: AbstractEquationPDAE{dType, tType}

    d::Int
    m::Int
    n::Int
    ϑ::ϑType
    f::fType
    g::gType
    g̅::g̅Type
    ϕ::ϕType
    ψ::ψType
    v::vType
    Ω::ΩType
    ∇H::∇HType
    t₀::tType
    q₀::Array{dType, N}
    p₀::Array{dType, N}
    λ₀::Array{dType, N}
    μ₀::Array{dType, N}
    parameters::pType
    periodicity::Vector{dType}

    function VDAE(DT::DataType, N::Int, d::Int, m::Int, n::Int, ϑ::ϑType, f::fType,
                  g::gType, g̅::g̅Type, ϕ::ϕType, ψ::ψType, t₀::tType,
                  q₀::AbstractArray{dType}, p₀::AbstractArray{dType},
                  λ₀::AbstractArray{dType}, μ₀::AbstractArray{dType};
                  v::vType=nothing, Ω::ΩType=nothing, ∇H::∇HType=nothing,
                  parameters=nothing, periodicity=zeros(DT,d)) where {
                        dType <: Number, tType <: Number, ϑType <: Function,
                        fType <: Function, gType <: Function, g̅Type <: Function,
                        ϕType <: Function, ψType <: Function, vType <: Function,
                        ΩType <: Union{Function,Nothing}, ∇HType <: Union{Function,Nothing}}

        @assert d == size(q₀,1) == size(p₀,1) == size(λ₀,1)
        @assert m == size(μ₀,1)
        @assert n == size(q₀,2) == size(p₀,2)
        @assert ndims(q₀) == ndims(p₀) == N ∈ (1,2)

        new{DT, tType, ϑType, fType, gType, g̅Type, ϕType, ψType, vType, ΩType, ∇HType, typeof(parameters), N}(d, m, n, ϑ, f, g, g̅, ϕ, ψ, v, Ω, ∇H, t₀,
                convert(Array{DT}, q₀), convert(Array{DT}, p₀), convert(Array{DT}, λ₀), convert(Array{DT}, μ₀), parameters, periodicity)
    end
end

function VDAE(ϑ, f, g, g̅, ϕ, ψ, t₀::Number, q₀::AbstractArray{DT}, p₀::AbstractArray{DT}, λ₀::AbstractArray{DT}=zero(q₀), μ₀::AbstractArray{DT}=zero(q₀); kwargs...) where {DT}
    VDAE(DT, ndims(q₀), size(q₀,1), size(μ₀,1), size(q₀,2), ϑ, f, g, g̅, ϕ, ψ, t₀, q₀, p₀, λ₀, μ₀; kwargs...)
end

function VDAE(ϑ, f, g, g̅, ϕ, ψ, q₀::AbstractArray, p₀::AbstractArray, λ₀::AbstractArray=zero(q₀), μ₀::AbstractArray=zero(q₀); kwargs...)
    VDAE(ϑ, f, g, g̅, ϕ, ψ, zero(eltype(q₀)), q₀, p₀, λ₀, μ₀; kwargs...)
end

Base.hash(ode::VDAE, h::UInt) = hash(ode.d, hash(ode.m, hash(ode.n,
          hash(ode.ϑ, hash(ode.f, hash(ode.g, hash(ode.g̅, hash(ode.ϕ,
          hash(ode.ψ, hash(ode.v, hash(ode.Ω, hash(ode.∇H,
          hash(ode.t₀, hash(ode.q₀, hash(ode.p₀, hash(ode.λ₀, hash(ode.μ₀,
          hash(ode.parameters, hash(ode.periodicity, h)))))))))))))))))))

Base.:(==)(ode1::VDAE, ode2::VDAE) = (
                                ode1.d == ode2.d
                             && ode1.m == ode2.m
                             && ode1.n == ode2.n
                             && ode1.ϑ == ode2.ϑ
                             && ode1.f == ode2.f
                             && ode1.g == ode2.g
                             && ode1.g̅ == ode2.g̅
                             && ode1.ϕ == ode2.ϕ
                             && ode1.ψ == ode2.ψ
                             && ode1.v == ode2.v
                             && ode1.Ω == ode2.Ω
                             && ode1.∇H == ode2.∇H
                             && ode1.t₀ == ode2.t₀
                             && ode1.q₀ == ode2.q₀
                             && ode1.p₀ == ode2.p₀
                             && ode1.λ₀ == ode2.λ₀
                             && ode1.μ₀ == ode2.μ₀
                             && ode1.parameters == ode2.parameters
                             && ode1.periodicity == ode2.periodicity)

function Base.similar(ode::VDAE, q₀, p₀, λ₀=get_λ₀(q₀, ode.λ₀), μ₀=get_λ₀(q₀, ode.μ₀); kwargs...)
    similar(ode, ode.t₀, q₀, p₀, λ₀, μ₀; kwargs...)
end

function Base.similar(ode::VDAE, t₀::TT, q₀::AbstractArray{DT}, p₀::AbstractArray{DT},
                      λ₀::AbstractArray{DT}=get_λ₀(q₀, ode.λ₀), μ₀=get_λ₀(q₀, ode.μ₀);
                      parameters=ode.parameters, periodicity=ode.periodicity) where {DT  <: Number, TT <: Number}
    @assert ode.d == size(q₀,1) == size(p₀,1) == size(λ₀,1)
    VDAE(ode.ϑ, ode.f, ode.g, ode.g̅, ode.ϕ, ode.ψ, t₀, q₀, p₀, λ₀, μ₀;
         v=ode.v, Ω=ode.Ω, ∇H=ode.∇H, parameters=parameters, periodicity=periodicity)
end

@inline Base.ndims(ode::VDAE) = ode.d

@inline periodicity(equation::VDAE) = equation.periodicity


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 = "finufft"
version = v"1.1.2"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/flatironinstitute/finufft.git", "5b92daf4244c92844dd2640ef80457c030e1bf25")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/finufft/
# The buggy makefile forces `-fopenmp` even if `OMP=OFF`
make lib -j${nproc} OMP="OFF" OMPFLAGS="" CFLAGS="-fPIC -O3 -funroll-loops -fcx-limited-range"
cp lib/libfinufft.so "${libdir}/libfinufft.${dlext}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="FFTW_jll", uuid="f5851436-0d7a-5f13-b9de-f02708fd171a"))
]

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


text:
"""
From: Theodorus Meuwissen
Sent: mandag 28. desember 2020 21:26

Hi Xijiang

It seems the DRP don’t work well.  Even when using EBV, the use of weights does
not work well.  Please try to use deregressed EBV Garrick et al.  (Genetics
Selection Evolution 2009 41:55)

Basically:

DRP = EBV/R2

DRP = deregressed proof which is to be analysed; R2 = reliability of the EBV = EDC / (EDC+lambda)

EDC = effective daughter contributions

Lambda = (4-h2)/h2 if EDC reflects number of daughters; h2 is heritability being used.

If EDC reflects effective number of records; lamda is from animal model.

The weight of the deregressed proofs is: lambda*R2/(1-R2)

Correction which does not affect relative weights but actual weights so that they are on the correct scale.
"""
function dr_ebv()

end


text:
function lazy_append(a::Triangulation,b::Triangulation)
  AppendedTriangulation(a,b)
end

struct AppendedTriangulation{Dc,Dp} <: Triangulation{Dc,Dp}
  a::Triangulation{Dc,Dp}
  b::Triangulation{Dc,Dp}
  function AppendedTriangulation(a::Triangulation{Dc,Dp},b::Triangulation{Dc,Dp}) where {Dc,Dp}
    @assert get_background_triangulation(a) === get_background_triangulation(b)
    new{Dc,Dp}(a,b)
  end
end

TriangulationStyle(::Type{<:AppendedTriangulation}) = SubTriangulation()

get_background_triangulation(trian::AppendedTriangulation) = get_background_triangulation(trian.a)

function get_cell_coordinates(trian::AppendedTriangulation)
  a = get_cell_coordinates(trian.a)
  b = get_cell_coordinates(trian.b)
  lazy_append(a,b)
end

function get_cell_ref_coordinates(trian::AppendedTriangulation)
  a = get_cell_ref_coordinates(trian.a)
  b = get_cell_ref_coordinates(trian.b)
  lazy_append(a,b)
end

function get_reffes(trian::AppendedTriangulation)
  vcat(get_reffes(trian.a),get_reffes(trian.b))
end

function get_cell_type(trian::AppendedTriangulation)
  a = get_cell_type(trian.a)
  b = get_cell_type(trian.b) .+ Int8(length(get_reffes(trian.a)))
  lazy_append(a,b)
end

#function reindex(f::AbstractArray,trian::AppendedTriangulation)
#  a = reindex(f,trian.a)
#  b = reindex(f,trian.b)
#  lazy_append(a,b)
#end

function get_cell_to_bgcell(trian::AppendedTriangulation)
  a = get_cell_to_bgcell(trian.a)
  b = get_cell_to_bgcell(trian.b)
  lazy_append(a,b)
end

function get_cell_ref_map(trian::AppendedTriangulation)
  a = get_cell_ref_map(trian.a)
  b = get_cell_ref_map(trian.b)
  lazy_append(a,b)
end

#function restrict(f::AbstractArray,trian::AppendedTriangulation)
#  a = restrict(f,trian.a)
#  b = restrict(f,trian.b)
#  lazy_append(a,b)
#end

# In this case, we do not want a lazy_append since it will become difficult to
# compress / expand the reffes.
#function get_cell_reffe(trian::AppendedTriangulation)
#  a = get_cell_reffe(trian.a)
#  b = get_cell_reffe(trian.b)
#  lazy_append(a,b)
#end

function get_cell_shapefuns(trian::AppendedTriangulation)
  a = get_cell_shapefuns(trian.a)
  b = get_cell_shapefuns(trian.b)
  lazy_append(a,b)
end

function get_cell_map(trian::AppendedTriangulation)
  a = get_cell_map(trian.a)
  b = get_cell_map(trian.b)
  lazy_append(a,b)
end

function get_facet_normal(trian::AppendedTriangulation)
  cm = get_cell_map(trian)
  a = get_facet_normal(trian.a)
  b = get_facet_normal(trian.b)
  lazy_append(a,b)
end

#function CellQuadrature(trian::AppendedTriangulation,degree1::Integer,degree2::Integer)
#  quad1 = CellQuadrature(trian.a,degree1)
#  quad2 = CellQuadrature(trian.b,degree2)
#  lazy_append(quad1,quad2)
#end
#
#function CellQuadrature(trian::AppendedTriangulation,degree::Integer)
#  CellQuadrature(trian,degree,degree)
#end

#function compress_contributions(cell_mat::AppendedArray,trian::AppendedTriangulation)
#  if length(cell_mat.a) == num_cells(trian.a) && length(cell_mat.b) == num_cells(trian.b)
#    a = compress_contributions(cell_mat.a,trian.a)
#    b = compress_contributions(cell_mat.b,trian.b)
#    return lazy_append(a,b)
#  else
#    return cell_mat
#  end
#end
#
#function compress_ids(cell_ids::AppendedArray,trian::AppendedTriangulation)
#  if length(cell_ids.a) == num_cells(trian.a) && length(cell_ids.b) == num_cells(trian.b)
#    a = compress_ids(cell_ids.a,trian.a)
#    b = compress_ids(cell_ids.b,trian.b)
#    return lazy_append(a,b)
#  else
#    return cell_ids
#  end
#end

function compress_contributions(cell_mat,trian::AppendedTriangulation)
  cell_to_bgcell = get_cell_to_bgcell(trian)
  ccell_mat = compress_contributions(cell_mat,cell_to_bgcell)
  ccell_mat
end

function compress_ids(cell_ids,trian::AppendedTriangulation)
  cell_to_bgcell = get_cell_to_bgcell(trian)
  compress_ids(cell_ids,cell_to_bgcell)
end


text:
using OutlierDetectionTrees
using OutlierDetectionTest

test_meta.(eval.(OutlierDetectionTrees.MODELS))

data = TestData()
run_test(detector) = test_detector(detector, data)

run_test(IForestDetector())


text:
# Julia wrapper for header: /usr/include/scip/scip_pricer.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludePricer(scip, name, desc, priority, delay, pricercopy, pricerfree, pricerinit, pricerexit, pricerinitsol, pricerexitsol, pricerredcost, pricerfarkas, pricerdata)
    ccall((:SCIPincludePricer, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Cint, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_PRICERDATA}), scip, name, desc, priority, delay, pricercopy, pricerfree, pricerinit, pricerexit, pricerinitsol, pricerexitsol, pricerredcost, pricerfarkas, pricerdata)
end

function SCIPincludePricerBasic(scip, pricerptr, name, desc, priority, delay, pricerredcost, pricerfarkas, pricerdata)
    ccall((:SCIPincludePricerBasic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_PRICER}}, Cstring, Cstring, Cint, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_PRICERDATA}), scip, pricerptr, name, desc, priority, delay, pricerredcost, pricerfarkas, pricerdata)
end

function SCIPsetPricerCopy(scip, pricer, pricercopy)
    ccall((:SCIPsetPricerCopy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PRICER}, Ptr{Cvoid}), scip, pricer, pricercopy)
end

function SCIPsetPricerFree(scip, pricer, pricerfree)
    ccall((:SCIPsetPricerFree, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PRICER}, Ptr{Cvoid}), scip, pricer, pricerfree)
end

function SCIPsetPricerInit(scip, pricer, pricerinit)
    ccall((:SCIPsetPricerInit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PRICER}, Ptr{Cvoid}), scip, pricer, pricerinit)
end

function SCIPsetPricerExit(scip, pricer, pricerexit)
    ccall((:SCIPsetPricerExit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PRICER}, Ptr{Cvoid}), scip, pricer, pricerexit)
end

function SCIPsetPricerInitsol(scip, pricer, pricerinitsol)
    ccall((:SCIPsetPricerInitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PRICER}, Ptr{Cvoid}), scip, pricer, pricerinitsol)
end

function SCIPsetPricerExitsol(scip, pricer, pricerexitsol)
    ccall((:SCIPsetPricerExitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PRICER}, Ptr{Cvoid}), scip, pricer, pricerexitsol)
end

function SCIPfindPricer(scip, name)
    ccall((:SCIPfindPricer, libscip), Ptr{SCIP_PRICER}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetPricers(scip)
    ccall((:SCIPgetPricers, libscip), Ptr{Ptr{SCIP_PRICER}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNPricers(scip)
    ccall((:SCIPgetNPricers, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNActivePricers(scip)
    ccall((:SCIPgetNActivePricers, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPsetPricerPriority(scip, pricer, priority)
    ccall((:SCIPsetPricerPriority, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PRICER}, Cint), scip, pricer, priority)
end

function SCIPactivatePricer(scip, pricer)
    ccall((:SCIPactivatePricer, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PRICER}), scip, pricer)
end

function SCIPdeactivatePricer(scip, pricer)
    ccall((:SCIPdeactivatePricer, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PRICER}), scip, pricer)
end


text:
"""
```
SlipSystem(;
    crystalStruct::AbstractCrystalStruct,
    slipPlane::AbstractArray,
    bVec::AbstractArray
)
```
Creates a [`SlipSystem`](@ref).

Ensures `slipPlane ⟂ bVec`.
"""
function SlipSystem(;
    crystalStruct::AbstractCrystalStruct,
    slipPlane::AbstractArray,
    bVec::AbstractArray,
)
    if sum(slipPlane .!= 0) != 0 && sum(bVec .!= 0) != 0
        if ndims(slipPlane) == 1
            @assert slipPlane ⋅ bVec ≈ 0 "SlipSystem: slip plane, n == $(slipPlane), and Burgers vector, b = $(bVec), must be orthogonal."
        else
            idx = findall(x -> !(x ≈ 0), vec(sum(slipPlane .* bVec, dims = 1)))
            @assert isempty(idx) "SlipSystem: entries of the slip plane, n[:, $idx] = $(slipPlane[:, idx]), and Burgers vector, b[:, $idx] = $(bVec[:, idx]), are not orthogonal."
        end
    end

    return SlipSystem(crystalStruct, slipPlane, bVec)
end

"""
```
DislocationParameters(;
    mobility::AbstractMobility,
    dragCoeffs = (edge = 1.0, screw = 2.0, climb = 1e9),
    coreRad = 1.0,
    coreRadMag = 1.0,
    coreEnergy = 1 / (4 * π) * log(coreRad / 0.1),
    minSegLen = 2 * coreRad,
    maxSegLen = 20 * coreRad,
    minArea = minSegLen^2 / sqrt(2),
    maxArea = 100 * minArea,
    collisionDist = minSegLen / 2,
    slipStepCritLen = maxSegLen / 2,
    slipStepCritArea = 0.5 * (slipStepCritLen^2) * sind(1),
    remesh = true,
    collision = true,
    separation = true,
    virtualRemesh = true,
    parCPU = false,
    parGPU = false,
)
```
Creates [`DislocationParameters`](@ref). Automatically calculates derived values, asserts values are reasonable and provides sensible default values.
"""
function DislocationParameters(;
    mobility::AbstractMobility,
    dragCoeffs = (edge = 1.0, screw = 2.0, climb = 1e9),
    coreRad = 1.0,
    coreRadMag = 1.0,
    coreEnergy = 1 / (4 * π) * log(coreRad / 0.1),
    minSegLen = 2 * coreRad,
    maxSegLen = 20 * coreRad,
    minArea = minSegLen^2 / sqrt(2),
    maxArea = 100 * minArea,
    collisionDist = minSegLen / 2,
    slipStepCritLen = maxSegLen / 2,
    slipStepCritArea = 0.5 * (slipStepCritLen^2) * sind(1),
    remesh = true,
    collision = true,
    separation = true,
    virtualRemesh = true,
    parCPU = false,
    parGPU = false,
)
    coreRad == minSegLen == maxSegLen == 0 ? nothing :
    @assert coreRad < minSegLen < maxSegLen
    minArea == maxArea == 0 ? nothing : @assert minArea < maxArea

    return DislocationParameters(
        mobility,
        dragCoeffs,
        coreRad,
        coreRad^2,
        coreRadMag,
        coreEnergy,
        minSegLen,
        maxSegLen,
        minSegLen * 2,
        minSegLen^2,
        minArea,
        maxArea,
        minArea^2,
        maxArea^2,
        collisionDist,
        collisionDist^2,
        slipStepCritLen,
        slipStepCritArea,
        remesh,
        collision,
        separation,
        virtualRemesh,
        parCPU,
        parGPU,
    )
end

"""
```
DislocationLoop(
    loopType::AbstractDlnStr,
    numSides,
    nodeSide,
    numLoops,
    segLen,
    slipSystem,
    _slipPlane,
    _bVec,
    label,
    buffer,
    range,
    dist,
)
```
Fallback for creating a generic [`DislocationLoop`](@ref).
"""
function DislocationLoop(
    loopType::AbstractDlnStr,
    numSides,
    nodeSide,
    numLoops,
    segLen,
    slipSystemIdx,
    slipSystem,
    label,
    buffer,
    range,
    dist,
)
    nodeTotal::Int = 0
    links = zeros(MMatrix{2, nodeTotal, Int})
    coord = zeros(MMatrix{3, nodeTotal})
    slipPlane = zeros(MMatrix{3, 0})
    bVec = zeros(MMatrix{3, 0})

    return DislocationLoop(
        loopType,
        numSides,
        nodeSide,
        numLoops,
        segLen,
        0,
        label,
        links,
        slipPlane,
        bVec,
        coord,
        buffer,
        range,
        dist,
    )
end
"""
```
DislocationLoop(
    loopType::loopPure,
    numSides,
    nodeSide,
    numLoops,
    segLen,
    slipSystem,
    _slipPlane::AbstractArray,
    _bVec::AbstractArray,
    label::AbstractVector{nodeTypeDln},
    buffer,
    range,
    dist::AbstractDistribution,
)
```
Constructor for `loopPure` [`DislocationLoop`](@ref)s.

# Inputs

- `loopType::loopPure`: can be either `loopPrism()` or `loopShear()` to make prismatic or shear loops.
- `numSides`: number of sides in the loop
- `nodeSide`: nodes per side of the loop
- `numLoops`: number of loops to generate when making the dislocation network
- `segLen`: length of each initial dislocation segment
- `slipSystem`: slip system from [`SlipSystem`](@ref) the loop belongs to
- `_slipPlane`: slip plane vector
- `_bVec`: Burgers vector
- `label`: node labels of type [`nodeTypeDln`](@ref)
- `buffer`: buffer for increasing the spread of the generated dislocation loops in the network
- `range`: range on which the loops will be distributed in the network
- `dist`: distribution used to generate the network
"""
function DislocationLoop(
    loopType::loopPure,
    numSides,
    nodeSide,
    numLoops,
    segLen,
    slipSystemIdx,
    slipSystem::SlipSystem,
    label::AbstractVector{nodeTypeDln},
    buffer,
    range,
    dist::AbstractDistribution,
)
    nodeTotal = numSides * nodeSide # Calculate total number of nodes for memory allocation.
    numSegLen = length(segLen) # Number of segment lengths.

    # Validate input.
    @assert length(label) == nodeTotal "DislocationLoop: All $nodeTotal nodes must be labelled. There are $(length(label)) labels currently defined."
    @assert numSegLen == nodeTotal "DislocationLoop: All $nodeTotal segments must have their lengths defined. There are $numSegLen lengths currently defined."

    _slipPlane = slipSystem.slipPlane[:, slipSystemIdx]
    _bVec = slipSystem.bVec[:, slipSystemIdx]
    # Normalise vectors.
    elemT = eltype(_slipPlane)
    _slipPlane = _slipPlane / norm(_slipPlane)
    _bVec = _bVec / norm(_bVec)

    # Pick rotation axis for segments.
    # Shear loops rotate around slip plane vector. They have screw, mixed and edge segments.
    if typeof(loopType) == loopShear
        rotAxis = SVector{3, elemT}(_slipPlane[1], _slipPlane[2], _slipPlane[3])
        # Prismatic loops rotate around Burgers vector. All segments are edge.
    else
        rotAxis = SVector{3, elemT}(_bVec[1], _bVec[2], _bVec[3])
        # Catch all.
    end

    # Allocate arrays.
    links = zeros(MMatrix{2, nodeTotal, Int})
    coord = zeros(MMatrix{3, nodeTotal})
    slipPlane = MMatrix{3, nodeTotal}(repeat(_slipPlane, inner = (1, numSegLen)))
    bVec = MMatrix{3, nodeTotal}(repeat(_bVec, inner = (1, numSegLen)))
    seg = zeros(MMatrix{3, numSegLen})

    # Create initial segments.
    staticSlipPlane = SVector{3, elemT}(_slipPlane[1], _slipPlane[2], _slipPlane[3])
    staticBVec = SVector{3, elemT}(_bVec[1], _bVec[2], _bVec[3])
    @inbounds @simd for i in eachindex(segLen)
        seg[:, i] = makeSegment(segEdge(), staticSlipPlane, staticBVec) * segLen[i]
    end

    θ = externalAngle(numSides)  # External angle of a regular polygon with numSides.

    # Loop over polygon's sides.
    origin = SVector{3, elemT}(0, 0, 0)
    @inbounds for i in 1:numSides
        # Index for side i.
        idx = (i - 1) * nodeSide
        # Rotate segments by external angle of polygon to make polygonal loop.
        modIdx = mod(i - 1, numSegLen) + 1
        staticSeg = SVector{3, elemT}(seg[1, modIdx], seg[2, modIdx], seg[3, modIdx])
        rseg = rot3D(staticSeg, rotAxis, origin, θ * (i - 1))
        # DO NOT add @simd, this loop works by adding rseg to the previous coordinate to make the loop. Loop over the nodes per side.
        for j in 1:nodeSide
            # Count first node once.
            if i == j == 1
                coord[:, 1] .= 0 # Initial coordinate is on the origin.
                continue
            end
            if idx + j <= nodeTotal
                # Add segment vector to previous coordinate.
                coord[:, idx + j] += @views coord[:, idx + j - 1] + rseg
            end
        end
    end

    # Find centre of the loop and make it zero.
    meanCoord = mean(coord, dims = 2)
    coord .-= meanCoord

    # Create links matrix.
    @inbounds @simd for j in 1:(nodeTotal - 1)
        links[:, j] .= (j, j + 1)
    end
    links[:, nodeTotal] .= (nodeTotal, 1)

    return DislocationLoop(
        loopType,
        numSides,
        nodeSide,
        numLoops,
        segLen,
        slipSystemIdx,
        label,
        links,
        slipPlane,
        bVec,
        coord,
        buffer,
        range,
        dist,
    )
end
"""
```
DislocationLoop(
    loopType::loopImpure,
    numSides,
    nodeSide,
    numLoops,
    segLen,
    slipSystem,
    _slipPlane::AbstractArray,
    _bVec::AbstractArray,
    label::AbstractVector{nodeTypeDln},
    buffer,
    range,
    dist::AbstractDistribution,
)
```
Fallback [`DislocationLoop`](@ref) constructor for other as of yet unimplemented `loopImpure`.
"""
function DislocationLoop(
    loopType::loopImpure,
    numSides,
    nodeSide,
    numLoops,
    segLen,
    slipSystemIdx,
    slipSystem::SlipSystem,
    label::AbstractVector{nodeTypeDln},
    buffer,
    range,
    dist::AbstractDistribution,
)
    @warn "DislocationLoop: Constructor for $(typeof(loopType)) not defined, defaulting to prismatic loop."
    return DislocationLoop(
        loopPrism(),
        numSides,
        nodeSide,
        numLoops,
        segLen,
        slipSystemIdx,
        slipSystem,
        label,
        buffer,
        range,
        dist,
    )
end
"""
```
DislocationLoop(;
    loopType::AbstractDlnStr,
    numSides,
    nodeSide,
    numLoops,
    segLen,
    slipSystem,
    _slipPlane,
    _bVec,
    label,
    buffer,
    range,
    dist,
)
```
Create a [`DislocationLoop`](@ref).
"""
function DislocationLoop(;
    loopType::AbstractDlnStr,
    numSides,
    nodeSide,
    numLoops,
    segLen,
    slipSystemIdx,
    slipSystem::SlipSystem,
    label,
    buffer,
    range,
    dist,
)
    return DislocationLoop(
        loopType,
        numSides,
        nodeSide,
        numLoops,
        segLen,
        slipSystemIdx,
        slipSystem,
        label,
        buffer,
        range,
        dist,
    )
end

"""
```
DislocationNetwork(;
    links::AbstractArray,
    slipPlane::AbstractArray,
    bVec::AbstractArray,
    coord::AbstractArray,
    label::AbstractVector{nodeTypeDln},
    nodeVel::AbstractArray,
    nodeForce::AbstractArray,
    numNode = length(label),
    numSeg = size(links, 2),
    maxConnect = 4,
    connectivity::AbstractArray = zeros(Int, 1 + 2 * maxConnect, length(label)),
    linksConnect::AbstractArray = zeros(Int, 2, size(links, 2)),
    segIdx::AbstractArray = zeros(Int, size(links, 2), 3),
    segForce::AbstractArray = zeros(3, 2, size(links, 2)),
)
```
Create a [`DislocationNetwork`](@ref). We recommend generating networks from [`DislocationLoop`](@ref) unless you want a special case.
"""
function DislocationNetwork(;
    links::AbstractArray,
    slipPlane::AbstractArray,
    bVec::AbstractArray,
    coord::AbstractArray,
    label::AbstractVector{nodeTypeDln},
    nodeVel::AbstractArray,
    nodeForce::AbstractArray,
    numNode = length(label),
    numSeg = size(links, 2),
    maxConnect = 4,
    connectivity::AbstractArray = zeros(Int, 1 + 2 * maxConnect, length(label)),
    linksConnect::AbstractArray = zeros(Int, 2, size(links, 2)),
    extSeg::AbstractArray = zeros(Bool, size(links, 2)),
    segIdx::AbstractArray = zeros(Int, size(links, 2), 3),
    segForce::AbstractArray = zeros(3, 2, size(links, 2)),
)
    @assert size(links, 1) == size(segForce, 2) == 2
    @assert size(bVec, 1) == size(slipPlane, 1) == size(coord, 1) size(segForce, 1) == 3
    @assert size(links, 2) == size(bVec, 2) == size(slipPlane, 2) == size(segForce, 3)
    @assert size(coord, 2) == length(label)
    @assert length(numNode) == length(numSeg) == 1

    typeof(numNode) <: AbstractVector ? numNodeArr = numNode : numNodeArr = [numNode]
    typeof(numSeg) <: AbstractVector ? numSegArr = numSeg : numSegArr = [numSeg]

    return DislocationNetwork(
        numNodeArr,
        numSegArr,
        maxConnect,
        label,
        links,
        connectivity,
        linksConnect,
        slipPlane,
        extSeg,
        segIdx,
        bVec,
        coord,
        nodeVel,
        nodeForce,
        segForce,
    )
end
"""
```
DislocationNetwork(
    sources::DislocationLoopCollection,
    maxConnect = 4,
    args...;
    memBuffer = nothing,
    checkConsistency = true,
    kw...,
)
```
Creates a [`DislocationNetwork`](@ref) out of a [`DislocationLoopCollection`](@ref).

# Inputs

- `args...` are optional arguments that will be passed on to the [`loopDistribution`](@ref) function which distributes the loops in `sources` according to the type of their `dist` variable.
- `kw...` are optional keyword arguments that will also be passed to `loopDistribution`.
- `memBuffer` is the numerical value for allocating memory in advance. The quantity, `memBuffer × N`, where `N` is the total number of nodes in `sources`, will be the initial number of entries allocated in the matrices that keep the network's data. If no `memBuffer` is provided, the number of entries allocated will be `round(N*log2(N)).
"""
function DislocationNetwork(
    sources::DislocationLoopCollection,
    maxConnect = 4,
    args...;
    memBuffer = nothing,
    checkConsistency = true,
    kw...,
)
    # Initialisation.
    nodeTotal::Int = 0
    lims = zeros(MMatrix{3, 2})
    # Calculate node total.
    for i in eachindex(sources)
        nodeTotal += sources[i].numLoops * length(sources[i].label)
    end
    # Memory buffer.
    isnothing(memBuffer) ? nodeBuffer = Int(round(nodeTotal * log2(nodeTotal))) :
    nodeBuffer = nodeTotal * Int(memBuffer)

    # Allocate memory.
    links = zeros(Int, 2, nodeBuffer)
    slipPlane = zeros(3, nodeBuffer)
    bVec = zeros(3, nodeBuffer)
    coord = zeros(3, nodeBuffer)
    label = zeros(nodeTypeDln, nodeBuffer)
    nodeVel = zeros(Float64, 3, nodeBuffer)
    nodeForce = zeros(Float64, 3, nodeBuffer)
    numNode = nodeTotal
    numSeg = nodeTotal
    segForce = zeros(Float64, 3, 2, nodeBuffer)

    initIdx = 1
    # Fill the matrices that will make up the network.
    makeNetwork!(
        links,
        slipPlane,
        bVec,
        coord,
        label,
        sources,
        lims,
        initIdx,
        args...;
        kw...,
    )

    # Calculate number of segments and indexing matrix.
    segIdx, extSeg = getSegmentIdx(links, label)
    # Generate connectivity and linksConnect matrix.
    connectivity, linksConnect = makeConnect(links, maxConnect)

    # Create network.
    network = DislocationNetwork(
        [numNode],
        [numSeg],
        maxConnect,
        label,
        links,
        connectivity,
        linksConnect,
        slipPlane,
        extSeg,
        segIdx,
        bVec,
        coord,
        nodeVel,
        nodeForce,
        segForce,
    )

    # Check that the network is generated properly.
    checkConsistency ? checkNetwork(network) : nothing

    return network
end
"""
```
DislocationNetwork!(
    network::DislocationNetwork,
    sources::DislocationLoopCollection,
    args...;
    memBuffer = nothing,
    checkConsistency = true,
    kw...,
)
```
Adds a [`DislocationLoopCollection`](@ref) to an existing [`DislocationNetwork`](@ref).
"""
function DislocationNetwork!(
    network::DislocationNetwork,
    sources::DislocationLoopCollection,
    args...;
    memBuffer = nothing,
    checkConsistency = true,
    kw...,
)
    # For comments see DislocationNetwork. It is a 1-to-1 translation except that this one modifies the network in-place.

    iszero(network) && return DislocationNetwork(
        sources,
        args...;
        checkConsistency = checkConsistency,
        kw...,
    )

    nodeTotal::Int = 0
    lims = zeros(MMatrix{3, 2})
    for i in eachindex(sources)
        nodeTotal += sources[i].numLoops * length(sources[i].label)
    end
    numNode = nodeTotal

    # Allocate memory.
    available = length(findall(x -> x == noneDln, network.label))
    if nodeTotal > available
        newEntries = Int(round(nodeTotal * log2(nodeTotal)))
        network = push!(network, newEntries)
    end

    links = network.links
    slipPlane = network.slipPlane
    bVec = network.bVec
    coord = network.coord
    label = network.label

    # Since the network has already been created, initIdx is the next available index to store new data.
    initIdx::Int = 1
    first = findfirst(x -> x == noneDln, label)
    isnothing(first) ? initIdx = 1 : initIdx = first
    makeNetwork!(
        links,
        slipPlane,
        bVec,
        coord,
        label,
        sources,
        lims,
        initIdx,
        args...;
        kw...,
    )
    network.numNode[1] += numNode
    network.numSeg[1] += numNode

    getSegmentIdx!(network)
    makeConnect!(network)

    checkConsistency ? checkNetwork(network) : nothing
    return network
end

"""
```
makeNetwork!(
    links,
    slipPlane,
    bVec,
    coord,
    label,
    sources,
    lims,
    initIdx,
    args...;
    kw...,
)
```
Internal function called by [`DislocationNetwork`](@ref) to fill the arrays that define the network.
"""
function makeNetwork!(
    links,
    slipPlane,
    bVec,
    coord,
    label,
    sources,
    lims,
    initIdx,
    args...;
    kw...,
)
    nodeTotal::Int = 0
    elemT = eltype(coord)
    @inbounds @simd for i in eachindex(sources)
        idx = initIdx + nodeTotal
        nodesLoop = length(sources[i].label)
        numLoops = sources[i].numLoops
        numNodes = numLoops * nodesLoop
        # Calculate the normalised displacements for all loops in sources[i] according to their distribution.
        disp = loopDistribution(sources[i].dist, numLoops, args...; kw...)
        # Calculate the real spatial limits of the distributions.
        limits!(lims, mean(sources[i].segLen), sources[i].range, sources[i].buffer)
        # Fill out the data for all loops specified in sources[i].
        for j in 1:numLoops
            # The number of nodes in the loop is nodesLoop, so that's our stride inside sources[i]
            idxi = idx + (j - 1) * nodesLoop
            idxf = idxi + nodesLoop - 1
            # Links are numbered sequentially in network so we have to account for previously assigned links.
            links[:, idxi:idxf] .=
                sources[i].links[:, 1:nodesLoop] .+ (nodeTotal + initIdx - 1)
            slipPlane[:, idxi:idxf] = sources[i].slipPlane[:, 1:nodesLoop]
            bVec[:, idxi:idxf] = sources[i].bVec[:, 1:nodesLoop]
            coord[:, idxi:idxf] = sources[i].coord[:, 1:nodesLoop]
            label[idxi:idxf] = sources[i].label[1:nodesLoop]
            # Map the normalised displacements to real space using the real limits and translate the nodes' coordinates accordingly.
            staticDisp = SVector{3, elemT}(disp[1, j], disp[2, j], disp[3, j])
            viewCoord = @view coord[:, idxi:idxf]
            translatePoints!(viewCoord, lims, staticDisp)
            nodeTotal += nodesLoop
        end
    end
    return nothing
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
module MathOptInterface

"""
    ModelLike

Abstract supertype for objects that implement the "Model" interface for defining
an optimization problem.
"""
abstract type ModelLike end

# This allows to use `ModelLike`s in broadcast calls without the need to
# embed it in a `Ref`
Base.broadcastable(model::ModelLike) = Ref(model)

function Base.show(io::IO, model::ModelLike)
    return Utilities.print_with_acronym(io, summary(model))
end

"""
    AbstractOptimizer <: ModelLike

Abstract supertype for objects representing an instance of an optimization
problem tied to a particular solver. This is typically a solver's in-memory
representation. In addition to `ModelLike`, `AbstractOptimizer` objects let you
solve the model and query the solution.
"""
abstract type AbstractOptimizer <: ModelLike end

"""
    optimize!(optimizer::AbstractOptimizer)

Optimize the problem contained in `optimizer`.

Before calling `optimize!`, the problem should first be constructed using the
incremental interface (see [`supports_incremental_interface`](@ref)) or [`copy_to`](@ref).
"""
function optimize! end

"""
    optimize!(dest::AbstractOptimizer, src::ModelLike)::Tuple{IndexMap,Bool}

A "one-shot" call that copies the problem from `src` into `dest` and then uses
`dest` to optimize the problem.

Returns a tuple of an [`IndexMap`](@ref) and a `Bool` `copied`.

 * The [`IndexMap`](@ref) object translates variable and constraint indices from
   the `src` model to the corresponding indices in the `dest` optimizer. See
   [`copy_to`](@ref) for details.
 * If `copied == true`, `src` was copied to `dest` and then cached, allowing
   incremental modification if supported by the solver.
 * If `copied == false`, a cache of the model was _not_ kept in `dest`.
   Therefore, only the solution information (attributes for which
   [`is_set_by_optimize`](@ref) is true) is available to query.

!!! note
    The main purpose of `optimize!` method with two arguments is for use in
    [`Utilities.CachingOptimizer`](@ref).

!!! warning
    The new `optimize!` method with two arguments is an experimental new feature
    of MOI v0.10.2 that may break in MOI v1.0.

## Relationship to the single-argument `optimize!`

The default fallback of `optimize!(dest::AbstractOptimizer, src::ModelLike)` is
```julia
function optimize!(dest::AbstractOptimizer, src::ModelLike)
    index_map = copy_to(dest, src)
    optimize!(dest)
    return index_map, true
end
```
Therefore, subtypes of [`AbstractOptimizer`](@ref) should either implement this
two-argument method, or implement both [`copy_to(::Optimizer, ::ModelLike)`](@ref copy_to)
and `optimize!(::Optimizer)`.
"""
function optimize!(dest, src)
    # The arguments above are untyped to avoid ambiguities.
    index_map = copy_to(dest, src)
    optimize!(dest)
    return index_map, true
end

"""
    compute_conflict!(optimizer::AbstractOptimizer)

Computes a minimal subset of constraints such that the model with the other
constraint removed is still infeasible.

Some solvers call a set of conflicting constraints an Irreducible Inconsistent
Subsystem (IIS).

See also [`ConflictStatus`](@ref) and [`ConstraintConflictStatus`](@ref).

### Note

If the model is modified after a call to `compute_conflict!`, the implementor
is not obliged to purge the conflict. Any calls to the above attributes may
return values for the original conflict without a warning. Similarly, when
modifying the model, the conflict can be discarded.
"""
function compute_conflict!(optimizer::AbstractOptimizer)
    return throw(
        ArgumentError(
            "The optimizer $(typeof(optimizer)) does not support " *
            "`compute_conflict!`",
        ),
    )
end

"""
    write_to_file(model::ModelLike, filename::String)

Writes the current model data to the given file.
Supported file types depend on the model type.
"""
function write_to_file end

"""
    read_from_file(model::ModelLike, filename::String)

Read the file `filename` into the model `model`. If `model` is non-empty, this may
throw an error.

Supported file types depend on the model type.

### Note

Once the contents of the file are loaded into the model, users can query the variables via
`get(model, ListOfVariableIndices())`. However, some filetypes, such as LP files, do not
maintain an explicit ordering of the variables. Therefore, the returned list may be in an
arbitrary order. To avoid depending on the order of the indices, users should look up each
variable index by name: `get(model, VariableIndex, "name")`.
"""
function read_from_file end

"""
    is_empty(model::ModelLike)

Returns `false` if the `model` has any model attribute set or has any variables
or constraints.

Note that an empty model can have optimizer attributes set.
"""
function is_empty end

"""
    empty!(model::ModelLike)

Empty the model, that is, remove all variables, constraints and model attributes
but not optimizer attributes.
"""
function empty! end

"""
    supports_incremental_interface(model::ModelLike)

Return a `Bool` indicating whether `model` supports building incrementally via
[`add_variable`](@ref) and [`add_constraint`](@ref).

The main purpose of this function is to determine whether a model can be loaded
into `model` incrementally or whether it should be cached and copied at once
instead.
"""
supports_incremental_interface(::ModelLike) = false

"""
    copy_to(dest::ModelLike, src::ModelLike)::IndexMap

Copy the model from `src` into `dest`.

The target `dest` is emptied, and all previous indices to variables and
constraints in `dest` are invalidated.

Returns an [`IndexMap`](@ref) object that translates variable and constraint
indices from the `src` model to the corresponding indices in the `dest` model.

## Notes

 * If a constraint that in `src` is not supported by `dest`, then an
   [`UnsupportedConstraint`](@ref) error is thrown.
 * If an [`AbstractModelAttribute`](@ref), [`AbstractVariableAttribute`](@ref),
   or [`AbstractConstraintAttribute`](@ref) is set in `src` but not supported by
   `dest`, then an [`UnsupportedAttribute`](@ref) error is thrown.

[`AbstractOptimizerAttribute`](@ref)s are _not_ copied  to the `dest` model.

## IndexMap

Implementations of `copy_to` must return an [`IndexMap`](@ref). For technical
reasons, this type is defined in the Utilities submodule as
`MOI.Utilities.IndexMap`. However, since it is an integral part of the MOI API,
we provide `MOI.IndexMap` as an alias.

### Example

```julia
# Given empty `ModelLike` objects `src` and `dest`.

x = add_variable(src)

is_valid(src, x)   # true
is_valid(dest, x)  # false (`dest` has no variables)

index_map = copy_to(dest, src)
is_valid(dest, x) # false (unless index_map[x] == x)
is_valid(dest, index_map[x]) # true
```
"""
function copy_to(dest, src)
    return error(
        "`copy_to` is not supported by the solver `$(typeof(dest))`. Did " *
        "you mean to call " *
        "`optimize!(dest::AbstractOptimizer, src::ModelLike)` instead?",
    )
end

include("error.jl")
include("indextypes.jl")
include("functions.jl")
include("sets.jl")
include("attributes.jl")
include("constraints.jl")
include("modifications.jl")
include("variables.jl")
include("nlp.jl")

include("precompile.jl")

# submodules
include("Utilities/Utilities.jl") # MOI.Utilities
include("Test/Test.jl")
include("Bridges/Bridges.jl")     # MOI.Bridges
include("Benchmarks/Benchmarks.jl")
include("FileFormats/FileFormats.jl")

include("instantiate.jl")

_precompile_()

"""
    IndexMap()

The dictionary-like object returned by [`copy_to`](@ref).

## IndexMap

Implementations of [`copy_to`](@ref) must return an [`IndexMap`](@ref). For
technical reasons, the `IndexMap` type is defined in the Utilties submodule as
`MOI.Utilities.IndexMap`. However, since it is an integral part of the MOI API,
we provide this `MOI.IndexMap` as an alias.
"""
const IndexMap = Utilities.IndexMap

end


text:
using LinearAlgebra

@testset "generators" begin
    @testset "erdos_renyi" begin
        @test_throws DomainError erdos_renyi_undirected(10,11.)
        @test_throws DomainError erdos_renyi_directed(10,11.)

        n = 100
        avgdegs = range(1., stop=2*log(n), length=100)
        compsizes = map( (dbar) ->
                maximum(scomponents(erdos_renyi_undirected(n,dbar)).sizes),
            avgdegs )

        @test is_undirected(erdos_renyi_undirected(10,2.))

        erdos_renyi_undirected(0, 0.)
        erdos_renyi_undirected(0, 0)
        erdos_renyi_undirected(5, 0)
        erdos_renyi_undirected(5, 0.)
        erdos_renyi_undirected(5, 1.)
        erdos_renyi_undirected(5, 1)
        @test all(diag(sparse_transpose(erdos_renyi_undirected(10, 0.5))) .== 0.)
        erdos_renyi_directed(0, 0)
        erdos_renyi_directed(0, 0)
        erdos_renyi_directed(5, 0.)
        erdos_renyi_directed(5, 1.)
        @test all(diag(sparse_transpose(erdos_renyi_directed(10, 0.5))) .== 0.)
    end

    @testset "chung_lu" begin
        @test_throws ArgumentError chung_lu_undirected([1,2],2)
        @test_throws ArgumentError chung_lu_undirected([1,3])
        @test size(chung_lu_undirected([1,2]),1) == 2
    end

    @testset "havel_hakimi" begin
        @test is_graphical_sequence([1,1]) == true
        @test is_graphical_sequence([0,1,1]) == true
        @test is_graphical_sequence([2,2,2]) == true
        @test is_graphical_sequence([1,1,2]) == true
        @test is_graphical_sequence([1,2,2]) == false
        @test is_graphical_sequence([1,1,3]) == false
        @test is_graphical_sequence([0]) == true
        @test is_graphical_sequence(Int[]) == true
        @test is_graphical_sequence([1,1,1,1]) == true

        @test Matrix(sparse_transpose(havel_hakimi_graph(Int[]))) == zeros(Bool,0,0)
        @test Matrix(sparse_transpose(havel_hakimi_graph([0]))) == zeros(Bool,1,1)
        @test Matrix(sparse_transpose(havel_hakimi_graph([1,1]))) == [0 1; 1 0]
        @test Matrix(sparse_transpose(havel_hakimi_graph([2,2,2]))) == [0 1 1; 1 0 1; 1 1 0]
        six_star = [0 1 1 1 1 1; 1 0 0 0 0 0; 1 0 0 0 0 0; 1 0 0 0 0 0; 1 0 0 0 0 0; 1 0 0 0 0 0]
        @test Matrix(sparse_transpose(havel_hakimi_graph([5,1,1,1,1,1]))) == six_star

        @test_throws ArgumentError is_graphical_sequence([-1])
        @test_throws ArgumentError havel_hakimi_graph([1,2,2])
    end

    @testset "pa_graph" begin
        @test typeof(pa_graph(10,5,5)) == MatrixNetwork{Bool}
        @test typeof(pa_edges!(2,1,[(1,1)])) == Vector{Tuple{Int,Int}}
        @test_throws ArgumentError pa_edges!(5,2,Vector{Tuple{Int,Int}}())
        @test is_empty(pa_graph(0,0,0)) == true
        @test all(diag(sparse_transpose(pa_graph(10, 2, 3))) .== 0.)
        @test is_undirected(pa_graph(10, 12, 3))
        @test_throws ArgumentError pa_graph(-1,10,3)
        @test_throws ArgumentError pa_graph(5,10,-3)
        @test maximum(map(first, pa_edges!(5,1,[(1,2),(2,1)],2))) == 7
        @test size(pa_graph(10,5,5),1) == 10
        @test size(pa_graph(10,-5,5),1) == 10
        @test nnz(sparse_transpose(pa_graph(10,-5,5))) == 20
    end

    @testset "gpa_graph" begin
        @test typeof(gpa_graph(10,0.5,0.3,2)) == MatrixNetwork{Bool}
        @test typeof(gpa_graph(10,0.5,0.3,2,Val{true})) == MatrixNetwork{Bool}
        @test typeof(gpa_edges!(4,.75,.25,[(1,2)],1)) == Vector{Tuple{Int,Int}}
        @test_throws ArgumentError gpa_edges!(4,.75,.25,[(1,1)],1)
        @test typeof(gpa_edges!(4,.75,.25,[(1,1)],1,Val{true})) == Vector{Tuple{Int,Int}}
        @test is_empty(gpa_graph(0,0.0,0.0,0)) == true
        @test all(diag(sparse_transpose(gpa_graph(10, .3,.4, 3))) .== 0.)
        @test is_undirected(gpa_graph(10, .2, .8,4))
        @test_throws ArgumentError gpa_graph(-1,.3,.4,5)
        @test_throws ArgumentError gpa_graph(8,.8,.4,2)
        @test_throws ArgumentError gpa_graph(10,.3,.4,-2)
        @test size(gpa_graph(12,.5,.4,2),1) == 12
        @test size(gpa_graph(21,.1,.1,3),1) == 21
    end

    @testset "roach_graph" begin
        @test_throws ArgumentError roach_graph(-1)
        @test is_empty(roach_graph(0))
        @test is_connected(roach_graph(5))
        @test is_connected(roach_graph(5, Val{true})[1])
        @test maximum(bfs(roach_graph(5, Val{false}),20)[1]) == 11
        @test maximum(bfs(roach_graph(10),40)[1]) == 21
        @test (bfs(roach_graph(10),40)[1])[1] == 20

        G,xy = roach_graph(5, Val{true})
        filt = vec(all(xy .>= 0,dims = 2))
        Asub = sparse_transpose(G)[filt,filt]
        # Asub should be an antennae...
        @test bfs(Asub,1)[1][5] == 4
        filt = vec(all(xy .<= 0,dims = 2))
        Asub = sparse_transpose(G)[filt,filt]
        # Asub should be an line...
        @test bfs(Asub,1)[1][5] == 4
    end

    @testset "lollipop_graph" begin
        @test_throws ArgumentError lollipop_graph(-1)
        @test_throws ArgumentError lollipop_graph(-1,-1)
        @test_throws ArgumentError lollipop_graph(5,-1,Val{true})

        @test is_connected(lollipop_graph(5))
        @test is_connected(lollipop_graph(5, Val{true})[1])
        @test is_connected(lollipop_graph(5, 10))
        @test is_connected(lollipop_graph(10, 5))

        @test maximum(bfs(lollipop_graph(5,10),1)[1]) == 6
        @test maximum(bfs(lollipop_graph(10,5),1)[1]) == 11

        G,xy = lollipop_graph(6, 5, Val{true})
        filt = vec(all(xy .<= 0, dims = 2))
        Asub = sparse_transpose(G)[filt,filt]
        # Asub should be an tail (or a line graph
        @test bfs(Asub,1)[1][6] == 5
        Asub = sparse_transpose(G)[.~filt,.~filt]
        @test nnz(Asub) == 5*4

    end

    @testset "partial_duplication" begin

        A = sprand(100,100,.2)
        A = max.(A,A')
        B = MatrixNetwork(A)
        @test_throws ArgumentError partial_duplication(B,-1,.1)
        @test_throws ArgumentError partial_duplication(B,1,-.1)
        @test_throws ArgumentError partial_duplication(B,1,1.1)
        @test_throws ArgumentError partial_duplication(MatrixNetwork(sprand(20,10,.1)),1,.1)

        #check edge list creation and conversion
        C,dup_vertices = partial_duplication(B,0,1.0)
        @test length(dup_vertices) == 0
        @test is_undirected(C)
        @test sparse(C) == A

        @inferred partial_duplication(B,100,.5)

        #check symmetry is preserved
        C,dup_vertices = partial_duplication(B,100,.5)
        @test is_undirected(C)
        @test length(dup_vertices) == 100

        C,_ = partial_duplication(B,100,.1)
        @test is_undirected(C)

        #check extremal proability case
        C,_ = partial_duplication(B,100,0.0)
        @test norm(sparse(C),2) == norm(A,2) # no new entries added

    end

    @testset "_get_outedges" begin

        A = sprand(100,100,.2)
        A = max.(A,A')
        B = MatrixNetwork(A)

        @test_throws ArgumentError MatrixNetworks._get_outedges(B,-1)
        @test_throws ArgumentError MatrixNetworks._get_outedges(B,101)
        @inferred MatrixNetworks._get_outedges(B,50)

        row_idx = rand(1:100)
        SA_Is, SA_Vs = findnz(A[row_idx,:])
        MN_Is, MN_Vs = MatrixNetworks._get_outedges(B,row_idx)

        @test MN_Is == SA_Is
        @test MN_Vs == SA_Vs

    end
end

max_stars_repo_name: nassarhuda/MatrixNetworks

max_issues_repo_name: nassarhuda/MatrixNetworks

max_forks_repo_name: nassarhuda/MatrixNetworks


text:
#
# misc.jl --
#
# Implement reading/writing of FITS data from/to FITS files.
#
#------------------------------------------------------------------------------
#
# This file is part of OIFITS.jl which is licensed under the MIT "Expat"
# License:
#
# Copyright (C) 2015: Éric Thiébaut.
#
#------------------------------------------------------------------------------

using FITSIO
using FITSIO.Libcfitsio

import FITSIO: TableHDU, ASCIITableHDU
import Base: read

const _EXTENSION = Dict("IMAGE" => :image_hdu,
                        "TABLE" => :ascii_table,
                        "BINTABLE" => :binary_table)

# Guess HDU type (no warranty to work for primary HDU nor for incomplete
# header).
function get_hdutype(hdr::FITSHeader)
    if haskey(hdr, "XTENSION")
       return get(_EXTENSION, uppercase(rstrip(hdr["XTENSION"])), :unknown)
    elseif haskey(hdr, "SIMPLE") && hdr["SIMPLE"] == true
        return :image_hdu
    else
        return :unknown
    end
end

# Convert low-level handle into high level HDU type.
function make_hdu(ff::FITSFile)
    fits_assert_open(ff)
    hdutype = fits_get_hdu_type(ff)
    n = fits_get_hdu_num(ff)
    hdutype == :image_hdu    ? ImageHDU(ff, n) :
    hdutype == :binary_table ? TableHDU(ff, n) :
    hdutype == :ascii_table  ? ASCIITableHDU(ff, n) :
    error("current FITS HDU is not a table")
end

# Low level version.
read_table(ff::FITSFile) = read_table(make_hdu(ff))

function read_table(hdu::Union{TableHDU,ASCIITableHDU})
    hdr = read_header(hdu)
    data = Dict{ASCIIString,Any}()
    ncols = get_integer(hdr, "TFIELDS", 0)
    for k in 1:ncols
        name = uppercase(strip(get_string(hdr, "TTYPE$k", "")))
        if haskey(data, name)
            warn("duplicate column name: \"$name\"")
            continue
        end
        data[name] = read_column(hdu.fitsfile, k)
        units = strip(get_string(hdr, "TUNIT$k", ""))
        if length(units) > 0
            data[name*".units"] = units
        end
    end
    return data
end

# Read the entire table from disk. (High level version.)
read(hdu::Union{TableHDU,ASCIITableHDU}) = read_table(hdu)


text:
import Cairo
using Compose
using Colors
using GeometryBasics
using StaticArrays
using FactoryRendering

# rendering practice
Compose.set_default_graphic_size(20cm,20cm)
# vtxs coordinates and shapes
vtx_shape = GeometryBasics.Ngon(SVector(map(i->0.6*Point(sin(i*π/2+π/4),cos(i*π/2+π/4),0.0),0:3)...))
vtxs = Vector{Point3{Float64}}()
s = [15,20]
for i in 1:s[1]
    for j in 1:s[2]
        push!(vtxs, Point3{Float64}(i,j,0.0))
    end
end
floor_model = EntityRenderModel(
    configs=vtxs,
    shapes = map(v->vtx_shape,vtxs),
    colors = HSVA(50,0,1,0.2),
    )
# robot coordinates and shapes
r_vtxs = vtxs[[1,28,137]]
r_shape = Circle(Point(0.0,0.0),0.3)
robot_model = FactoryRendering.EntityRenderModel(
    configs=r_vtxs,
    shapes = map(v->r_shape,r_vtxs),
    colors = map(v->HSVA(200,1,1,0.9),r_vtxs)
)
o_vtxs=vtxs[[28,137,41,267]]
object_model = FactoryRendering.EntityRenderModel(
    configs=o_vtxs,
    shapes = Circle(Point(0.0,0.0),0.15),
    colors = HSVA(0,0,0,0.9)
)
# OrthographicProjection
theta = π/3
azimuth = π/8
proj = FactoryRendering.OrthographicProjection(theta,azimuth)
# proj = FactoryRendering.BirdsEyeProjection()

unit_box = FactoryRendering.construct_unit_box(floor_model,proj,ones(2))
Compose.compose(context(units=unit_box),
    get_render_layer(object_model,proj),
    get_render_layer(robot_model,proj),
    get_render_layer(floor_model,proj),
    )


text:
module TestSystem

using Test, PIMD

β, m, V, n = 1.4, 2.3, harmonicPotential, 3
ω2 = n / β^2
quantum = OneDimensionalQuantumSystem(β, m, V)

@testset "Primitive" begin
    primitive = Primitive(quantum, n)
    @test primitive.A ≈ [2 -1 -1; -1 2 -1; -1 -1 2] * m * n / β^2
end

@testset "Staging" begin
    staging = Staging(quantum, n)
    @test staging.A ≈ [0 0 0; 0 ω2 0; 0 0 ω2] * staging.M
end

@testset "Normal" begin
    normal = Normal(quantum, n)
    @test normal.A ≈ [0 0 0; 0 ω2 0; 0 0 ω2] * normal.M
end

end


text:
using SortingLab
using Test
using CategoricalArrays
using DataFrames, Random
using SortingAlgorithms

N = 204900
K = 100

include("fsort.jl")
include("fsort-missing.jl")
include("fsortperm-missing.jl")

# test fsortperm
# a = rand(1:Int(N/K), N);
# ca = copy(a);
#
# @time ab = fsortperm(a, rev = true)
# @test ca == a
# @test issorted(a[ab], rev = true)
#
# @time ab = fsortperm(a);
# @test ca == a
# @test issorted(a[ab])
#
# a_2048 = rand(1:2048, N)
# @time ab = fsortperm(a_2048);
# @test issorted(a_2048[ab])
#
# @time ab = fsortperm(a_2048, rev = true);
# @test issorted(a_2048[ab], rev = true)

# categorical sort

pools = "id".*lpad.(1:100, 3, '0');
byvec = CategoricalArray{String, 1}(rand(UInt32(1):UInt32(length(pools)), N), CategoricalPool(pools, false));
# @benchmark byvec_sorted = fsort($byvec)
# @benchmark byvec_sorted = fsort($byvec)
byvec = compress(byvec);

# @benchmark fsort($byvec) samples = 50 seconds = 120
# BenchmarkTools.Trial:
#   memory estimate:  1002.61 KiB
#   allocs estimate:  325
#   --------------
#   minimum time:     1.148 ms (0.00% GC)
#   median time:      1.588 ms (0.00% GC)
#   mean time:        1.829 ms (7.72% GC)
#   maximum time:     9.459 ms (67.31% GC)
#   --------------
#   samples:          2724
#   evals/sample:     1
# @benchmark SortingLab.fsort2($byvec) samples = 50 seconds = 120

byvec_sorted = fsort(byvec);
@test issorted(byvec_sorted)

fsort!(byvec)
@test issorted(byvec)

# String sort

# const M=1000; const K=100;
svec1 = rand([randstring(rand(1:4)) for k in 1:N÷K], N);
@time res1 = radixsort(svec1)
@test issorted(res1)

# test sortperm
idx = fsortperm(svec1)
@test issorted(svec1[idx])

svec1 = rand([randstring(rand(1:8)) for k in 1:N÷K], N);
@time res1 = radixsort(svec1)
@test issorted(res1)

# test sortperm
idx = fsortperm(svec1)
@test issorted(svec1[idx])

svec1 = rand([randstring(rand(1:32)) for k in 1:N÷K], N);
@time res1 = radixsort(svec1)
@test issorted(res1)

@time res1 = radixsort(svec1, true)
@test issorted(res1, rev = true)

# test sortperm
idx = fsortperm(svec1)
@test issorted(svec1[idx])

@time sort!(svec1);
@test issorted(svec1)

svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:N÷K], N);
@time radixsort!(svec1);
@test issorted(svec1)

svec1 = rand([string(rand(Char.(32:126), rand(1:16))...) for k in 1:N÷K], N);
@time radixsort!(svec1);
@test issorted(svec1)

svec1 = rand([string(rand(Char.(32:126), rand(1:24))...) for k in 1:N÷K], N);
@time radixsort!(svec1);
@test issorted(svec1)

svec1 = rand([string(rand(Char.(32:126), rand(1:32))...) for k in 1:N÷K], N);
@time radixsort!(svec1);
@test issorted(svec1)


text:
using GadgetUnits
using GadgetIO

function readSingleCRShockDataFromOutputFile(file::String)

        # read file into memory
        f = open(file)
        lines = readlines(f)
        close(f)

        # filter only relevant lines. Output of every line starts with "CR DATA".
        g = occursin.("CR DATA", lines)
        lines = lines[g]

        # get number of output lines
        N = length(lines)

        # init datatype for analyze
        cr = CRShockData(N)

        for i ∈ 1:N
                line_content            = split(lines[i])
                cr.dt[i]                = parse(Float64,line_content[3])
                cr.Mach[i]              = parse(Float64,line_content[4])
                cr.Shock_Speed[i]       = parse(Float64,line_content[5])
                cr.Shock_Compress[i]    = parse(Float64,line_content[6])
                cr.Shock_Energy_In[i]   = parse(Float64,line_content[7])
                cr.Shock_Energy_Real[i] = parse(Float64,line_content[8])
                cr.Energy_P[i]          = parse(Float64,line_content[9])
                cr.Energy_e[i]          = parse(Float64,line_content[10])
        end

        return cr
end

"""
    getCRMomentumDistributionFromPartID( snap_file::String, ID::Integer;
                                         pmin::Real=1.0, pmax::Real=1.0e6,
                                         Nbins::Integer=0, mode::Int64=3)

Reads the spectra from a single SPH particle via the particle ID.
"""
function getCRMomentumDistributionFromPartID(snap_file::String, ID::Integer;
                                             pmin::Real=1.0, pmax::Real=1.0e6,
                                             Nbins::Integer=0, mode::Int64=3,
                                             protons::Bool=true, electrons::Bool=true)

    h = head_to_obj(snap_file)
    info = read_info(snap_file)

    if info == 1
        if Nbins == 0
            error("Can't read spectrum! No info block present!\nSupply number of momentum bins to proceed!")
        else
            info = Array{InfoLine,1}(undef,7)
            info[1] = InfoLine("ID",    UInt32, Int32(1),     [1, 0, 0, 0, 0, 0])
            info[2] = InfoLine("CRpN", Float32, Int32(Nbins), [1, 0, 0, 0, 0, 0])
            info[3] = InfoLine("CRpS", Float32, Int32(Nbins), [1, 0, 0, 0, 0, 0])
            info[4] = InfoLine("CRpC", Float32, Int32(1),     [1, 0, 0, 0, 0, 0])
            info[5] = InfoLine("CReN", Float32, Int32(Nbins), [1, 0, 0, 0, 0, 0])
            info[6] = InfoLine("CReS", Float32, Int32(Nbins), [1, 0, 0, 0, 0, 0])
            info[7] = InfoLine("CReC", Float32, Int32(1),     [1, 0, 0, 0, 0, 0])
        end
    end

    # read block positions to speed up IO
    block_positions = GadgetIO.get_block_positions(snap_file)

    id = read_block(snap_file, "ID",
                    info=info[getfield.(info, :block_name) .== "ID"][1],
                    parttype=0, block_position=block_positions["ID"])

    # select the position of the requested ID
    part = findfirst( id .== UInt32(ID) )[1]

    # protons
    if protons
        CRpN = Float64.(
            read_block(snap_file, "CRpN",
                        info=info[getfield.(info, :block_name) .== "CRpN"][1],
                        parttype=0, block_position=block_positions["CRpN"])[:,part]
                )
        CRpS = read_block(snap_file, "CRpS",
                        info=info[getfield.(info, :block_name) .== "CRpS"][1],
                        parttype=0, block_position=block_positions["CRpS"])[:,part]
        CRpC = read_block(snap_file, "CRpC",
                        info=info[getfield.(info, :block_name) .== "CRpC"][1],
                        parttype=0, block_position=block_positions["CRpC"])[part]

        Nbins = size(CRpS,1)

    end

    # electrons
    if electrons
        CReN = Float64.(
            read_block(snap_file, "CReN",
                        info=info[getfield.(info, :block_name) .== "CReN"][1],
                        parttype=0, block_position=block_positions["CReN"])[:,part]
                )
        CReS = read_block(snap_file, "CReS",
                        info=info[getfield.(info, :block_name) .== "CReS"][1],
                        parttype=0, block_position=block_positions["CReS"])[:,part]
        CReC = read_block(snap_file, "CReC",
                        info=info[getfield.(info, :block_name) .== "CReC"][1],
                        parttype=0, block_position=block_positions["CReC"])[part]

        Nbins = size(CReS,1)
    end

    par = CRMomentumDistributionConfig(pmin, pmax, Nbins, mode)

    if protons && electrons
        return CRMomentumDistribution( CRpN, CRpS, CRpC, par.pmin, par.pmax, par.mc_e ),
               CRMomentumDistribution( CReN, CReS, CReC, par.pmin, par.pmax, par.mc_p )
    elseif protons
        return CRMomentumDistribution( CRpN, CRpS, CRpC, par.pmin, par.pmax, par.mc_e )
    elseif electrons
        return CRMomentumDistribution( CReN, CReS, CReC, par.pmin, par.pmax, par.mc_p )
    end
end

"""
    write_crp_cre_to_txt( t::Vector{<:Real}, CRp::CRMomentumDistribution, CRe::CRMomentumDistribution,
                          output_file::String )

Write CR Proton and Electron spectra for a series of time steps `t` to a txt file.
"""
function write_crp_cre_to_txt(t::Vector{<:Real}, CRp::CRMomentumDistribution, CRe::CRMomentumDistribution,
                              output_file::String)

    data = Matrix{Float64}(undef, length(t), 1+4*length(CRp[1].norm))

    for i = 1:length(t)
        data[i,:] = [t[i] CRp[i].bound[1:end-1]' CRp[i].norm' CRe[i].bound[1:end-1]' CRe[i].norm' ]
    end

    writedlm(output_file, data)
end

"""
    read_crp_cre_from_txt(filename::String)

Read CR Proton and Electron spectra from a txt file.
"""
function read_crp_cre_from_txt(filename::String)

    data = readdlm(filename)

    N = size(data,1)
    Nbins = Int64((size(data,2) - 1) / 8)

    CRp = Array{CRMomentumDistribution,1}(undef,N)
    CRe = Array{CRMomentumDistribution,1}(undef,N)

    t = data[:,1]

    for i = 1:N
        CRp[i] = CRMomentumDistribution([data[i,2:2Nbins+1]; data[i,2Nbins+1]], data[i,2Nbins+2:4Nbins+1])
        CRe[i] = CRMomentumDistribution([data[i,4Nbins+2:6Nbins+1]; data[i,6Nbins+1]], data[i,6Nbins+2:8Nbins+1])
    end

    return t, CRp, CRe
end

function write_cr_to_txt(t::Vector{<:Real}, CR::Vector{CRMomentumDistribution}, output_file::String)

    data = Matrix{Float64}(undef, length(t), 1+2*length(CR[1].norm))

    for i = 1:length(t)
        data[i,:] = [t[i] CR[i].bound[1:end-1]' CR[i].norm' ]
    end

    writedlm(output_file, data)
end

function read_cr_from_txt(fi)

    data = readdlm(fi)

    N = size(data,1)
    Nbins = Int64((size(data,2) - 1) / 4)

    CR = Array{CRMomentumDistribution,1}(undef,N)

    t = data[:,1]

    for i = 1:N
        CR[i] = CRMomentumDistribution([data[i,2:2Nbins+1]; data[i,2Nbins+1]], data[i,2Nbins+2:4Nbins+1])
    end

    return t, CR
end


text:
export addhook!, delhook!, gethook, sethook

struct AgentInterface
  agent::Agent
  hooks::Dict{String,Tuple{String,NTuple{N,String} where N}}
end
AgentInterface(agent) = AgentInterface(agent, Dict{String,Tuple{String,NTuple{N,String} where N}}())

# Adds a hook to an interface
# TODO: Validate that path exists
function addhook!(interface::AgentInterface, name::String, obj::String, path::Tuple)
  @assert !haskey(interface.hooks, name) "Hook $name already exists"
  interface.hooks[name] = (name, map(param->(param isa Val ? param : Val{param}()), path))
end

# Deletes a hook from an interface
function delhook!(interface::AgentInterface, name::String)
  @assert haskey(interface.hooks, name) "Hook $name does not exist"
  Base.delete!(interface.hooks, name)
end

# Gets the value of an interface hook
# FIXME: Actually lookup the parameter
function gethook(interface::AgentInterface, name::String)
  interface.nodes[interface.hooks[name]]
end

# Sets the value of an interface hook
# FIXME: Actually lookup the parameter
sethook!(interface::AgentInterface, name::String, value) = (interface.hooks[name] = value)


text:
# BondBasisSelector(Bsel::ACE.SparseBasis;
#                   isym=:be, bond_weight = 1.0, env_weight = 1.0) =
#    ACE.CategorySparseBasis(isym, [:bond, :env];
#             maxorder = ACE.maxorder(Bsel),
#             p = Bsel.p,
#             weight = Bsel.weight,
#             maxlevels = Bsel.maxlevels,
#             minorder_dict = Dict( :bond => 1),
#             maxorder_dict = Dict( :bond => 1),
#             weight_cat = Dict(:bond => bond_weight, :env=> env_weight)
#          )

# function SymmetricBond_basis(ϕ::ACE.AbstractProperty, env::ACE.BondEnvelope, Bsel::ACE.SparseBasis; RnYlm = nothing, bondsymmetry=nothing, kwargs...)
#    BondSelector =  BondBasisSelector(Bsel; kwargs...)
#    if RnYlm === nothing
#        RnYlm = RnYlm_1pbasis(;   r0 = ACE.cutoff_radialbasis(env),
#                                            rin = 0.0,
#                                            trans = PolyTransform(2, ACE.cutoff_radialbasis(env)),
#                                            pcut = 2,
#                                            pin = 0,
#                                            kwargs...
#                                        )
#    end
#    filterfun = _->true
#    if bondsymmetry == "Invariant"
#       filterfun = ACE.EvenL(:be, [:bond])
#    end
#    if bondsymmetry == "Covariant"
#       filterfun = x -> !(ACE.EvenL(:be, [:bond])(x))
#    end
#    Bc = ACE.Categorical1pBasis([:bond, :env]; varsym = :be, idxsym = :be )
#    B1p =  Bc * RnYlm * env
#    return ACE.SymmetricBasis(ϕ, B1p, BondSelector; filterfun = filterfun)
# end

import ACE: AbstractSparseBasis, maxorder, Prodb, Onepb, OneParticleBasis,
            degree
import Base: filter

struct SparseBondBasis <: AbstractSparseBasis
   maxorder::Int
   weight::Dict
   maxlevels::Dict
   p::Float64
   besym::Symbol
   bondsym::Symbol
   envsym::Symbol
   ksym::Symbol
   lsym::Symbol
   weight_cat::Dict
end

@noinline function SparseBondBasis(;
               besym = :be, bondsym = :bond, envsym = :env, ksym = :k,
               lsym = :l,
               maxorder = nothing,
               p = 1,
               weight = Dict{Symbol, Float64}(),
               default_maxdeg = nothing,
               maxlevels = nothing,
               weight_cat = Dict(bondsym => 1.0, envsym => 1.0),
               )
   @show maxorder, default_maxdeg
   if (default_maxdeg != nothing) && (maxlevels == nothing)
      return SparseBondBasis(maxorder, weight,
                          Dict("default" => default_maxdeg),
                          p, besym, bondsym, envsym, ksym, lsym, weight_cat)
   elseif (default_maxdeg == nothing) && (maxlevels != nothing)
      return SparseBondBasis(maxorder, weight, maxlevels,
                          p, besym, bondsym, envsym, ksym, lsym, weight_cat)
   else
      @error """Either both or neither optional arguments `maxlevels` and
                `default_maxdeg` were provided. To avoid ambiguity ensure that
                exactly one of these arguments is provided."""
   end
end

function Base.filter(b::Onepb, Bsel::SparseBondBasis, basis::OneParticleBasis)
   d = degree(b, basis)
   k = b[Bsel.ksym]
   return (k == 1 && b[Bsel.besym] == Bsel.envsym) ||
            (d == k-1 && b[Bsel.besym] == Bsel.bondsym)
end

function Base.filter(bb::Prodb, Bsel::SparseBondBasis, basis::OneParticleBasis)
   if length(bb) == 0; return true; end
   has1bond = count((b[Bsel.besym] == Bsel.bondsym) for b in bb) == 1
   isinvariant = sum( b[Bsel.lsym] for b in bb ) == 0
   return has1bond && isinvariant
end

# maxorder and maxlevel are inherited from the abstract interface

level(b::Union{Prodb, Onepb}, Bsel::SparseBondBasis, basis::OneParticleBasis) =
      cat_weighted_degree(b, Bsel, basis)

# Category-weighted degree function
cat_weighted_degree(b::Onepb, Bsel::SparseBondBasis, basis::OneParticleBasis) =
      degree(b, basis, Bsel.weight) * Bsel.weight_cat[getproperty(b, Bsel.isym)]

cat_weighted_degree(bb::Prodb, Bsel::SparseBondBasis, basis::OneParticleBasis) = (
      length(bb) == 0 ? 0.0
                      : norm(cat_weighted_degree.(bb, Ref(Bsel), Ref(basis)), Bsel.p)
      )


text:
# Julia wrapper for header: /usr/include/curl/curl.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@c Int32 curl_strequal (Ptr{Uint8}, Ptr{Uint8}) libcurl
@c Int32 curl_strnequal (Ptr{Uint8}, Ptr{Uint8}, size_t) libcurl
@c CURLFORMcode curl_formadd (Ptr{Ptr{Void}}, Ptr{Ptr{Void}}) libcurl
@c Int32 curl_formget (Ptr{Void}, Ptr{None}, curl_formget_callback) libcurl
@c None curl_formfree (Ptr{Void},) libcurl
@c Ptr{Uint8} curl_getenv (Ptr{Uint8},) libcurl
@c Ptr{Uint8} curl_version () libcurl
@c Ptr{Uint8} curl_easy_escape (Ptr{CURL}, Ptr{Uint8}, Int32) libcurl
@c Ptr{Uint8} curl_escape (Ptr{Uint8}, Int32) libcurl
@c Ptr{Uint8} curl_easy_unescape (Ptr{CURL}, Ptr{Uint8}, Int32, Ptr{Int32}) libcurl
@c Ptr{Uint8} curl_unescape (Ptr{Uint8}, Int32) libcurl
@c None curl_free (Ptr{None},) libcurl
@c CURLcode curl_global_init (Int32,) libcurl
@c CURLcode curl_global_init_mem (Int32, curl_malloc_callback, curl_free_callback, curl_realloc_callback, curl_strdup_callback, curl_calloc_callback) libcurl
@c None curl_global_cleanup () libcurl
@c Ptr{Void} curl_slist_append (Ptr{Void}, Ptr{Uint8}) libcurl
@c None curl_slist_free_all (Ptr{Void},) libcurl
@c time_t curl_getdate (Ptr{Uint8}, Ptr{time_t}) libcurl
@c Ptr{CURLSH} curl_share_init () libcurl
@c CURLSHcode curl_share_setopt (Ptr{CURLSH}, CURLSHoption) libcurl
@c CURLSHcode curl_share_cleanup (Ptr{CURLSH},) libcurl
@c Ptr{curl_version_info_data} curl_version_info (CURLversion,) libcurl
@c Ptr{Uint8} curl_easy_strerror (CURLcode,) libcurl
@c Ptr{Uint8} curl_share_strerror (CURLSHcode,) libcurl
@c CURLcode curl_easy_pause (Ptr{CURL}, Int32) libcurl
@c Ptr{CURL} curl_easy_init () libcurl
@c CURLcode curl_easy_setopt (Ptr{CURL}, CURLoption) libcurl
@c CURLcode curl_easy_perform (Ptr{CURL},) libcurl
@c None curl_easy_cleanup (Ptr{CURL},) libcurl
@c CURLcode curl_easy_getinfo (Ptr{CURL}, CURLINFO) libcurl
@c Ptr{CURL} curl_easy_duphandle (Ptr{CURL},) libcurl
@c None curl_easy_reset (Ptr{CURL},) libcurl
@c CURLcode curl_easy_recv (Ptr{CURL}, Ptr{None}, size_t, Ptr{size_t}) libcurl
@c CURLcode curl_easy_send (Ptr{CURL}, Ptr{None}, size_t, Ptr{size_t}) libcurl
@c Ptr{CURLM} curl_multi_init () libcurl
@c CURLMcode curl_multi_add_handle (Ptr{CURLM}, Ptr{CURL}) libcurl
@c CURLMcode curl_multi_remove_handle (Ptr{CURLM}, Ptr{CURL}) libcurl
@c CURLMcode curl_multi_fdset (Ptr{CURLM}, Ptr{fd_set}, Ptr{fd_set}, Ptr{fd_set}, Ptr{Int32}) libcurl
@c CURLMcode curl_multi_perform (Ptr{CURLM}, Ptr{Int32}) libcurl
@c CURLMcode curl_multi_cleanup (Ptr{CURLM},) libcurl
@c Ptr{CURLMsg} curl_multi_info_read (Ptr{CURLM}, Ptr{Int32}) libcurl
@c Ptr{Uint8} curl_multi_strerror (CURLMcode,) libcurl
@c CURLMcode curl_multi_socket (Ptr{CURLM}, curl_socket_t, Ptr{Int32}) libcurl
@c CURLMcode curl_multi_socket_action (Ptr{CURLM}, curl_socket_t, Int32, Ptr{Int32}) libcurl
@c CURLMcode curl_multi_socket_all (Ptr{CURLM}, Ptr{Int32}) libcurl
@c CURLMcode curl_multi_timeout (Ptr{CURLM}, Ptr{Int32}) libcurl
@c CURLMcode curl_multi_setopt (Ptr{CURLM}, CURLMoption) libcurl
@c CURLMcode curl_multi_assign (Ptr{CURLM}, curl_socket_t, Ptr{None}) libcurl


text:
using StatsBase

include("functions_align.jl")

function bootstrap(path_to_fasta, n_rep, out_prefix)
    # get the dimensions of the alignment
    array_dim = get_alignment_dimensions(path_to_fasta)
    l = array_dim[1] # number of sites in each sequence
    n = array_dim[2] # number of sequences in the fasta file

    for bootstrap in 1:n_rep
        open(out_prefix * string(bootstrap) * ".fasta", "w") do io

            # V is the sample of sites:
            V = sample(collect(1:l), (l,1), replace = true)

            records = Channel{fasta_record}((channel_arg) -> read_fasta_alignment(path_to_fasta, channel_arg))

            for record in records
                id = record.description
                sequence = record.seq

                # A is an array of one fasta entry's nucleotides
                A = Array{Char,1}(undef, l)
                # BS is the same array, but sampled with replacement (i.e., one bootstrap)
                BS = Array{Char,1}(undef, l)

                counter = 1
                for nuc in sequence
                    A[counter] = nuc
                    counter+=1
                end

                for i in 1:length(V)
                    BS[i] = A[V[i]]
                end

                new_seq = join(BS)

                println(io, ">" * id)
                println(io, new_seq)
            end

        close(io)
        end

    end
end

max_stars_repo_name: cov-ert/julialign

max_issues_repo_name: cov-ert/julialign

max_forks_repo_name: cov-ert/julialign


text:
# TODO: two separate servers should work okay
module TestServer
  using FactCheck
  using Connection
  using Server
  using Client
  using Config
  using Event
  include("Helper.jl")

  facts("Tests server creation and stop without run") do
    con = Server.create(IP, PORT)
    Server.stop(con)
  end
  facts("Tests two servers creation and stop without run") do
    con1 = Server.create(IP, PORT)
    con2 = Server.create(IP, PORT + 1)
    Server.stop(con2)
    Server.stop(con1)
  end
  facts("Tests two servers creation and stop one without run") do
    con1 = Server.create(IP, PORT)
    con2 = Server.create(IP, PORT + 1)
    Server.stop(con2)
    @fact Server.isOk(con1) --> true
    Server.stop(con1)
  end
  facts("Tests server creation and start") do
    con = Server.create(IP, PORT)
    Server.run(con)
    Server.stop(con)
  end
  facts("Tests two servers on the same port creation") do
    local answer::Command = Command(0, [])
    msg  = "Hello"
    con1 = Server.create(IP, PORT)
    con2 = Server.create(IP, PORT)
    con3 = Client.create(IP, PORT)
    Event.on(con1.observer, Connection.EVENT_BEFORE_RESPONSE, (sock, cmd, ans)->answer = Command(cmd.fn, deepcopy(cmd.args)))

    @fact Server.isOk(con1) --> true
    @fact Server.isOk(con2) --> false

    Server.run(con1)
    Server.run(con2)
    Client.request(con3, 1, msg)
    wait(()->answer.fn !== 0)

    @fact answer.fn      --> 1
    @fact answer.args[1] --> msg
    Client.stop(con3)
    Server.stop(con2)
    Server.stop(con1)
  end
  facts("Tests stopping the server and create again") do
    con = Server.create(IP, PORT)
    @fact Server.isOk(con) --> true
    Server.stop(con)

    con = Server.create(IP, PORT)
    @fact Server.isOk(con) --> true
    Server.stop(con)
  end
  facts("Tests server receiving command") do
    #
    # server and client creation
    #
    local answer::Command        = Command(0, Array{Any,1}())
    local scon::ServerConnection = Server.create(IP, PORT)
    local ccon::ClientConnection = Client.create(IP, PORT)
    Event.on(scon.observer, Connection.EVENT_BEFORE_RESPONSE, (sock, cmd, ans)->answer = Command(cmd.fn, deepcopy(cmd.args)))
    Server.run(scon)
    Client.request(ccon, 1, 10) # function - 1, parameter - 10
    wait(()->answer.fn !== 0)
    @fact answer.fn      --> 1
    @fact answer.args[1] --> 10

    Client.stop(ccon)
    Server.stop(scon)
  end
  facts("Tests server receiving commands from two clients") do
    local answer::Array{Command, 1} = Command[]
    local scon::ServerConnection = Server.create(IP, PORT)
    local ccon1::ClientConnection = Client.create(IP, PORT)
    local ccon2::ClientConnection = Client.create(IP, PORT)

    Event.on(scon.observer, Connection.EVENT_BEFORE_RESPONSE, (sock, cmd, ans)->push!(answer, Command(cmd.fn, deepcopy(cmd.args))))
    Server.run(scon)
    Client.request(ccon1, 1, 10) # function - 1, parameter - 10
    Client.request(ccon2, 2, 20) # function - 2, parameter - 20
    #
    # This line is very important. yield() runs the interconnection
    # loop between client and server.
    #
    wait(()->length(answer) > 0)
    @fact length(answer) --> 2
    @fact (answer[1].fn === 1 || answer[1].fn === 2) --> true
    @fact (answer[1].args[1] === 10 || answer[1].args[1] === 20) --> true
    @fact (answer[2].fn === 1 || answer[2].fn === 2) --> true
    @fact (answer[2].args[1] === 10 || answer[2].args[1] === 20) --> true

    Client.stop(ccon2)
    Client.stop(ccon1)
    Server.stop(scon)
  end
  facts("Tests sending the message to stopped server") do
    local answer::Command        = Command(0, Array{Any,1}())
    local scon::ServerConnection = Server.create(IP, PORT)
    local ccon::ClientConnection = Client.create(IP, PORT)

    Event.on(scon.observer, Connection.EVENT_BEFORE_RESPONSE, (sock, cmd, ans)->answer = Command(cmd.fn, deepcopy(cmd.args)))
    Server.run(scon)
    Client.request(ccon, 1, 10) # function - 1, parameter - 10
    Server.stop(scon)
    Client.request(ccon, 2, 20) # function - 2, parameter - 20

    wait(()->!Server.isOk(scon))
    @fact answer.fn      --> 1
    @fact answer.args[1] --> 10

    Client.stop(ccon)
  end
  facts("Tests client recreation and send a message") do
    local answer::Command        = Command(0, Array{Any,1}())
    local scon::ServerConnection = Server.create(IP, PORT)
    local ccon::ClientConnection = Client.create(IP, PORT)

    Event.on(scon.observer, Connection.EVENT_BEFORE_RESPONSE, (sock, cmd, ans)->answer = Command(cmd.fn, deepcopy(cmd.args)))
    Server.run(scon)
    Client.request(ccon, 1, 10)
    Client.stop(ccon)
    ccon = Client.create(IP, PORT)
    Client.request(ccon, 2, 20)
    wait(()->answer.fn === 2)

    @fact answer.fn      --> 2
    @fact answer.args[1] --> 20

    Client.stop(ccon)
    Server.stop(scon)
  end
  facts("Tests big data send") do
    local answer::Command        = Command(0, Array{Any,1}())
    local scon::ServerConnection = Server.create(IP, PORT)
    local ccon::ClientConnection = Client.create(IP, PORT)
    local data::Array{Int, 1}    = zeros(Int, 1000000)

    Event.on(scon.observer, Connection.EVENT_BEFORE_RESPONSE, (sock, cmd, ans)->answer = Command(cmd.fn, deepcopy(cmd.args)))
    Server.run(scon)
    Client.request(ccon, 1, data)
    wait(()->answer.fn !== 0)

    @fact answer.fn --> 1
    @fact answer.args[1] == data --> true

    Client.stop(ccon)
    Server.stop(scon)
  end
  facts("Tests client without server creation") do
    con = Client.create(IP, PORT)
    @fact Client.request(con, 1, 10) --> false
    Client.stop(con)
  end
  facts("Tests pooling") do
    local scon::ServerConnection = Server.create(IP, PORT)
    local ccon::ClientConnection = Client.create(IP, PORT)
    local result::Int            = 0

    Event.on(scon.observer, Connection.EVENT_BEFORE_RESPONSE, (sock, cmd, ans)->result += cmd.args[1])
    Server.run(scon)
    Client.request(ccon, 1, 1)
    wait(()->result !== 0)
    Client.request(ccon, 1, 2)
    wait(()->result > 1)

    Client.request(ccon, 1, 3)
    Client.request(ccon, 1, 4)
    wait(()->result > 5)
    wait(()->result > 9)

    @fact result --> 10

    Client.stop(ccon)
    Server.stop(scon)
  end
  facts("Tests pooling from many clients") do
    local scon::ServerConnection  = Server.create(IP, PORT)
    local ccon1::ClientConnection = Client.create(IP, PORT)
    local ccon2::ClientConnection = Client.create(IP, PORT)
    local ccon3::ClientConnection = Client.create(IP, PORT)
    local ccon4::ClientConnection = Client.create(IP, PORT)
    local ccon5::ClientConnection = Client.create(IP, PORT)
    local result::Int             = 0

    Event.on(scon.observer, Connection.EVENT_BEFORE_RESPONSE, (sock, cmd, ans)->result += cmd.args[1])
    Server.run(scon)
    Client.request(ccon1, 1, 1)
    wait(()->result > 0)
    Client.request(ccon2, 1, 2)
    wait(()->result > 2)
    Client.request(ccon3, 1, 3)
    wait(()->result > 5)
    Client.request(ccon4, 1, 1)
    wait(()->result > 6)
    Client.request(ccon5, 1, 3)
    wait(()->result > 9)

    @fact result --> 10

    result = 0
    Client.request(ccon3, 1, 3)
    Client.request(ccon5, 1, 4)
    Client.request(ccon1, 1, 2)
    Client.request(ccon4, 1, 1)
    Client.request(ccon2, 1, 1)
    wait(()->result > 10)

    @fact result --> 11

    Client.stop(ccon5)
    Client.stop(ccon4)
    Client.stop(ccon3)
    Client.stop(ccon2)
    Client.stop(ccon1)
    Server.stop(scon)
  end
  facts("Tests rending a request after client has stopped") do
    local result::Int = 0
    local res::Bool
    scon = Server.create(IP, PORT)
    ccon = Client.create(IP, PORT)

    Event.on(scon.observer, Connection.EVENT_BEFORE_RESPONSE, (sock, cmd, ans)->result = cmd.args[1])
    Server.run(scon)
    Client.stop(ccon)
    res = Client.request(ccon, 1, 10)

    @fact res    --> false
    @fact result --> 0

    Server.stop(scon)
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
include("../gcc_common.jl")
build_and_upload_gcc(v"5.2.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil