text:
# This file contains functions for the Revenue Efficiency DEA model
"""
    RevenueDEAModel
An data structure representing a revenue DEA model.
"""
struct RevenueDEAModel <: AbstractRevenueDEAModel
    n::Int64
    m::Int64
    s::Int64
    rts::Symbol
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
end

"""
    dearevenue(X, Y, P)
Compute revenue efficiency using data envelopment analysis for
inputs `X`, outputs `Y` and price of outputs `P`.

# Optional Arguments
- `rts=:VRS`: chooses variable returns to scale. For constant returns to scale choose `:CRS`.

# Examples
```jldoctest
julia> X = [5 3; 2 4; 4 2; 4 8; 7 9.0];

julia> Y = [7 4; 10 8; 8 10; 5 4; 3 6.0];

julia> P = [3 2; 3 2; 3 2; 3 2; 3 2.0];

julia> dearevenue(X, Y, P)
Revenue DEA Model
DMUs = 5; Inputs = 2; Outputs = 2
Orientation = Output; Returns to Scale = VRS
──────────────────────────────────
    Revenue  Technical  Allocative
──────────────────────────────────
1  0.644444   0.777778    0.828571
2  1.0        1.0         1.0
3  1.0        1.0         1.0
4  0.5        0.5         1.0
5  0.456522   0.6         0.76087
──────────────────────────────────
```
"""
function dearevenue(X::Matrix, Y::Matrix, P::Matrix; rts::Symbol = :VRS)::RevenueDEAModel
    # Check parameters
    nx, m = size(X)
    ny, s = size(Y)

    np, sp = size(P)

    if nx != ny
        error("number of observations is different in inputs and outputs")
    end
    if np != ny
        error("number of observations is different in output prices and outputs")
    end
    if sp != s
        error("number of output prices and outputs is different")
    end

    # Compute efficiency for each DMU
    n = nx

    Yefficient = zeros(n,s)
    refficiency = zeros(n)
    rlambdaeff = spzeros(n, n)

    for i=1:n
        # Value of inputs and outputs to evaluate
        x0 = X[i,:]
        p0 = P[i,:]

        # Create the optimization model
        deamodel = Model(with_optimizer(GLPK.Optimizer))
        @variable(deamodel, Yeff[1:s])
        @variable(deamodel, lambda[1:n] >= 0)

        @objective(deamodel, Max, sum(p0[j] .* Yeff[j] for j in 1:s))

        @constraint(deamodel, [j in 1:m], sum(X[t,j] * lambda[t] for t in 1:n) <= x0[j])
        @constraint(deamodel, [j in 1:s], sum(Y[t,j] * lambda[t] for t in 1:n) >= Yeff[j])

        # Add return to scale constraints
        if rts == :CRS
            # No contraint to add for constant returns to scale
        elseif rts == :VRS
            @constraint(deamodel, sum(lambda) == 1)
        else
            error("Invalid returns to scale $rts. Returns to scale should be :CRS or :VRS")
        end

        # Optimize and return results
        JuMP.optimize!(deamodel)

        Yefficient[i,:]  = JuMP.value.(Yeff)
        rlambdaeff[i,:] = JuMP.value.(lambda)

    end

    # Revenue, technical and allocative efficiency
    refficiency  = vec( sum(P .* Y, dims = 2) ./ sum(P .* Yefficient, dims = 2) )
    techefficiency = 1 ./ efficiency(dea(X, Y, orient = :Output, rts = rts, slack = false))
    allocefficiency = refficiency ./ techefficiency
    return RevenueDEAModel(n, m, s, rts, refficiency, rlambdaeff, techefficiency, allocefficiency)

end

function dearevenue(X::Vector, Y::Matrix, P::Matrix, rts::Symbol = :VRS)::RevenueDEAModel
    X = X[:,:]
    return dearevenue(X, Y, P, rts = rts)
end

function dearevenue(X::Matrix, Y::Vector, P::Vector; rts::Symbol = :VRS)::RevenueDEAModel
    Y = Y[:,:]
    P = P[:,:]
    return dearevenue(X, Y, P, rts = rts)
end

function dearevenue(X::Vector, Y::Vector, P::Vector; rts::Symbol = :VRS)::RevenueDEAModel
    X = X[:,:]
    Y = Y[:,:]
    P = P[:,:]
    return dearevenue(X, Y, P, rts = rts)
end

function Base.show(io::IO, x::RevenueDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    eff = efficiency(x)
    techeff = efficiency(x, :Technical)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "Revenue DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = Output")
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Revenue", "Technical", "Allocative"], ["$i" for i in 1:n]))

    else

    end
end


text:
struct NamedVector{K,L,T,V} <: FieldVector{L,T}
    nt::NamedTuple{K,V}
    NamedVector{K,L,T,V}(nt::NT) where {K,L,T,V<:Tuple,NT<:NamedTuple{K,Tuple{Vararg{T,L}}}} = new{K,L,T,V}(nt)
end
NamedVector{K,L,T,V}(t::Tuple{Vararg{T,L}}) where {K,L,T,V<:Tuple} = NamedVector{K,L,T,V}(NamedTuple{K}(t))
NamedVector{K,L,T,V}(x::T, xs::T...) where {K,L,T,V} = NamedVector{K,L,T,V}((x, xs...))

NamedVector{K,L,T}(t::V) where {K,L,T,V<:Tuple} = NamedVector{K,L,T,V}(t)
NamedVector{K,L,T}(t::V) where {K,L,T,V<:NamedTuple{K,V1}} where V1 = NamedVector{K,L,T,V1}(t)
NamedVector(nt::NamedTuple{K,V}) where {K,V<:Tuple{Vararg{Any,L}}} where L = NamedVector{K,L}(nt)
NamedVector{K}(t::V) where {K,V<:Tuple{Vararg{<:Any,L}}} where L = NamedVector{K,L}(t)
# Type promotion step...
function NamedVector{K,L}(t::V) where {K,L,V<:Tuple{Vararg{<:Any,L}}}
    t_uniform_type = map(x -> convert(promote_type(map(typeof, t)...), x), t)
    T = promote_type(map(typeof, t)...)
    NamedVector{K,L,T,NTuple{L,T}}(NamedTuple{K,NTuple{L,T}}(t_uniform_type))
end
function NamedVector{K,L}(nt::NT) where {K,L,NT<:NamedTuple{K,<:Tuple{Vararg{<:Any,L}}}}
    # TODO: why doesn't this compile away...
    NamedVector{K,L}(Tuple(nt))
end
function NamedVector{K,L}(v::SVector{L,T}) where {K,L,T}
    NamedVector{K,L,T,NTuple{L,T}}(NamedTuple{K,NTuple{L,T}}(Tuple(v)))
end
NamedVector(; kw...) = NamedVector{keys(kw),length(kw)}(values(kw))

Base.parent(a::NamedVector) = getfield(a, :nt)
Base.values(a::NamedVector) = values(parent(a))
Base.getproperty(a::NamedVector, x::Symbol) = getproperty(parent(a), x)
Base.@propagate_inbounds Base.getindex(a::NamedVector, i::Int) = getindex(parent(a), i)
Base.@propagate_inbounds Base.setindex!(a::NamedVector, x, i::Int) = (setindex!(parent(a), i, x); a)
# Base.map(f, as::Union{NamedVector{K},NamedTuple{K}}...) where K = begin
Base.map(f, a1::NamedVector{K,L}, as::NamedVector{K}...) where {K,L} = begin
    NamedVector{K,L}(map(f, map(SVector, (a1, as...))...))
end

_maybeparent(a::NamedVector) = parent(a)
_maybeparent(nt::NamedTuple) = nt

Base.tail(a::NamedVector{K,L,T}) where {K,L,T} = NamedVector(Base.tail(parent(a)))

ext: jl


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

@doc raw"""FlexPersistentVolumeSource represents a generic persistent volume resource that is provisioned/attached using an exec based plugin.

    IoK8sApiCoreV1FlexPersistentVolumeSource(;
        driver=nothing,
        fsType=nothing,
        options=nothing,
        readOnly=nothing,
        secretRef=nothing,
    )

    - driver::String : Driver is the name of the driver to use for this volume.
    - fsType::String : Filesystem type to mount. Must be a filesystem type supported by the host operating system. Ex. \"ext4\", \"xfs\", \"ntfs\". The default filesystem depends on FlexVolume script.
    - options::Dict{String, String} : Optional: Extra command options if any.
    - readOnly::Bool : Optional: Defaults to false (read/write). ReadOnly here will force the ReadOnly setting in VolumeMounts.
    - secretRef::IoK8sApiCoreV1SecretReference : Optional: SecretRef is reference to the secret object containing sensitive information to pass to the plugin scripts. This may be empty if no secret object is specified. If the secret object contains more than one secret, all secrets are passed to the plugin scripts.
"""
mutable struct IoK8sApiCoreV1FlexPersistentVolumeSource <: SwaggerModel
    driver::Any # spec type: Union{ Nothing, String } # spec name: driver
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    options::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: options
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    secretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretReference } # spec name: secretRef

    function IoK8sApiCoreV1FlexPersistentVolumeSource(;driver=nothing, fsType=nothing, options=nothing, readOnly=nothing, secretRef=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1FlexPersistentVolumeSource, Symbol("driver"), driver)
        setfield!(o, Symbol("driver"), driver)
        validate_property(IoK8sApiCoreV1FlexPersistentVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1FlexPersistentVolumeSource, Symbol("options"), options)
        setfield!(o, Symbol("options"), options)
        validate_property(IoK8sApiCoreV1FlexPersistentVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1FlexPersistentVolumeSource, Symbol("secretRef"), secretRef)
        setfield!(o, Symbol("secretRef"), secretRef)
        o
    end
end # type IoK8sApiCoreV1FlexPersistentVolumeSource

const _property_map_IoK8sApiCoreV1FlexPersistentVolumeSource = Dict{Symbol,Symbol}(Symbol("driver")=>Symbol("driver"), Symbol("fsType")=>Symbol("fsType"), Symbol("options")=>Symbol("options"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("secretRef")=>Symbol("secretRef"))
const _property_types_IoK8sApiCoreV1FlexPersistentVolumeSource = Dict{Symbol,String}(Symbol("driver")=>"String", Symbol("fsType")=>"String", Symbol("options")=>"Dict{String, String}", Symbol("readOnly")=>"Bool", Symbol("secretRef")=>"IoK8sApiCoreV1SecretReference")
Base.propertynames(::Type{ IoK8sApiCoreV1FlexPersistentVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1FlexPersistentVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1FlexPersistentVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1FlexPersistentVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1FlexPersistentVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1FlexPersistentVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1FlexPersistentVolumeSource)
    (getproperty(o, Symbol("driver")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1FlexPersistentVolumeSource }, name::Symbol, val)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Base.Test

srand(123)

@testset "basic functions" begin
    @test I[1,1] == 1 # getindex
    @test I[1,2] == 0 # getindex
    @test I === I' # transpose
    @test ndims(I) == 2
    @test one(UniformScaling{Float32}) == UniformScaling(one(Float32))
    @test zero(UniformScaling{Float32}) == UniformScaling(zero(Float32))
    @test eltype(one(UniformScaling{Float32})) == Float32
    @test zero(UniformScaling(rand(Complex128))) == zero(UniformScaling{Complex128})
    @test one(UniformScaling(rand(Complex128))) == one(UniformScaling{Complex128})
    @test eltype(one(UniformScaling(rand(Complex128)))) == Complex128
    @test -one(UniformScaling(2)) == UniformScaling(-1)
    @test sparse(3I,4,5) == spdiagm(fill(3,4),0,4,5)
    @test sparse(3I,5,4) == spdiagm(fill(3,4),0,5,4)
end

@testset "istriu, istril, issymmetric, ishermitian, isapprox" begin
    @test istriu(I)
    @test istril(I)
    @test issymmetric(I)
    @test issymmetric(UniformScaling(complex(1.0,1.0)))
    @test ishermitian(I)
    @test !ishermitian(UniformScaling(complex(1.0,1.0)))
    @test UniformScaling(4.00000000000001) ≈ UniformScaling(4.0)
    @test UniformScaling(4.32) ≈ UniformScaling(4.3) rtol=0.1 atol=0.01
end

@testset "* and / with number" begin
    α = randn()
    @test α .* UniformScaling(1.0) == UniformScaling(1.0) .* α
    @test UniformScaling(α)./α == UniformScaling(1.0)
end

@test copy(UniformScaling(one(Float64))) == UniformScaling(one(Float64))
@test sprint(show,UniformScaling(one(Float32))) == "UniformScaling{Float32}\n1.0*I"

λ = complex(randn(),randn())
J = UniformScaling(λ)
@testset "transpose, conj, inv" begin
    @test ndims(J) == 2
    @test transpose(J) == J
    @test J*eye(2) == conj(J'eye(2)) # ctranpose (and A(c)_mul_B)
    @test I + I === UniformScaling(2) # +
    @test inv(I) == I
    @test inv(J) == UniformScaling(inv(λ))
    @test cond(I) == 1
    @test cond(J) == (λ ≠ zero(λ) ? one(real(λ)) : oftype(real(λ), Inf))
end

B = bitrand(2,2)
@test B + I == B + eye(B)
@test I + B == B + eye(B)

@testset "binary ops with matrices" begin
    let AA = randn(2, 2)
        for SS in (sprandn(3,3, 0.5), speye(Int, 3))
            @testset for atype in ("Array", "SubArray")
                if atype == "Array"
                    A = AA
                    S = SS
                else
                    A = view(AA, 1:2, 1:2)
                    S = view(SS, 1:3, 1:3)
                end

                @test @inferred(A + I) == A + eye(A)
                @test @inferred(I + A) == A + eye(A)
                @test @inferred(I - I) === UniformScaling(0)
                @test @inferred(B - I) == B - eye(B)
                @test @inferred(I - B) == eye(B) - B
                @test @inferred(A - I) == A - eye(A)
                @test @inferred(I - A) == eye(A) - A
                @test @inferred(I*J) === UniformScaling(λ)
                @test @inferred(B*J) == B*λ
                @test @inferred(J*B) == B*λ
                @test @inferred(I*A) !== A # Don't alias
                @test @inferred(I*S) !== S # Don't alias
                @test @inferred(A*I) !== A # Don't alias
                @test @inferred(S*I) !== S # Don't alias

                @test @inferred(S*J) == S*λ
                @test @inferred(J*S) == S*λ
                @test @inferred(A*J) == A*λ
                @test @inferred(J*A) == A*λ
                @test @inferred(J*ones(3)) == ones(3)*λ
                @test @inferred(λ*J) === UniformScaling(λ*J.λ)
                @test @inferred(J*λ) === UniformScaling(λ*J.λ)
                @test @inferred(J/I) === J
                @test @inferred(I/A) == inv(A)
                @test @inferred(A/I) == A
                @test @inferred(I/λ) === UniformScaling(1/λ)
                @test @inferred(I\J) === J

                if atype == "Array"
                    T = LowerTriangular(randn(3,3))
                else
                    T = LowerTriangular(view(randn(3,3), 1:3, 1:3))
                end
                @test @inferred(T + J) == full(T) + J
                @test @inferred(J + T) == J + full(T)
                @test @inferred(T - J) == full(T) - J
                @test @inferred(J - T) == J - full(T)
                @test @inferred(T\I) == inv(T)

                if atype == "Array"
                    T = LinAlg.UnitLowerTriangular(randn(3,3))
                else
                    T = LinAlg.UnitLowerTriangular(view(randn(3,3), 1:3, 1:3))
                end
                @test @inferred(T + J) == full(T) + J
                @test @inferred(J + T) == J + full(T)
                @test @inferred(T - J) == full(T) - J
                @test @inferred(J - T) == J - full(T)
                @test @inferred(T\I) == inv(T)

                if atype == "Array"
                    T = UpperTriangular(randn(3,3))
                else
                    T = UpperTriangular(view(randn(3,3), 1:3, 1:3))
                end
                @test @inferred(T + J) == full(T) + J
                @test @inferred(J + T) == J + full(T)
                @test @inferred(T - J) == full(T) - J
                @test @inferred(J - T) == J - full(T)
                @test @inferred(T\I) == inv(T)

                if atype == "Array"
                    T = LinAlg.UnitUpperTriangular(randn(3,3))
                else
                    T = LinAlg.UnitUpperTriangular(view(randn(3,3), 1:3, 1:3))
                end
                @test @inferred(T + J) == full(T) + J
                @test @inferred(J + T) == J + full(T)
                @test @inferred(T - J) == full(T) - J
                @test @inferred(J - T) == J - full(T)
                @test @inferred(T\I) == inv(T)

                @test @inferred(I\A) == A
                @test @inferred(A\I) == inv(A)
                @test @inferred(λ\I) === UniformScaling(1/λ)
            end
        end
    end
end

@testset "hcat and vcat" begin
    @test_throws ArgumentError hcat(I)
    @test_throws ArgumentError [I I]
    @test_throws ArgumentError vcat(I)
    @test_throws ArgumentError [I; I]
    @test_throws ArgumentError [I I; I]
    for T in (Matrix, SparseMatrixCSC)
        A = T(rand(3,4))
        B = T(rand(3,3))
        @test (hcat(A,2I))::T == hcat(A,2eye(3,3))
        @test (vcat(A,2I))::T == vcat(A,2eye(4,4))
        @test (hcat(I,3I,A,2I))::T == hcat(eye(3,3),3eye(3,3),A,2eye(3,3))
        @test (vcat(I,3I,A,2I))::T == vcat(eye(4,4),3eye(4,4),A,2eye(4,4))
        @test (hvcat((2,1,2),B,2I,I,3I,4I))::T ==
            hvcat((2,1,2),B,2eye(3,3),eye(6,6),3eye(3,3),4eye(3,3))
    end
end

max_stars_repo_name: odow/julia

max_issues_repo_name: odow/julia

max_forks_repo_name: odow/julia


text:
using RealNeuralNetworks
using RealNeuralNetworks.NodeNets.DBFs
using Test

seg = zeros(UInt32,(100,100,100))
seg[48:52,48:52, 1:end ] .= one(UInt32)
seg[49:52, 49:52, 48:52] .= one(UInt32)
seg[47:54, 47:54, 71:78] .= one(UInt32)

points = RealNeuralNetworks.NodeNets.PointArrays.from_seg(seg)
@testset "test dbf computation" begin
    dbf1 = DBFs.compute_DBF(points)
    #dbf2 = DBFs.compute_DBF(points, boundary_point_indexes)
    bin_im = DBFs.create_binary_image( seg )
    bin_im2 = DBFs.create_binary_image( points )
    @show size(bin_im)
    @show size(bin_im2)
    # the size of bin_im2 is smaller than bin_im
    #@test all(bin_im .== bin_im2)
    dbf3 = DBFs.compute_DBF(points, bin_im )
    #@show dbf1
    #@show dbf3
    # it seems that the dbf from binary image is not correct
    # map((x,y) -> @test_approx_eq_eps(x,y,1), dbf1, dbf3)
end


text:
using Gadfly, RDatasets

set_default_plot_size(5inch, 3inch)

iris = dataset("datasets", "iris")

fs = [ x->0.3x+1, x->0.3x+1.1 ]
cs = ["versicolor", "virginica"]
plot(iris,
  layer(x=:SepalLength, y=:SepalWidth, color=:Species),
  layer(fs, 4, 8, color=cs, order=2)
)


text:
abstract type AbstractContext end
abstract type AbstractNormaliseDataContext <: AbstractContext end
abstract type AbstractNormalisationMethod end
abstract type AbstractCovarianceMatrix end
abstract type AbstractFittingMethod end
abstract type AbstractObservations end
abstract type AbstractUncertainObservations <: AbstractObservations end
abstract type AbstractGauge end
abstract type AbstractOptimisationScheme end
abstract type AbstractLatentEllipseParametrisation end
abstract type AbstractEllipseParameterConversion end
abstract type AbstractOptimisationResult end

@with_kw struct OptimisationResult{T₁ <: AbstractVector} <: AbstractOptimisationResult
    minimiser::T₁ = SVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
    minimum::Float64 = 0.0
    info_str::String = ""
end

struct AlgebraicToGeometric <: AbstractEllipseParameterConversion end
struct GeometricToAlgebraic <: AbstractEllipseParameterConversion end

struct FirstLatentEllipseParametrisation <: AbstractLatentEllipseParametrisation end
struct SecondLatentEllipseParametrisation <: AbstractLatentEllipseParametrisation end

@with_kw struct ManualEstimation{T <: Union{AbstractVecOrMat, Nothing}} <: AbstractOptimisationScheme
    𝛉::T = nothing
end

struct IsotropicScalingTranslation <: AbstractNormalisationMethod end
struct IsotropicScaling <: AbstractNormalisationMethod end

struct DirectEllipseFit <: AbstractFittingMethod end

@with_kw struct LevenbergMarquardt{T <:  AbstractOptimisationScheme} <: AbstractOptimisationScheme
    max_iter::Integer = 10000
    max_func_eval::Integer = 500
    γ::Float64 = 1.2
    γ_inc::Float64 = 10
    γ_dec::Float64 = 0.1
    λ::Float64 = 0.01
    λ_lower_bound::Float64 = 1e-10
    tol_𝛉::Float64 = 1e-8             # change in parameters tolerance
    tol_cost::Float64 = 1e-12         # change in cost tolerance
    tol_Δ::Float64 = 1e-8             # change in update tolerance
    tol_∇::Float64 = 1e-12            # change in gradient tolerance
    seed::T = ManualEstimation()
end

"""
```
    GuaranteedEllipseFit(; parametrisation = SecondLatentEllipseParametrisation(), optimisation_scheme = LevenbergMarquardt())
```

Configures the choice of implicit parametrisation and the optimisation algorithm; currently only the default values are supported.
You can also choose to display summary information about the optimisation process in the console by setting `verbose = true`.

"""
@with_kw struct GuaranteedEllipseFit <: AbstractFittingMethod
    parametrisation::AbstractLatentEllipseParametrisation = SecondLatentEllipseParametrisation()
    optimisation_scheme::LevenbergMarquardt = LevenbergMarquardt()
    verbose::Bool = false
end

GuaranteedEllipseFit(𝛉::AbstractVector) = GuaranteedEllipseFit(optimisation_scheme = LevenbergMarquardt(seed = ManualEstimation(𝛉)))

struct UnitNorm <: AbstractGauge end
struct ToNormalisedSpace end
struct FromNormalisedSpace end

struct JacobianMatrix{T₁ <: Union{Nothing, AbstractObservations}, T₂ <: AbstractMatrix}
    observations::T₁
    jacobian::T₂
end


text:
file = open("input.dat", "w")
for i = 1:800
    write(file, rand(Float32, 1024, 1024))
end
close(file)


text:
export CuDeviceTexture

abstract type TextureInterpolationMode end
struct NearestNeighbour      <: TextureInterpolationMode end
struct LinearInterpolation   <: TextureInterpolationMode end
struct CubicInterpolation    <: TextureInterpolationMode end

abstract type TextureMemorySource end
struct ArrayMemory   <: TextureMemorySource end
struct LinearMemory  <: TextureMemorySource end

"""
    CuDeviceTexture{T,N,M,NC,I}

`N`-dimensional device texture with elements of type `T`. This type is the device-side
counterpart of [`CuTexture{T,N,P}`](@ref), and can be used to access textures using regular
indexing notation. If `NC` is true, indices used by these accesses should be normalized,
i.e., fall into the `[0,1)` domain. The `I` type parameter indicates the kind of
interpolation that happens when indexing into this texture. The source memory of the
texture is specified by the `M` parameter, either linear memory or a texture array.

Device-side texture objects cannot be created directly, but should be created host-side
using [`CuTexture{T,N,P}`](@ref) and passed to the kernal as an argument.

!!! warning
    Experimental API. Subject to change without deprecation.
"""
struct CuDeviceTexture{T,N,M<:TextureMemorySource,NC,I<:TextureInterpolationMode} <: AbstractArray{T,N}
    dims::Dims{N}
    handle::CUtexObject
end

Base.convert(::Type{CUtexObject}, t::CuDeviceTexture) = t.handle

## array interface

Base.elsize(::Type{<:CuDeviceTexture{T}}) where {T} = sizeof(T)

Base.size(tm::CuDeviceTexture) = tm.dims
Base.sizeof(tm::CuDeviceTexture) = Base.elsize(x) * length(x)

## low-level operations

# Source: NVVM IR specification 1.4
# NOTE: tex1Dfetch (integer coordinates) is unsupported, as it can be easily done using ldg

struct Vec4{T}
    a::T
    b::T
    c::T
    d::T
end

Base.Tuple(x::Vec4) = tuple(x.a, x.b, x.c, x.d)

for (dispatch_rettyp, julia_rettyp, llvm_rettyp) in
        ((Signed,        Vec4{UInt32},  :v4u32),
         (Unsigned,      Vec4{Int32},   :v4s32),
         (AbstractFloat, Vec4{Float32}, :v4f32))

    eltyp = Union{dispatch_rettyp, NTuple{<:Any,dispatch_rettyp}}

    # tex1D only supports array memory
    @eval tex(texObject::CuDeviceTexture{<:$eltyp,1,ArrayMemory}, x::Number) =
        Tuple(ccall($("llvm.nvvm.tex.unified.1d.$llvm_rettyp.f32"), llvmcall,
                    $julia_rettyp, (CUtexObject, Float32), texObject, x))

    # tex2D and tex3D supports all memories
    for dims in 2:3
        llvm_dim = "$(dims)d"
        julia_args = (:x, :y, :z)[1:dims]
        julia_sig = ntuple(_->Float32, dims)
        julia_params = ntuple(i->:($(julia_args[i])::Number), dims)

        @eval tex(texObject::CuDeviceTexture{<:$eltyp,$dims,}, $(julia_params...)) =
            Tuple(ccall($("llvm.nvvm.tex.unified.$llvm_dim.$llvm_rettyp.f32"), llvmcall,
                        $julia_rettyp, (CUtexObject, $(julia_sig...)), texObject, $(julia_args...)))
    end
end

## hardware-supported indexing

# we only support Float32 indices
@inline Base.getindex(t::CuDeviceTexture, idx::Vararg{Real,N}) where {N} =
    Base.getindex(t, ntuple(i->Float32(idx[i]), N)...)

@inline function Base.getindex(t::CuDeviceTexture{T,N,<:Any,true,I}, idx::Vararg{Float32,N}) where
                              {T,N,I<:Union{NearestNeighbour,LinearInterpolation}}
    # normalized coordinates range between 0 and 1, and can be used as-is
    vals = tex(t, idx...)
    return unpack(T, vals)
end

@inline function Base.getindex(t::CuDeviceTexture{T,N,<:Any,false,I}, idx::Vararg{Float32,N}) where
                              {T,N,I<:Union{NearestNeighbour,LinearInterpolation}}
    # non-normalized coordinates should be adjusted for 1-based indexing
    vals = tex(t, ntuple(i->idx[i]-0.5, N)...)
    return unpack(T, vals)
end

# unpack single-channel texture fetches as values, tuples otherwise
@inline unpack(::Type{T},           vals::NTuple) where T = unpack(T, vals[1])
@inline unpack(::Type{NTuple{1,T}}, vals::NTuple) where T = unpack(T, vals[1])
@inline unpack(::Type{NTuple{C,T}}, vals::NTuple) where {C,T} = ntuple(i->unpack(T, vals[i]), C)

@inline unpack(::Type{T}, val::T) where {T} = val
@inline unpack(::Type{T}, val::Real) where {T <: Integer} = unsafe_trunc(T, val)
@inline unpack(::Type{Float16}, val::Float32) = convert(Float16, val)

## cubic indexing (building on linear filtering)

# Source: GPU Gems 2, Chapter 20: Fast Third-Order Texture Filtering
#         CUDA sample: bicubicTextures

# cubic B-spline basis functions
w0(a::Float32) = (1.0f0/6.0f0)*(a*(a*(-a + 3.0f0) - 3.0f0) + 1.0f0)
w1(a::Float32) = (1.0f0/6.0f0)*(a*a*(3.0f0*a - 6.0f0) + 4.0f0)
w2(a::Float32) = (1.0f0/6.0f0)*(a*(a*(-3.0f0*a + 3.0f0) + 3.0f0) + 1.0f0)
w3(a::Float32) = (1.0f0/6.0f0)*(a*a*a)

# amplitude functions
g0(a::Float32) = w0(a) + w1(a)
g1(a::Float32) = w2(a) + w3(a)

# offset functions
# NOTE: +0.5 offset to compensate for CUDA linear filtering convention
h0(a::Float32) = -1.0f0 + w1(a) / (w0(a) + w1(a)) + 0.5f0
h1(a::Float32) = 1.0f0 + w3(a) / (w2(a) + w3(a)) + 0.5f0

@inline function Base.getindex(t::CuDeviceTexture{T,1,<:Any,false,CubicInterpolation},
                               x::Float32) where {T}
    x -= 1.0f0
    px = floor(x)   # integer position
    fx = x - px     # fractional position

    g0x = g0(fx)
    g1x = g1(fx)
    h0x = h0(fx)
    h1x = h1(fx)

    vals = g0x .* tex(t, px + h0x) .+ g1x .* tex(t, px + h1x)
    return (unpack(T, vals))
end

@inline function Base.getindex(t::CuDeviceTexture{T,2,<:Any,false,CubicInterpolation},
                               x::Float32, y::Float32) where {T}
    x -= 1.0f0
    y -= 1.0f0
    px = floor(x)   # integer position
    py = floor(y)
    fx = x - px     # fractional position
    fy = y - py

    g0x = g0(fx)
    g1x = g1(fx)
    h0x = h0(fx)
    h1x = h1(fx)
    h0y = h0(fy)
    h1y = h1(fy)

    vals = g0(fy) .* (g0x .* tex(t, px + h0x, py + h0y) .+
                      g1x .* tex(t, px + h1x, py + h0y)) .+
           g1(fy) .* (g0x .* tex(t, px + h0x, py + h1y) .+
                      g1x .* tex(t, px + h1x, py + h1y))
    return (unpack(T, vals))
end


text:
# create ground state
# imprint vortex
# evolve in real time
using Revise, ProjectedGPE, VortexDistributions

using PyPlot

function groundstate(γ,tf,Ncut)
#== start parameters ==#
sinfo = Params()
#fundamental constants/units
  ħ = 1.0545718e-34; kB = 1.38064852e-23; amu = 1.660339040e-27; Bohr = 5.29e-11
#Rb87 mass and scattering length
  m = 86.909180527*amu; as = 100*Bohr
#trap frequencies
  ωx = 2π; ωy = ωx; ωz = 0.
#choice of time, length, energy units
  t0 = 1.0/ωy; x0 = sqrt(ħ*t0/m); E0 = ħ/t0
#interactions and chemical potential, e.g., g  = (4*pi*ħ^2*as/m)*x0^3/E0 (dimensionless 3D)
  g = 0.1; μ  = 10.0
#damping parameters (dimensionless)
  #γ = 0.5
  ℳ  = 0.0
#time evolution parameters
  ti = 0.0; #tf = 10.0/γ  #evolve for 2 damping times
  Nt = 50; t = collect(linspace(ti,tf,Nt)); dt = 0.01π/μ #integrate step size [ - should have dt ≪ 2π/μ]
#
@pack sinfo = ωx,ωy,ωz,γ,ℳ,g,t0,x0,E0,μ,ti,tf,Nt,t,dt
#== end parameters ==#

#= Initialize cfield (dimensionless) =#
  basis = "Hermite"
  ecut = Ncut*ħ*ωy/E0
  Ω = [ωx; ωy]*t0
  cinfo = makecinfo(ecut,Ω,basis)
  tinfo = maketinfo(cinfo,4)
  @unpack en,P,M = cinfo

#= create initial state and extra fields =#
  c0  = randn(M...)+im*randn(M...); c0=P.*c0
  ψ0  = c2x(c0,tinfo) #initial condition
  ψ = similar(ψ0)  #a field to write to in place

#= PGPE (in place)       =#
#add spectral terms into Lgp!, or x-space terms into nlin!
  function Lgp!(dc,c,p,t)
      nlin!(dc,c,tinfo)
      dc .= P.*(-im*(1-im*γ)*((en - μ).*c .+ g*dc))
  end
  function nlin!(dc,c,tinfo)
      c2x!(ψ,c,tinfo)
      x2c!(dc,abs2.(ψ).*ψ,tinfo)
  end
#========================#

tspan = (t[1],t[end])
prob = ODEProblem(Lgp!,c0,tspan)
alg =    DP5() #abstol = 1e-3, reltol = 1e-3

println("Started evolution ...")
#@time sol = solve(prob,alg,dt=dt,saveat=t,alg_hints=[:stiff],save_everystep=false,dense=false);
@time sol = solve(prob,alg,dt=dt,saveat=t);
println("... Finished.")
return sinfo,cinfo,sol
end

#EVOLVE
sinfo,cinfo,sol = groundstate(0.3,10/γ,30.5)

function addvortex()
@unpack ħ,m,ωx,ωy,ωz,γ,ℳ,g,x0,t0,E0,μ,ti,tf,Nt,t,dt = sinfo
#@unpack M,Ω,ecut,P,en = cinfo

#= Plot solution =#
#lets be explicit about units:
Rx = sqrt(2μ*E0/m/ωx^2)/x0
Ry = sqrt(2μ*E0/m/ωy^2)/x0
yMax=1.5Ry
xMax=1.5Rx

Nx = 600
Ny = Nx
x = collect(linspace(-xMax,xMax,Nx))
y = collect(linspace(-yMax,yMax,Ny))
tinfop = maketinfoplot(cinfo,x,y)
ψ = c2x(sol[end],tinfop)

#= imprint vortex =#
xv = 0.2*Rx; yv = 0.0; σv = 1
testvort = [xv yv σv]
dx = x[2]-x[1];dy = y[2]-y[1]

ix = find(abs.(x-xv).<=dx);ix = ix[1]
iy = find(abs.(y-yv).<=dy);iy = iy[1]
#ξ = x0/(sqrt(μ))
ξ = 1/(sqrt(g*abs2.(ψ)[ix,iy]))

makevortex!(ψ,testvort,x,y,ξ)
#project back to cfield
return xgrid2c(ψ,x,y,cinfo)
end

c0 = addvortex()

function realtime(tf,Ncut)

#== start parameters ==#
sinfo = Params()
#fundamental constants/units
  ħ = 1.0545718e-34; kB = 1.38064852e-23; amu = 1.660339040e-27; Bohr = 5.29e-11
#Rb87 mass and scattering length
  m = 86.909180527*amu; as = 100*Bohr
#trap frequencies
  ωx = 2π; ωy = ωx; ωz = 0.
#choice of time, length, energy units
  t0 = 1.0/ωy; x0 = sqrt(ħ*t0/m); E0 = ħ/t0
#interactions and chemical potential, e.g., g  = (4*pi*ħ^2*as/m)*x0^3/E0 (dimensionless 3D)
  g = 0.1; μ  = 10.0
#damping parameters (dimensionless)
  γ = 0.0
  ℳ  = 0.0
#time evolution parameters
  ti = 0.0; #tf = 10.0/γ  #evolve for 2 damping times
  Nt = 50; t = collect(linspace(ti,tf,Nt)); dt = 0.01π/μ #integrate step size [ - should have dt ≪ 2π/μ]
#
@pack sinfo = ωx,ωy,ωz,γ,ℳ,g,t0,x0,E0,μ,ti,tf,Nt,t,dt
#== end parameters ==#

#= Initialize cfield (dimensionless) =#
  basis = "Hermite"
  ecut = Ncut*ħ*ωy/E0
  Ω = [ωx; ωy]*t0
  cinfo = makecinfo(ecut,Ω,basis)
  tinfo = maketinfo(cinfo,4)
  @unpack en,P,M = cinfo

#= create initial state and extra fields =#
  ψ0  = c2x(c0,tinfo) #initial condition
  ψ = similar(ψ0)  #a field to write to in place

#= PGPE (in place)       =#
#add spectral terms into Lgp!, or x-space terms into nlin!
  function Lgp!(dc,c,p,t)
      nlin!(dc,c,tinfo)
      dc .= P.*(-im*(1-im*γ)*((en - μ).*c .+ g*dc))
  end
  function nlin!(dc,c,tinfo)
      c2x!(ψ,c,tinfo)
      x2c!(dc,abs2.(ψ).*ψ,tinfo)
  end
#========================#

tspan = (t[1],t[end])
prob = ODEProblem(Lgp!,c0,tspan)
alg =    DP5() #abstol = 1e-3, reltol = 1e-3

println("Started evolution ...")
#@time sol = solve(prob,alg,dt=dt,saveat=t,alg_hints=[:stiff],save_everystep=false,dense=false);
@time sol = solve(prob,alg,dt=dt,saveat=t);
println("... Finished.")
return sinfo,cinfo,sol
end

#should probably be careful not to write over initial state here
sinfo2,cinfo2,sol2 = realtime(60.,30.5)

function findvortex()
## Transform to cartesian grid
@unpack ħ,m,ωx,ωy,ωz,γ,ℳ,g,x0,t0,E0,μ,ti,tf,Nt,t,dt = sinfo2
#@unpack M,Ω,ecut,P,en = cinfo

#lets be explicit about units:
Rx = sqrt(2μ*E0/m/ωx^2)/x0
Ry = sqrt(2μ*E0/m/ωy^2)/x0
yMax=1.5Ry
xMax=1.5Rx

Nx = 600
Ny = Nx
x = collect(linspace(-xMax,xMax,Nx))
y = collect(linspace(-yMax,yMax,Ny))
tinfop = maketinfoplot(cinfo,x,y)
ψ = c2x(sol2[1],tinfop)
mask = complex(Float64.(sqrt.(x.^2.+y'.^2).<0.7*Rx*ones(real(ψ))))

xt = zeros(t)
yt = zeros(t)
for j in eachindex(t)
    c2x!(ψ,sol2[j],tinfop)
    xt[j], yt[j], _ = findvortices(x,y,ψ.*mask)
end
return xt, yt
end

xt,yt = findvortex()

@unpack ħ,m,ωx,ωy,ωz,γ,ℳ,g,x0,t0,E0,μ,ti,tf,Nt,t,dt = sinfo2
plot(t,xt)
plot(t,yt)
#=
figure(figsize=(15,5))
subplot(1,2,1)
pcolormesh(x,y,g*abs2.(ψ'))
xlabel("x/x0")
ylabel("y/x0")
colorbar()
subplot(1,2,2)
pcolormesh(x,y,angle.(ψ'))
xlabel("x/x0")
ylabel("y/x0")
colorbar()

figure(figsize=(15,5))
subplot(1,2,1)
pcolormesh(x,y,g*abs2.(ψ'))
xlabel("x/x0")
ylabel("y/x0")
colorbar()
subplot(1,2,2)
pcolormesh(x,y,angle.(ψ'))
xlabel("x/x0")
ylabel("y/x0")
colorbar()

figure()
plot(x,g*abs2.(ψ[:,300]))
plot(x,μ*(1-x.^2/Rx^2).*vortexcore(x-xv,ξ).^2)
=#


text:
include("compute_ui.jl")
include("compute_z.jl")
include("compute_b.jl")
include("compute_yi.jl")
include("compute_dij.jl")

function get_num_coeffs(twojmax :: Int, n_elems :: Int, chem_flag :: Bool)
    J = twojmax
    if J % 2 == 0
        m = J/2 + 1
        num_coeffs = Int( m * (m+1) * (2*m+1) / 6 )
    elseif J % 2 == 1
        m = (J+1)/2
        num_coeffs = Int( m * (m+1) * (m+2) / 3 )
    else
        AssertionError("twojmax must be an integer!")
    end
    if chem_flag
        num_coeffs = num_coeffs * n_elems * n_elems * n_elems
    else
        num_coeffs = num_coeffs
    end
    return num_coeffs
end

function neighborlist(A::AbstractSystem, snap::SNAPParams)

    # Construct cutoff
    num_elements = length(snap.elements)
    cutsq = zeros(num_elements, num_elements)
    cutmax = 0.0
    for i = 1:num_elements
        cut = 2.0*snap.radii[i]*snap.rcutfac
        if cut > cutmax
            cutmax = cut
        end
        cutsq[i, i] = cut*cut

        for j = (i+1):num_elements
            cut = (snap.radii[i] + snap.radii[j]) * snap.rcutfac
            cutsq[i, j] = cut*cut
            cutsq[j, i] = cut*cut
        end
    end

    # Convert Positions to Matrix for Ball tree
    X = to_array(A)

    # Create Metric for Periodic Boundary Conditions
    periodic = periodicity(A)
    L = @SVector [ periodic[i] ? ustrip(bounding_box(A)[i][i]) : Inf for i = 1:3]
    d = Distances.PeriodicEuclidean(L)

    # Build Ball tree
    tree = BallTree(X, d)

    # elements = unique([element(ai) for ai in A])
    elements = unique(atomic_symbol(A))

    # Intialize empty vectors
    i = Vector{Int64}[] # i
    j = Vector{Int64}[] # j
    R = Vector{Float64}[] # Distances
    r = Vector{SVector{3, Float64}}[] # Positions

    # Fill vectors
    for n = 1:length(X)
        n_element = findall(x->x==Symbol(A[n].atomic_symbol), elements)[1]
        neighbors = inrange(tree, X[n], sqrt(cutmax), true)
        itemp = Int64[]
        jtemp = Int64[]
        Rtemp = Float64[]
        rtemp = SVector{3, Float64}[]
        for m in 1:length(neighbors)
            m_element = findall(x->x==Symbol(A[m].atomic_symbol), elements)[1]
            rr = get_distance(L, X[n], X[neighbors[m]])
            rsq = dot(rr,rr)
            if (neighbors[m] != n) & (rsq <= cutsq[n_element, m_element])
                push!(itemp, n)
                push!(jtemp, neighbors[m])
                push!(rtemp, rr)
                push!(Rtemp, sqrt(rsq))
            end
        end

        push!(i, itemp)
        push!(j, jtemp)
        push!(R, Rtemp)
        push!(r, rtemp)
    end
    return NeighborList(i, j, R, r)
end


text:
function rand(::AbstractRNG, ::SamplerTrivial{CloseOpen01{T}}) where {P, T<:ArbFloat{P}}
    oldprec = precision(BigFloat)
    setprecision(BigFloat, P+9)
    denom = BigInt(2)^P
    bigfloat = BigFloat(rand(BigInt(0):denom))
    bigfloat = bigfloat / BigFloat(denom)
    result = ArbFloat{P}(bigfloat)
    setprecision(BigFloat, oldprec)
    return result
end

@inline function rand(::AbstractRNG, ::SamplerTrivial{CloseOpen01{ArbFloat}})
    P = workingprecision(ArbFloat)
    return rand(ArbFloat{P})
end

function rand(::AbstractRNG, ::SamplerType{T}) where {P, T<:ArbReal{P}}
    arbfloat = rand(ArbFloat{P})
    return ArbReal{P}(arbfloat)
end

@inline function rand(::AbstractRNG, ::SamplerType{ArbReal})
    P = workingprecision(ArbReal)
    return rand(ArbReal{P})
end

function rand(::AbstractRNG, ::SamplerType{T}) where {P, T<:ArbComplex{P}}
    re = ArbReal{P}(rand(ArbFloat{P}))
    im = ArbReal{P}(rand(ArbFloat{P}))
    return ArbComplex{P}(re, im)
end

@inline function rand(::AbstractRNG, ::SamplerType{ArbComplex})
    P = workingprecision(ArbFloat)
    return rand(ArbComplex{P})
end

function randn(rng::AbstractRNG, ::Type{ArbFloat{N}}) where {N}
    urand1, urand2 = rand(ArbFloat{N}, 2)
    urand1 = urand1 + urand1 - 1
    urand2 = urand2 + urand2 - 1
    s = urand1*urand1 + urand2*urand2

    while s >= 1 || s === 0
        urand1, urand2 = rand(ArbFloat{N}, 2)
        urand1 = urand1 + urand1 - 1
        urand2 = urand2 + urand2 - 1
        s = urand1*urand1 + urand2*urand2
    end

    s = sqrt( -log(s) / s )
    return (urand1 + urand2) * s
end

randn(rng::AbstractRNG, ::Type{ArbFloat})  = randn(rng, ArbFloat{workingprecision(ArbFloat)})

function randn(rng::AbstractRNG, ::Type{ArbReal{N}}) where {N}
    urand1, urand2 = rand(ArbReal{N}, 2)
    urand1 = urand1 + urand1 - 1
    urand2 = urand2 + urand2 - 1
    s = urand1*urand1 + urand2*urand2

    while s >= 1 || s === 0
        urand1, urand2 = rand(ArbReal{N}, 2)
        urand1 = urand1 + urand1 - 1
        urand2 = urand2 + urand2 - 1
        s = urand1*urand1 + urand2*urand2
    end

    s = sqrt( -log(s) / s )
    return (urand1 + urand2) * s
end

randn(rng::AbstractRNG, ::Type{ArbReal})  = randn(rng, ArbReal{workingprecision(ArbReal)})


text:
"""
State machine used for these tests:

@startuml
state TestParentChildEventHandlingMachine {
    state TestParentChildEventHandlingChild
    TestParentChildEventHandlingChild : ChildHandledEvent / set handled_child_handled_event
}

TestParentChildEventHandlingMachine : ChildHandledEvent / set handled_child_handled_event
TestParentChildEventHandlingMachine : ChildUnhandledEvent / set_handled_child_unhandled_event
@enduml
"""

using Logging, Test, HSM

struct ChildHandledEvent <: HSM.AbstractHsmEvent
end

struct ChildUnhandledEvent <: HSM.AbstractHsmEvent
end

macro test_parent_child_event_state(name)
    return :(
        mutable struct $name <: HSM.AbstractHsmState
            state_info::HSM.HsmStateInfo
            handled_child_handled_event::Bool
            handled_child_unhandled_event::Bool

            $name(parent_state) = new(HSM.HsmStateInfo(parent_state), false, false)
        end
    )
end

@test_parent_child_event_state(TestParentChildEventHandlingMachine)
@test_parent_child_event_state(TestParentChildEventHandlingChild)

function HSM.on_event!(state::TestParentChildEventHandlingMachine, event::ChildHandledEvent)
    @debug "on_event!(TestParentChildEventHandlingMachine, ChildHandledEvent)"
    state.handled_child_handled_event = true
    return true
end

function HSM.on_event!(state::TestParentChildEventHandlingMachine, event::ChildUnhandledEvent)
    @debug "on_event!(TestParentChildEventHandlingMachine, ChildUnhandledEvent)"
    state.handled_child_unhandled_event = true
    return true
end

function HSM.on_event!(state::TestParentChildEventHandlingChild, event::ChildHandledEvent)
    @debug "on_event!(TestParentChildEventHandlingMachine, ChildHandledEvent)"
    state.handled_child_handled_event = true
    return true
end

@testset "handle_event!() -- Parent state handler not called when child handles event" begin
    machine = TestParentChildEventHandlingMachine(nothing)
    child = TestParentChildEventHandlingChild(machine)
    HSM.transition_to_state!(machine, child)

    HSM.handle_event!(child, ChildHandledEvent())
    @test machine.handled_child_handled_event == false
    @test machine.handled_child_unhandled_event == false
    @test child.handled_child_handled_event == true
    @test child.handled_child_unhandled_event == false
end

@testset "handle_event!() -- Parent state handler called when child does not handle event" begin
    machine = TestParentChildEventHandlingMachine(nothing)
    child = TestParentChildEventHandlingChild(machine)
    HSM.transition_to_state!(machine, child)

    HSM.handle_event!(child, ChildUnhandledEvent())
    @test machine.handled_child_handled_event == false
    @test machine.handled_child_unhandled_event == true
    @test child.handled_child_handled_event == false
    @test child.handled_child_unhandled_event == false
end


text:
module CircGeometry

include("struct_defs.jl")
include("is_inside_functions.jl")
include("shuffling_functions.jl")
include("main_functions.jl")
include("io.jl")

export MaterialParameters,
       OutlineCircle,
       OutlineRectangle,
       OutlinePolygon,
       generate_porous_structure,
       compute_between_buffer,
       compute_volume_fraction,
       Point,
       translate!

export write_circ,
       csv_to_polygon,
       save_image

end # module

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
# Low-level C API for libCEED

module C

using Libdl, libCEED_jll

const UINT_MAX = typemax(Cuint)
include(joinpath(@__DIR__, "generated", "libceed_bindings.jl"))

const CEED_STRIDES_BACKEND = Ref{Ptr{CeedInt}}()
const CEED_BASIS_COLLOCATED = Ref{CeedBasis}()
const CEED_VECTOR_ACTIVE = Ref{CeedVector}()
const CEED_VECTOR_NONE = Ref{CeedVector}()
const CEED_ELEMRESTRICTION_NONE = Ref{CeedElemRestriction}()
const CEED_QFUNCTION_NONE = Ref{CeedQFunction}()
const CEED_REQUEST_IMMEDIATE = Ref{CeedRequest}()
const CEED_REQUEST_ORDERED = Ref{CeedRequest}()

function __init__()
    global libceed_handle = dlopen(libceed)
    # some global variables
    CEED_STRIDES_BACKEND[] = cglobal((:CEED_STRIDES_BACKEND, libceed))
    CEED_BASIS_COLLOCATED[] =
        unsafe_load(cglobal((:CEED_BASIS_COLLOCATED, libceed), CeedBasis))
    CEED_VECTOR_ACTIVE[] = unsafe_load(cglobal((:CEED_VECTOR_ACTIVE, libceed), CeedVector))
    CEED_VECTOR_NONE[] = unsafe_load(cglobal((:CEED_VECTOR_NONE, libceed), CeedVector))
    CEED_ELEMRESTRICTION_NONE[] =
        unsafe_load(cglobal((:CEED_ELEMRESTRICTION_NONE, libceed), CeedElemRestriction))
    CEED_QFUNCTION_NONE[] =
        unsafe_load(cglobal((:CEED_QFUNCTION_NONE, libceed), CeedQFunction))
    CEED_REQUEST_IMMEDIATE[] =
        unsafe_load(cglobal((:CEED_REQUEST_IMMEDIATE, libceed), Ptr{CeedRequest}))
    CEED_REQUEST_ORDERED[] =
        unsafe_load(cglobal((:CEED_REQUEST_ORDERED, libceed), Ptr{CeedRequest}))
end

end # module

max_stars_repo_name: jedbrown/libCEED

max_issues_repo_name: jedbrown/libCEED

max_forks_repo_name: jedbrown/libCEED


text:
# --------------------------------------------------------------------------- #
# Inspired by model "FCRN-A" model from (Xie, Noble, Zisserman; 2015)

struct FCRNA{I <: Image, BN} <: FCModel{I, BN}
  weights :: Vector{Any}
  bnmoments :: Vector{Any}
end

# --------------------------------------------------------------------------- #
# Network architecture (no batch normalization)

function density{I <: Image}(w, s, x, ::Type{FCRNA{I, false}})
  y = rconv(w[2],  pool(rconv(w[1], x)))
  y = rconv(w[4],  pool(rconv(w[3], pool(y))))
  y = rconv(w[6],  uconv(w[5], y))
  y = rconv(w[8],  uconv(w[7], y))
  y = rconv(w[10], uconv(w[9], y))

  y = conv4(w[11], y, padding=pad(w[11]))
  return reshape(y, size(y)[[1,2,4]]...)
end

# --------------------------------------------------------------------------- #
# Initialization (no batch normalization)

function init_weights{I<:Image}(::Type{FCRNA{I, false}},
                                wr :: Function = wr,
                                wu :: Function = wu;
                                potency :: Integer = 6,
                                seed :: Integer = rand(1:10000))
  Knet.setseed(seed)
  c = imgchannels(I)
  s = 2^(potency-1)

  return Any[
    wr(c, s), wr(s, 2s), wr(2s, 4s), wr(4s, 16s),
    wu(16s, 16s), wr(16s, 4s), wu(4s, 4s), wr(4s, 2s),
    wu(2s, 2s), wr(2s, s), wr(s, 1) ]
end

init_state{I<:Image}(::Type{FCRNA{I, false}}) = Any[]

# --------------------------------------------------------------------------- #
# Network architecture (batch normalization)

function density{I <: Image}(w, s, x, ::Type{FCRNA{I, true}})
  y = rbconv(w[4],  s[2],  w[3],  pool(rbconv(w[2], s[1], w[1], x)))
  y = rbconv(w[8],  s[4],  w[7],  pool(rbconv(w[6], s[3], w[5], pool(y))))
  y = rbconv(w[12], s[6],  w[11], ubconv(w[10], s[5], w[9], y))
  y = rbconv(w[16], s[8],  w[15], ubconv(w[14], s[7], w[13], y))
  y = rbconv(w[20], s[10], w[19], ubconv(w[18], s[9], w[17], y))

  y = conv4(w[21], y, padding=pad(w[21]))
  return reshape(y, size(y)[[1,2,4]]...)
end

# --------------------------------------------------------------------------- #
# Initialization (batch normalization)

function init_weights{I<:Image}(::Type{FCRNA{I, true}},
                                wr :: Function = wr,
                                wu :: Function = wu;
                                potency :: Integer = 6,
                                seed :: Integer = rand(1:10000))

  # Need the same convolutional weights as in the case without batch
  # normalization. However, also need to insert the trainable bn-values at
  # the right positions of the weights vector
  s = 2^(potency-1)

  bm = [s, 2s, 4s, 16s, 16s, 4s, 4s, 2s, 2s, s]
  weights = init_weights(FCRNA{I, false}, wr, wu,
                         potency = potency, seed = seed)
  return Any[ odd(i) ? weights[ceil(Int, i/2)] : bn(bm[div(i,2)]) for i in 1:21 ]
end

function init_state{I<:Image}(::Type{FCRNA{I, true}})
  return Any[ bnmoments() for i in 1:10 ]
end

# --------------------------------------------------------------------------- #
# Convenience-constructor

function FCRNA(I,
               wr :: Function = wr,
               wu :: Function = wu;
               bn :: Bool = false,
               potency :: Integer = 6,
               seed :: Integer = rand(1:10000))

  weights = init_weights(FCRNA{I, bn}, wr, wu,
                         potency = potency, seed = seed)
  state = init_state(FCRNA{I, bn})

  return FCRNA{I, bn}(weights, state)
end

max_stars_repo_name: tscode/DCellC

max_issues_repo_name: tscode/DCellC


text:
module UnwrappingTests

using Unwrapping, Test, Statistics, GZip
using Unwrapping: unwrap, fft_unwrap, wrap, arc

const TESTDIR = @__DIR__

function read_test_data(filename::AbstractString)
    name = basename(filename)
    gzipped = endswith(name, ".gz")
    last = length(name)
    if gzipped
        last -= 3
    end
    i3 = findprev(isequal('.'), name, last)
    (i3 === nothing || SubString(name, i3:i3+3) != ".dat") &&
        error("invalid file extension")
    i2 = findprev(isequal('x'), name, i3-1)
    i2 === nothing && error("invalid file name (missing 'x' in dimensions)")
    i1 = findprev(isequal('-'), name, i2-1)
    i1 === nothing && error("invalid file name (missing dimensions)")
    n1 = tryparse(Int, SubString(name, i1+1:i2-1); base=10)
    (n1 === nothing || n1 < 1) &&
        error("invalid file name (bad first dimension)")
    n2 = tryparse(Int, SubString(name, i2+1:i3-1); base=10)
    (n2 === nothing || n2 < 1) &&
        error("invalid file name (bad second dimension)")
    A = Array{UInt8}(undef, n1, n2)
    if gzipped
        GZip.open(filename, "r") do io
            read!(io, A)
        end
    else
        open(filename, "r") do io
            read!(io, A)
        end
    end
    return A
end

function make_screen(T::Type{<:AbstractFloat}, dims::Dims{2})
    A = Array{T}(undef, dims)
    n1, n2 = dims
    p = Unwrapping.twopi(T) # period
    c1 = T(n1+1)/2
    c2 = T(n2+1)/2
    # keep |ai|*ni + |bi| << p to avoid excessive wrapping
    a1 = p/(7*n1)
    a2 = p/(9*n2)
    b1 = T(-0.4)
    b2 = T(+0.3)
    for i2 in 1:n2
        x2 = i2 - c2
        for i1 in 1:n1
            x1 = i1 - c1
            A[i1,i2] = (a1*x1 + b1)*x1 + (a2*x2 + b2)*x2
        end
    end
    return A
end

gauss = read_test_data(joinpath(TESTDIR, "gaussian-434x342.dat.gz"))
peaks = read_test_data(joinpath(TESTDIR, "peaks-434x342.dat.gz"))

#dims = (384, 288)

@testset "DCT and FFT-based methods (T = $T)" for T in (Float32, Float64)
    dct_tol = 200*eps(T)
    fft_tol = 0.005
    #A = make_screen(T, dims)
    for A in (T(0.2)*gauss, T(0.2)*peaks)
        Aw = arc.(A)
        A1 = unwrap(Aw; debiased=false)
        A2 = unwrap(Aw; debiased=true)
        B1 = fft_unwrap(Aw; debiased=false)
        B2 = fft_unwrap(Aw; debiased=true)

        # std() method is insensitive to a bias
        @test std(A - A1) ≤ dct_tol*std(A - Aw)
        @test std(A - A2) ≤ dct_tol*std(A - Aw)
        @test std(A - B1) ≤ fft_tol*std(A - Aw)
        @test std(A - B2) ≤ fft_tol*std(A - Aw)
        # A2 and B2 are debiased (modulo 2π)
        @test abs(arc(mean(A - A2))) ≤ 2e-3
        @test abs(arc(mean(A - B2))) ≤ 2e-3
    end
end

end # module


text:
import Base.exp

abstract type Function end

mutable struct Variable
    data::Array
    grad::Union{Array,Nothing}
    creator::Union{Function,Nothing}
    Variable(data) = new(data, nothing, nothing)
end
function backward(v::Variable)::Nothing
    if isnothing(v.grad)
        v.grad = ones(size(v.data))
    end
    funcs::Array{Function} = [v.creator]
    while !isempty(funcs)
        f::Function = pop!(funcs) # No need to check for nothingness at the first time
        x::Variable, y::Variable = f.input, f.output
        x.grad = backward(f, y.grad)
        if !isnothing(x.creator)
            push!(funcs, x.creator)
        end
    end
end

function _function(f::Function, forward::Core.Function, input::Variable)::Variable
    x::Array = input.data
    y::Array = forward(x)
    output = Variable(y)
    output.creator = f
    f.input = input
    f.output = output
    return output
end

mutable struct Square <: Function
    input::Union{Variable,Nothing}
    output::Union{Variable,Nothing}
    Square() = new()
end
(s::Square)(input::Variable)::Variable = _function(s, x -> x .^ 2, input)
backward(s::Square, gy::Array)::Array = 2 * s.input.data .* gy
square(x::Variable)::Variable = Square()(x)

mutable struct Exp <: Function
    input::Union{Variable,Nothing}
    output::Union{Variable,Nothing}
    Exp() = new()
end
(s::Exp)(input::Variable)::Variable = _function(s, x -> exp.(x), input)
backward(e::Exp, gy::Array)::Array = exp.(e.input.data) .* gy
exp(x::Variable)::Variable = Exp()(x)


text:
"""
A Julia suite for chaos and nonlinear dynamics
"""
module DynamicalSystems

using StaticArrays

export SVector, SMatrix, @SVector, @SMatrix, Size

using Reexport

@reexport using DelayEmbeddings
@reexport using DynamicalSystemsBase
@reexport using Entropies
@reexport using ChaosTools
@reexport using RecurrenceAnalysis

# Update messages:
using Scratch
display_update = true
version_number = "2.2.0"
update_name = "update_v$(version_number)"

if display_update
    # Get scratch space for this package
    versions_dir = @get_scratch!("versions")
    if !isfile(joinpath(versions_dir, update_name))
        printstyled(
            stdout,
            """
            \nUpdate message: DynamicalSystems v$(version_number)\n
            Direct keyword acceptance of DifferentialEquations.jl related keywords is
            no longer supported in DynamicalSystems.jl. This means that this is deprecated:
            ```
            lyapunov(ds, T; alg = Tsit5(), abstol = 1e-9)
            ```
            Instead, use the explicit keyword `diffeq`, and give it a named tuple with
            the keyword arguments, i.e.:
            ```
            lyapunov(ds, T; diffeq = (alg = Tsit5(), abstol = 1e-9))
            ```

            Also, the online documentation now uses Makie.jl,
            references InteractiveDynamics.jl and showcases some of its functions.
            """;
            color = :light_magenta,
        )
        touch(joinpath(versions_dir, update_name))
    end
end

end # module


text:
module FastArrayOps
import Base.LinAlg: BlasReal, BlasComplex, BlasFloat, BlasInt, BlasChar
const libblas = Base.libblas_name
export fast_scale!,     unsafe_fast_scale!,
       fast_add!,       unsafe_fast_add!,
       fast_addscal!,   unsafe_fast_addscal!,
       fast_copy!,      unsafe_fast_copy!,
       fast_fill!,      unsafe_fast_fill!
export fast_check1, fast_check2, fast_check3, nmax2nel, nel2nmax, fast_args2range, fast_range2args

# WARNING: FastArrayOps.jl gets overwritten by FastArrayOps_src.jl when running make.jl

## CONSTANTS

const NLIM_SCALE = 13
const NLIM_SCALE_OOP1 = 80
const NLIM_SCALE_OOP2 = 100000
const NLIM_SCALEARR = typemax(Int)
const NLIM_SCALEARR_OOP = typemax(Int)
const NLIM_ADD = typemax(Int)
const NLIM_ADD_OOP = typemax(Int)
const NLIM_ADDARR = 13
const NLIM_ADDARR_OOP1 = 30
const NLIM_ADDARR_OOP2 = 100000
const NLIM_ADDARRSCAL = 13
const NLIM_ADDARRSCAL_OOP1 = 30
const NLIM_ADDARRSCAL_OOP2 = 100000
const NLIM_COPY1 = 80
const NLIM_COPY2 = 100000
const NLIM_FILL = 13

## UTILS

function nmax2nel(i::Int, inc::Int, nmax::Int)
    @assert 0 < i
    nmax < i && return 0
    return div(nmax - i, abs(inc)) + 1
end
function nel2nmax(i::Int, inc::Int, nel::Int)
    @assert 0 < i
    nel < 0 && return i - 1
    return i + (nel- 1)*abs(inc)
end
function fast_args2range(i::Int, inc::Int, n::Int)
    @assert 0 < i
    r = i:abs(inc):nel2nmax(i, abs(inc), n)
    if inc < 0
        r = reverse(r)
    end
    return r
end
function fast_range2args(r::Range)
    inc = step(r)
    if inc > 0
        return (first(r), inc, length(r))
    else
        return (last(r), inc, length(r))
    end
end

function fast_check1(x, ix, incx, n)
    ##!fast_check1
    return 0
end
function fast_check2(x, ix, incx, y, iy, incy, n)
    ##!fast_check2
    return 0
end
function fast_check3(x, ix, incx, y, iy, incy, z, iz, incz, n)
    ##!fast_check3
    return 0
end

# utils for 0 fill
inttype(::Type{Float32}) = Int32
inttype(::Type{Float64}) = Int64
inttype(::Type{Complex64}) = Int64
inttype(::Type{Complex128}) = Int128
function fast_reinterpret1{T<:Number}(::Type{T}, a::Array)
    @assert length(a) == 1
    ccall(:jl_reshape_array, Array{T,1}, (Any, Any, Any), Array{T,1}, a, (1,))
end
# are float 0s zero bits?
function zerobits()
    v = reinterpret(inttype(Float32), convert(Float32, 0))
    v += reinterpret(inttype(Float64), convert(Float64, 0))
    v += fast_reinterpret1(inttype(Complex64), [convert(Complex64, 0)])[1]
    v += fast_reinterpret1(inttype(Complex128), [convert(Complex128, 0)])[1]
    if v == 0
        return true
    else
        return false
    end
end
const ZEROFLOAT = zerobits()

for (fscal, fcopy, faxpy, ftbmv, fsbmv, elty) in (
                        (:dscal_, :dcopy_, :daxpy_, :dtbmv_, :dsbmv_, :Float64),
                        (:sscal_, :scopy_, :saxpy_, :stbmv_, :ssbmv_, :Float32),
                        (:zscal_, :zcopy_, :zaxpy_, :ztbmv_, :zsbmv_, :Complex128),
                        (:cscal_, :ccopy_, :caxpy_, :ctbmv_, :csbmv_, :Complex64))

## SCALE METHODS

for (f, isunsafe) in ( (:fast_scale!, false), (:unsafe_fast_scale!, true) )
@eval begin

# x = a*x
# =======
# general
function ($f)(x::Array{$elty}, ix::Int, incx::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!fast_check1
    end
    if n < $NLIM_SCALE*incx
        ##!scalarr1_for#*
    else
        ##!blas_scale
    end
    return x
end
# inc1
function ($f)(x::Array{$elty}, ix::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!set1_inc1
        ##!fast_check1
    end
    if n < $NLIM_SCALE #*incx
        ##!scalarr1_for_inc1#*
    else
        ##!set1_inc1
        ##!blas_scale
    end
    return x
end

# x = a*y
# =======
# general
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, incy::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!fast_check2
    end
    mul = max(abs(incx), abs(incy))
    if n < $NLIM_SCALE_OOP1*mul || n*mul > $NLIM_SCALE_OOP2
        ##!scalarr1_foroop#*
    else
        ##!blas_copy
        ##!blas_scale
    end
    return x
end
# inceq
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!set2_inceq
        ##!fast_check2
    end
    mul = abs(incx)
    if n < $NLIM_SCALE_OOP1*mul || n*mul > $NLIM_SCALE_OOP2
        ##!scalarr1_foroop_inceq#*
    else
        ##!set2_inceq
        ##!blas_copy
        ##!blas_scale
    end
    return x
end
# inc1
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, iy::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!set2_inc1
        ##!fast_check2
    end
    if n < $NLIM_SCALE_OOP1 || n > $NLIM_SCALE_OOP2
        ##!scalarr1_foroop_inc1#*
    else
        ##!set2_inc1
        ##!c_memcpy
        ##!blas_scale
    end
    return x
end
# inc1ieq
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, a::$elty, n::Int)
    $isunsafe || begin
        ##!set2_inc1ieq
        ##!fast_check2
    end
    if n < $NLIM_SCALE_OOP1 || n > $NLIM_SCALE_OOP2
        ##!scalarr1_foroop_inc1ieq#*
    else
        ##!set2_inc1ieq
        ##!c_memcpy
        ##!blas_scale
    end
    return x
end

# x = x.*y
# ========
# general
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, incy::Int, n::Int)
    $isunsafe || begin
        ##!fast_check2
    end
    ##!arr2xy_for#*
    return x
end
# inceq
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, n::Int)
    $isunsafe || begin
        ##!set2_inceq
        ##!fast_check2
    end
    ##!arr2xy_for_inceq#*
    return x
end
# inc1
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, iy::Int, n::Int)
    $isunsafe || begin
        ##!set2_inc1
        ##!fast_check2
    end
    ##!arr2xy_for_inc1#*
    return x
end
# inc1ieq
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, n::Int)
    $isunsafe || begin
        ##!set2_inc1ieq
        ##!fast_check2
    end
    ##!arr2xy_for_inc1ieq#*
    return x
end

# x = y.*z
# ========

# general
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, incy::Int, z::Array{$elty}, iz::Int, incz::Int, n::Int)
    $isunsafe || begin
        ##!fast_check3
    end
    ##!arr2yz_foroop#*
    return x
end
# inceq
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, z::Array{$elty}, iz::Int, n::Int)
    $isunsafe || begin
        ##!set3_inceq
        ##!fast_check3
    end
    ##!arr2yz_foroop_inceq#*
    return x
end
# inc1
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, iy::Int, z::Array{$elty}, iz::Int, n::Int)
    $isunsafe || begin
        ##!set3_inc1
        ##!fast_check3
    end
    ##!arr2yz_foroop_inc1#*
    return x
end
# inc1ieq
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, z::Array{$elty}, n::Int)
    $isunsafe || begin
        ##!set3_inc1ieq
        ##!fast_check3
    end
    ##!arr2yz_foroop_inc1ieq#*
    return x
end

end # eval begin
end # for

## ADD METHODS

for (f, isunsafe) in ( (:fast_add!, false), (:unsafe_fast_add!, true) )
@eval begin

# x = x + a
# =======
# general
function ($f)(x::Array{$elty}, ix::Int, incx::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!fast_check1
    end
    ##!scalarr1_for#+
    return x
end

# inc1
function ($f)(x::Array{$elty}, ix::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!set1_inc1
        ##!fast_check1
    end
    ##!scalarr1_for_inc1#+
    return x
end

# x = y + a
# =======
# general
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, incy::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!fast_check2
    end
    ##!scalarr1_foroop#+
    return x
end
# inceq
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!set2_inceq
        ##!fast_check2
    end
    ##!scalarr1_foroop_inceq#+
    return x
end
# inc1
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, iy::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!set2_inc1
        ##!fast_check2
    end
    ##!scalarr1_foroop_inc1#+
    return x
end
# inc1ieq
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, a::$elty, n::Int)
    $isunsafe || begin
        ##!set2_inc1ieq
        ##!fast_check2
    end
    ##!scalarr1_foroop_inc1ieq#+
    return x
end

# x = x + y
# =========
# general
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, incy::Int, n::Int)
    $isunsafe || begin
        ##!fast_check2
    end
    if n < $NLIM_ADDARR #*mul # || n*mul > $NLIM_SCALEARR
        ##!arr2xy_for#+
    else
        a = 1
        ##!blas_axpy
    end
    return x
end
# inceq
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, n::Int)
    $isunsafe || begin
        ##!set2_inceq
        ##!fast_check2
    end
    if n < $NLIM_ADDARR
        ##!arr2xy_for_inceq#+
    else
        a = 1
        ##!set2_inceq
        ##!blas_axpy
    end
    return x
end
# inc1
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, iy::Int, n::Int)
    $isunsafe || begin
        ##!set2_inc1
        ##!fast_check2
    end
    if n < $NLIM_ADDARR
        ##!arr2xy_for_inc1#+
    else
        a = 1
        ##!set2_inc1
        ##!blas_axpy
    end
    return x
end
# inc1ieq
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, n::Int)
    $isunsafe || begin
        ##!set2_inc1ieq
        ##!fast_check2
    end
    if n < $NLIM_ADDARR
        ##!arr2xy_for_inc1ieq#+
    else
        a = 1
        ##!set2_inc1ieq
        ##!blas_axpy
    end
    return x
end

# x = y + z
# =========
# general
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, incy::Int, z::Array{$elty}, iz::Int, incz::Int, n::Int)
    $isunsafe || begin
        ##!fast_check3
    end
    if n < $NLIM_ADDARR_OOP1 || n > $NLIM_ADDARR_OOP2 #*mul # || n*mul > $NLIM_SCALEARR
        ##!arr2yz_foroop#+
    else
        a = 1
        ##!blas_copy
        y, iy, incy = z, iz, incz
        ##!blas_axpy
    end
    return x
end
# inceq
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, z::Array{$elty}, iz::Int, n::Int)
    $isunsafe || begin
        ##!set3_inceq
        ##!fast_check3
    end
    if n < $NLIM_ADDARR_OOP1 || n > $NLIM_ADDARR_OOP2
        ##!arr2yz_foroop_inceq#+
    else
        a = 1
        ##!set3_inceq
        ##!blas_copy
        y, iy, incy = z, iz, incz
        ##!blas_axpy
    end
    return x
end
# inc1
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, iy::Int, z::Array{$elty}, iz::Int, n::Int)
    $isunsafe || begin
        ##!set3_inc1
        ##!fast_check3
    end
    if n < $NLIM_ADDARR_OOP1 || n > $NLIM_ADDARR_OOP2
        ##!arr2yz_foroop_inc1#+
    else
        a = 1
        ##!set3_inc1
        ##!c_memcpy
        y, iy, incy = z, iz, incz
        ##!blas_axpy
    end
    return x
end
# inc1ieq
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, z::Array{$elty}, n::Int)
    $isunsafe || begin
        ##!set3_inc1ieq
        ##!fast_check3
    end
    if n < $NLIM_ADDARR_OOP1 || n > $NLIM_ADDARR_OOP2
        ##!arr2yz_foroop_inc1ieq#+
    else
        a = 1
        ##!set3_inc1ieq
        ##!c_memcpy
        y, iy, incy = z, iz, incz
        ##!blas_axpy
    end
    return x
end

end # eval begin
end # for

## ADDSCAL METHODS

for (f, isunsafe) in ( (:fast_addscal!, false), (:unsafe_fast_addscal!, true) )
@eval begin

# x = x + a*y
# =========
# general
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, incy::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!fast_check2
    end
    if n < $NLIM_ADDARRSCAL #*mul # || n*mul > $NLIM_SCALEARR
        ##!addarrscal_for
    else
        ##!blas_axpy
    end
    return x
end
# inceq
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!set2_inceq
        ##!fast_check2
    end
    if n < $NLIM_ADDARRSCAL
        ##!addarrscal_for_inceq
    else
        ##!set2_inceq
        ##!blas_axpy
    end
    return x
end
# inc1
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, iy::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!set2_inc1
        ##!fast_check2
    end
    if n < $NLIM_ADDARRSCAL
        ##!addarrscal_for_inc1
    else
        ##!set2_inc1
        ##!blas_axpy
    end
    return x
end
# inc1ieq
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, a::$elty, n::Int)
    $isunsafe || begin
        ##!set2_inc1ieq
        ##!fast_check2
    end
    if n < $NLIM_ADDARRSCAL
        ##!addarrscal_for_inc1ieq
    else
        ##!set2_inc1ieq
        ##!blas_axpy
    end
    return x
end

# x = y + a*z
# =========
# general
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, incy::Int, z::Array{$elty}, iz::Int, incz::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!fast_check3
    end
    if n < $NLIM_ADDARRSCAL_OOP1 || n > $NLIM_ADDARRSCAL_OOP2 #*mul # || n*mul > $NLIM_SCALEARR
        ##!addarrscal_foroop
    else
        ##!blas_copy
        y, iy, incy = z, iz, incz
        ##!blas_axpy
    end
    return x
end
# inceq
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, z::Array{$elty}, iz::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!set3_inceq
        ##!fast_check3
    end
    if n < $NLIM_ADDARRSCAL_OOP1 || n > $NLIM_ADDARRSCAL_OOP2
        ##!addarrscal_foroop_inceq
    else
        ##!set3_inceq
        ##!blas_copy
        y, iy, incy = z, iz, incz
        ##!blas_axpy
    end
    return x
end
# inc1
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, iy::Int, z::Array{$elty}, iz::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!set3_inc1
        ##!fast_check3
    end
    if n < $NLIM_ADDARRSCAL_OOP1 || n > $NLIM_ADDARRSCAL_OOP2
        ##!addarrscal_foroop_inc1
    else
        ##!set3_inc1
        ##!c_memcpy
        y, iy, incy = z, iz, incz
        ##!blas_axpy
    end
    return x
end
# inc1ieq
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, z::Array{$elty}, a::$elty, n::Int)
    $isunsafe || begin
        ##!set3_inc1ieq
        ##!fast_check3
    end
    if n < $NLIM_ADDARRSCAL_OOP1 || n > $NLIM_ADDARRSCAL_OOP2
        ##!addarrscal_foroop_inc1ieq
    else
        ##!set3_inc1ieq
        ##!c_memcpy
        y, iy, incy = z, iz, incz
        ##!blas_axpy
    end
    return x
end

end # eval begin
end # for

## COPY METHODS

for (f, isunsafe) in ( (:fast_copy!, false), (:unsafe_fast_copy!, true) )
@eval begin

# x = y
# =====
# general
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, incy::Int, n::Int)
    $isunsafe || begin
        ##!fast_check2
    end
    mul = max(abs(incx), abs(incy))
    if n < $NLIM_COPY1*mul || n*mul > $NLIM_COPY2
        ##!copy_foroop
    else
        ##!blas_copy
    end
    return x
end
# inceq
function ($f)(x::Array{$elty}, ix::Int, incx::Int, y::Array{$elty}, iy::Int, n::Int)
    $isunsafe || begin
        ##!set2_inceq
        ##!fast_check2
    end
    mul = abs(incx)
    if n < $NLIM_COPY1*mul || n*mul > $NLIM_COPY2
        ##!copy_foroop_inceq
    else
        ##!set2_inceq
        ##!blas_copy
    end
    return x
end
# inc1
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, iy::Int, n::Int)
    $isunsafe || begin
        ##!set2_inc1
        ##!fast_check2
    end
    ##!set2_inc1
    ##!c_memcpy
    return x
end
# inc1ieq
function ($f)(x::Array{$elty}, ix::Int, y::Array{$elty}, n::Int)
    $isunsafe || begin
        ##!set2_inc1ieq
        ##!fast_check2
    end
    ##!set2_inc1ieq
    ##!c_memcpy
    return x
end

end # eval begin
end # for

## FILL METHODS

for (f, isunsafe) in ( (:fast_fill!, false), (:unsafe_fast_fill!, true) )
@eval begin

# x = a
# =======
# general
function ($f)(x::Array{$elty}, ix::Int, incx::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!fast_check1
    end
    ##!fill_foroop
    return x
end
# inc1
function ($f)(x::Array{$elty}, ix::Int, a::$elty, n::Int)
    $isunsafe || begin
        ##!set1_inc1
        ##!fast_check1
    end
    if a == 0 && n > $NLIM_FILL && ZEROFLOAT
        ##!set1_inc1
        ##!c_memset
    else
        ##!fill_foroop_inc1
    end
    return x
end

end # eval begin
end # for

end # for

end # module


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

@with_kw struct ProposalCovTunerConfig <: AbstractMCMCTunerConfig
    λ::Float64 = 0.5
    α::IntervalSets.ClosedInterval{Float64} = ClosedInterval(0.15, 0.35)
    β::Float64 = 1.5
    c::IntervalSets.ClosedInterval{Float64} = ClosedInterval(1e-4, 1e2)
end

export ProposalCovTunerConfig

AbstractMCMCTunerConfig(algorithm::MetropolisHastings) = ProposalCovTunerConfig()

(config::ProposalCovTunerConfig)(chain::MCMCIterator{<:MetropolisHastings}; init_proposal::Bool = true) =
    ProposalCovTuner(config, chain, init_proposal)

mutable struct ProposalCovTuner{
    C<:MCMCIterator{<:MetropolisHastings},
    S<:MCMCBasicStats
} <: AbstractMCMCTuner
    config::ProposalCovTunerConfig
    chain::C
    stats::S
    iteration::Int
    scale::Float64
end

export ProposalCovTuner

function ProposalCovTuner(
    config::ProposalCovTunerConfig,
    chain::MCMCIterator{<:MetropolisHastings},
    init_proposal::Bool = true
)
    m = nparams(chain)
    scale = 2.38^2 / m
    tuner = ProposalCovTuner(config, chain, MCMCBasicStats(chain), 1, scale)

    if init_proposal
        tuning_init_proposal!(tuner)
    end

    tuner
end

isviable(tuner::ProposalCovTuner) = nsamples(tuner.chain.state) >= 2

function tuning_init_proposal!(tuner::ProposalCovTuner)
    chain = tuner.chain

    # ToDo: Generalize, currently limited to HyperRectBounds
    bounds = param_bounds(chain.target)
    vol = spatialvolume(bounds)
    flat_var = (vol.hi - vol.lo).^2 / 12

    m = length(flat_var)
    Σ_unscaled = full(PDiagMat(flat_var))
    Σ = Σ_unscaled * tuner.scale

    next_cycle!(chain)
    chain.state.pdist = set_cov(chain.state.pdist, Σ)

    chain
end

function tuning_update!(tuner::ProposalCovTuner; ll::LogLevel = LOG_NONE)
    chain = tuner.chain
    stats = tuner.stats
    config = tuner.config

    α_min = minimum(config.α)
    α_max = maximum(config.α)

    c_min = minimum(config.c)
    c_max = maximum(config.c)

    β = config.β

    state = chain.state

    t = tuner.iteration
    λ = config.λ
    c = tuner.scale
    Σ_old = full(get_cov(state.pdist))

    S = convert(Array, stats.param_stats.cov)
    a_t = 1 / t^λ
    new_Σ_unscal = (1 - a_t) * (Σ_old/c) + a_t * S

    α = eff_acceptance_ratio(state)

    if α_min <= α <= α_max
        chain.tuned = true
        @log_msg ll "MCMC chain $(chain.id) tuned, acceptance ratio = $α"
    else
        chain.tuned = false
        @log_msg ll "MCMC chain $(chain.id) *not* tuned, acceptance ratio = $α"

        if α > α_max && c < c_max
            tuner.scale = c * β
        elseif α < α_min && c > c_min
            tuner.scale = c / β
        end
    end

    Σ_new = full(Hermitian(new_Σ_unscal * tuner.scale))

    if !isposdef(Σ_new)
        ev = eig(Σ_new)[1]
        ev_min = minimum(ev)
        Σ_11 = Σ_new[1,1]
        Σ_new += max(10 * eps(typeof(ev_min)), -ev_min) * eye(length(ev))
        Σ_new ./= (Σ_new[1,1] / Σ_11)
    end

    next_cycle!(chain)
    state.pdist = set_cov(state.pdist, Σ_new)
    tuner.iteration += 1

    chain
end

function run_tuning_cycle!(
    callbacks,
    tuners::AbstractVector{<:ProposalCovTuner},
    exec_context::ExecContext = ExecContext();
    ll::LogLevel = LOG_NONE,
    kwargs...
)
    run_tuning_iterations!(callbacks, tuners, exec_context; ll=ll, kwargs...)
    tuning_update!.(tuners; ll = ll)
    nothing
end

function run_tuning_iterations!(
    callbacks,
    tuners::AbstractVector{<:ProposalCovTuner},
    exec_context::ExecContext;
    max_nsamples::Int64 = Int64(1000),
    max_nsteps::Int64 = Int64(10000),
    max_time::Float64 = Inf,
    ll::LogLevel = LOG_NONE
)
    chains = map(x -> x.chain, tuners)
    user_callbacks = mcmc_callback_vector(callbacks, chains)

    combined_callbacks = broadcast(tuners, user_callbacks) do tuner, user_callback
        (level, chain) -> begin
            if level == 1
                get_samples!(tuner.stats, chain, true)
            end
            user_callback(level, chain)
        end
    end

    mcmc_iterate!(combined_callbacks, chains, exec_context, max_nsamples = max_nsamples, max_nsteps = max_nsteps, max_time = max_time)
    nothing
end


text:
using Test
include("lib/lib.jl")
include("main/main.jl")
include("dev.jl")


text:
function parse_kw(ps::ParseState, ::Type{Val{Tokens.FUNCTION}})
    startbyte = ps.t.startbyte
    start_col = ps.t.startpos[2] + 4

    # Parsing
    kw = INSTANCE(ps)
    format_kw(ps)
    # signature

    if isoperator(ps.nt.kind) && ps.nt.kind != Tokens.EX_OR && ps.nnt.kind == Tokens.LPAREN
        start1 = ps.nt.startbyte
        next(ps)
        op = OPERATOR(ps)
        next(ps)
        if issyntaxunarycall(op)
            sig = EXPR{UnarySyntaxOpCall}(EXPR[op, INSTANCE(ps)], 0, Variable[], "")
        else
            sig = EXPR{Call}(EXPR[op, INSTANCE(ps)], 0, Variable[], "")
        end
        @catcherror ps startbyte @default ps @closer ps paren parse_comma_sep(ps, sig)
        next(ps)
        push!(sig.args, INSTANCE(ps))
        sig.span = ps.nt.startbyte - start1
        @default ps @closer ps inwhere @closer ps ws @closer ps block while !closer(ps)
            @catcherror ps startbyte sig = parse_compound(ps, sig)
        end
    else
        @catcherror ps startbyte sig = @default ps @closer ps inwhere @closer ps block @closer ps ws parse_expression(ps)
    end

    while ps.nt.kind == Tokens.WHERE
        @catcherror ps startbyte sig = @default ps @closer ps inwhere @closer ps block @closer ps ws parse_compound(ps, sig)
    end

    if sig isa EXPR{InvisBrackets} && !(sig.args[2] isa EXPR{TupleH})
        sig = EXPR{TupleH}(sig.args, sig.span, Variable[], "")
    end

    _get_sig_defs!(sig)

    block = EXPR{Block}(EXPR[], 0, Variable[], "")
    @catcherror ps startbyte @default ps @scope ps Scope{Tokens.FUNCTION} parse_block(ps, block, start_col)

    # Construction
    if isempty(block.args)
        if sig isa EXPR{Call} || sig isa EXPR{BinarySyntaxOpCall} && !(sig.args[1] isa EXPR{OPERATOR{PlusOp,Tokens.EX_OR,false}})
            args = EXPR[sig, block]
        else
            args = EXPR[sig]
        end
    else
        args = EXPR[sig, block]
    end

    next(ps)

    ret = EXPR{FunctionDef}(EXPR[kw], ps.nt.startbyte - startbyte, Variable[], "")
    for a in args
        push!(ret.args, a)
    end
    push!(ret.args, INSTANCE(ps))

    ret.defs = [Variable(Expr(_get_fname(sig)), :Function, ret)]
    return ret
end

"""
    parse_call(ps, ret)

Parses a function call. Expects to start before the opening parentheses and is passed the expression declaring the function name, `ret`.
"""
function parse_call(ps::ParseState, ret::EXPR{OPERATOR{PlusOp,Tokens.EX_OR,false}})
    startbyte = ps.t.startbyte
    arg = @precedence ps 20 parse_expression(ps)
    ret = EXPR{UnarySyntaxOpCall}(EXPR[ret, arg], ret.span + arg.span, Variable[], "")
    return ret
end
function parse_call(ps::ParseState, ret::EXPR{OPERATOR{DeclarationOp,Tokens.DECLARATION,false}})
    startbyte = ps.t.startbyte
    arg = @precedence ps 20 parse_expression(ps)
    ret = EXPR{UnarySyntaxOpCall}(EXPR[ret, arg], ret.span + arg.span, Variable[], "")
    return ret
end
function parse_call(ps::ParseState, ret::EXPR{OP}) where OP <: OPERATOR{TimesOp,Tokens.AND}
    startbyte = ps.t.startbyte
    arg = @precedence ps 20 parse_expression(ps)
    ret = EXPR{UnarySyntaxOpCall}(EXPR[ret, arg], ret.span + arg.span, Variable[], "")
    return ret
end
function parse_call(ps::ParseState, ret::EXPR{OPERATOR{ComparisonOp,Tokens.ISSUBTYPE,false}})
    startbyte = ps.t.startbyte
    arg = @precedence ps 13 parse_expression(ps)
    ret = EXPR{Call}(EXPR[ret; arg.args], ret.span + arg.span, Variable[], "")
    return ret
end

function parse_call(ps::ParseState, ret::EXPR{OPERATOR{ComparisonOp,Tokens.ISSUPERTYPE,false}})
    startbyte = ps.t.startbyte
    arg = @precedence ps 13 parse_expression(ps)
    ret = EXPR{Call}(EXPR[ret; arg.args], ret.span + arg.span, Variable[], "")
    return ret
end

function parse_call(ps::ParseState, ret::EXPR{OP}) where OP <: OPERATOR{20,Tokens.NOT}
    startbyte = ps.t.startbyte
    arg = @precedence ps 13 parse_expression(ps)
    if arg isa EXPR{TupleH}
        ret = EXPR{Call}(EXPR[ret; arg.args], ret.span + arg.span, Variable[], "")
    else
        ret = EXPR{UnaryOpCall}(EXPR[ret, arg], ret.span + arg.span, Variable[], "")
    end
    return ret
end

function parse_call(ps::ParseState, ret::EXPR{OP}) where OP <: OPERATOR{PlusOp,Tokens.PLUS}
    startbyte = ps.t.startbyte
    arg = @precedence ps 13 parse_expression(ps)
    if arg isa EXPR{TupleH}
        ret = EXPR{Call}(EXPR[ret; arg.args], ret.span + arg.span, Variable[], "")
    else
        ret = EXPR{UnaryOpCall}(EXPR[ret, arg], ret.span + arg.span, Variable[], "")
    end
    return ret
end

function parse_call(ps::ParseState, ret::EXPR{OP}) where OP <: OPERATOR{PlusOp,Tokens.MINUS}
    startbyte = ps.t.startbyte
    arg = @precedence ps 13 parse_expression(ps)
    if arg isa EXPR{TupleH}
        ret = EXPR{Call}(EXPR[ret; arg.args], ret.span + arg.span, Variable[], "")
    else
        ret = EXPR{UnaryOpCall}(EXPR[ret, arg], ret.span + arg.span, Variable[], "")
    end
    return ret
end

function parse_call(ps::ParseState, ret)
    startbyte = ps.t.startbyte

    next(ps)
    ret = EXPR{Call}(EXPR[ret, INSTANCE(ps)], ret.span - ps.t.startbyte, Variable[], "")
    format_lbracket(ps)
    @default ps @closer ps paren parse_comma_sep(ps, ret)
    next(ps)
    push!(ret.args, INSTANCE(ps))
    format_rbracket(ps)
    ret.span += ps.nt.startbyte

    # if length(ret.args) > 0 && ismacro(ret.args[1])
    #     ret.head = MACROCALL
    # end
    # if ret.head isa HEAD{Tokens.CCALL} && length(ret.args) > 1 && ret.args[2] isa IDENTIFIER && (ret.args[2].val == :stdcall || ret.args[2].val == :fastcall || ret.args[2].val == :cdecl || ret.args[2].val == :thiscall)
    #     arg = splice!(ret.args, 2)
    #     push!(ret.args, EXPR(arg, [], arg.span))
    # end

    return ret
end

function parse_comma_sep(ps::ParseState, ret::EXPR, kw = true, block = false, formatcomma = true)
    startbyte = ps.nt.startbyte

    @catcherror ps startbyte @nocloser ps inwhere @noscope ps @nocloser ps newline @closer ps comma while !closer(ps)
        block && (ps.trackscope = true)
        a = parse_expression(ps)

        if kw && !ps.closer.brace && a isa EXPR{BinarySyntaxOpCall} && a.args[2] isa EXPR{OPERATOR{AssignmentOp,Tokens.EQ,false}}
            a = EXPR{Kw}(a.args, a.span, Variable[], "")
            # remove format message for kw args
            if !isempty(ps.diagnostics) && ps.nt.startbyte - a.args[3].span - a.args[2].span <= last(last(ps.diagnostics).loc) <= ps.nt.startbyte - a.args[3].span
                pop!(ps.diagnostics)
            end
        end
        push!(ret.args, a)
        if ps.nt.kind == Tokens.COMMA
            next(ps)
            push!(ret.args, INSTANCE(ps))
            if formatcomma
                format_comma(ps)
            else
                format_no_rws(ps)
            end
        end
        if ps.ws.kind == SemiColonWS
            break
        end
    end

    if ps.ws.kind == SemiColonWS
        if block
            body = EXPR{Block}(EXPR[pop!(ret.args)], 0, Variable[], "")
            body.span = body.args[1].span
            # if last(body.args) isa EXPR{BinarySyntaxOpCall} && last(body.args).args[2] isa EXPR{OP} where OP <: OPERATOR{AssignmentOp,Tokens.EQ}
            #     _track_assignment(ps, last(body.args).args[1], last(body.args).args[3], last(body.args).defs)
            # end
            @nocloser ps newline @closer ps comma while @nocloser ps semicolon !closer(ps)
                @catcherror ps startbyte a = parse_expression(ps)
                push!(body.args, a)
                body.span += a.span
            end
            push!(ret.args, body)
            return body
        else
            ps.nt.kind == Tokens.RPAREN && return
            paras = EXPR{Parameters}(EXPR[], -ps.nt.startbyte, Variable[], "")
            @nocloser ps inwhere @nocloser ps newline @nocloser ps semicolon @closer ps comma while !closer(ps)
                @catcherror ps startbyte a = parse_expression(ps)
                if kw && !ps.closer.brace && a isa EXPR{BinarySyntaxOpCall} && a.args[2] isa EXPR{OPERATOR{AssignmentOp,Tokens.EQ,false}}
                    a = EXPR{Kw}(a.args, a.span, Variable[], "")
                    # remove format message for kw args
                    if !isempty(ps.diagnostics) && ps.nt.startbyte - a.args[3].span - a.args[2].span <= last(last(ps.diagnostics).loc) <= ps.nt.startbyte - a.args[3].span
                        pop!(ps.diagnostics)
                    end
                end
                push!(paras.args, a)
                if ps.nt.kind == Tokens.COMMA
                    next(ps)
                    push!(paras.args, INSTANCE(ps))
                    format_comma(ps)
                end
            end
            paras.span += ps.nt.startbyte
            push!(ret.args, paras)
        end
    end
end

function _get_sig_defs!(sig1)
    params = _get_fparams(sig1)
    sig1.defs = Variable[Variable(p, :DataType, sig1) for p in params]

    sig = sig1
    while sig isa EXPR{BinarySyntaxOpCall} && (sig.args[2] isa EXPR{OPERATOR{DeclarationOp,Tokens.DECLARATION,false}} || sig.args[2] isa EXPR{OPERATOR{WhereOp,Tokens.WHERE,false}})
        if sig.args[2] isa EXPR{OPERATOR{WhereOp,Tokens.WHERE,false}}
            haswhere = true
        end
        sig = sig.args[1]
    end
    for i = 2:length(sig.args)
        arg = sig.args[i]
        if !(arg isa EXPR{P} where P <: PUNCTUATION)
            a = _arg_id(arg)
            !(a isa EXPR{IDENTIFIER}) && continue
            t = get_t(arg)
            push!(sig1.defs, Variable(Symbol(a.val), t, sig1))
        end
    end
end

# NEEDS FIX
_arg_id(x) = x
_arg_id(x::EXPR{IDENTIFIER}) = x
_arg_id(x::EXPR{Quotenode}) = x.val
_arg_id(x::EXPR{Curly}) = _arg_id(x.args[1])
_arg_id(x::EXPR{Kw}) = _arg_id(x.args[1])

function _arg_id(x::EXPR{UnarySyntaxOpCall})
    if x.args[2] isa EXPR{OPERATOR{7,Tokens.DDDOT,false}}
        return _arg_id(x.args[1])
    else
        return x
    end
end

function _arg_id(x::EXPR{BinarySyntaxOpCall})
    if x.args[2] isa EXPR{OPERATOR{DeclarationOp,Tokens.DECLARATION,false}} || x.args[2] isa EXPR{OPERATOR{WhereOp,Tokens.WHERE,false}}
        return _arg_id(x.args[1])
    else
        return x
    end
end

_get_fparams(x::EXPR, args = Symbol[]) = args

function _get_fparams(x::EXPR{Call}, args = Symbol[])
    if x.args[1] isa EXPR{Curly}
       _get_fparams(x.args[1], args)
    end
    unique(args)
end

function _get_fparams(x::EXPR{Curly}, args = Symbol[])
    for i = 3:length(x.args)
        a = x.args[i]
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            if a isa EXPR{IDENTIFIER}
                push!(args, Expr(a))
            elseif a isa EXPR{BinarySyntaxOpCall} && a.args[2] isa EXPR{OPERATOR{ComparisonOp,Tokens.ISSUBTYPE,false}}
                push!(args, Expr(a).args[1])
            end
        end
    end
    unique(args)
end

function _get_fparams(x::EXPR{BinarySyntaxOpCall}, args = Symbol[])
    if x.args[2] isa EXPR{OPERATOR{WhereOp,Tokens.WHERE,false}}
        if x.args[1] isa EXPR{BinarySyntaxOpCall} && x.args[1].args[2] isa EXPR{OPERATOR{WhereOp,Tokens.WHERE,false}}
            _get_fparams(x.args[1], args)
        end
        for i = 3:length(x.args)
            a = x.args[i]
            if !(a isa EXPR{P} where P <: PUNCTUATION)
                if a isa EXPR{IDENTIFIER}
                    push!(args, Expr(a))
                elseif a isa EXPR{BinarySyntaxOpCall} && a.args[2] isa EXPR{OPERATOR{ComparisonOp,Tokens.ISSUBTYPE,false}}
                    push!(args, Expr(a).args[1])
                end
            end
        end
    end
    return unique(args)
end

_get_fname(sig::EXPR{FunctionDef}) = _get_fname(sig.args[2])
_get_fname(sig::EXPR{IDENTIFIER}) = sig
_get_fname(sig::EXPR{Tuple}) = NOTHING
function _get_fname(sig::EXPR{BinarySyntaxOpCall})
    if sig.args[2] isa EXPR{OPERATOR{DeclarationOp,Tokens.DECLARATION,false}} || sig.args[2] isa EXPR{OPERATOR{WhereOp,Tokens.WHERE,false}}
        return _get_fname(sig.args[1])
    else
        return get_id(sig.args[1])
    end
end
_get_fname(sig) = get_id(sig.args[1])

_get_fsig(fdecl::EXPR{FunctionDef}) = fdecl.args[2]
_get_fsig(fdecl::EXPR{BinarySyntaxOpCall}) = fdecl.args[1]

declares_function(x) = false
declares_function(x::EXPR{FunctionDef}) = true
declares_function(x::EXPR{BinarySyntaxOpCall}) = x.args[2] isa EXPR{OPERATOR{AssignmentOp,Tokens.EQ,false}} && x.args[1] isa EXPR{Call}


text:
module SecureIO

using Nettle

struct SecureSocket <: IO
    socket::IO
    size::Union{Integer,Nothing}
    ch::Channel{UInt8}
    enc::Encryptor
    dec::Decryptor
end

function SecureSocket(socket::IO,key;size=nothing)
    size!=nothing && @assert mod(size,16)==0

    key32 = hexdigest("sha256", "$key")[1:32]
    enc = Encryptor("AES256", key32)
    dec = Decryptor("AES256", key32)

    SecureSocket(socket,size,Channel{UInt8}(Inf),enc,dec)
end

import Base.isopen
isopen(s::SecureSocket) = isopen(s.socket)

import Base.close
close(s::SecureSocket) = close(s.socket)

function addpadding(text::Vector{UInt8},size)
    # Only last two bytes are used to encode the padding boundary. That limits the possible size.
    @assert length(text) + 2 <= size <= 2^16

    endbytes = reinterpret(UInt8, Int16[length(text)])
    paddedtext = [text; UInt8[0 for i in 1:(size - length(text) - 2)]; endbytes]

    return paddedtext
end

function trimpadding(text::Vector{UInt8})
    endbytes = text[end-1:end]
    n = reinterpret(Int16,endbytes)[1]
    return text[1:n]
end

function stack(io::IO,msg::Vector{UInt8})
    frontbytes = reinterpret(UInt8,Int16[length(msg)])
    item = UInt8[frontbytes...,msg...]
    write(io,item)
end

function unstack(io::IO)
    sizebytes = [read(io,UInt8),read(io,UInt8)]
    size = reinterpret(Int16,sizebytes)[1]

    msg = UInt8[]
    for i in 1:size
        push!(msg,read(io,UInt8))
    end
    return msg
end

function unstack(io::IOBuffer)
    bytes = take!(io)
    size = reinterpret(Int16,bytes[1:2])[1]
    msg = bytes[3:size+2]
    if length(bytes)>size+2
        write(io,bytes[size+3:end])
    end
    return msg
end

import Base.write
function _write_fixed(s::SecureSocket,plaintext::Vector{UInt8},size)
    if size - 2 < length(plaintext)
        error("Message with length $(length(plaintext)) does not fit in $size - 2 bytes")
    end

    #paddedtext = add_padding_PKCS5(Vector{UInt8}(plaintext), size)
    paddedtext = addpadding(plaintext, size)

    msgenc = encrypt(s.enc,paddedtext)
    stack(s.socket,msgenc)
end

function _write_elastic(s::SecureSocket,plaintext::Vector{UInt8})
    n = length(plaintext)

    if mod(n+2,16)==0
        size = n+2
    else
        size = (div(n+2,16) + 1)*16
    end

    paddedtext = addpadding(plaintext, size)
    msgenc = encrypt(s.enc,paddedtext)
    stack(s.socket,msgenc)
end

function _write(s::SecureSocket,msg::Vector{UInt8})
    if s.size==nothing
        _write_elastic(s,msg)
    else
        _write_fixed(s,msg,s.size)
    end
end

_write(s::SecureSocket,msg::UInt8) = _write(s,UInt8[msg])
_write(s::SecureSocket,msg::String) = _write(s,Vector{UInt8}(msg))

write(s::SecureSocket,msg::UInt8) = _write(s,msg)
write(s::SecureSocket,msg::String) = _write(s,msg)
write(s::SecureSocket,msg::Vector{UInt8}) = _write(s,msg)

function getbytes(s::SecureSocket)
    ciphertext = unstack(s.socket) ### If s.socket is a buffer.

    deciphertext = decrypt(s.dec,ciphertext)

    #str = trim_padding_PKCS5(deciphertext)
    msg = trimpadding(deciphertext)

    if typeof(msg)==UInt8
        return msg
    else
        return take!(IOBuffer(msg))
    end
end

import Base.read
function read(s::SecureSocket,x::Type{UInt8})
    if isready(s.ch)==false
        bytes = getbytes(s)
        if typeof(bytes)==UInt8
            put!(s.ch,bytes)
        else
            for byte in bytes
                put!(s.ch,byte)
            end
        end
    end

    return take!(s.ch)
end

export SecureSocket, read, write

end # module


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module HighlyAbundant

using IterTools, NumberTheory, RecordSearch

export ModuleHighlyAbundant
export I002093, F002093, L002093, V002093
export I034885, F034885, L034885, V034885

"""

* I002093, F002093, L002093, V002093, I034885, F034885, L034885, V034885
"""
const ModuleHighlyAbundant = ""

# -----------------------------------------------------------------
# Indices!

"""

Iterate over the first ``n`` highly abundant numbers.
"""
I002093(n) = Records(σ, n, true, true)

"""

Iterate over the highly abundant numbers which do not exceed ``n`` (``1 ≤ i ≤ n``).
"""
F002093(n) = Records(σ, n, false, true)

"""

Return the first ``n`` highly abundant numbers as an array.
"""
L002093(n) = collect(I002093(n))

"""

Return the ``n``-th highly abundant number. Indexing starts with 1.
"""
V002093(n) = nth(I002093(n), n)
# V002093(n) = nth(I002093(n+1), n+1) Indexing starts at 0

# -----------------------------------------------------------------
# Values!

"""

Iterate over the first ``n`` record values of sigma.
"""
I034885(n) = Records(σ, n, true, false)

"""

Iterate over the record values of sigma the indices of which do not exceed ``n`` (``1 ≤ r ≤ n``).
"""
F034885(n) = Records(σ, n, false, false)

"""

Return the first ``n`` record values of sigma as an array.
"""
L034885(n) = collect(I034885(n))

"""

Return the ``n``-th record of sigma. Indexing starts with 1.
"""
V034885(n) = nth(I034885(n), n)
# V034885(n) = nth(I034885(n+1), n+1) Indexing starts at 0.

#START-TEST-########################################################

function test() end

function demo()
    # 1, 2, 3, 4, 6, 8, 10, 12, 16, 18, 20, 24, 30, 36, 42, 48, 60, ...

    println("\nIterate over the first 20 highly abundant numbers.")
    for r ∈ I002093(20)
        print(r, ", ")
    end
    println("...")

    println("\nReturn the first 20 highly abundant numbers as an array.")
    println(L002093(20))

    println("\nIterate over the highly abundant numbers which do not exceed 20.")
    for r ∈ F002093(20)
        print(r, ", ")
    end
    println("...")

    println("\nReturn the 9-th highly abundant number.")
    println("9 -> ", V002093(9))
    println()

    # 1, 3, 4, 7, 12, 15, 18, 28, 31, 39, 42, 60, 72, 91, 96, 124, 168, ...

    println("\nIterate over the first 20 record values of sigma.")
    for r ∈ I034885(20)
        print(r, ", ")
    end
    println("...")

    println("\nReturn the first 20 record values of sigma as an array.")
    println(L034885(20))

    println("\nIterate over the record values of sigma the indices of which do not exceed 20.")
    for r ∈ F034885(20)
        print(r, ", ")
    end
    println("...")

    println("\nReturn the 9-th record of sigma.")
    println("9 -> ", V034885(9))
    println()
end

"""

for n ∈ 1:100 V034885(n)
    0.433958 seconds (802.79 k allocations: 12.632 MiB)
"""
function perf()
    GC.gc()
    @time (for n ∈ 1:100 V034885(n) end)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=
Iterate over the first 20 highly abundant numbers.
1, 2, 3, 4, 6, 8, 10, 12, 16, 18, 20, 24, 30, 36, 42, 48, 60, 72, 84, 90, ...

Return the first 20 highly abundant numbers as an array.
fmpz[1, 2, 3, 4, 6, 8, 10, 12, 16, 18, 20, 24, 30, 36, 42, 48, 60, 72, 84, 90]

Iterate over the highly abundant numbers which do not exceed 20.
1, 2, 3, 4, 6, 8, 10, 12, 16, 18, 20, ...

Return the 9-th highly abundant number. (Indexing starts at 1.)
9 -> 16

Iterate over the first 20 record values of sigma.
1, 3, 4, 7, 12, 15, 18, 28, 31, 39, 42, 60, 72, 91, 96, 124, 168, 195, 224, 234, ...

Return the first 20 record values of sigma as an array.
fmpz[1, 3, 4, 7, 12, 15, 18, 28, 31, 39, 42, 60, 72, 91, 96, 124, 168, 195, 224, 234]

Iterate over the record values of sigma the indices of which do not
exceed 20.
1, 3, 4, 7, 12, 15, 18, 28, 31, 39, 42, ...

Return the 9-th record of sigma. (Indexing starts at 1.)
9 -> 31

0.418483 seconds (842.91 k allocations: 13.243 MiB)
 =#

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
using NewtonsMethod
using Test

@testset "newtonroot" begin
    f(x) = x^2
    fp(x) = 2*x
    @test newtonroot(f, fp; x0=1.0) ≈ 0.0 atol=1.0e-7
    @test newtonroot(f; x0=1.0) ≈ 0.0 atol=1.0e-7
end


text:
# Example showing the use of waveform triggering and segmented capture
using AcqSynth, PlotlyJS
const acq = AcqSynth

# As in the first example, we setup the board
boardnum = 0
acq.init_board(boardnum)

# Let's change to single channel mode
chan = acq.AIN0
acq.set_channels(boardnum,chan)
# We can verify that this worked properly
acq.get_channels(boardnum) == chan

# Check the frequency, in single channel mode, this will be twice the normal
acq.get_frequency(boardnum)

# Configure the board for waveform triggering
hysteresis = 32
threshold = 2^11-220 # middle value for a 12-bit board
ttype = acq.WAVEFORM_TRIGGER
# This function configures the waveform parameters
acq.set_waveform_trigger_params(boardnum,threshold,hysteresis)
# This function sets the trigger type
acq.set_trigger(boardnum,ttype,acq.RISING_EDGE,chan)

# Configure segmented capture
count = 16 # number of segments to acquire
depth = 2^8 # number of samples per segment
acq.set_segmented_capture(boardnum,count,depth)
acq.get_segmented_capture(boardnum) == (count,depth)

acq.set_pretrigger_mem(boardnum,0)

# Setup acquisition, the card will wait until a trigger is detected
numblocks = 100 # let's grab 2 blocks
acq.setup_acquire(boardnum,numblocks)

# Grab the data from the card
data = acq.get_samples_12(boardnum,numblocks) # transfer data as 12-bit samples

# Plot the stuff
t = (1:length(data))./4E9 # time vector, we are sampling at 4GSPS
plot(scatter(x=t[1:2000],y=data[1:2000],name="AIN0"))

# # We can configure other things
# # ECL trigger
# acq.get_ECL_trigger_delay(boardnum)
# acq.get_ECL_trigger(boardnum)
# acq.set_ECL_trigger(boardnum,0)
# acq.set_ECL_trigger_delay(boardnum,512)
# # Configure averager
# count = 2^6 # number of segments to acquire
# depth = 2^12 # number of samples per segment
# acq.set_averager(boardnum,count,depth)
# acq.get_averager(boardnum) == (count,depth)