text:
using KernelAbstractions
using Test

import KernelAbstractions.NDIteration: NDRange, StaticSize, DynamicSize

@kernel function index(A)
    I = @index(Global, NTuple)
    @show A[I...]
end

@kernel function literal_pow(A)
    A[1] = 2^11
end

@kernel function checked(A, a, b)
    A[1] = Base.Checked.checked_add(a, b)
end

function compiler_testsuite()
    kernel = index(CPU(), DynamicSize(), DynamicSize())
    iterspace = NDRange{1, StaticSize{(128,)}, StaticSize{(8,)}}();
    ctx = KernelAbstractions.mkcontext(kernel, 1, nothing, iterspace, Val(KernelAbstractions.NoDynamicCheck()))

    @test KernelAbstractions.Cassette.overdub(ctx, KernelAbstractions.__index_Global_NTuple, CartesianIndex(1)) == (1,)

    let (CI, rt) = @ka_code_typed literal_pow(CPU())(zeros(Int,1), ndrange=1)
        # test that there is no invoke of overdub
        @test !any(stmt->(stmt isa Expr) && stmt.head == :invoke, CI.code)
    end

    if VERSION >= v"1.5"
        let (CI, rt) = @ka_code_typed checked(CPU())(zeros(Int,1), 1, 2, ndrange=1)
            # test that there is no invoke of overdub
            @test !any(stmt->(stmt isa Expr) && stmt.head == :invoke, CI.code)
        end
    end
end


text: module TestModuleA end

ext: jl


text:
using covid19abm
using Test

@testset "demographics" begin
    initialize()
end


text:
function IntegerBootstrapping(arg0::jdouble)
    return IntegerBootstrapping((jdouble,), arg0)
end


text:
"""
Methods for the query generator: specifically, to

1. generate sparsity coefficients b and subsampling matrices M
2. get the indices of a signal subsample
3. compute a subsampled and delayed Walsh-Hadamard transform.
"""

# potential TODO figure out how to load in parts of the signal bit by bit
# see how Amirali did it

include("input_signal.jl")
using StatsBase
using Random

"""
A semi-arbitrary fixed choice of the sparsity coefficient.
See get_b for full signature.
"""
function get_b_simple(signal::Signal)
    # I don't know why, but changing this to signal.n \div 3 makes it break
    if signal.n == 4
        return 2
    elseif (signal.n == 8) || (signal.n == 16)
        return 4
    elseif (signal.n < 16)
        return signal.n ÷ 2
    else
        return signal.n ÷ 4
    end
end

get_b_lookup = Dict(
    :simple => get_b_simple,
    :random => get_b_simple
)

"""
Get the sparsity coefficient for the signal.

Arguments
---------
signal : TestSignal
The signal whose WHT we want.

method : Symbol
The method to use. All methods referenced must use this signature (minus "method".)

Returns
-------
b : Int64
The sparsity coefficient.
"""
function get_b(signal::Signal; method=:simple)::Int64
    return get_b_lookup[method](signal)
end

function get_partial_I(n, b, k)
    M = falses(n, b)
    for j in 1:b
        M[k+j, j] = 1
    end
    return M
end

"""
A semi-arbitrary fixed choice of the subsampling matrices. See get_Ms for full signature.
"""
function get_Ms_simple(n::Int64, b::Int64)
    # @assert n % b == 0 "n must be exactly divisible by b"
    num_to_get = Int64(ceil(n / b))

    Ms = BitArray[]
    push!(Ms, get_partial_I(n, b, n-b))
    for i in num_to_get-2:-1:0
        push!(Ms, get_partial_I(n, b, b*i))
    end
    return Ms
end

function get_Ms_random(n::Int64, b::Int64)
    num_to_get = n ÷ b
    Ms = BitArray[]

    for i in 1:num_to_get
        push!(Ms, bitrand(n, b))
    end
    return Ms
end

get_Ms_lookup = Dict(
    :simple => get_Ms_simple,
    :random => get_Ms_random
)

function get_Ms(n::Int64, b::Int64; method=:simple)::Array{BitArray,1}
    return get_Ms_lookup[method](n, b)
end

"""
Gets the delays matrix [0; I], of dimension (n+1, n). See get_D for full signature.
"""
function get_D_identity_like(n::Int64; kwargs...)
    return BitArray(vcat(zeros(1, n), Matrix{Bool}(I, n, n)))
end

"""
Gets a random delays matrix of dimension (num_delays, n). See get_D for full signature.
"""
function get_D_random(n::Int64; kwargs...)
    num_delays = kwargs[:num_delays]
    choices = sample(0:2^n-1, num_delays, replace=false)
    return @pipe dec_to_bin.(choices, n) |> hcat(_...) |> transpose |> BitArray
end

"""
Get a repetition code based (NSO-SPRIGHT) delays matrix. See get_D for full signature.
Not sure of correctness based on the paper: this is not (num_delays, n).
"""
function get_D_nso(n::Int64; kwargs...)::BitArray{2}
    num_delays = kwargs[:num_delays]
    p1 = num_delays ÷ n
    random_offsets = get_D_random(n; num_delays=p1)
    D = BitArray(undef, 0, n)
    identity_like = get_D_identity_like(n)
    for row in eachrow(random_offsets)
        modulated_offsets = @pipe [(row .⊻ r) for r in eachrow(identity_like)] |> hcat(_...) |> transpose
        D = vcat(D, modulated_offsets)
    end
    return D
end

function get_generator_matrix(codelen::Int64, n::Int64, code::Symbol)
    return zeros((codelen, n))
end

"""
Get a channel code based (SO-SPRIGHT) delays matrix. See get_D for full signature.
Returns a matrix of dimension (p1+p2+p3, n).
"""
function get_D_so(n::Int64; kwargs...)::BitArray{2}
    p1, p2, p3 = kwargs[:num_delays]
    code = kwargs[:code]
    D1 = get_D_random(n, num_delays=p1)
    D2 = zeros((p2, n))
    D3 = get_generator_matrix(p3, n, code)
    return vcat(D1, D2, D3)
end

get_D_lookup = Dict(
    :identity_like => get_D_identity_like,
    :random => get_D_random,
    :nso => get_D_nso,
    :so => get_D_so
)

"""
Delay generator: gets a delays matrix.

Arguments
---------
n : Int64

number of bits: log2 of the signal length.

Returns
-------
D : num_delays×n Array{Bool,2}

The delays matrix; if num_delays is not specified in kwargs, see the relevant sub-function for a default.
"""
function get_D(n::Int64; method=:identity_like, kwargs...)
    return get_D_lookup[method](n; kwargs...)
end

"""
Query generator: creates indices for signal subsamples.
Arguments
---------
M : n×b Array{Bool,2}
The subsampling matrix.

d : n-element Array{Bool,1}
The subsampling offset.

Returns
-------
indices : B-element Array{Int64,1}
The (decimal) subsample indices.
"""
function subsample_indices(M, d)
    L = binary_ints(size(M, 2))
    inds_binary = Bool.((M*L .⊻ d) .% 2)
    return @pipe inds_binary |> eachcol |> collect |> bin_to_dec.(_) .+ 1
end

"""
Creates random delays, subsamples according to M and the random delays,
and returns the subsample WHT along with the delays.

Arguments
---------
signal : Signal
The signal to subsample, delay, and compute the WHT of.

M : n×b Array{Bool,2}
The subsampling matrix.

num_delays : Int64
The number of delays to apply; or, the number of rows in the delays matrix.

force_identity_like : Bool
Whether to make D = [0; I] like in the noiseless case; for debugging.
"""
function compute_delayed_wht(signal::Signal, M, D)
    return compute_delayed_subtransform(signal, M, D, fwht)
end

function compute_delayed_fft(signal::Signal, M, D)
    return compute_delayed_subtransform(signal, M, D, fft)
end

function compute_delayed_subtransform(signal::Signal, M, D, transform::Function)
    inds = map(d -> subsample_indices(M, d), D |> eachrow |> collect)
    used_inds = reduce(union, inds)
    samples_to_transform = map(x -> signal[x], inds)
    return hcat(transform.(samples_to_transform)...), used_inds
end


text:
# # ::: MAIN FUNCTIONS
# function DenseLayer(input_size, output_size, batch_size)
#     # w init = glorot_uniform from https://github.com/FluxML/Flux.jl/blob/master/src/utils.jl
#     # b init = 0
#     w = [(CuArrays.rand(output_size, input_size) .- 0.5) .* sqrt(24/(output_size+input_size)) for _ in 1:batch_size]
#     b = [CuArrays.zeros(output_size) for _ in 1:batch_size]
#
#     dlocaly_dlocalX = [CuArrays.zeros(output_size, input_size) for _ in 1:batch_size]
#     dlocaly_dw = [CuArrays.zeros(input_size) for _ in 1:batch_size]
#
#     function forward!(x)
#         y = [CuArrays.zeros(output_size) for i in 1:batch_size]
#
#         for i in 1:batch_size
#
#             t = time()
#
#             y[i] = w[i] * x[i] + b[i]
#
#             println("post lin: ", time()-t)
#
#             dlocaly_dlocalX[i] = w[i]
#             dlocaly_dw[i] = x[i]
#         end
#         return y
#     end
#
#     function backward!(Δ, lr)
#         dprediction_dlocalX = [CuArrays.zeros(input_size) for _ in 1:batch_size]
#
#         for i in 1:batch_size
#             dprediction_dlocalX[i] .= dlocaly_dlocalX[i]' * Δ[i]
#
#             w[i] .-= lr .* (Δ[i] * dlocaly_dw[i]')
#             b[i] .-= lr .* Δ[i]
#
#             dlocaly_dw[i] .= 0.
#             dlocaly_dlocalX[i] .= 0.
#         end
#
#         return dprediction_dlocalX
#     end
#
#     return forward!, backward!
# end
#
# function ActivationLayer(in_out_size, batch_size, act_f, dact_f)
#     dlocaly_dlocalX = [CuArrays.zeros(in_out_size, in_out_size) for _ in 1:batch_size]
#
#     function forward!(x)
#         y = [similar(x[i]) for i in 1:batch_size]
#
#         for i in 1:batch_size
#             dlocaly_dlocalX[i] = dact_f(x[i])
#
#             y[i] = act_f(x[i])
#         end
#         return y
#     end
#
#     function backward!(Δ, lr)
#         dprediction_dlocalX = [similar(Δ[i]) for i in 1:batch_size]
#
#         for i in 1:batch_size
#             dprediction_dlocalX[i] = dlocaly_dlocalX[i] * Δ[i]
#         end
#
#         return dprediction_dlocalX
#     end
#
#     return forward!, backward!
# end
#
# function BatchNormLayer(in_out_size, batch_size)
#     α = [CuArrays.rand(in_out_size) for _ in 1:batch_size]
#     β = [rand() for _ in 1:batch_size]
#     ϵ = 0.00001 # to prevent devision by zero
#
#     dlocaly_dα = [CuArrays.zeros(in_out_size) for _ in 1:batch_size]
#     dlocaly_dlocalX = [CuArrays.zeros(in_out_size) for _ in 1:batch_size]
#
#     function forward!(x) # ϵ to avoid dividing by zero
#         y = [similar(x[i]) for i in 1:batch_size]
#
#         N = sum([length(x[i]) for i in 1:batch_size])
#         μ = sum([sum(x[i]) for i in 1:batch_size]) / N
#         σ = sum([sum((x[i] .- μ) .^ 2) for i in 1:batch_size]) / N
#
#         for i in 1:batch_size
#             bn = (x[i] .- μ) ./ sqrt(σ + ϵ)
#             y[i] = α[i] .* bn .+ β[i]
#
#             dlocaly_dlocalX[i] = α[i] .* sqrt(σ + ϵ) # d_BN(x[i])
#             dlocaly_dα[i] = bn
#         end
#
#         return y
#     end
#
#     function backward!(Δ, lr)
#         dprediction_dlocalX = [similar(Δ[i]) for i in 1:batch_size]
#
#         for i in 1:batch_size
#             dprediction_dlocalX[i] = Δ[i] .* dlocaly_dlocalX[i]
#
#             α[i] .-= lr .* (Δ[i] .* dlocaly_dα[i])
#             β[i] -= lr * sum(Δ[i])
#         end
#
#         return dprediction_dlocalX
#     end
#
#     return forward!, backward!
# end
#
#
# function feed_forward!(x, forward_layer_functions::Array)
#     for ff_layer in forward_layer_functions
#         x = ff_layer(x)
#     end
#     return x
# end
#
# function feed_backwards!(d_pred, backward_layer_functions::Array; lr=0.006)
#     for fb_layer in reverse(backward_layer_functions)
#         d_pred = fb_layer(d_pred, lr)
#     end
#     return d_pred
# end


text:
# Implements following functions
# - descent path: monotonically decreasing states visited (best solution throughout search)
# - acceptance path (current solution throughout search)

@recipe function descent_path(df::DataFrame, type::Type{SearchTrace}; xax = :iteration, yax = :value)
    ldf = df[:, :accepted .== true & :delta .<= 0]
    x = ldf[:,xax]
    y = ldf[:,yax]

    seriestype --> :scatter
    return x,y
end

max_stars_repo_name: Michiel-VL/PolyaViz

max_issues_repo_name: Michiel-VL/PolyaViz


text:
const WEIGHT          = 1
const INCOME          = 2
const WEIGHTED_INCOME = 3
const POPN_ACCUM      = 4
const INCOME_ACCUM    = 5
const DEFAULT_FGT_ALPHAS = [ 0.0, 0.50, 1.0, 1.50, 2.0, 2.5 ];
const DEFAULT_ATKINSON_ES = [ 0.25, 0.50, 0.75, 1.0, 1.25, 1.50, 1.75, 2.0, 2.25 ];
const DEFAULT_ENTROPIES = [ 1.25, 1.50, 1.75, 2.0, 2.25, 2.5 ];

const OutputDict = Dict{ Symbol, Any }
const OutputDictArray = Array{ OutputDict, 1 }

function sortAndAccumulate( aug :: Array{<:Real,2}, sortdata :: Bool, nrows :: Integer ) :: Array{<:Real,2}
    if sortdata
        aug = sortslices( aug, alg=QuickSort, dims=1,lt=((x,y)->isless(x[INCOME],y[INCOME])))
    end
    cumulative_weight :: Float64 = 0.0
    cumulative_income :: Float64 = 0.0
    for row in 1:nrows
            cumulative_weight += aug[row,WEIGHT]
            cumulative_income += aug[row,WEIGHTED_INCOME]
            aug[row,POPN_ACCUM] = cumulative_weight
            aug[row,INCOME_ACCUM] = cumulative_income
    end
    return aug
end # sortAndAccumulate

"
internal function that makes a sorted array
with cumulative income and population added
"
function makeaugmented(
    data,
    weightcol :: Symbol,
    incomecol :: Symbol,
    sortdata  :: Bool = true,
    deletenegatives :: Bool = true ) :: Array{Float64,2}
    @assert TableTraits.isiterabletable( data ) "data needs to implement IterableTables"
    iter = IteratorInterfaceExtensions.getiterator(data)
    nrows = length(iter)
    aug = zeros( nrows, 5 )
    r = 0
    for row in iter
        if ( ismissing( row[weightcol])) || ( ismissing( row[incomecol])) ||
            ( deletenegatives && row[incomecol] < 0.0 )
            ;
        else
            r += 1
            aug[r,WEIGHT] = get(row[weightcol]) ## this is the datavalue thing; see:
            aug[r,INCOME] = get(row[incomecol])
            aug[r,WEIGHTED_INCOME] = aug[r,WEIGHT]*aug[r,INCOME]
        end # not missing or negative
    end
    aug = aug[1:r,:]
    aug = sortAndAccumulate( aug, sortdata, r )
    return aug
end

"
internal function that makes a sorted array
with cumulative income and population added
"
function makeaugmented(
    data            :: Array{<:Real,2},
    weightpos       :: Integer = 1,
    incomepos       :: Integer = 2,
    sortdata        :: Bool = true,
    deletenegatives :: Bool = true ) :: Array{Float64,2}

    nrows = size( data )[1]
    aug = zeros( nrows, 5 )
    r = 0
    for row in 1:nrows
        if( deletenegatives && data[row,incomepos] < 0.0 )
            ;
        else
            r += 1
            aug[r,WEIGHT] = data[row,weightpos]
            aug[r,INCOME] = data[row,incomepos]
            aug[r,WEIGHTED_INCOME] = data[row,incomepos]*data[row,weightpos]
        end # not neg or negs accepted
    end # row loop
    aug = sortAndAccumulate( aug, sortdata, r )
    return aug[1:r,:]
end

"""
calculate a Gini coefficient on one of our sorted arrays
"""
function makegini( data :: Array{Float64, 2} ) :: Float64
    lorenz :: Float64 = 0.0

    nrows = size( data )[1]
    if nrows == 0
        return 0.0
    end
    lastr = data[nrows,:]
    for row in 1:nrows
        lorenz += data[row, WEIGHT]*((2.0*data[row,INCOME_ACCUM]) - data[row,WEIGHTED_INCOME])
    end
    return 1.0-(lorenz/lastr[INCOME_ACCUM])/lastr[POPN_ACCUM]
end

"""
generate a subset of one of our datasets with just the elements whose incomes
are below the line. Probably possible in 1 line, once I get the hang of this
a bit more.
"""
function makeallbelowline( data :: Array{Float64, 2}, line :: Float64 ) :: Array{Float64, 2 }
    nrows = size( data )[1]
    ncols = size( data )[2]
    outa = zeros( Float64, nrows, ncols ) # Array{Float64}( undef, 0, 5 )
    @assert ncols == 5 "data should have 5 cols"
    nout = 0
    for row in 1:nrows
        if data[row,INCOME] < line
            nout += 1
            outa[ nout, : ] .= data[row,:]
        end
    end
    return outa[1:nout,:]
end

"""
Create a dictionary of poverty measures.

This is based on the [World Bank's Poverty Handbook](http://documents.worldbank.org/curated/en/488081468157174849/Handbook-on-poverty-and-inequality)
by  Haughton and Khandker.

Arguments:
* `rawdata` - an nxm array of real nunbers; each row is an observation; one col should be a weight, another is income;
   positions assumed to be 1 and 2 unless weight and incomepos are supplied
* `line` - a poverty line, assumed same for all obs (so income is assumed to be already equivalised)
* `foster_greer_thorndyke_alphas` - coefficients for FGT poverty measures; note that FGT(0)
   corresponds to headcount and FGT(1) to gap; count and gap are computed directly anyway
   but it's worth checking one against the other.
* `growth` is (e.g.) 0.01 for 1% per period, and is used for 'time to exit' measure.
"""
function makepoverty(
    rawdata                       :: Array{<:Real, 2},
    line                          :: Real,
    growth                        :: Real,
    weightpos                     :: Integer = 1,
    incomepos                     :: Integer = 2,
    foster_greer_thorndyke_alphas :: AbstractArray{<:Real, 1} = DEFAULT_FGT_ALPHAS ) :: OutputDict

    data = makeaugmented( rawdata, weightpos, incomepos )
    makepovertyinternal(
        data = data,
        line = line,
        growth = growth,
        foster_greer_thorndyke_alphas = foster_greer_thorndyke_alphas )
end

"""
As above, but using the QueryVerse IterableTables interface
rawdata - basically anything resembling a dataframe; see: [https://github.com/queryverse/IterableTables.jl]
throws an exception if rawdata doesn't support iterabletable interface
"""
function makepoverty(
    rawdata,
    line                          :: Real,
    growth                        :: Real,
    weightcol                     :: Symbol,
    incomecol                     :: Symbol,
    foster_greer_thorndyke_alphas :: AbstractArray{<:Real, 1} = DEFAULT_FGT_ALPHAS,
     ) :: OutputDict
    @assert TableTraits.isiterabletable( rawdata ) "data needs to implement IterableTables"
    data = makeaugmented( rawdata, weightcol, incomecol )
    makepovertyinternal(
        data = data,
        line = line,
        growth = growth,
        foster_greer_thorndyke_alphas = foster_greer_thorndyke_alphas )
end

"""
Internal version, once we have our datatset
"""
function makepovertyinternal(
    ;
    data                          :: Array{Float64, 2},
    line                          :: Real,
    growth                        :: Real = 0.0,
    foster_greer_thorndyke_alphas :: AbstractArray{<:Real, 1} = DEFAULT_FGT_ALPHAS,) :: OutputDict

    pv = Dict{ Symbol, Any}()
    nrows = size( data )[1]
    ncols = size( data )[2]
    population = data[ nrows, POPN_ACCUM ]
    total_income = data[ nrows, INCOME_ACCUM ]

    nfgs = size( foster_greer_thorndyke_alphas )[1]
    @assert ncols == 5 "data should have 5 cols"
    pv[:fgt_alphas] = foster_greer_thorndyke_alphas
    pv[:headcount] = 0.0
    pv[:gap] = 0.0
    pv[:watts] = 0.0
    pv[:foster_greer_thorndyke] = zeros( Float64, nfgs )
    pv[:time_to_exit] = 0.0

    belowline = makeallbelowline( data, line )
    nbrrows = size( belowline )[1]

    pv[:gini_amongst_poor] = makegini( belowline )
    for row in 1:nbrrows
        inc :: Float64= belowline[row,INCOME]
        weight :: Float64 = belowline[row,WEIGHT]
        gap :: Float64 = line - inc
        @assert gap >= 0 "poverty gap must be postive"
        pv[:headcount] += weight
        pv[:gap] += weight*gap/line
        if belowline[row,INCOME ] > 0
            pv[:watts] += weight*log(line/inc)
        end
        for p in 1:nfgs
            fg = foster_greer_thorndyke_alphas[p]
            pv[:foster_greer_thorndyke][p] += weight*((gap/line)^fg)
        end
    end # main loop
    pv[:watts] /= population
    if growth > 0.0
        pv[:time_to_exit] = pv[:watts]/growth
    end
    pv[:gap] /= population
    pv[:headcount] /= population
    pv[:foster_greer_thorndyke] ./= population
    #
    # Gini of poverty gaps; see: WB pp 74-5
    #
    # create a 'Gini of the Gaps'
    # the sort routine in makeaugmented does a really
    # bad job here either because the data
    # is mostly zeros or because it's reverse sorted
    # (smallest income -> biggest gap)
    # we we can just create the dataset in reverse
    # and use that
    gdata = zeros( Float64, nrows, 5 )
    for row in 1:nrows
        gap = max( 0.0, line - data[row,INCOME] )
        gpos = nrows - row + 1
        gdata[gpos,INCOME] = gap;
        gdata[gpos,WEIGHT] = data[row,WEIGHT]
    end
    gdata = makeaugmented( gdata, 1, 2, false )
    pv[:poverty_gap_gini] = makegini( gdata )

    pv[:sen] = pv[:headcount]*pv[:gini_amongst_poor]+pv[:gap]*(1.0-pv[:gini_amongst_poor])
    pv[:shorrocks] = pv[:headcount]*pv[:gap]*(1.0+pv[:poverty_gap_gini])
    return pv
end # makepoverty

"""
Make a wee dict with :theil_between and :theil_within
See WB eqns 6.7/6.8.
TODO
1. there are some papers on decomposing Atkinson, but I
don't understand them ..
2. the over time 3-part version

 popindic : Inequal for the population as a whole
 subindices : an array of dics, one for each subgroup of interest
"""
function adddecomposedtheil( popindic :: OutputDict, subindices :: OutputDictArray ) :: OutputDict
    popn = popindic[:total_population]
    income = popindic[:total_income]
    avinc = popindic[:average_income]
    within = zeros(2)
    between = zeros(2)
    totalpop = 0.0
    totalinc = 0.0
    for ind in subindices
        popshare = ind[:total_population]/popn
        incshare = ind[:total_income]/income
        totalpop += popshare
        totalinc += incshare

        within[1] += ind[:theil][1]*popshare
        between[1] += popshare*log(avinc/ind[:average_income])

        within[2] += ind[:theil][2]*incshare
        between[2] += incshare*log(incshare/popshare)

    end
    overall1 = popindic[:theil][1]
    overall2 = popindic[:theil][2]

    @assert totalpop ≈ 1.0
    @assert totalinc ≈ 1.0
    @assert within[1]+between[1] ≈ popindic[:theil][1]
    @assert within[2]+between[2] ≈ popindic[:theil][2]
    md = OutputDict()
    md[:theil_between] = between
    md[:theil_within] = within
    md
end

"""
Make a dictionary of inequality measures.
This is mainly taken from chs 5 and 6 of the World Bank book.

1. `rawdata` a matrix with cols with weights and incomes
2. `atkinson_es` inequality aversion values for the Atkinson indexes
3. `generalised_entropy_alphas`
4. `weightpos` - column with weights
5. `incomepos` - column with incomes
"""
function makeinequality(
    rawdata                    :: Array{<:Real, 2 },
    weightpos                  :: Integer = 1,
    incomepos                  :: Integer = 2,
    atkinson_es                :: AbstractArray{<:Real, 1} = DEFAULT_ATKINSON_ES,
    generalised_entropy_alphas :: AbstractArray{<:Real, 1} = DEFAULT_ENTROPIES ) :: OutputDict
    data = makeaugmented( rawdata, weightpos, incomepos )
    return makeinequalityinternal(
        data = data,
        atkinson_es = atkinson_es,
        generalised_entropy_alphas = generalised_entropy_alphas
    )
end

"""
As above but using the iterable table interface; see: https://github.com/queryverse/IterableTables.jl
weightcol and incomecol are names of columns
"""
function makeinequality(
    rawdata,
    weightcol                  :: Symbol,
    incomecol                  :: Symbol,
    atkinson_es                :: AbstractArray{<:Real, 1} = DEFAULT_ATKINSON_ES,
    generalised_entropy_alphas :: AbstractArray{<:Real, 1} = DEFAULT_ENTROPIES ) :: OutputDict
    @assert TableTraits.isiterabletable( rawdata ) "data needs to implement IterableTables"
    data = makeaugmented( rawdata, weightcol, incomecol )
    return makeinequalityinternal(
        data = data,
        atkinson_es = atkinson_es,
        generalised_entropy_alphas = generalised_entropy_alphas
    )

end

"""
Chop a dataset with populations and incomes
into numbins groups in a form suitable for
e.g. a Gini curve.
* col1 is cumulative population,
*    2 cumulative income/whatever,
*    3 threshold income level.
"""
function binify(
    rawdata   :: Array{<:Real, 2 },
    numbins   :: Integer,
    weightpos :: Integer = 1,
    incomepos :: Integer = 2 ) :: AbstractArray{<:Real, 2}
    data = makeaugmented( rawdata, weightpos, incomepos, true, true )
    return binifyinternal( data, numbins )
end

"""
As above, but using any DataFrame like thing that supports the TableTraits.isiterabletable interface
"""
function binify(
    rawdata,
    numbins   :: Integer,
    weightcol :: Symbol,
    incomecol :: Symbol ) :: AbstractArray{<:Real, 2}
    @assert TableTraits.isiterabletable( rawdata ) "data needs to implement IterableTables"
    data = makeaugmented( rawdata, weightcol, incomecol, true, true )
    return binifyinternal( data, numbins )
end

## FIXME this can go horribly wrong with small unbalanced numbers of rows -
# try pop=[22,1,1,1,1,1,1,1,1,1 inc==[1,2,2,2,2,2,2,2,2,2]
function binifyinternal(
    data      :: Array{<:Real, 2 },
    numbins   :: Integer ) :: AbstractArray{<:Real, 2}
    nrows = size( data )[1]
    ncols = size( data )[2]
    out = zeros( Float64, numbins, 3 )
    total_population = data[ nrows, POPN_ACCUM ]
    total_income = data[ nrows, INCOME_ACCUM ]
    bin_size :: Float64 = 1.0/numbins
    bno = 0
    thresh = bin_size
    popsharelast = 0.0
    incomelast = 0.0
    incomeaccumlast = 0.0
    incomeaccum = 0.0
    # print( "total_population $total_population total_income $total_income \n")
    for row in 1:nrows
        income = data[row,INCOME]
        incomeaccum = data[row,INCOME_ACCUM]/total_income
        popshare = data[row,POPN_ACCUM]/total_population
        if popshare ≈ thresh
            bno += 1
            out[bno,1] = popshare
            out[bno,2] = incomeaccum
            out[bno,3] = income
            # print( "row $row popshare $popshare ≈ thresh $thresh incomeaccum $incomeaccum \n")
            thresh += bin_size
        elseif( popsharelast < thresh ) && ( popshare > thresh)
            bno += 1
            # print( "row $row popsharelast $popsharelast thresh $thresh popshare $popshare\n" )
            pgap = popshare - popsharelast
            p1 = (thresh - popsharelast)
            p2 = (popshare - thresh)
            out[bno,1] = ((popshare*p2)+(popsharelast*p1))/pgap
            out[bno,2] = ((incomeaccum*p2)+(incomeaccumlast*p1))/pgap
            out[bno,3] = ((income*p2)+(incomelast*p1))/pgap
            thresh += bin_size
        end
        popsharelast = popshare
        incomelast = income
        incomeaccumlast = incomeaccum
    end
    return out
end

"""
Make a dictionary of inequality measures.
This is mainly taken from chs 5 and 6 of the World Bank book.

1. `rawdata` a matrix with cols with weights and incomes
2. `atkinson_es` inequality aversion values for the Atkinson indexes
3. `generalised_entropy_alphas`
4. `weightpos` - column with weights
5. `incomepos` - column with incomes
"""
function makeinequalityinternal(
    ;
    data                       :: Array{Float64, 2},
    atkinson_es                :: AbstractArray{<:Real, 1} = DEFAULT_ATKINSON_ES,
    generalised_entropy_alphas :: AbstractArray{<:Real, 1} = DEFAULT_ENTROPIES ) :: OutputDict
    nrows = size( data )[1]
    ncols = size( data )[2]
    @assert ncols == 5 "data should have 5 cols but has $ncols"

    nats = size( atkinson_es )[1]
    neps = size( generalised_entropy_alphas )[1]
    iq = OutputDict()
    # initialise atkinsons; 1 for e = 1 0 otherwise
    iq[:atkinson_es] = atkinson_es
    iq[:atkinson] = zeros( Float64, nats )
    iq[:atkinson] = [i == 1.0 ? 1.0 : 0.0 for i in atkinson_es]
    iq[:generalised_entropy_alphas] = generalised_entropy_alphas
    iq[:generalised_entropy] = zeros( Float64, neps )
    iq[:negative_or_zero_income_count] = 0
    iq[:hoover] = 0.0
    iq[:theil] = zeros(Float64,2)
    iq[:gini] = makegini( data )

    total_income = data[nrows,INCOME_ACCUM]
    total_population = data[nrows,POPN_ACCUM]
    y_bar = total_income/total_population
    iq[:total_income] = total_income
    iq[:total_population] = total_population
    iq[:average_income] = y_bar
    popsharelast = 0.0
    incomelast = 0.0
    for row in 1:nrows
        income = data[row,INCOME]
        weight = data[row,WEIGHT]
        if income > 0.0
            y_yb  :: Float64 = income/y_bar
            yb_y  :: Float64 = y_bar/income
            ln_y_yb :: Float64 = log( y_yb )
            ln_yb_y :: Float64 = log( yb_y )
            iq[:hoover] += weight*abs( income - y_bar )
            iq[:theil][1] += weight*ln_yb_y
            iq[:theil][2] += weight*y_yb*ln_y_yb
            for i in 1:nats
                    es :: Float64 = iq[:atkinson_es][i]
                    if es != 1.0
                        iq[:atkinson][i] += weight*(y_yb^(1.0-es))
                    else
                        iq[:atkinson][i] *= (income)^(weight/total_population)
                    end # e = 1 case
            end # atkinsons
            for i in 1:neps
                alpha :: Float64 = iq[:generalised_entropy_alphas][i]
                iq[:generalised_entropy][i] += weight*(y_yb^alpha)
            end # entropies
        else
            iq[:negative_or_zero_income_count] += 1
        end # positive income
    end # main loop
    deciles = binify( data, 10 )
    iq[:median] = deciles[5,3]
    # top 10/bottom 40
    iq[:palma] = (1.0-deciles[9,2])/deciles[4,2]
    iq[:deciles] = deciles
    iq[:hoover] /= 2.0*total_income
    for i in 1:neps
        alpha :: Float64 = iq[:generalised_entropy_alphas][i]
        aq :: Float64 = (1.0/(alpha*(alpha-1.0)))
        iq[:generalised_entropy][i] =
            aq*((iq[:generalised_entropy][i]/total_population)-1.0)
    end # entropies
    for i in 1:nats
        es :: Float64 = iq[:atkinson_es][i]
        if es != 1.0
            iq[:atkinson][i] = 1.0 - ( iq[:atkinson][i]/total_population )^(1.0/(1.0-es))
        else
            iq[:atkinson][i] = 1.0 - ( iq[:atkinson][i]/y_bar )
        end # e = 1
    end
    iq[:theil] ./= total_population
    return iq
end # makeinequalityinternal


text:
#=
test_rules_visualization.jl:
- Julia version:
- Author: bramb
- Date: 2019-01-24
=#

using Test

include("util.jl")

@testset "rules_visualization" begin
    using HRGVisualization
    rules = Dict{String, HyperGraph{String}}(
        "S" => HyperGraph([HyperEdge("JOHN", ["_"], ["_"])]),
        "JOHN" => HyperGraph([HyperEdge("John", ["_"], ["3"]), HyperEdge("LOVES",["3"], ["_"])]),
        "LOVES" => HyperGraph([HyperEdge("loves", ["_"], ["4"]), HyperEdge("MARY", ["4"], ["_"])]),
        "MARY" => HyperGraph([HyperEdge("Mary", ["_"], ["_"])])
    )

    @test HRGVisualization._is_nonterminal("S") == true
    @test HRGVisualization._is_nonterminal("JOHN") == true
    @test HRGVisualization._is_nonterminal("loves") == false

    nonterminals = HRGVisualization._get_nonterminals(rules["S"])
    @test nonterminals == ["JOHN"]

    nonterminals = HRGVisualization._get_nonterminals(rules["JOHN"])
    @test nonterminals == ["LOVES"]

    nonterminals = HRGVisualization._get_nonterminals(rules["LOVES"])
    @test nonterminals == ["MARY"]

    nonterminals = HRGVisualization._get_nonterminals(rules["MARY"])
    @test nonterminals == []

    @test HRGVisualization._order_keys(rules) == ["S", "JOHN", "LOVES", "MARY"]

    rules_dot = to_dot(rules)
    expected = """
    digraph HRG_Rules {
        rankdir=LR
        edge [arrowsize=0.5]
        subgraph cluster_1 {
            label="MARY ⇒"
            rule1_n1_1[shape=circle;width=0.05;label=""]
            rule1_n1_2[shape=circle;width=0.05;label=""]
            rule1_he1[shape=box;label="Mary"]
            {rule1_n1_1} -> rule1_he1 -> {rule1_n1_2}
        }
        subgraph cluster_2 {
            label="LOVES ⇒"
            rule2_n1[shape=point;label="";xlabel="4";fontsize=10]
            rule2_n2_1[shape=circle;width=0.05;label=""]
            rule2_n2_2[shape=circle;width=0.05;label=""]
            rule2_he1[shape=box;label="loves"]
            rule2_he2[shape=box;label="MARY"]
            {rule2_n2_1} -> rule2_he1 -> {rule2_n1}
            {rule2_n1} -> rule2_he2 -> {rule2_n2_2}
        }
        subgraph cluster_3 {
            label="JOHN ⇒"
            rule3_n1[shape=point;label="";xlabel="3";fontsize=10]
            rule3_n2_1[shape=circle;width=0.05;label=""]
            rule3_n2_2[shape=circle;width=0.05;label=""]
            rule3_he1[shape=box;label="John"]
            rule3_he2[shape=box;label="LOVES"]
            {rule3_n2_1} -> rule3_he1 -> {rule3_n1}
            {rule3_n1} -> rule3_he2 -> {rule3_n2_2}
        }
        subgraph cluster_4 {
            label="S ⇒"
            rule4_n1_1[shape=circle;width=0.05;label=""]
            rule4_n1_2[shape=circle;width=0.05;label=""]
            rule4_he1[shape=box;label="JOHN"]
            {rule4_n1_1} -> rule4_he1 -> {rule4_n1_2}
        }
    }
    """
    _test_normalized_strings_are_equal(rules_dot,expected)
    _print_dot(rules_dot)
end


text:
struct CustomPotentialParameters <: PotentialParameters
    a::AbstractFloat
end

import NBodySimulator.get_accelerating_function
function get_accelerating_function(p::CustomPotentialParameters, simulation::NBodySimulation)
    (dv, u, v, t, i) -> begin
        custom_accel = SVector(p.a, 0.0, 0.0); dv .= custom_accel
    end
end

@testset "A sysem with custom potential" begin
    m1 = 1.0
    m2 = 1.0

    t1 = 0.0
    t2 = 1.0
    τ = (t2 - t1) / 100

    p1 = MassBody(SVector(1, 0.0, 0.0), SVector(0.0, 0.0, 0.0), m1)
    p2 = MassBody(SVector(-1, 0.0, 0.0), SVector(0.0, 0.0, 0.0), m2)

    parameters = CustomPotentialParameters(1.5)
    system = PotentialNBodySystem([p1, p2], Dict(:custom_potential_params => parameters))
    simulation = NBodySimulation(system, (t1, t2))
    simResult = run_simulation(simulation, VelocityVerlet(), dt=τ)

    v2 = get_velocity(simResult, t2, 1)
    r2 = get_position(simResult, t2, 1)

    ε = 1e-6
    @test 1.5 ≈ v2[1] atol = ε
    @test 1.75 ≈ r2[1] atol = ε

    potential_system = PotentialNBodySystem([p1, p2]; potentials=[:lennard_jones, :electrostatic, :gravitational, :magnetostatic])
    @test sprint(io -> show(io, potential_system)) ==
"Potentials: \nLennard-Jones:\n\tϵ:1.0\n\tσ:1.0\n\tR:2.5\nElectrostatic:\n\tk:9.0e9\nMagnetostatic:\n\tμ/4π:1.0e-7\nGravitational:\n\tG:6.67408e-11\n"

end


text:
# ---
# title: 1437. Check If All 1's Are at Least Length K Places Away
# id: problem1437
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/check-if-all-1s-are-at-least-length-k-places-away/description/>
# hidden: true
# ---
#
# Given an array `nums` of 0s and 1s and an integer `k`, return `True` if all
# 1's are at least `k` places away from each other, otherwise return `False`.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/04/15/sample_1_1791.png)**
#
#
#
#     Input: nums = [1,0,0,0,1,0,0,1], k = 2
#     Output: true
#     Explanation: Each of the 1s are at least 2 places away from each other.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/04/15/sample_2_1791.png)**
#
#
#
#     Input: nums = [1,0,0,1,0,1], k = 2
#     Output: false
#     Explanation: The second 1 and third 1 are only one apart from each other.
#
# **Example 3:**
#
#
#
#     Input: nums = [1,1,1,1,1], k = 0
#     Output: true
#
#
# **Example 4:**
#
#
#
#     Input: nums = [0,1,0,1], k = 1
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10^5`
#   * `0 <= k <= nums.length`
#   * `nums[i]` is `0` or `1`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Base.Threads: Atomic, SpinLock, atomic_add!

# ListSegment represents a segment from a GrowingList; it has a list `data` to hold the elements, filled up to `currentlength`, and possibly a reference to the next segment, if it is not the final segment.
mutable struct ListSegment{T}
    data::Vector{T}
    currentlength::Int
    next::Base.RefValue{ListSegment{T}}
end
ListSegment{T}(data::Vector{T}, currentlength::Int) where T =
    ListSegment{T}(data, currentlength, Ref{ListSegment{T}}())

# getindex, assumes that index is inbounds, traverses the linked list
function _unsafe_getindex(l::ListSegment, i::Int)
    if i <= l.currentlength
        getindex(l.data, i)
    else
        _unsafe_getindex(l.next[], i - l.currentlength)
    end
end

# get or push a new element at the end; in itself not thread safe, should be protected by the lock in the parent GrowingList
function _unsafe_get!(l::ListSegment{T}, n::Int, default, newlength) where T
    N = length(l.data)
    if n > N
        if isassigned(l.next)
            return _unsafe_get!(l.next[], n - N, default, newlength)
        else
            newsegment = Vector{T}(undef, newlength)
            newsegment[1] = default
            l.next = Ref(ListSegment{T}(newsegment, 1))
            return default
        end
    else
        if n <= l.currentlength
            @inbounds return getindex(l.data, n)
        else
            @assert n == l.currentlength+1
            l.data[n] = default
            l.currentlength += 1
            return default
        end
    end
end

"""
    GrowingList{T} <: AbstractVector{T}

    GrowingList([iter,]; sizehint = max(16, length(iter)), growthfactor = 2.)

A thread safe vector / list data structure where new elements can be added at the back.
Once an element is set, it cannot be changed or removed. This ensures thread safe
`getindex` of that element without requiring a lock. The `length` of a `GrowingList`
instance can also be probed without a lock, but the return value will be a lower bound,
i.e. the list can already have increased in length at the same time.

New elements can be added using the syntax
`get!(l::GrowingList, i::Int, value)`
`get!(value_generator::Callable, l::GrowingList, i::Int)`

where the new element `value` or `value_generator()` will only be added if `i` is
`length(l)+1`. If multiple tasks or threads try to `get!` the same index `i`, only one of
them will actually be adding that element. The `value` or `value_generator()` produced by
the different threads should be the same to avoid unpredictable results.

The list is grown by adding new segments using a linked list data structure. This guarantees that existing data does never have to move in memory, which is required in order to make `getindex` threadsafe without lock.
"""
mutable struct GrowingList{T} <: AbstractVector{T}
    first::ListSegment{T}
    totallength::Atomic{Int}
    growthfactor::Float64
    lock::SpinLock
    function GrowingList{T}(iter;
                            sizehint = max(16, length(iter)),
                            growthfactor = 2.) where {T}
        firstsegment = Vector{T}(undef, sizehint)
        i = 0
        next = iterate(iter)
        @inbounds while i < sizehint && next !== nothing
            i += 1
            val, state = next
            firstsegment[i] = val
            next = iterate(iter, state)
        end
        first = ListSegment{T}(firstsegment, i)
        while next !== nothing
            i += 1
            val, state = next
            _unsafe_getindex(first, i, val, ceil(Int, (i-1)*growthfactor))
            next = iterate(iter, state)
        end
        return new{T}(first, Atomic{Int}(i), growthfactor, SpinLock())
    end
end
GrowingList(v::Vector{T}; sizehint = max(16, length(v)), growthfactor = 2.) where {T} =
    GrowingList{T}(v; sizehint = sizehint, growthfactor = growthfactor)

GrowingList{T}(; sizehint = 16, growthfactor = 2.) where {T} =
    GrowingList{T}((); sizehint = sizehint, growthfactor = growthfactor)

GrowingList(; sizehint = 16, growthfactor = 2.) =
    GrowingList{Any}((); sizehint = sizehint, growthfactor = growthfactor)

Base.length(l::GrowingList) = l.totallength[]
function _raise_length!(l::GrowingList)
    atomic_add!(l.totallength, 1)
end

Base.size(l::GrowingList) = (length(l),)

@inline function Base.getindex(l::GrowingList, n::Int)
    @boundscheck checkbounds(l, n)
    return _unsafe_getindex(l.first, n)
end

function Base.get!(l::GrowingList, n::Int, default)
    if n <= length(l)
        return _unsafe_getindex(l.first, n)
    else
        ll = l.lock
        lock(ll)
        len = length(l)
        nextlen = len + 1
        if n <= len # try again, maybe already ok now
            unlock(ll)
            return _unsafe_getindex(l.first, n)
        elseif n == nextlen
            _unsafe_get!(l.first, n, default, ceil(Int, (l.growthfactor-1)*len))
            _raise_length!(l)
            unlock(ll)
            return default
        else
            unlock(ll)
            _inserterror(nextlen)
        end
    end
end

@noinline _inserterror(len::Int) =
    throw(ArgumentError("can only insert new element at next index: " * string(len)))

function Base.get!(default::Base.Callable, l::GrowingList, n::Int)
    if n <= l.totallength[]
        return _unsafe_getindex(l.first, n)
    else
        v = default()
        ll = l.lock
        lock(ll)
        len = length(l)
        nextlen = len + 1
        if n <= len # try again, maybe already ok now
            unlock(ll)
            return _unsafe_getindex(l.first, n)
        elseif n == nextlen
            _unsafe_get!(l.first, n, v, ceil(Int, (l.growthfactor-1)*len))
            _raise_length!(l)
            unlock(ll)
            return v
        else
            unlock(ll)
            _inserterror(nextlen)
        end
    end
end


text: println( lu([1 3 5 ; 2 4 7 ; 1 1 0]) )

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
module BasicFWI
using jInv.Mesh
using jInv.Utils
using Distributed
using SparseArrays
using LinearAlgebra
using DelimitedFiles
using Statistics

export BasicFWIparam
import jInv.ForwardShare.ForwardProbType
"""
	BasicFWIparam <: ForwardProbType

	type defining Full Waveform Inversion (FWI) problem.

	Constructor: getBasicFWIparam(omega,gamma,Q,P,Mesh)

	Fields:
		omega     - frequencies
		gamma     - attenuation
		Sources   - matrix of sources
		Receivers - receiver matrix
		Mesh      - forward mesh
		Fields    -
		Ainv      - for linear solver, LU factorizations
"""
mutable struct BasicFWIparam <: ForwardProbType
    omega     # frequencies
    gamma     # attenuation
    Sources   # Sources
    Receivers # Receivers
    Mesh      # Mesh
	Fields    # Fields
	# Helmholtz # Helmholtz operators
	Ainv      # LU factorization
	originalSources	  # sources before trace estimation
	ExtendedSources
	# HinvtPs
end

export getBasicFWIparam
"""
	pFor = getBasicFWIparam(omega,gamma,Q,P,Mesh,doDistribute=false)

	constructs BasicFWI param

"""
function getBasicFWIparam(omega,gamma,Q,P,Mesh,doDistribute=false)
    nfreq = length(omega)

    if doDistribute && (nfreq > 1)
        pFor     = Array{RemoteChannel}(undef, nfreq)
        probsMax = ceil(Integer,nfreq/nworkers())
        nprobs   = zeros(maximum(workers()))

        i=1; nextidx() = (idx = i; i+=1; idx)

        @sync begin
            for p=workers()
                @async begin
                    while true
                        if nprobs[p]>=probsMax
                            break
                        end
                        idx = nextidx()
                        if (idx>nfreq)
                            break
                        end
                        pFor[idx] = initRemoteChannel(getBasicFWIparam,p,omega[idx],gamma,Q,P,Mesh)
                        nprobs[p] +=1
                        wait(pFor[idx])
                    end
                end
            end
        end
    else
        pFor = BasicFWIparam(omega,gamma,Q,P,Mesh,[],[],copy(Q),copy(Q))
    end
    return pFor
end

include("getData.jl")
include("getSensMatVec.jl")
include("getSensTMatVec.jl")
include("getMassMatrix.jl")

include("getHelmholtzOperator.jl")
include("freqCont.jl")
include("solvers.jl")

end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#

function build_tree!(scene::Scene, world::Object3D)::Nothing
    tree                = scene.tree
    stack               = scene.stack
    allCollisionElements = scene.allCollisionElements
    empty!(tree)
    empty!(stack)
    empty!(scene.allVisuElements)
    empty!(allCollisionElements)

    fillVisuElements!(scene, world, world)
    append!(stack, world.children)
    while length(stack) > 0
        obj = pop!(stack)
        fillVisuElements!(scene, obj, world)
        if canCollide(obj)
            push!(allCollisionElements, obj)
        end
        push!(tree, obj)
        # Visit children of frame
        append!(stack, obj.children)
    end
    if length(allCollisionElements) > 1
        scene.collide = true
    else
        scene.collide = false
    end
    return nothing
end

insert_and_dedup!(v::Vector, x) = (splice!(v, searchsorted(v,x), x); v)

function addIndicesOfCutJointsToSuperObj(scene::Scene)
    tmp = collect(values(scene.noCPairsHelp))
    for i=1:length(tmp)
        if length(tmp[i]) == 2
            insert_and_dedup!(scene.noCPairs[minimum(tmp[i])], maximum(tmp[i]))
        else
            error("...from addIndicesOfCutJointsToSuperObj: problems with amount of cut joints")
        end
    end
    return nothing
end

function createAABB_noCPairs(scene::Scene, superObjsRow::SuperObjsRow)
    if length(superObjsRow.superObjCollision.superObj) > 0 && !isempty(superObjsRow.noCPair)
        push!(scene.noCPairs, superObjsRow.noCPair)
    else
        push!(scene.noCPairs, [0])
    end

    if length(superObjsRow.superObjCollision.superObj) > 0
        AABBrow = [Basics.BoundingBox() for i = 1:length(superObjsRow.superObjCollision.superObj)]
        push!(scene.AABB, AABBrow)
    else
        push!(scene.AABB, [])
    end
    return nothing
end

function changeParentToRootObj(newParent::Object3D, obj::Object3D)
    if !isRootObject(newParent)
        error("from changeParentToRootObj: new parent is not a root object!")
    end

    child_r_rel  = obj.r_rel
    child_R_rel  = obj.R_rel
    parent_r_rel = obj.parent.r_rel
    parent_R_rel = obj.parent.R_rel

    obj.r_rel = parent_r_rel + parent_R_rel' * child_r_rel
    obj.R_rel = child_R_rel * parent_R_rel

    # Basics.deleteItem(obj.parent.children, obj) # klappt an dieser Stelle nicht
    # Reverse obj, so that newParent is the new parent
    obj.parent = newParent
    push!(newParent.children, obj)
    return nothing
end

function changeParent(newParent::Object3D, obj::Object3D)
    if isRootObject(newParent)
        child_r_abs  = obj.r_abs
        child_R_abs  = obj.R_abs
        parent_r_abs_new = newParent.r_abs
        parent_R_abs_new = newParent.R_abs

        obj.r_rel = parent_R_abs_new * (child_r_abs -parent_r_abs_new)
        obj.R_rel = child_R_abs * parent_R_abs_new'

        Basics.deleteItem(obj.parent.children, obj)
        obj.parent = newParent
        push!(newParent.children, obj)
    elseif isRootObject(newParent.parent)
        changeParent(newParent.parent, obj)
    else
        error("from changeParent: this should not happen")
    end
end

function getRootObj(obj::Object3D)
    if isRootObject(obj)
        return obj
    else
        getRootObj(obj.parent)
end; end

function getWorldObj(obj::Object3D)
    root = getRootObj(obj)
    if isWorld(root)
        return root
    else
        getWorldObj(root.parent)
    end
end

# the indices of super objects, which can't collide, are stored in a list
function fillStackOrBuffer!(scene::Scene, superObj::SuperObjsRow, obj::Object3D, rootSuperObj::Object3D)::Nothing
    n_children =  length(obj.children)
    help = fill(false, n_children)

    for i = 1:n_children
        child = obj.children[i]
        if isNotWorld(child)
            if isNotFixed(child)
                push!(scene.buffer, child)
                child.isRootObj = true
                child.computeAcceleration = true
                if !child.canCollide    # !isFree(child) &&  !( typeof(child.joint) <: Modia3D.AbstractPrismatic )
                    push!(superObj.noCPair, length(scene.buffer))
                end
            else
            if isMovable(child)
                push!(scene.buffer, child)
                child.isRootObj = true
                child.computeAcceleration = true
            else
                push!(scene.stack, child)
                assignAccVelo(scene.treeAccVelo, child)
                if !(obj == rootSuperObj)
                    changeParentToRootObj(rootSuperObj, child)
                    help[i] = true
    end; end; end; end; end

    if !isempty(obj.children)
        deleteat!(obj.children,help)
    end
    return nothing
end

# it builds the elements which may collide
# to reduce the amount of collision pairs, some assumptions are made:
#   elements which are rigidly attached, can't collide
#     these elements form together a super object
#   elements which are directly connected with a joint can't collide
#     these elements are excluded from the collision list
function build_superObjs!(scene::Scene, world::Object3D)::Nothing
    if !scene.initSuperObj
    stack = scene.stack
    buffer = scene.buffer
    treeAccVelo  = scene.treeAccVelo
    empty!(stack)
    empty!(buffer)
    empty!(scene.allVisuElements)
    empty!(treeAccVelo)

    world.computeAcceleration = true
    world.isRootObj = true # all objs stored in buffer are root objs
    push!(buffer, world)
    actPos = 1
    nPos   = 1

    hasOneCollisionSuperObj  = false
    hasMoreCollisionSuperObj = false

    while actPos <= nPos
        superObjsRow = SuperObjsRow()
        AABBrow      = Vector{Basics.BoundingBox}[]
        rootSuperObj = buffer[actPos]

        fillVisuElements!(scene, rootSuperObj, world)
        if rootSuperObj != world
            assignAll(scene, superObjsRow, rootSuperObj, world, actPos)
            push!(treeAccVelo, rootSuperObj)
        end
        fillStackOrBuffer!(scene, superObjsRow, rootSuperObj, rootSuperObj)

        if isMovable(rootSuperObj)
            push!(scene.allowedToMove, true)
        else
            push!(scene.allowedToMove, nothing)
        end

        while length(stack) > 0
            frameChild = pop!(stack)
            fillVisuElements!(scene, frameChild, world)
            assignAll(scene, superObjsRow, frameChild, world, actPos)
            fillStackOrBuffer!(scene, superObjsRow, frameChild, rootSuperObj)
        end

        if length(superObjsRow.superObjCollision.superObj) > 0 && hasOneCollisionSuperObj == true
            hasMoreCollisionSuperObj = true
        elseif length(superObjsRow.superObjCollision.superObj) > 0
            hasOneCollisionSuperObj = true
        end

        createAABB_noCPairs(scene, superObjsRow)
        push!(scene.superObjs, superObjsRow)
        nPos = length(buffer)
        actPos += 1
    end
    addIndicesOfCutJointsToSuperObj(scene)
    @assert(length(scene.noCPairs) == length(scene.superObjs))

    #=
    println("first elements of tree")
    println("[")
    for a in scene.treeAccVelo
        println("parent = ", Modia3D.fullName(a.parent))
        println(Modia3D.fullName(a))
        println("[")
        for b in a.children
        println("children = ", Modia3D.fullName(b))
        end
        println("]")
    end
    println("]")
    =#
    #=
    println("superObjRow.superObjCollision.superObj")
    for superObjRow in scene.superObjs
    println("[")
    for a in superObjRow.superObjCollision.superObj
        println(Modia3D.fullName(a))
    end
    println("]")
    println(" ")
    end
    =#

    #=
    println("allowedToMove")
    for a in scene.allowedToMove
        show(a)
        println(" ")
    end
    =#

    #=
    println("superObjRow.superObjMovable.superObj")
    for superObjRow in scene.superObjs
    println("[")
    for a in superObjRow.superObjMovable.superObj
        println(Modia3D.fullName(a))
    end
    println("]")
    println(" ")
    end
    =#
    #=
    println("superObjRow.superObjMass.superObj")
    for superObjRow in scene.superObjs
    println("[")
    for a in superObjRow.superObjMass.superObj
        println(Modia3D.fullName(a))
    end
    println("]")
    println(" ")
    end
    =#

    if hasMoreCollisionSuperObj
        scene.collide = true
    else
        scene.collide = false
    end
    scene.initSuperObj = true
    end
    return nothing
end

function visualizeWorld!(world::Object3D; sceneOptions = SceneOptions())::Nothing
    scene = Scene(sceneOptions)
    scene.analysis = Modia3D.KinematicAnalysis
    initAnalysis!(world, scene)
    updatePosition!(world)
    visualize!(Modia3D.renderer[1], 0.0)
    visualize!(Modia3D.renderer[1], 1.0)
    closeAnalysis!(scene)
    return nothing
end

function makeTreeAvailable(scene::Scene)
    if scene.options.useOptimizedStructure
        scene.treeForComputation = scene.treeAccVelo
    else
        for obj in scene.tree
            if featureHasMass(obj.feature)
                # Copy feature mass properties to object mass properties
                obj.m    = obj.feature.massProperties.m
                obj.r_CM = obj.feature.massProperties.rCM
                obj.I_CM = obj.feature.massProperties.I
            end
        end
        scene.treeForComputation = scene.tree
    end
end

function makeJointsAvailable(scene::Scene)
    tree = scene.treeForComputation
    empty!(scene.revolute)
    empty!(scene.prismatic)
    empty!(scene.freeMotion)

    for obj in tree
        jointKind = obj.jointKind

        if jointKind == FixKind
            if obj.R_rel === Modia3D.NullRotation
                obj.jointKind = FixTranslationKind
            end

        elseif jointKind == FixTranslationKind
            continue

        elseif jointKind == RevoluteKind
            push!(scene.revolute, obj.joint)
            obj.jointIndex = length(scene.revolute)

        elseif jointKind == PrismaticKind
            push!(scene.prismatic, obj.joint)
            obj.jointIndex = length(scene.prismatic)

        elseif jointKind == FreeMotionKind
            push!(scene.freeMotion, obj.joint)
            obj.jointIndex = length(scene.freeMotion)

            if hasNoParent(obj.parent)
                obj.jointKind = AbsoluteFreeMotionKind
            end

        elseif jointKind == AbsoluteFreeMotionKind
            push!(scene.freeMotion, obj.joint)
            obj.jointIndex = length(scene.freeMotion)

        else
            error("Bug in Modia3D/src/Composition/handler.jl: jointKind = $jointKind is not known:\njoint = ", obj.joint)
        end
    end
end

function chooseAndBuildUpTree(world::Object3D, scene::Scene)
    # Build tree for optimized structure or standard structure
    # collision handling is only available for optimized structure
    if scene.options.useOptimizedStructure
        build_superObjs!(scene, world)
        if scene.options.enableContactDetection && scene.collide
            initializeContactDetection!(world, scene)
            if length(world.contactVisuObj1) > 0
                append!(scene.allVisuElements, world.contactVisuObj1)
                append!(scene.allVisuElements, world.contactVisuObj2)
            end
            if length(world.supportVisuObj1A) > 0
                append!(scene.allVisuElements, world.supportVisuObj1A)
                append!(scene.allVisuElements, world.supportVisuObj2A)
                append!(scene.allVisuElements, world.supportVisuObj3A)
                append!(scene.allVisuElements, world.supportVisuObj1B)
                append!(scene.allVisuElements, world.supportVisuObj2B)
                append!(scene.allVisuElements, world.supportVisuObj3B)
            end
            if length(world.AABBVisu) > 0
                append!(scene.allVisuElements, world.AABBVisu)
            end
        end
        initializeMassComputation!(scene)
    else
        build_tree!(scene, world)
        if scene.options.enableContactDetection
            error("Collision handling is only possible with the optimized structure. Please set useOptimizedStructure = true in Modia3D.SceneOptions.")
        end
    end
    if length(scene.allVisuElements) > 0
        scene.visualize = scene.options.enableVisualization
    else
        scene.visualize = false
        scene.exportAnimation = false
    end
    makeTreeAvailable(scene)
    makeJointsAvailable(scene)
    scene.initAnalysis = true
    return nothing
end

function initAnalysis!(world::Object3D, scene::Scene)
    # Initialize spanning tree and visualization (if visualization desired and visual elements present)
    chooseAndBuildUpTree(world, scene)
    if scene.visualize
        initializeVisualization(Modia3D.renderer[1], scene.allVisuElements)
    end
    return nothing
end

function errorMessageCollision(functionName::String)
    error("\nError message from ", functionName, ":\n",
            "There are too less solids which can collide.")
end

function visualize!(scene::Scene, time::Float64)
    if scene.visualize
        visualize!(Modia3D.renderer[1], time)
    end
end

function emptyScene!(scene)
    empty!(scene.stack)
    empty!(scene.buffer)

    empty!(scene.superObjs)
    empty!(scene.updateVisuElements)
    empty!(scene.treeAccVelo)
    empty!(scene.tree)
    empty!(scene.allVisuElements)
    empty!(scene.allCollisionElements)
    empty!(scene.noCPairs)
    empty!(scene.noCPairsHelp)
    empty!(scene.allowedToMove)
    empty!(scene.AABB)

    scene.collide = false
    scene.initSuperObj = false

    # Close Analysis
    scene.initAnalysis = false
end

function closeAnalysis!(scene::Scene)
    # Close Visualisation
    closeVisualization(Modia3D.renderer[1])
    empty!(scene.allVisuElements)
    scene.visualize = false
    # Close Collision detection
    if scene.collide
        closeContactDetection!(scene.options.contactDetection)
    end
    emptyScene!(scene)
end


text:
PREAMBLE = BitArray([1,1,0,1,0,0,0,1])

function make_subframe_ninv(id_data)
    return BitArray(vcat(falses(240),falses(8), reverse(id_data), falses(19), falses(22), PREAMBLE))
end
function make_subframe_inv(id_data)
    return BitArray(vcat(falses(240),falses(8), reverse(id_data), falses(19), falses(22), .!PREAMBLE))
end

test_buffer_ninv = Vector{BitArray}(undef, 0)
push!(test_buffer_ninv, make_subframe_ninv([0,0,1]))
push!(test_buffer_ninv, make_subframe_ninv([0,1,0]))
push!(test_buffer_ninv, make_subframe_ninv([0,1,1]))
push!(test_buffer_ninv, make_subframe_ninv([1,0,0]))
push!(test_buffer_ninv, make_subframe_ninv([1,0,1]))

test_buffer_inv = Vector{BitArray}(undef, 0)
push!(test_buffer_inv, make_subframe_inv([0,0,1]))
push!(test_buffer_inv, make_subframe_inv([0,1,0]))
push!(test_buffer_inv, make_subframe_inv([0,1,1]))
push!(test_buffer_inv, make_subframe_inv([1,0,0]))
push!(test_buffer_inv, make_subframe_inv([1,0,1]))

false_frame = zeros(300)

@testset "Test Not Inverted Preamble" begin
    buffer = vcat(test_buffer_ninv[5], test_buffer_ninv[4], test_buffer_ninv[3], test_buffer_ninv[2], test_buffer_ninv[1])
    out = GNSSDecoder.find_preamble(buffer)
    @test out == true

    buffer = vcat(false_frame, test_buffer_ninv[4], test_buffer_ninv[3], test_buffer_ninv[2], test_buffer_ninv[1])
    out = GNSSDecoder.find_preamble(buffer)
    @test out == false

    buffer = vcat(test_buffer_ninv[1], test_buffer_ninv[2], test_buffer_ninv[3], test_buffer_ninv[4], test_buffer_ninv[5])
    out = GNSSDecoder.find_preamble(buffer)
    @test out == false

end

@testset "Test Inverted Preamble" begin
    buffer = vcat(test_buffer_inv[5], test_buffer_inv[4], test_buffer_inv[3], test_buffer_inv[2], test_buffer_inv[1])
    out = GNSSDecoder.find_preamble(buffer)
    @test out == true

    buffer = vcat(false_frame, test_buffer_inv[4], test_buffer_inv[3], test_buffer_inv[2], test_buffer_inv[1])
    out = GNSSDecoder.find_preamble(buffer)
    @test out == false

    buffer = vcat(test_buffer_inv[1], test_buffer_inv[2], test_buffer_inv[3], test_buffer_inv[4], test_buffer_inv[5])
    out = GNSSDecoder.find_preamble(buffer)
    @test out == false
end

@testset "Test Mixed Preamble" begin
    buffer = vcat(test_buffer_inv[5], test_buffer_ninv[4], test_buffer_inv[3], test_buffer_ninv[2], test_buffer_inv[1])
    out = GNSSDecoder.find_preamble(buffer)
    @test out == true

    buffer = vcat(false_frame, test_buffer_ninv[4], test_buffer_ninv[3], test_buffer_inv[2], test_buffer_inv[1])
    out = GNSSDecoder.find_preamble(buffer)
    @test out == false

    buffer = vcat(test_buffer_ninv[1], test_buffer_inv[2], test_buffer_inv[3], test_buffer_inv[4], test_buffer_ninv[5])
    out = GNSSDecoder.find_preamble(buffer)
    @test out == false
end


text:
ClassicalOrthogonalPolynomials.checkpoints(d::UnitDisk{T}) where T = [SVector{2,T}(0.1,0.2), SVector{2,T}(0.2,0.3)]

"""
    ZernikeWeight(a, b)

is a quasi-vector representing `r^(2a) * (1-r^2)^b`
"""
struct ZernikeWeight{T} <: Weight{T}
    a::T
    b::T
end

"""
    ZernikeWeight(b)

is a quasi-vector representing `(1-r^2)^b`
"""

ZernikeWeight(b) = ZernikeWeight(zero(b), b)
ZernikeWeight{T}(b) where T = ZernikeWeight{T}(zero(T), b)
ZernikeWeight{T}() where T = ZernikeWeight{T}(zero(T))
ZernikeWeight() = ZernikeWeight{Float64}()

copy(w::ZernikeWeight) = w

axes(::ZernikeWeight{T}) where T = (Inclusion(UnitDisk{T}()),)

==(w::ZernikeWeight, v::ZernikeWeight) = w.a == v.a && w.b == v.b

function getindex(w::ZernikeWeight, xy::StaticVector{2})
    r = norm(xy)
    r^(2w.a) * (1-r^2)^w.b
end

"""
    Zernike(a, b)

is a quasi-matrix orthogonal `r^(2a) * (1-r^2)^b`
"""
struct Zernike{T} <: BivariateOrthogonalPolynomial{T}
    a::T
    b::T
    Zernike{T}(a::T, b::T) where T = new{T}(a, b)
end
Zernike{T}(a, b) where T = Zernike{T}(convert(T,a), convert(T,b))
Zernike(a::T, b::V) where {T,V} = Zernike{float(promote_type(T,V))}(a, b)
Zernike{T}(b) where T = Zernike{T}(zero(b), b)
Zernike{T}() where T = Zernike{T}(zero(T))

"""
    Zernike(b)

is a quasi-matrix orthogonal `(1-r^2)^b`
"""
Zernike(b) = Zernike(zero(b), b)
Zernike() = Zernike{Float64}()

axes(P::Zernike{T}) where T = (Inclusion(UnitDisk{T}()),blockedrange(oneto(∞)))

==(w::Zernike, v::Zernike) = w.a == v.a && w.b == v.b

copy(A::Zernike) = A

Base.summary(io::IO, P::Zernike) = print(io, "Zernike($(P.a), $(P.b))")

orthogonalityweight(Z::Zernike) = ZernikeWeight(Z.a, Z.b)

zerniker(ℓ, m, a, b, r::T) where T = sqrt(convert(T,2)^(m+a+b+2-iszero(m))/π) * r^m * normalizedjacobip((ℓ-m) ÷ 2, b, m+a, 2r^2-1)
zerniker(ℓ, m, b, r) = zerniker(ℓ, m, zero(b), b, r)
zerniker(ℓ, m, r) = zerniker(ℓ, m, zero(r), r)

function zernikez(ℓ, ms, a, b, rθ::RadialCoordinate{T}) where T
    r,θ = rθ.r,rθ.θ
    m = abs(ms)
    zerniker(ℓ, m, a, b, r) * (signbit(ms) ? sin(m*θ) : cos(m*θ))
end

zernikez(ℓ, ms, a, b, xy::StaticVector{2}) = zernikez(ℓ, ms, a, b, RadialCoordinate(xy))
zernikez(ℓ, ms, b, xy::StaticVector{2}) = zernikez(ℓ, ms, zero(b), b, xy)
zernikez(ℓ, ms, xy::StaticVector{2,T}) where T = zernikez(ℓ, ms, zero(T), xy)

function getindex(Z::Zernike{T}, rθ::RadialCoordinate, B::BlockIndex{1}) where T
    ℓ = Int(block(B))-1
    k = blockindex(B)
    m = iseven(ℓ) ? k-isodd(k) : k-iseven(k)
    zernikez(ℓ, (isodd(k+ℓ) ? 1 : -1) * m, Z.a, Z.b, rθ)
end

getindex(Z::Zernike, xy::StaticVector{2}, B::BlockIndex{1}) = Z[RadialCoordinate(xy), B]
getindex(Z::Zernike, xy::StaticVector{2}, B::Block{1}) = [Z[xy, B[j]] for j=1:Int(B)]
getindex(Z::Zernike, xy::StaticVector{2}, JR::BlockOneTo) = mortar([Z[xy,Block(J)] for J = 1:Int(JR[end])])

###
# Jacobi matrices
###
function jacobimatrix(::Val{1}, Z::Zernike{T}) where T
    if iszero(Z.a)
        α = Z.b     # extract second basis parameter

    k = mortar(Base.OneTo.(oneto(∞)))     # k counts the the angular mode (+1)
    n = mortar(Fill.(oneto(∞),oneto(∞)))  # n counts the block number which corresponds to the order (+1)

    # repeatedly used for sorting
    keven = iseven.(k)
    kodd = isodd.(k)
    neven = iseven.(n)
    nodd = isodd.(n)

    # h1-h5 are helpers for our different sorting scheme

    ## Compute super diagonal of super diagonal blocks.
    dufirst = neven .* (k .== 2) .* n .* (n .+ 2*α) ./ 2

    ## Compute even-block entries in super diagonal of super diagonal blocks
    h1 = n .- (n .+ 2 .- k .- keven .* (n .> 2)) .÷ 2
    dueven = ((nodd .* k) .>= 2) .* h1 .* (h1 .+ α)

    ## Compute odd-block entries in super diagonal of super diagonal blocks
    h2 = (n .- 2 .+ k .+ kodd .* (n .> 2)) .÷ 2
    duodd = (((neven .* k) .>= 2)  .- (neven .* (k .== 2))) .* h2 .* (h2 .+ α)

    ## Compute even-block sub diagonal elements of super diagonal blocks
    h3 = n .- (k .+ 1 .+ kodd .+ n) .÷ 2
    dleven = ((k .<= (n .- 2)) .- ((k .< 2) .* nodd)) .* neven .* h3 .* (h3 .+ α)

    ## Compute odd-block sub diagonal elements of super diagonal blocks
    h4 = n .- (k .+ 1 .+ keven .+ n) .÷ 2
    dlodd = (k .> 1) .* (n .> 3) .* nodd .* h4 .* (h4 .+ α)

    ## Compute and add in special case odd-block sub diagonal elements of super diagonal blocks
    h5 = n .- 2 .+ k .+ keven
    dlspecial = nodd .* (k .== 1) .* h5 .* (h5 .+ 2*α) ./ 2

    # finalize bands with explicit formula
    quotient = 4 .* (n .+ (α-1)) .* (n .+ α)
    du = sqrt.( (dufirst .+ dueven .+ duodd ) ./ quotient)
    dl = sqrt.( (dleven .+ dlodd .+ dlspecial)  ./ quotient)

    return Symmetric(BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, du, Zeros((axes(n,1),)), dl)', axes(n,1), (-1,1), (1,1)))
    else
        error("Implement for non-zero first basis parameter.")
    end
end

function jacobimatrix(::Val{2}, Z::Zernike{T}) where T
    if iszero(Z.a)
        α = Z.b     # extract second basis parameter

        k = mortar(Base.OneTo.(oneto(∞)))     # k counts the the angular mode (+1)
        n = mortar(Fill.(oneto(∞),oneto(∞)))  # n counts the block number which corresponds to the order (+1)

        # repeatedly used for sorting
        keven = iseven.(k)
        kodd = isodd.(k)
        neven = iseven.(n)
        nodd = isodd.(n)

        # h1-h4 are helpers for our different sorting scheme

        # first entries for all blocks
        h1 = (n .- nodd)
        l1 = (k .== 1) .* (h1 .* (h1 .+ 2*α) ./ 2)

        # Even blocks
        h0 = (kodd .* ((k .÷ 2) .+ 1) .- (keven .* ((k .÷ 2) .- 1)))
        h2 =  (k .>= 2) .* ((n .÷ 2 .- 1) .+ h0)
        l2 = neven .* (h2 .* (h2 .+ α))

        # Odd blocks
        h3 = (n .> k .>= 2) .* (((n .+ 1) .÷ 2) .- h0)
        l3 = nodd .* (h3 .* (h3 .+ α))
        # Combine for diagonal of super diagonal block
        d = sqrt.((l1 .+ l2 .+ l3) ./ (4 .* (n .+ (α-1)) .* (n .+ α)))

        # The off-diagonals of the super diagonal block are negative, shifted versions of the diagonal with some entries skipped
        dl = (-1) .* (nodd .* kodd .+ neven .* keven) .* Vcat(0 , d)
        du = (-1) .* (nodd .* keven .+ neven .* kodd) .* view(d,2:∞)

        # generate and return bands
        return Symmetric(BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, dl, Zeros((axes(n,1),)), d, Zeros((axes(n,1),)), du)', axes(n,1), (-1,1), (2,2)))
    else
        error("Implement for non-zero first basis parameter.")
    end
end

###
# Transforms
###

const FiniteZernike{T} = SubQuasiArray{T,2,Zernike{T},<:Tuple{<:Inclusion,<:BlockSlice{BlockRange1{OneTo{Int}}}}}

function grid(S::FiniteZernike{T}) where T
    N = blocksize(S,2) ÷ 2 + 1 # polynomial degree
    M = 4N-3

    r = sinpi.((N .-(0:N-1) .- one(T)/2) ./ (2N))

    # The angular grid:
    θ = (0:M-1)*convert(T,2)/M
    RadialCoordinate.(r, π*θ')
end

_angle(rθ::RadialCoordinate) = rθ.θ

function plotgrid(S::FiniteZernike{T}) where T
    N = blocksize(S,2) ÷ 2 + 1 # polynomial degree
    g = grid(parent(S)[:,Block.(OneTo(2N))]) # double sampling
    θ = [map(_angle,g[1,:]); 0]
    [permutedims(RadialCoordinate.(1,θ)); g g[:,1]; permutedims(RadialCoordinate.(0,θ))]
end

function plotgrid(S::SubQuasiArray{<:Any,2,<:Zernike})
    kr,jr = parentindices(S)
    Z = parent(S)
    plotgrid(Z[kr,Block.(OneTo(Int(findblock(axes(Z,2),maximum(jr)))))])
end

function plotvalues(u::ApplyQuasiVector{T,typeof(*),<:Tuple{Zernike, AbstractVector}}, x) where T
    Z,c = u.args
    CS = blockcolsupport(c)
    N = Int(last(CS)) ÷ 2 + 1 # polynomial degree
    F = ZernikeITransform{T}(2N, Z.a, Z.b)
    C = F * c[Block.(OneTo(2N))] # transform to grid
    [permutedims(u[x[1,:]]); # evaluate on edge of disk
     C C[:,1];
     fill(u[x[end,1]], 1, size(x,2))] # evaluate at origin and repeat
end

function plotvalues(u::ApplyQuasiVector{T,typeof(*),<:Tuple{Weighted{<:Any,<:Zernike}, AbstractVector}}, x) where T
    U = plotvalues(unweighted(u), x)
    w = weight(u.args[1])
    w[x] .* U
end

struct ZernikeTransform{T} <: Plan{T}
    N::Int
    disk2cxf::FastTransforms.FTPlan{T,2,FastTransforms.DISK}
    analysis::FastTransforms.FTPlan{T,2,FastTransforms.DISKANALYSIS}
end

struct ZernikeITransform{T} <: Plan{T}
    N::Int
    disk2cxf::FastTransforms.FTPlan{T,2,FastTransforms.DISK}
    synthesis::FastTransforms.FTPlan{T,2,FastTransforms.DISKSYNTHESIS}
end

function ZernikeTransform{T}(N::Int, a::Number, b::Number) where T<:Real
    Ñ = N ÷ 2 + 1
    ZernikeTransform{T}(N, plan_disk2cxf(T, Ñ, a, b), plan_disk_analysis(T, Ñ, 4Ñ-3))
end
function ZernikeITransform{T}(N::Int, a::Number, b::Number) where T<:Real
    Ñ = N ÷ 2 + 1
    ZernikeITransform{T}(N, plan_disk2cxf(T, Ñ, a, b), plan_disk_synthesis(T, Ñ, 4Ñ-3))
end

*(P::ZernikeTransform{T}, f::AbstractArray) where T = P * convert(Matrix{T}, f)
*(P::ZernikeTransform{T}, f::Matrix{T}) where T = ModalTrav(P.disk2cxf \ (P.analysis * f))
*(P::ZernikeITransform, f::AbstractVector) = P.synthesis * (P.disk2cxf * ModalTrav(f).matrix)

factorize(S::FiniteZernike{T}) where T = TransformFactorization(grid(S), ZernikeTransform{T}(blocksize(S,2), parent(S).a, parent(S).b))

@simplify function *(Δ::Laplacian, WZ::Weighted{<:Any,<:Zernike})
    @assert WZ.P.a == 0 && WZ.P.b == 1
    T = eltype(eltype(WZ))
    WZ.P * ModalInterlace{T}(broadcast(k ->  Diagonal(-cumsum(k:8:∞)), 4:4:∞), (ℵ₀,ℵ₀), (0,0))
end

@simplify function *(Δ::Laplacian, Z::Zernike)
    a,b = Z.a,Z.b
    @assert a == 0
    T = promote_type(eltype(eltype(Δ)),eltype(Z)) # TODO: remove extra eltype
    D = Derivative(Inclusion(ChebyshevInterval{T}()))
    Δs = BroadcastVector{AbstractMatrix{T}}((C,B,A) -> 4(HalfWeighted{:b}(C)\(D*HalfWeighted{:b}(B)))*(B\(D*A)), Normalized.(Jacobi.(b+2,a:∞)), Normalized.(Jacobi.(b+1,(a+1):∞)), Normalized.(Jacobi.(b,a:∞)))
    Δ = ModalInterlace(Δs, (ℵ₀,ℵ₀), (-2,2))
    Zernike(a,b+2) * Δ
end

###
# Fractional Laplacian
###

function *(L::AbsLaplacianPower, WZ::Weighted{<:Any,<:Zernike{<:Any}})
    @assert axes(L,1) == axes(WZ,1) && WZ.P.a == 0 && WZ.P.b == L.α
    WZ.P * Diagonal(WeightedZernikeFractionalLaplacianDiag{typeof(L.α)}(L.α))
end

# gives the entries for the (negative!) fractional Laplacian (-Δ)^(α) times (1-r^2)^α * Zernike(α)
struct WeightedZernikeFractionalLaplacianDiag{T} <: AbstractBlockVector{T}
    α::T
end

axes(::WeightedZernikeFractionalLaplacianDiag) = (blockedrange(oneto(∞)),)
copy(R::WeightedZernikeFractionalLaplacianDiag) = R

MemoryLayout(::Type{<:WeightedZernikeFractionalLaplacianDiag}) = LazyLayout()
Base.BroadcastStyle(::Type{<:Diagonal{<:Any,<:WeightedZernikeFractionalLaplacianDiag}}) = LazyArrayStyle{2}()

getindex(W::WeightedZernikeFractionalLaplacianDiag, k::Integer) = W[findblockindex(axes(W,1),k)]

function Base.view(W::WeightedZernikeFractionalLaplacianDiag{T}, K::Block{1}) where T
    l = Int(K)
    if isodd(l)
        m = Vcat(0,interlace(Array(2:2:l),Array(2:2:l)))
    else #if iseven(l)
        m = Vcat(interlace(Array(1:2:l),Array(1:2:l)))
    end
    return convert(AbstractVector{T}, 2^(2*W.α)*fractionalcfs2d.(l-1,m,W.α))
end

# generic d-dimensional ball fractional coefficients without the 2^(2*β) factor. m is assumed to be entered as abs(m)
function fractionalcfs(l::Integer, m::Integer, α::T, d::Integer) where T
    n = (l-m)÷2
    return exp(loggamma(α+n+1)+loggamma((2*α+2*n+d+2*m)/2)-loggamma(one(T)+n)-loggamma((2*one(T)*m+2*n+d)/2))
end
# 2 dimensional special case, again without the 2^(2*β) factor
fractionalcfs2d(l::Integer, m::Integer, β) = fractionalcfs(l,m,β,2)

function \(A::Zernike{T}, B::Zernike{V}) where {T,V}
    TV = promote_type(T,V)
    A.a == B.a && A.b == B.b && return Eye{TV}((axes(A,2),))
    st = Int(A.a - B.a + A.b - B.b)
    ModalInterlace{TV}((Normalized.(Jacobi{TV}.(A.b,A.a:∞)) .\ Normalized.(Jacobi{TV}.(B.b,B.a:∞))) .* convert(TV, 2)^(-st/2), (ℵ₀,ℵ₀), (0,2st))
end

function \(A::Zernike{T}, wB::Weighted{V,Zernike{V}}) where {T,V}
    TV = promote_type(T,V)
    B = wB.P
    A.a == B.a == A.b == B.b == 0 && return Eye{TV}((axes(A,2),))
    c = Int(B.a - A.a + B.b - A.b)
    @assert iszero(B.a)
    ModalInterlace{TV}((Normalized.(Jacobi{TV}.(A.b, A.a:∞)) .\ HalfWeighted{:a}.(Normalized.(Jacobi{TV}.(B.b, B.a:∞)))) .* convert(TV, 2)^(-c/2), (ℵ₀,ℵ₀), (2Int(B.b), 2Int(A.a+A.b)))
end


text:
module DimmWitted

immutable DW
  _dw::Ptr{Void}
  modelrepl::Cint
  datarepl::Cint
  accessmethod::Cint
  shareddatatype
end

_libpath = ""
_data_type = ""
_model_type = ""

_nogc = {0.0,0.0}

MR_PERCORE = 0
MR_PERNODE = 1
MR_PERMACHINE = 2
MR_SINGLETHREAD_DEBUG = 3

DR_FULL = 0
DR_SHARDING = 1

AC_ROW = 0
AC_COL = 1
AC_C2R = 2

function set_libpath(path)
	global _libpath
	_libpath = path
end

function get_libpath()
	global _libpath
	return _libpath
end

function hi()
	@eval ccall(($(string("Hello")), $(_libpath)), Void, ())
	return nothing
end

function open{DATATYPE, MODELTYPE}(examples::SparseMatrixCSC{DATATYPE,Int64}, model::Array{MODELTYPE,1}, _modelrepl, _datarepl, _acmethod, _shared_data=[1,2])

	modelrepl = convert(Cint, _modelrepl)
	datarepl = convert(Cint, _datarepl)
	acmethod = convert(Cint, _acmethod)

	global _libpath, _dw, _data_type, _model_type, _nogc

	_model = model

	nrows = examples.m
	ncols = examples.n
	colptr = examples.colptr
	rowptr = examples.rowval
	nnz = examples.nzval

	@eval immutable tmptype
		idx::Clonglong
		data::$(DATATYPE)
	end

	append!(_nogc, {_model, examples, _shared_data, tmptype})

	nmodelel = size(model, 1)

	_data_type = tmptype
	_model_type = MODELTYPE
	_shared_tupe = typeof(_shared_data)

	t = _shared_tupe
	n = length(_shared_data)

	_dw = @eval ccall( ($(string("SparseDimmWitted_Open2")), $(_libpath)), Ptr{Void}, (Any, Any, Any, Clonglong, Clonglong, Clonglong, Clonglong, Ptr{Void}, Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Void}, Cint, Cint, Cint, Any, Cint, Ptr{Void}), $(Array{DATATYPE}), $(Array{tmptype}), $(Array{MODELTYPE}), $(nrows), $(ncols), $(length(nnz)), $(nmodelel), $(nnz), $(rowptr), $(colptr) ,$(_model), $(modelrepl), $(datarepl), $(acmethod), $(t), $(n), $(_shared_data))

	dw = DW(_dw, modelrepl, datarepl, acmethod, typeof(_shared_data))

	println("[JULIA-DW] Created Sparse DimmWitted Object: ", dw._dw)

	return dw
end

function open{DATATYPE, MODELTYPE}(examples::Array{DATATYPE,2}, model::Array{MODELTYPE,1}, _modelrepl, _datarepl, _acmethod, _shared_data=[1,2])

	modelrepl = convert(Cint, _modelrepl)
	datarepl = convert(Cint, _datarepl)
	acmethod = convert(Cint, _acmethod)

	global _libpath, _dw, _data_type, _model_type, _nogc

	_examples = examples
	_model = model
	_examples_c = examples.'

	append!(_nogc, {_examples, _model, _examples_c, _shared_data})

	nrows = size(examples, 1)
	ncols = size(examples, 2)
	nmodelel = size(model, 1)

	_data_type = DATATYPE
	_model_type = MODELTYPE
	_shared_tupe = typeof(_shared_data)

	t = _shared_tupe
	n = length(_shared_data)

	_dw = @eval ccall( ($(string("DenseDimmWitted_Open2")), $(_libpath)), Ptr{Void}, (Any, Any, Clonglong, Clonglong, Clonglong, Ptr{Void}, Ptr{Void}, Cint, Cint, Cint, Any, Cint, Ptr{Void}), $(Array{DATATYPE}), $(Array{MODELTYPE}), $(nrows), $(ncols), $(nmodelel), $(_examples_c), $(_model), $(modelrepl), $(datarepl), $(acmethod), $(t), $(n), $(_shared_data))

	dw = DW(_dw, modelrepl, datarepl, acmethod, typeof(_shared_data))

	println("[JULIA-DW] Created Dense DimmWitted Object: ", dw._dw)

	return dw
end

function check_is_safe(func, ret, parameter)

	const stdout = STDOUT
	const rd, wr = redirect_stdout()
	code_llvm(func, parameter)
	str = readavailable(rd)
	close(rd)
	redirect_stdout(STDERR)

	if contains(str, "alloc") || contains(
			replace(
				replace(str, string("julia_",func), ""),
				"julia_type", "")
		, "julia_")
		println(str)
		return false
	else
		return true
	end
end

function set_n_numa_nodes(_dw, n_numa_nodes)
	@eval ccall(($(string("set_n_numa_node")), $(_libpath)), Cuint, (Ptr{Void}, Cint, Cint, Cint, Cint), $(_dw._dw), $(n_numa_nodes), $(_dw.modelrepl), $(_dw.datarepl), $(_dw.accessmethod))
end

function set_n_threads_per_node(_dw, n_thread_per_node)
	@eval ccall(($(string("set_n_thread_per_node")), $(_libpath)), Cuint, (Ptr{Void}, Cint, Cint, Cint, Cint), $(_dw._dw), $(n_thread_per_node), $(_dw.modelrepl), $(_dw.datarepl), $(_dw.accessmethod))
end

function register_row2(_dw, func, supress=false)

	global _data_type, _model_type, _libpath, _nogc

	is_safe = check_is_safe(func, Cdouble, (Array{_data_type,1}, Array{_model_type,1}, _dw.shareddatatype))
	if is_safe == false && supress==false
		error("Your function contains LLVM LR `alloc` or `call` other julia functions. We cannot register this function because it protentially is not thread-safe. Use register_row(_dw",",",func,",true) to register this function AT YOUR OWN RISK!")
	end

	const func_c = cfunction(func, Cdouble, (Array{_data_type,1}, Array{_model_type,1}, _dw.shareddatatype))

	append!(_nogc, {func_c, func})

	handle = @eval ccall(($(string("DenseDimmWitted_Register_Row2")), $(_libpath)), Cuint, (Ptr{Void}, Ptr{Void}, Cint, Cint, Cint), $(_dw._dw), $(func_c), $(_dw.modelrepl), $(_dw.datarepl), $(_dw.accessmethod))

	println("[JULIA-DW] Registered Row Function ", func, " Handle=", handle)

	return handle
end

function register_row(_dw, func, supress=false)

	global _data_type, _model_type, _libpath, _nogc

	is_safe = check_is_safe(func, Cdouble, (Array{_data_type,1}, Array{_model_type,1}))
	if is_safe == false && supress==false
		error("Your function contains LLVM LR `alloc` or `call` other julia functions. We cannot register this function because it protentially is not thread-safe. Use register_row(_dw",",",func,",true) to register this function AT YOUR OWN RISK!")
	end

	const func_c = cfunction(func, Cdouble, (Array{_data_type,1}, Array{_model_type,1}))

	append!(_nogc, {func_c, func})

	handle = @eval ccall(($(string("DenseDimmWitted_Register_Row2")), $(_libpath)), Cuint, (Ptr{Void}, Ptr{Void}, Cint, Cint, Cint), $(_dw._dw), $(func_c), $(_dw.modelrepl), $(_dw.datarepl), $(_dw.accessmethod))

	println("[JULIA-DW] Registered Row Function ", func, " Handle=", handle)

	return handle
end

function register_c2r(_dw, func, supress=false)

	is_safe = check_is_safe(func, Cdouble, (Array{_data_type,1}, Cint, Array{Array{Cdouble, 1},1}, Array{_model_type,1}))
	if is_safe == false && supress==false
		error("Your function contains LLVM LR `alloc` or `call` other julia functions. We cannot register this function because it protentially is not thread-safe. Use register_c2r(_dw",",",func,",true) to register this function AT YOUR OWN RISK!")
	end

	global _data_type, _model_type, _libpath, _nogc

	const func_c = cfunction(func, Cdouble, (Array{_data_type,1}, Cint, Array{Array{Cdouble, 1},1}, Array{_model_type,1}))

	append!(_nogc, {func_c, func})

	handle = @eval ccall(($(string("DenseDimmWitted_Register_C2R2")), $(_libpath)), Cuint, (Ptr{Void}, Ptr{Void}, Cint, Cint, Cint), $(_dw._dw), $(func_c), $(_dw.modelrepl), $(_dw.datarepl), $(_dw.accessmethod))

	println("[JULIA-DW] Registered Column-to-row Function ", func, " Handle=", handle)

	return handle
end

function register_c2r2(_dw, func, supress=false)

	is_safe = check_is_safe(func, Cdouble, (Array{_data_type,1}, Cint, Array{Array{Cdouble, 1},1}, Array{_model_type,1}, _dw.shareddatatype))
	if is_safe == false && supress==false
		error("Your function contains LLVM LR `alloc` or `call` other julia functions. We cannot register this function because it protentially is not thread-safe. Use register_c2r(_dw",",",func,",true) to register this function AT YOUR OWN RISK!")
	end

	global _data_type, _model_type, _libpath, _nogc

	const func_c = cfunction(func, Cdouble, (Array{_data_type,1}, Cint, Array{Array{Cdouble, 1},1}, Array{_model_type,1}, _dw.shareddatatype))

	append!(_nogc, {func_c, func})

	handle = @eval ccall(($(string("DenseDimmWitted_Register_C2R2")), $(_libpath)), Cuint, (Ptr{Void}, Ptr{Void}, Cint, Cint, Cint), $(_dw._dw), $(func_c), $(_dw.modelrepl), $(_dw.datarepl), $(_dw.accessmethod))

	println("[JULIA-DW] Registered Column-to-row Function ", func, " Handle=", handle)

	return handle
end

function register_col(_dw, func, supress=false)

	is_safe = check_is_safe(func, Cdouble, (Array{_data_type,1}, Array{_model_type,1}, _dw.shareddatatype))
	if is_safe == false && supress==false
		error("Your function contains LLVM LR `alloc` or `call` other julia functions. We cannot register this function because it protentially is not thread-safe. Use register_row(_dw",",",func,",true) to register this function AT YOUR OWN RISK!")
	end

	global _data_type, _model_type, _libpath, _nogc

	const func_c = cfunction(func, Cdouble, (Array{_data_type,1}, Array{_model_type,1}, _dw.shareddatatype))

	append!(_nogc, {func_c, func})

	handle = @eval ccall(($(string("DenseDimmWitted_Register_Col2")), $(_libpath)), Cuint, (Ptr{Void}, Ptr{Void}, Cint, Cint, Cint), $(_dw._dw), $(func_c), $(_dw.modelrepl), $(_dw.datarepl), $(_dw.accessmethod))

	println("[JULIA-DW] Registered Col Function ", func, " Handle=", handle)

	return handle
end

function register_col(_dw, func, supress=false)

	is_safe = check_is_safe(func, Cdouble, (Array{_data_type,1}, Array{_model_type,1}))
	if is_safe == false && supress==false
		error("Your function contains LLVM LR `alloc` or `call` other julia functions. We cannot register this function because it protentially is not thread-safe. Use register_row(_dw",",",func,",true) to register this function AT YOUR OWN RISK!")
	end

	global _data_type, _model_type, _libpath, _nogc

	const func_c = cfunction(func, Cdouble, (Array{_data_type,1}, Array{_model_type,1}))

	append!(_nogc, {func_c, func})

	handle = @eval ccall(($(string("DenseDimmWitted_Register_Col2")), $(_libpath)), Cuint, (Ptr{Void}, Ptr{Void}, Cint, Cint, Cint), $(_dw._dw), $(func_c), $(_dw.modelrepl), $(_dw.datarepl), $(_dw.accessmethod))

	println("[JULIA-DW] Registered Col Function ", func, " Handle=", handle)

	return handle
end

function register_model_avg(_dw, handle, func, supress=false)

	is_safe = check_is_safe(func, Cdouble, (Array{Array{_model_type,1},1}, Cint, Cint))
	if is_safe == false && supress==false
		error("Your function contains LLVM LR `alloc` or `call` other julia functions. We cannot register this function because it protentially is not thread-safe. Use register_model_avg(_dw",",",func,",true) to register this function AT YOUR OWN RISK!")
	end

	global _data_type, _model_type, _libpath, _nogc

	const func_c = cfunction(func, Cdouble, (Array{Array{Cdouble,1},1}, Cint, Cint))

	append!(_nogc, {func_c, func})

	@eval ccall(($(string("DenseDimmWitted_Register_ModelAvg2")), $(_libpath)), Void, (Ptr{Void}, Cuint, Ptr{Void}, Cint, Cint, Cint), $(_dw._dw), $(handle), $(func_c), $(_dw.modelrepl), $(_dw.datarepl), $(_dw.accessmethod))

	println("[JULIA-DW] Registered Avg Function ", func, " for Func Handle=", handle)

	return nothing

end

function exec(_dw, func_handle)
	global _libpath

	rs = @eval ccall(($(string("DenseDimmWitted_Exec2")), $(_libpath)), Cdouble, (Ptr{Void}, Cuint, Cint, Cint, Cint), $(_dw._dw), $(func_handle), $(_dw.modelrepl), $(_dw.datarepl), $(_dw.accessmethod))

	return rs
end

export libpath, get_libpath, open, register_row, exec

end

max_stars_repo_name: GloriaGo/dw

max_issues_repo_name: GloriaGo/dw

max_forks_repo_name: GloriaGo/dw


text:
type IntegerTriangle{T<:Integer}
    a::T
    b::T
    c::T
    p::T
    σ::T
end

function IntegerTriangle{T<:Integer}(a::T, b::T, c::T)
    p = a + b + c
    s = div(p, 2)
    σ = isqrt(s*(s-a)*(s-b)*(s-c))
    (x, y, z) = sort([a, b, c])
    IntegerTriangle(x, y, z, p, σ)
end

function isprimheronian{T<:Integer}(a::T, b::T, c::T)
    p = a + b + c
    iseven(p) || return false
    gcd(a, b, c) == 1 || return false
    s = div(p, 2)
    t = s*(s-a)*(s-b)*(s-c)
    0 < t || return false
    σ = isqrt(t)
    σ^2 == t
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Test SeisBlock type and methods

@testset "SeisBlock" begin

    @testset "Constructor" begin
        b = SeisBlock(rand(Float32,10,10))
        @test b.fileheader.bfh.ns == 10
        @test b.fileheader.bfh.DataSampleFormat == 5
        vec_ = rand(Float32,10)
        c = SeisBlock(vec_)
        mat_ = reshape(vec_,10,1)
        d = SeisBlock(mat_)
        @test c.data == d.data
    end

    @testset "Methods" begin
        b = SeisBlock(rand(Float32,10,10))
        @test size(b) == (10,10)
    end

    @testset "set_traceheader" begin
        b = SeisBlock(rand(Float32,10,10))
        set_traceheader!(b.traceheaders, :SourceX, 102*ones(Int32, 10))
        @test b.traceheaders[3].SourceX == 102
    end

    @testset "set_fileheader" begin
        b = SeisBlock(rand(Float32,10,10))
        set_fileheader!(b.fileheader.bfh, :ns, 12)
        @test b.fileheader.bfh.ns == 12
    end

    @testset "set_header" begin
        b = SeisBlock(rand(Float32,10,10))
        set_header!(b, :ns, 12)
        @test b.fileheader.bfh.ns == 12
        @test b.traceheaders[1].ns == 12

        set_header!(b, "ns", 13)
        @test b.fileheader.bfh.ns == 13
        @test b.traceheaders[1].ns == 13

        set_header!(b, :SourceX, 12)
        @test b.traceheaders[1].SourceX == 12

        set_header!(b, "SourceY", 13)
        @test b.traceheaders[1].SourceY == 13

        set_header!(b, :GroupY, 14*ones(Int32, 10))
        @test b.traceheaders[1].GroupY == 14

        set_header!(b, "GroupX", 15*ones(Int32, 10))
        @test b.traceheaders[1].GroupX == 15
    end
end


text:
@testset "Testing parameters.jl" begin

    @testset "Testing parameters creation and update" begin

        # Create
        parameters = AbrkgaObop.OBOPParameters(0.99)

        @test parameters.gama == 0.99
        @test parameters.max_generations == 137
        @test parameters.population_size == 100
        @test parameters.elite_size == 10
        @test parameters.mutation_size == 20

        # Update
        AbrkgaObop.update_parameters!(100,parameters)

        @test parameters.max_generations == 137
        @test parameters.population_size == 99
        @test parameters.elite_size == 21
        @test parameters.mutation_size == 9
    end;

end

max_stars_repo_name: LuizHNLorena/AbrkgaObop

max_issues_repo_name: LuizHNLorena/AbrkgaObop

max_forks_repo_name: LuizHNLorena/AbrkgaObop


text:
@testset "encoder/sixel" begin
    # Force sixel encoding
    old_encoder = ImageInTerminal.encoder_backend[1]
    ImageInTerminal.encoder_backend[1] = :Sixel

    @testset "lena" begin
        img = imresize(lena, 128, 128)
        io = IOBuffer()
        imshow256(io, img)
        seek(io, 0)
        # we can't reference test it right now, because ReferenceTests calls
        # FileIO.load, which returns an array for sixel file.
        @test assess_psnr(sixel_decode(io), img) > 30
    end
    @testset "small images" begin
        # small images still use ImageInTerminal's fallback encoding
        img = imresize(lena, 10, 10)
        io = IOBuffer()
        ensurecolor(imshow256, io, img)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))[1]
        @test !startswith(res, "\ePq\"") # not sixel encoding
    end
    @testset "vector" begin
        # vectors, no matter how large it is, does not use sixel
        img = rgb_line
        io = IOBuffer()
        ensurecolor(imshow256, io, img)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))[1]
        @test !startswith(res, "\ePq\"") # not sixel encoding

        io = IOBuffer()
        img = repeat(img, 10)
        ensurecolor(imshow256, io, img)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))[1]
        @test !startswith(res, "\ePq\"") # not sixel encoding
    end

    # restore encoder to previous one
    ImageInTerminal.encoder_backend[1] = old_encoder
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: apigatewaymanagementapi
using AWS.Compat
using AWS.UUIDs

"""
    delete_connection(connection_id)
    delete_connection(connection_id, params::Dict{String,<:Any})

Delete the connection with the provided id.

# Arguments
- `connection_id`:

"""
delete_connection(connectionId; aws_config::AbstractAWSConfig=global_aws_config()) = apigatewaymanagementapi("DELETE", "/@connections/$(connectionId)"; aws_config=aws_config)
delete_connection(connectionId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = apigatewaymanagementapi("DELETE", "/@connections/$(connectionId)", params; aws_config=aws_config)

"""
    get_connection(connection_id)
    get_connection(connection_id, params::Dict{String,<:Any})

Get information about the connection with the provided id.

# Arguments
- `connection_id`:

"""
get_connection(connectionId; aws_config::AbstractAWSConfig=global_aws_config()) = apigatewaymanagementapi("GET", "/@connections/$(connectionId)"; aws_config=aws_config)
get_connection(connectionId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = apigatewaymanagementapi("GET", "/@connections/$(connectionId)", params; aws_config=aws_config)

"""
    post_to_connection(data, connection_id)
    post_to_connection(data, connection_id, params::Dict{String,<:Any})

Sends the provided data to the specified connection.

# Arguments
- `data`: The data to be sent to the client specified by its connection id.
- `connection_id`: The identifier of the connection that a specific client is using.

"""
post_to_connection(Data, connectionId; aws_config::AbstractAWSConfig=global_aws_config()) = apigatewaymanagementapi("POST", "/@connections/$(connectionId)", Dict{String, Any}("Data"=>Data); aws_config=aws_config)
post_to_connection(Data, connectionId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = apigatewaymanagementapi("POST", "/@connections/$(connectionId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Data"=>Data), params)); aws_config=aws_config)