text:
using Clustering,
    LightGraphs,
    SparseArrays,
    LinearAlgebra,
    Arpack
export NystromMethod,
       LandmarkBasedRepresentation,
       DNCuts

"""
Large Scale Spectral Clustering with Landmark-Based Representation
Xinl ei Chen Deng Cai

# Members
- `landmark_selector::{T <: AbstractLandmarkSelection}` Method for extracting landmarks
- `number_of_landmarks::Integer` Number of landmarks to obtain
- `n_neighbors::Integer` Number of nearest neighbors
- `nev::Integer` Number of eigenvectors
- `w::Function` Number of clusters to obtain
- `normalize::Bool`
"""
struct LandmarkBasedRepresentation{T <: AbstractLandmarkSelection}
    landmark_selector::T
    number_of_landmarks::Integer
    n_neighbors::Integer
    nev::Integer
    w::Function
    normalize::Bool
end

"""
```
embedding(cfg::LandmarkBasedRepresentation, X)
```

"""
function embedding(cfg::LandmarkBasedRepresentation, X)

    n = number_of_patterns(X)
    p = cfg.number_of_landmarks
    landmark_indices = select_landmarks(cfg.landmark_selector, p, X)
    landmarks = get_element(X, landmark_indices)
    neighbors_cfg = KNNNeighborhood(landmarks, cfg.n_neighbors)

    I = zeros(Integer, n*cfg.n_neighbors)
    J = zeros(Integer, n*cfg.n_neighbors)
    V = zeros(n*cfg.n_neighbors)
    for i = 1:n
        i_neighbors =  neighbors(neighbors_cfg, get_element(X, i))
        weights = cfg.w(i, [], get_element(X, i), get_element(landmarks, i_neighbors))
        weights ./= sum(weights)
        for (j, (neigh, w)) in enumerate(zip(i_neighbors, weights))
            I[(i-1)*cfg.n_neighbors + j] = neigh
            J[(i-1)*cfg.n_neighbors + j] = i
            V[(i-1)*cfg.n_neighbors + j] = w
        end
    end
    Z = sparse(I, J, V, p, n)

    (svd) = svds(Z*Z', nsv=cfg.nev)[1]
    v = svd.S
    S = spdiagm(0=>1 ./ (sqrt.(v)))[1]
    B = S * svd.U' * Z
    if cfg.normalize
        B = normalize_cols(B)
    end
    return B'
end

"""
```julia

type NystromMethod{T<:AbstractLandmarkSelection}
landmarks_selector::T
number_of_landmarks::Integer
w::Function
nvec::Integer
end
```
The type ```NystromMethod``` proposed in  Spectral Grouping Using the Nystrom Method by Charless
Fowlkes, Serge Belongie, Fan Chung, and Jitendra Malik. It has to be defined:

- `landmarks_selector::T<:AbstractLandmarkSelection`. A mechanism to select the sampled
points.
- `number_of_landmarks::Integer`. The number of points to sample
- `w::Function`. The weight function for compute the similiarity. The signature of the weight function has to be `weight(i, j, e1,e2)`. Where `e1` and `e2` ara the data elements i-th and j-th respectivily, obtained via `get_element`, usually is a vector.
- `nvec::Integer`. The number of eigenvector to obtain.
- `threaded::Bool`. Default: True. Specifies whether the threaded version is used.
"""
struct NystromMethod{T <: AbstractLandmarkSelection} <: EigenvectorEmbedder
    landmarks_selector::T
    number_of_landmarks::Integer
    w::Function
    nvec::Integer
    threaded::Bool
end
function NystromMethod(landmarks_selector::T,
                       number_of_landmarks::Integer,
                       w::Function,
                       nvec::Integer) where T <: AbstractLandmarkSelection
    return NystromMethod(landmarks_selector,
                         number_of_landmarks,
                         w,
                         nvec,
                         true)
end

"""
create_A_B(cfg::NystromMethod, landmarks::Vector{Int},X)

Arguments:
- `cfg::NystromMethod`. The method configuration.
- `landmarks::Vector{T}`. A vector of integer that containts the \$n\$ indexes sampled from the data.
- `X` is the data that containt \$ N \$ patterns.

Let \$ W \\in \\mathbb{R}^{N \\times N}, W = \\begin{bmatrix} A & B^T \\\\ B & C \\end{bmatrix}, A
\\in \\mathbb{R}^{ n \\times n }, B \\in \\mathbb{R}^{(N-n) \\times n}, C \\in
\\mathbb{R}^{(N-n)\\times (N-n)} \$ . \$A\$ represents the subblock of weights among the random
samples, \$B\$ contains the weights from the random samples to the rest of the pixels, and
\$C\$ contains the weights between all of the remaining pixels.
The function computes \$A\$ and \$B\$ from the data ```X``` using the weight function defined in
```cfg```.
"""
function create_A_B(cfg::NystromMethod, landmarks::Vector{<:Integer}, X)
    n = number_of_patterns(X)
    p = length(landmarks)
    indexes_b = setdiff(collect(1:n), landmarks)
    m = length(indexes_b)
    n = p
    A = zeros(Float32, p, p)
    B = zeros(Float32, p, m)
    qq = length(get_element(X, 1))
    landmarks_m = zeros(Float32, length(get_element(X, 1)), length(landmarks))
    # Get a copy of the landamrks
    for j = 1:length(landmarks)
        get_element!(view(landmarks_m, :, j), X, landmarks[j])
    end
    Threads.@threads for j = 1:length(landmarks)
        A[:,j] = cfg.w(landmarks[j], landmarks, view(landmarks_m, :, j), landmarks_m)
    end

    vec_k = zeros(Float32, length(get_element(X, 1)), Threads.nthreads())
    Threads.@threads for k = 1:length(indexes_b)
       thread_id = Threads.threadid()
       get_element!(view(vec_k, :, thread_id), X, indexes_b[k])
       B[:,k] = cfg.w(indexes_b[k], landmarks, view(vec_k, :, thread_id), landmarks_m)
    end
    return (A, B)
end

function create_A_B_single_thread(cfg::NystromMethod, landmarks::Vector{<:Integer}, X)
    n = number_of_patterns(X)
    p = length(landmarks)
    indexes_b = setdiff(collect(1:n), landmarks)
    m  = length(indexes_b)
    n = p
    A = zeros(Float32, p, p)
    B = zeros(Float32, p, m)
    qq = length(get_element(X, 1))
    landmarks_m = zeros(Float32, length(get_element(X, 1)), length(landmarks))
    for j = 1:length(landmarks)
        get_element!(view(landmarks_m, :, j), X, landmarks[j])
    end
    for j = 1:length(landmarks)
        A[:,j] = cfg.w(landmarks[j], landmarks, view(landmarks_m, :, j), landmarks_m)
    end
    vec_k = zeros(Float32, length(get_element(X, 1)))
    for k = 1:length(indexes_b)
       get_element!(vec_k, X, indexes_b[k])
       B[:,k] = cfg.w(indexes_b[k], landmarks, vec_k, landmarks_m)
    end
    return (A, B)
end
"""
create_A_B(cfg::NystromMethod, X)

# Arguments:
- `cfg::NystromMethod`
- `X`

#Return values
- Sub-matrix A
- Sub-matrix B
- `Vector{Int}`. The sampled points used build the sub-matrices

This is an overloaded method. Computes the submatrix A and B according to
[`create_A_B(::NystromMethod, ::Vector{Int}, ::Any)`](@ref).
Returns the two submatrices and the sampled points used to calcluate it
"""
function create_A_B(cfg::NystromMethod, X)
  landmarks = select_landmarks(cfg.landmarks_selector, cfg.number_of_landmarks, X)
  if (cfg.threaded)
      (A, B) = create_A_B(cfg::NystromMethod, landmarks, X)
  else
      (A, B) = create_A_B_single_thread(cfg::NystromMethod, landmarks, X)
  end
  return (A, B, landmarks)
end
"""
embedding(cfg::NystromMethod, X)

This is an overloaded function
"""
function embedding(cfg::NystromMethod, X)
    n = number_of_patterns(X)
    landmarks = select_landmarks(cfg.landmarks_selector,
                                 cfg.number_of_landmarks,
                                 X)
    return embedding(cfg, landmarks, X)
end
"""
embedding(cfg::NystromMethod, landmarks::Vector{Int}, X)
# Arguments
- `cfg::[NystromMethod](@ref)`
- `landmarks::Vector{Int}`
- `x::Any`

# Return values
- `(E, L)`: The approximated eigenvectors, the aprooximated eigenvalues
Performs the eigenvector embedding according to
"""
function embedding(cfg::NystromMethod, landmarks::Vector{<:Integer}, X)
    if (cfg.threaded)
        (A, B) = create_A_B(cfg, landmarks, X)
    else
        (A, B) = create_A_B_single_thread(cfg, landmarks, X)
    end
    return embedding(cfg, A, B, landmarks)
end

function compute_dhat(AA::Matrix{T}, BB::Matrix{T}) where T
    n = size(AA, 1)
    m = size(BB, 2)
    dhat = zeros(T, n + m)
    dhat[1:n] = sum(AA, dims = 1) + sum(BB, dims = 2)'
    dhat[n + 1:end] = sum(BB, dims = 1) + sum(BB, dims = 2)' * pinv(AA) * BB
    dhat[dhat .< 0] .= 0
    return 1 ./ (sqrt.(dhat) .+ eps())
end
function compute_V(AA::Matrix{T}, BB::Matrix{T}, nvec::Integer) where T <: Number
    n = size(AA, 1)
    m = size(BB, 2)
    Asi = real(sqrt(Symmetric(pinv(AA))))
    F = svd(AA + ((Asi * (BB * BB')) * Asi))
    V_1 = (Asi * F.U) .* vec((1 ./ (sqrt.(F.S) .+ eps())))'
    VA = AA * V_1[:,1:nvec + 1]
    VB = BB' * V_1[:,1:nvec + 1]
    return vcat(VA, VB)
end
function normalize_A_and_B!(AA::Matrix, BB::Matrix)
    n    = size(AA, 1)
    m    = size(BB, 2)
    dhat = compute_dhat(AA, BB)
    vv   = view(dhat, 1:n)
    vb   = view(dhat, n .+ (1:m))
    for I in CartesianIndices(size(AA))
        @inbounds AA[I] *= vv[I[1]] * vv[I[2]]
    end
    for I in CartesianIndices(size(BB))
        @inbounds BB[I] *= vv[I[1]] * vb[I[2]]
    end
end
"""
embedding(cfg::NystromMethod, A::Matrix, B::Matrix, landmarks::Vector{Int})

Performs the eigenvector approximation given the two submatrices A and B.
"""
function embedding(cfg::NystromMethod, AA::Matrix, BB::Matrix, landmarks::Vector{<:Integer})
    n = size(AA, 1)
    m = size(BB, 2)
    normalize_A_and_B!(AA, BB)
    V = compute_V(AA, BB, cfg.nvec)
    indexes_b = setdiff(collect(1:(n + m)), landmarks)
    indexes   = sortperm(vcat(landmarks, indexes_b))
    for i = 2:cfg.nvec + 1
        V[:,i] = V[:,i] ./ V[:,1]
    end
    return V[indexes,2:cfg.nvec + 1]

end
"""
```@julia
type DNCuts
```
# Multiscale Combinatorial Grouping for Image Segmentation and Object Proposal Generation
## Jordi Pont-Tuset, Pablo Arbeláez, Jonathan T. Barron, Member, Ferran Marques, Jitendra Malik
"""
struct DNCuts <: EigenvectorEmbedder
    scales::Integer
    nev::Integer
    img_size
end
function pixel_decimate(img_size::Tuple{Int,Int}, steps)
    (nr, nc) = img_size
    decimated = CartesianIndices(img_size)[1:steps:nr, 1:steps:nc]
    return (LinearIndices(img_size)[decimated][:], size(decimated))
end

embedding(d::DNCuts, g::Graph) = embedding(d, adjacency_matrix(g))

"""
```
embedding(d::DNCuts, L)
```
"""
function embedding(d::DNCuts, W::AbstractMatrix{T}) where T<:Number
    matrices = []
    img_size = d.img_size
    for j = 1:d.scales
        (idx, img_size) = pixel_decimate(img_size,  2)
        B = W[:, idx]
        C = Diagonal(vec(1 ./ sum(B, dims=2))) * B
        push!(matrices, C)
        W = C' * B
    end
    ss = ShiMalikLaplacian(d.nev)
    V  = embedding(ss, NormalizedLaplacian(NormalizedAdjacency(CombinatorialAdjacency(W))))
    for s = d.scales:-1:1
        V = matrices[s] * V
    end
    return svd_whiten(V)
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api.concept.thing;
#
# import typedb.client.api.TypeDBTransaction;
# import typedb.client.api.concept.type.EntityType;
#
# import javax.annotation.CheckReturnValue;
#
# public interface Entity extends Thing {
#
#     @Override
#     @CheckReturnValue
#     default boolean isEntity() {
#         return true;
#     }
#
#     @Override
#     @CheckReturnValue
#     EntityType getType();
#
#     @Override
#     @CheckReturnValue
#     Entity.Remote asRemote(TypeDBTransaction transaction);
#
#     interface Remote extends Thing.Remote, Entity {}
# }

ext: jl


text:
################################################################################
# MIT License
#
# Copyright (c) 2021 Hajime Nakagami<nakagami@gmail.com>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
################################################################################
module Firebird
using Dates, TimeZones, DBInterface, Tables, Parsers, DecFP

export DBInterface

include("consts.jl")
include("errmsgs.jl")
include("tz_map.jl")
include("xsqlvar.jl")
include("utils.jl")
include("srp.jl")
include("arc4.jl")
include("chacha20.jl")
include("decfloat.jl")
include("wireprotocol.jl")
include("transaction.jl")
include("connection.jl")
include("statement.jl")
include("cursor.jl")
include("execute.jl")

Base.close(conn::Connection) = DBInterface.close!(conn)
Base.isopen(conn::Connection) = Firebird.isopen(conn)

end # module


text:
using Documenter, NextGenSeqUtils
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/MurrellGroup/NextGenSeqUtils.jl.git",
    julia = "nightly",
    osx = "osx"

makedocs()


text:
using Base: @get!, promote_eltype
using Base.Broadcast: bitcache_chunks, bitcache_size, dumpbitcache

_broadcast_shape(x...) = Base.to_shape(Base.Broadcast.broadcast_indices(x...))

# Get ref for value for a PooledDataArray, adding to the pool if
# necessary
_unsafe_pdaref!(Bpool, Brefdict::Dict, val::Missing) = 0
function _unsafe_pdaref!(Bpool, Brefdict::Dict, val)
    @get! Brefdict val begin
        push!(Bpool, val)
        length(Bpool)
    end
end

# Generate a branch for each possible combination of missing/not missing. This
# gives good performance at the cost of 2^narrays branches.
function gen_na_conds(f, nd, arrtype, outtype,
    daidx=findall(t -> t <: DataArray || t <: PooledDataArray, arrtype), pos=1, ismissing=())

    if pos > length(daidx)
        args = Any[Symbol("v_$(k)") for k = 1:length(arrtype)]
        for i = 1:length(daidx)
            if ismissing[i]
                args[daidx[i]] = missing
            end
        end

        # Needs to be gensymmed so that the compiler won't box it
        val = gensym("val")
        quote
            $val = $(Expr(:call, f, args...))
            $(if outtype <: DataArray
                :(@inbounds unsafe_dasetindex!(Bdata, Bc, $val, ind))
            elseif outtype <: PooledDataArray
                :(@inbounds (@nref $nd Brefs i) = _unsafe_pdaref!(Bpool, Brefdict, $val))
            end)
        end
    else
        k = daidx[pos]
        quote
            if $(Symbol("ismissing_$(k)"))
                $(gen_na_conds(f, nd, arrtype, outtype, daidx, pos+1, tuple(ismissing..., true)))
            else
                $(if arrtype[k] <: DataArray
                    :(@inbounds $(Symbol("v_$(k)")) = $(Symbol("data_$(k)"))[$(Symbol("state_$(k)_0"))])
                else
                    :(@inbounds $(Symbol("v_$(k)")) = $(Symbol("pool_$(k)"))[$(Symbol("r_$(k)"))])
                end)
                $(gen_na_conds(f, nd, arrtype, outtype, daidx, pos+1, tuple(ismissing..., false)))
            end
        end
    end
end

# Broadcast implementation for DataArrays
#
# TODO: Fall back on faster implementation for same-sized inputs when
# it is safe to do so.
Base.map!(f::F, B::Union{DataArray, PooledDataArray}, A0::AbstractArray, As::AbstractArray...) where {F} =
        broadcast!(f, B, A0, As...)
Base.map!(f::F, B::Union{DataArray, PooledDataArray}, A0, As...) where {F} =
        broadcast!(f, B, A0, As...)

@generated function _broadcast!(f, B::Union{DataArray, PooledDataArray}, As...)

    F  = :(f)
    nd = ndims(B)
    N  = length(As)

    dataarrays = find(t -> t <: DataArray, As)

    quote
        @boundscheck Base.Broadcast.check_broadcast_indices(indices(B), As...)
        # check_broadcast_shape(size(B), As...)
        @nexprs $N i->(A_i = As[i])

        @assert ndims(B) == $nd

        # Set up input DataArray/PooledDataArrays
        $(Expr(:block, [
            As[k] <: DataArray ? quote
                $(Symbol("na_$(k)")) = $(Symbol("A_$(k)")).na.chunks
                $(Symbol("data_$(k)")) = $(Symbol("A_$(k)")).data
                $(Symbol("state_$(k)_0")) = $(Symbol("state_$(k)_$(nd)")) = 1
                @nexprs $nd d->($(Symbol("skip_$(k)_d")) = size($(Symbol("data_$(k)")), d) == 1)
            end : As[k] <: PooledDataArray ? quote
                $(Symbol("refs_$(k)")) = $(Symbol("A_$(k)")).refs
                $(Symbol("pool_$(k)")) = $(Symbol("A_$(k)")).pool
            end : nothing
        for k = 1:N]...))

        # Set up output DataArray/PooledDataArray
        $(if B <: DataArray
            quote
                Bdata = B.data
                # Copy in case aliased
                # TODO: check for aliasing?
                Bna = falses(size(Bdata))
                Bc = Bna.chunks
                ind = 1
            end
        elseif B <: PooledDataArray
            quote
                Bpool = B.pool = similar(B.pool, 0)
                Brefs = B.refs
                Brefdict = Dict{eltype(Bpool),eltype(Brefs)}()
            end
        end)

        @nloops($nd, i, $(B <: DataArray ? (:Bdata) : (:Brefs)),
            # pre
            d->($(Expr(:block, [
                As[k] <: DataArray ? quote
                    $(Symbol("state_$(k)_")){d-1} = $(Symbol("state_$(k)_d"));
                    $(Symbol("j_$(k)_d")) = $(Symbol("skip_$(k)_d")) ? 1 : i_d
                end : (As[k] <: AbstractArray ? quote
                    $(Symbol("j_$(k)_d")) = size($(Symbol("A_$(k)")), d) == 1 ? 1 : i_d
                end : quote
                    $(Symbol("j_$(k)_d")) = 1
                end)
            for k = 1:N]...))),

            # post
            d->($(Expr(:block, [quote
                $(Symbol("skip_$(k)_d")) || ($(Symbol("state_$(k)_d")) = $(Symbol("state_$(k)_0")))
            end for k in dataarrays]...))),

            # body
            begin
                # Advance iterators for DataArray and determine missing status
                $(Expr(:block, [
                    As[k] <: DataArray ? quote
                        @inbounds $(Symbol("ismissing_$(k)")) = Base.unsafe_bitgetindex($(Symbol("na_$(k)")), $(Symbol("state_$(k)_0")))
                    end : As[k] <: PooledDataArray ? quote
                        @inbounds $(Symbol("r_$(k)")) = @nref $nd $(Symbol("refs_$(k)")) d->$(Symbol("j_$(k)_d"))
                        $(Symbol("ismissing_$(k)")) = $(Symbol("r_$(k)")) == 0
                    end : nothing
                for k = 1:N]...))

                # Extract values for other type
                $(Expr(:block, [
                    As[k] <: AbstractArray  && !(As[k] <: AbstractDataArray) ? quote
                        # ordinary AbstractArrays
                        @inbounds $(Symbol("v_$(k)")) = @nref $nd $(Symbol("A_$(k)")) d->$(Symbol("j_$(k)_d"))
                    end : quote
                        # non AbstractArrays (e.g. Strings and Numbers)
                        @inbounds $(Symbol("v_$(k)")) = $(Symbol("A_$(k)"))
                    end
                for k = 1:N]...))

                # Compute and store return value
                $(gen_na_conds(F, nd, As, B))

                # Increment state
                $(Expr(:block, [:($(Symbol("state_$(k)_0")) += 1) for k in dataarrays]...))
                $(if B <: DataArray
                    :(ind += 1)
                end)
            end)

        $(if B <: DataArray
            :(B.na = Bna)
        end)

        return B
    end
end
Base.Broadcast.broadcast!(f, B::Union{DataArray, PooledDataArray}, ::Type{T}, As...) where T =
    _broadcast!((t...) -> f(T, t...), B, As...)
Base.Broadcast.broadcast!(f, B::Union{DataArray, PooledDataArray}, A0::Number, As::Number...) =
    _broadcast!(f, B, A0, As...)
Base.Broadcast.broadcast!(f, B::Union{DataArray, PooledDataArray}, A0, As...) =
    _broadcast!(f, B, A0, As...)

Base.Broadcast.promote_containertype(::Type{DataArray}, ::Type{DataArray})             = DataArray
Base.Broadcast.promote_containertype(::Type{PooledDataArray}, ::Type{PooledDataArray}) = PooledDataArray
Base.Broadcast.promote_containertype(::Type{DataArray}, ::Type{Array})                 = DataArray
Base.Broadcast.promote_containertype(::Type{PooledDataArray}, ::Type{Array})           = PooledDataArray
Base.Broadcast.promote_containertype(::Type{Array}, ::Type{DataArray})                 = DataArray
Base.Broadcast.promote_containertype(::Type{Array}, ::Type{PooledDataArray})           = PooledDataArray
Base.Broadcast.promote_containertype(::Type{DataArray}, ::Type{PooledDataArray})       = DataArray
Base.Broadcast.promote_containertype(::Type{PooledDataArray}, ::Type{DataArray})       = DataArray
Base.Broadcast.promote_containertype(::Type{DataArray}, ct)                            = DataArray
Base.Broadcast.promote_containertype(::Type{PooledDataArray}, ct)                      = PooledDataArray
Base.Broadcast.promote_containertype(ct, ::Type{DataArray})                            = DataArray
Base.Broadcast.promote_containertype(ct, ::Type{PooledDataArray})                      = PooledDataArray
Base.Broadcast._containertype(::Type{T}) where T<:DataArray           = DataArray
Base.Broadcast._containertype(::Type{T}) where T<:PooledDataArray     = PooledDataArray
Base.Broadcast.broadcast_indices(::Type{T}, A) where T<:AbstractDataArray = indices(A)

@inline function broadcast_t(f, ::Type{T}, shape, A, Bs...) where {T}
    dest = Base.Broadcast.containertype(A, Bs...)(Missings.T(T), Base.index_lengths(shape...))
    return broadcast!(f, dest, A, Bs...)
end

# This is mainly to handle ismissing.(x) since ismissing is probably the only
# function that can guarantee that missings will never propagate
@inline function broadcast_t(f, ::Type{Bool}, shape, A, Bs...)
    dest = similar(BitArray, shape)
    return broadcast!(f, dest, A, Bs...)
end

# This one is almost identical to the version in Base and can hopefully be
# removed at some point. The main issue in Base is that it tests for
# isleaftype(T) which is false for Union{T,Missing}. If the test in Base
# can be modified to cover simple unions of leaftypes then this method
# can probably be deleted and the two _t methods adjusted to match the Base
# invokation from Base.Broadcast.broadcast_c
@inline function Base.Broadcast.broadcast_c(f, ::Type{S}, A, Bs...) where {S<:AbstractDataArray}
    T     = Base.Broadcast._broadcast_eltype(f, A, Bs...)
    shape = Base.Broadcast.broadcast_indices(A, Bs...)
    return broadcast_t(f, T, shape, A, Bs...)
end

# This one is much faster than normal broadcasting but the method won't get called
# in fusing operations like (!).(ismissing.(x))
Base.broadcast(::typeof(ismissing), da::DataArray) = copy(da.na)


text:
intensities = [][]
for every x_loc:
    for every y_loc:
        total_dists = []
        ind = 0
        for every plane_pos in plane_locations:
            # plane_pos = (x_plane, y_plane, z_plane)
            # z_loc = 0
            total_dists[ind] = sqrt(dX^2 + dY^2 + dZ^2)
            ind+=1
        # plane_data is a list of (range[], intensity[], timestamp)
        LIST_SHIFT_PER_TIME = scale # defined configuration constant
        ind = 0
        for every dist in total_dists:
            time_shift = dist_to_time(dist-total_dists[ind-1])
            index_shift_val = time_shift/LIST_SHIFT_PER_TIME
            plane_data[ind] = shift_list( plane_data[ind] , index_shift_val )
            ind+=1
        total_data_sum = []
        for every data in plane_data :
            total_data_sum+=data
        intensity = total_data_sum #UNFINISHED
        intensities[x_loc][y_loc] = peak_intensity
#done

max_stars_repo_name: jjimmykang/bwsi-backprojection

max_issues_repo_name: jjimmykang/bwsi-backprojection

max_forks_repo_name: jjimmykang/bwsi-backprojection


text:
using Test
include(joinpath(dirname(@__FILE__), "../../src/includeall.jl"))

@testset "Initializing one-dimensional derivatives" begin
    m = Li2020()
    delete!(get_derivatives(m), :∂p_∂w)
    delete!(get_derivatives(m), :∂²p_∂w²)
    init_derivatives!(m, Dict{Symbol, Vector{Int}}(:p => [1]))
    @test haskey(get_derivatives(m), :∂p_∂w)
    @test !haskey(get_derivatives(m), :∂²p_∂w²)
    m = Li2020()
    delete!(get_derivatives(m), :∂p_∂w)
    delete!(get_derivatives(m), :∂²p_∂w²)
    init_derivatives!(m, Dict{Symbol, Vector{Int}}(:p => [1, 2]))
    @test haskey(get_derivatives(m), :∂p_∂w)
    @test haskey(get_derivatives(m), :∂²p_∂w²)
    m = Li2020()
    delete!(get_derivatives(m), :∂p_∂w)
    delete!(get_derivatives(m), :∂²p_∂w²)
    init_derivatives!(m, Dict{Symbol, Vector{Int}}(:p => [2]))
    @test !haskey(get_derivatives(m), :∂p_∂w)
    @test haskey(get_derivatives(m), :∂²p_∂w²)
end

@testset "Initializing two-dimensional derivatives" begin
    m = Li2020()
    m.state_variables[:x] = 2
    delete!(get_derivatives(m), :∂p_∂w)
    delete!(get_derivatives(m), :∂²p_∂w²)

    # Just request 1 derivative
    init_derivatives!(m, Dict{Symbol, Vector{Tuple{Int,Int}}}(:p => [(1, 0)]))
    @test haskey(get_derivatives(m), :∂p_∂w)
    @test !haskey(get_derivatives(m), :∂²p_∂w²)
    delete!(get_derivatives(m), :∂p_∂w)
    init_derivatives!(m, Dict{Symbol, Vector{Tuple{Int,Int}}}(:p => [(1, 0), (2, 0)]))
    @test haskey(get_derivatives(m), :∂p_∂w)
    @test haskey(get_derivatives(m), :∂²p_∂w²)
    delete!(get_derivatives(m), :∂p_∂w)
    delete!(get_derivatives(m), :∂²p_∂w²)
    init_derivatives!(m, Dict{Symbol, Vector{Tuple{Int,Int}}}(:p => [(2, 0)]))
    @test !haskey(get_derivatives(m), :∂p_∂w)
    @test haskey(get_derivatives(m), :∂²p_∂w²)
    delete!(get_derivatives(m), :∂²p_∂w²)

    # Just request 1 derivative
    init_derivatives!(m, Dict{Symbol, Vector{Tuple{Int,Int}}}(:p => [(0, 1)]))
    @test haskey(get_derivatives(m), :∂p_∂x)
    @test !haskey(get_derivatives(m), :∂²p_∂x²)
    @test length(get_derivatives(m)) == 1
    delete!(get_derivatives(m), :∂p_∂x)
    init_derivatives!(m, Dict{Symbol, Vector{Tuple{Int,Int}}}(:p => [(0, 1), (0, 2)]))
    @test haskey(get_derivatives(m), :∂p_∂x)
    @test haskey(get_derivatives(m), :∂²p_∂x²)
    @test length(get_derivatives(m)) == 2
    delete!(get_derivatives(m), :∂p_∂x)
    delete!(get_derivatives(m), :∂²p_∂x²)
    init_derivatives!(m, Dict{Symbol, Vector{Tuple{Int,Int}}}(:p => [(0, 2)]))
    @test haskey(get_derivatives(m), :∂²p_∂x²)
    @test length(get_derivatives(m)) == 1
    delete!(get_derivatives(m), :∂²p_∂x²)

    # Now do multiple variables
    init_derivatives!(m, Dict{Symbol, Vector{Tuple{Int,Int}}}(:p => [(1, 1)]))
    @test haskey(get_derivatives(m), :∂²p_∂w∂x)
    @test length(get_derivatives(m)) == 1
    delete!(get_derivatives(m), :∂²p_∂w∂x)
    init_derivatives!(m, Dict{Symbol, Vector{Tuple{Int,Int}}}(:p => [(1, 0), (0, 1)]))
    @test haskey(get_derivatives(m), :∂p_∂w)
    @test haskey(get_derivatives(m), :∂p_∂x)
    @test length(get_derivatives(m)) == 2
    delete!(get_derivatives(m), :∂p_∂w)
    delete!(get_derivatives(m), :∂p_∂x)
    init_derivatives!(m, Dict{Symbol, Vector{Tuple{Int,Int}}}(:p => [(2, 0), (0, 2)]))
    @test haskey(get_derivatives(m), :∂²p_∂w²)
    @test haskey(get_derivatives(m), :∂²p_∂x²)
    @test length(get_derivatives(m)) == 2
    delete!(get_derivatives(m), :∂²p_∂w²)
    delete!(get_derivatives(m), :∂²p_∂x²)
end

@testset "Standard derivatives for continuous-time models" begin
    @test standard_derivs(1) == Vector{Int}([1, 2])
    @test standard_derivs(2) == Vector{Tuple{Int, Int}}([(1, 0), (0, 1), (2, 0), (1, 1), (0, 2)])
end


text:
using Fzf
using Test

@testset "Fzf.jl" begin
    # Write your tests here.
end


text:
@testset "389.find-the-difference.jl" begin
    @test find_the_difference("abcd", "abcde") == 'e'
    @test find_the_difference("", "y") == 'y'
    # rd_chars = rand('a':'z', 10000)
    # s = *(rd_chars...);
    # ch = rand('a':'z')
    # push!(rd_chars, ch)
    # t = *(shuffle!(rd_chars)...)
    # @test find_the_difference(s, t) == ch
end


text:
using SphericalFunctions
using Test, Random, ProgressMeter
using FastTransforms, Quaternionic

include("test_utilities.jl")

enabled_tests = lowercase.(ARGS)

help = ("help" ∈ enabled_tests || "--help" ∈ enabled_tests)
helptests = []

# This block is cribbed from StaticArrays.jl/test/runtests.jl
function addtests(fname)
    key = lowercase(splitext(fname)[1])
    if help
        push!(helptests, key)
    else
        if isempty(enabled_tests) || key in enabled_tests
            println("Running $key.jl")
            Random.seed!(42)
            include(fname)
        end
    end
end

@testset verbose=true "All tests" begin
    addtests("complex_powers.jl")
    addtests("indexing.jl")
    addtests("iterators.jl")
    addtests("associated_legendre.jl")
    addtests("wigner_matrices.jl")
    addtests("weights.jl")
    addtests("map2salm.jl")
end

if help
    println()
    println("Pass no args to run all tests, or select one or more of the following:")
    for helptest in helptests
        println("    ", helptest)
    end
end


text:
"""
    `est_marg_model(calls)`

    Function to estimate marginal model from without accounting for error.

    # Examples
    ```julia-repl
    julia>
    ```
"""
function est_marg_model(calls::Vector{Vector{MethCallCpgGrp}})::Vector{Float64}

    # Get total number of CpG sites
    N = length(calls[1])

    # For each call
    pvec = fill(0.0, N)
    nvec = fill(0.0, N)
    for call in calls
        for i = 1:length(call)
            if call[i].obs
                pvec[i] += call[i].log_pyx_m - call[i].log_pyx_u >= 0 ? 1 : 0
                nvec[i] += 1
            end
        end
    end

    # Return phat vector
    return pvec ./ nvec

end
"""
    `comp_mml_marg_model(pvec,𝒜s)`

    Function to compute average MML within analysis regions.

    # Examples
    ```julia-repl
    julia>
    ```
"""
function comp_mml_marg_model(pvec::Vector{Float64}, 𝒜s::Vector{UnitRange{Int64}})::Vector{Float64}

    # Return MML
    return [sum(pvec[x]) / length(x) for x in 𝒜s]

end
"""
    `comp_nme_marg_model(pvec,𝒜s)`

    Function to compute average NME within analysis regions.

    # Examples
    ```julia-repl
    julia>
    ```
"""
function comp_nme_marg_model(pvec::Vector{Float64}, 𝒜s::Vector{UnitRange{Int64}})::Vector{Float64}

    # Compute H(X) = -p⋅log(p) + (1-p)⋅log(p)
    cpg_nme = - 1.0 / (LOG2) * [x * log(max(x, eps(Float64))) + (1.0 - x) * log(max(1.0 - x, eps(Float64))) for x in pvec]

    # Return NME
    return [sum(cpg_nme[x]) / length(x) for x in 𝒜s]

end


text:
using DelimitedFiles
using StatsBase

input = joinpath(@__DIR__, "input")
entries = readlines(input)

struct Entry
    low::Int
    high::Int
    char::Char
    pass::String
end

function v1(entry)
    e = parse_entry(entry)
    pass = e.pass
    freq = countmap([c for c in pass])
    if haskey(freq, e.char)
        return freq[e.char] >= e.low && freq[e.char] <= e.high
    end
    return false
end

function v2(entry)
    e = parse_entry(entry)
    pass = e.pass
    return (pass[e.low] == e.char && pass[e.high] != e.char) || (pass[e.low] != e.char && pass[e.high] == e.char)
end

function parse_entry(entry::String)
    parts = split(entry, " ")
    limits = parts[1]
    low = parse(Int, split(limits, "-")[1])
    high = parse(Int, split(limits, "-")[2])
    rule = parts[2]
    char = rule[1]
    pass = parts[3]

    return Entry(low, high, char, pass)
end

p1 = count(v1, entries)
p2 = count(v2, entries)

@assert(p1 == 383)
@assert(p2 == 272)

println("-----------------------------------------------------------------------")
println("password philosophy -- part one :: $p1")
println("password philosophy -- part two :: $p2")
println("-----------------------------------------------------------------------")


text:
# TODO for all users, may need to distinguish between `atom` and `atomic`
# TODO predicate `var` is needed for general
isvar(t) = false
isvar(::Var) = true
atomic(t) = true
atomic(::Compound) = false

bothvar(a, b) = false
bothvar(a::Var, b::Var) = true
hasvar(v::Vector{<:Term}) = reduce(&, isvar.(v))

"""
functor(Term, Functor, Arity)
functor(f(a, b), f, 2)
"""
function functor(term; options...)
    length(term.args) == 3 || return false
    term, func, arity = term.args
    functor(term, func, arity; options...)
end

function functor(term::Var, func, arity; options...)
    hasvar([func, arity]) && return false
    ins = Compound(func.name, [Var(gensym()) for i in 1:arity.name])
    unifies(term, ins; options...)
end

function functor(term::Const, func, arity; options...)
    if bothvar(func, arity)
        unified = unifies(func, arity; options...)
        unified && return unifies(Const(0), arity; options...)
        return false
    elseif func isa Var
        arity.name == 0 || return false
        return unifies(term, func; options...)
    elseif arity isa Var
        return unifies(Const(0), arity; options...)
    else
        return term == func && arity.name == 0
    end
    false
end

function functor(term::Compound, func, arity; options...)
    if bothvar(func, arity)
        if is_term_list(term)
            li = to_term_list(term)
            unified = unifies(Const(:cons), func; options...)
            unified && return unified(Const(length(li)), arity; options...)
            return false
        else
            unified = unifies(Const(term.name), func; options...)
            unified || return false
            return unifies(Const(length(term.args)), arity; options...)
        end
    elseif func isa Var
        length(term.args) == arity.name || return false
        return unifies(Const(term.name), func; options...)
    elseif arity isa Var
        (func isa Const && term.name == func.name) || return false
        return unifies(Const(length(term.args)), arity; options...)
    else
        (func isa Const && term.name == func.name) || return false
        arity.name == length(term.args) || return false
        return true
    end
    false
end

"""
arg(N, Term, Arg)
arg(2, foo(a, b), b).
"""
function arg(term; options...)
    length(term.args) == 3 || return false
    n, vterm, varg = term.args
    hasvar([n, vterm]) && return false
    (n isa Const && n.name isa Integer) || return false
    arg(n.name, vterm, varg; options...)
end

arg(a, b, c; options...) = false
function arg(n::Int, vterm::Compound, varg; options...)
    if is_term_list(vterm)
        li = to_julia_list(vterm)
        n > 0 || return false
        res = n > length(li) ? Compound(:cend, []) : to_term_list(li[n:end])
        isvar(varg) && return unifies(res, varg; options...)
        return varg == res
    else
        0 < n <= length(vterm.args) || return false
        res = @inbounds vterm.args[n]
        isvar(varg) && return unifies(res, varg; options)
        return res == varg
    end
    false
end

"""
univ(Term, List)
univ(foo(hello, X), List) |- List = [foo, hello, X]
univ(Term, [baz, foo(1)]) |- Term = baz(foo(1))
"""
univ(a, b; options...) = false

# univ(foo(hello, X), List) |- List = [foo, hello, X]
function univ(term::Compound, list::Var; options...)
    head = Const(term.name)
    length(term.args) == 0 && return false
    li = [head; term.args] |> to_term_list
    return unifies(li, list; options...)
end

# univ(Term, [baz, foo(1)]) |- Term = baz(foo(1))
function univ(term::Var, list::Compound; options...)
    list = to_julia_list(list)
    length(list) >= 2 || return
    func = Compound(first(list).name, list[2:end])
    return unifies(func, term; options...)
end


text:
using Test
using LogicCircuits
using Random: bitrand, rand
using DataFrames: DataFrame

include("helper/gpu.jl")

@testset "Binary flows test" begin

    r = fully_factorized_circuit(PlainLogicCircuit, 10)
    input = [1 0 1 0 1 0 1 0 1 0;
            1 1 1 1 1 1 1 1 1 1;
            0 0 0 0 0 0 0 0 0 0;
            0 1 1 0 1 0 0 1 0 1]

    input = DataFrame(BitArray(input))
    @test r(input) == BitVector([1,1,1,1])

    vtree = PlainVtree(10, :balanced)
    r = fully_factorized_circuit(StructLogicCircuit, vtree)
    @test r(input) == BitVector([1,1,1,1])

    v, f = satisfies_flows(r, input)
    foreach(literal_nodes(r)) do n
        id = n.data.node_id
        @test v[:,id] == f[:,id] # invariant of logically valid circuits
    end

    cpu_gpu_agree(input) do d
        satisfies_flows(r, d)[1] # same value
    end

    cpu_gpu_agree(input) do d
        satisfies_flows(r, d)[2][3:end,:] # same flows (except constants)
    end

    num_vars = 7
    mgr = SddMgr(num_vars, :balanced)
    v = Dict([(i => compile(mgr, Lit(i))) for i=1:num_vars])
    c = (v[1] | !v[2] | v[3]) &
        (v[2] | !v[7] | v[6]) &
        (v[3] | !v[4] | v[5]) &
        (v[1] | !v[4] | v[6])
    input = DataFrame(bitrand(25,num_vars))

    r = smooth(PlainLogicCircuit(c)) # flows don't make sense unless the circuit is smooth; cannot smooth trimmed SDDs

    v, f = satisfies_flows(r, input)
    foreach(literal_nodes(r)) do n
        id = n.data.node_id
        @test v[:,id] .& v[:,end] == f[:,id] # invariant of all circuits
    end

    cpu_gpu_agree(input) do d
        satisfies_flows(r, d)[1] # same value
    end

    cpu_gpu_agree(input) do d
        satisfies_flows(r, d)[2][3:end,:] # same flows (except constants)
    end

    l1 = LogicCircuit(Lit(1))
    l2 = LogicCircuit(Lit(2))
    l3 = LogicCircuit(Lit(-1))
    l4 = LogicCircuit(Lit(-2))
    r = (l1 & l2) | (l3 & l4)
    input = DataFrame(bitrand(4,2))

    v, f = satisfies_flows(r, input)
    foreach(literal_nodes(r)) do n
        id = n.data.node_id
        @test v[:,id] .& v[:,end] == f[:,id] # invariant of all circuits
    end

    cpu_gpu_agree(input) do d
        satisfies_flows(r, d)[1] # same value
    end

    cpu_gpu_agree(input) do d
        satisfies_flows(r, d)[2][3:end,:] # same flows (except constants)
    end

end

@testset "Weighted binary flows test" begin

    r = fully_factorized_circuit(PlainLogicCircuit, 10)
    input = [1 0 1 0 1 0 1 0 1 0;
            1 1 1 1 1 1 1 1 1 1;
            0 0 0 0 0 0 0 0 0 0;
            0 1 1 0 1 0 0 1 0 1]

    input = DataFrame(BitArray(input))
    @test r(input) == BitVector([1,1,1,1])

    weights = [0.6, 0.6, 0.6, 0.6]
    weights = Array{Float32, 1}(weights)

    vtree = PlainVtree(10, :balanced)
    r = fully_factorized_circuit(StructLogicCircuit, vtree)
    @test r(input) == BitVector([1,1,1,1])

    v, f = satisfies_flows(r, input; weights = weights)
    foreach(literal_nodes(r)) do n
        id = n.data.node_id
        @test v[:,id] == f[:,id] # invariant of logically valid circuits
    end

    if CUDA.functional()
        @test all(satisfies_flows(r, input; weights = weights)[1] .≈ to_cpu(satisfies_flows(r, to_gpu(input); weights = to_gpu(weights))[1]))

        @test all(satisfies_flows(r, input; weights)[2][3:end,:] .≈ to_cpu(satisfies_flows(r, to_gpu(input); weights = to_gpu(weights))[2][3:end,:]))
    end

    num_vars = 7
    mgr = SddMgr(num_vars, :balanced)
    v = Dict([(i => compile(mgr, Lit(i))) for i=1:num_vars])
    c = (v[1] | !v[2] | v[3]) &
        (v[2] | !v[7] | v[6]) &
        (v[3] | !v[4] | v[5]) &
        (v[1] | !v[4] | v[6])
    input = DataFrame(bitrand(25,num_vars))

    weights = Array{Float32, 1}(ones(25) * 0.6)

    r = smooth(PlainLogicCircuit(c)) # flows don't make sense unless the circuit is smooth; cannot smooth trimmed SDDs

    v, f = satisfies_flows(r, input; weights = weights)
    foreach(literal_nodes(r)) do n
        id = n.data.node_id
        @test v[:,id] .& v[:,end] == f[:,id] # invariant of all circuits
    end

    if CUDA.functional()
        @test all(satisfies_flows(r, input; weights = weights)[1] .≈ to_cpu(satisfies_flows(r, to_gpu(input); weights = to_gpu(weights))[1]))

        @test all(satisfies_flows(r, input; weights)[2][3:end,:] .≈ to_cpu(satisfies_flows(r, to_gpu(input); weights = to_gpu(weights))[2][3:end,:]))
    end

    l1 = LogicCircuit(Lit(1))
    l2 = LogicCircuit(Lit(2))
    l3 = LogicCircuit(Lit(-1))
    l4 = LogicCircuit(Lit(-2))
    r = (l1 & l2) | (l3 & l4)
    input = DataFrame(bitrand(4,2))

    v, f = satisfies_flows(r, input; weights = weights)
    foreach(literal_nodes(r)) do n
        id = n.data.node_id
        @test v[:,id] .& v[:,end] == f[:,id] # invariant of all circuits
    end

    if CUDA.functional()
        @test all(satisfies_flows(r, input; weights = weights)[1] .≈ to_cpu(satisfies_flows(r, to_gpu(input); weights = to_gpu(weights))[1]))

        @test all(satisfies_flows(r, input; weights = weights)[2][3:end,:] .≈ to_cpu(satisfies_flows(r, to_gpu(input); weights = to_gpu(weights))[2][3:end,:]))
    end

end

@testset "Soft Flows test" begin
    l_a = LogicCircuit(Lit(1))
    l_na = LogicCircuit(Lit(-1))
    l_b = LogicCircuit(Lit(2))
    l_nb = LogicCircuit(Lit(-2))

    o_c = (l_a & l_b) | (l_na & l_b)
    o_d = (l_a & l_b) | (l_a & l_nb)
    r = o_c | o_d

    weights = [0.6, 0.6, 0.6]
    weights = Array{Float32, 1}(weights)

    df = DataFrame(BitMatrix([true true; false true; false false]))
    sdf = soften(df, 0.001; scale_by_marginal = false)

    v, f = satisfies_flows(r, sdf)
    @test all(v[:,o_c.data.node_id] .≈ sdf[:, 2])
    @test all(v[:,o_c.data.node_id] .≈ f[:,o_c.data.node_id])
    @test all(v[:,o_d.data.node_id] .≈ sdf[:, 1])
    @test all(v[:,o_d.data.node_id] .≈ f[:,o_d.data.node_id])

    v, f = satisfies_flows(r, sdf; weights = weights)
    @test all(v[:,o_c.data.node_id] .≈ sdf[:, 2])
    @test all(v[:,o_c.data.node_id] .≈ f[:,o_c.data.node_id])
    @test all(v[:,o_d.data.node_id] .≈ sdf[:, 1])
    @test all(v[:,o_d.data.node_id] .≈ f[:,o_d.data.node_id])

    if CUDA.functional()
        @test all(satisfies_flows(r, sdf; weights = weights)[1] .≈ to_cpu(satisfies_flows(r, to_gpu(sdf); weights = to_gpu(weights))[1]))
    end
end

@testset "Probabilistic Flows test" begin

    r = fully_factorized_circuit(PlainLogicCircuit, 10)

    input = DataFrame(Float64[1 0 1 0 1 0 1 0 1 0;
                    1 1 1 1 1 1 1 1 1 1;
                    0 0 0 0 0 0 0 0 0 0;
                    0 1 1 0 1 0 0 1 0 1])

    @test r(input) ≈ [1.0, 1.0, 1.0, 1.0]

    v, f = satisfies_flows(r, input)
    foreach(literal_nodes(r)) do n
        id = n.data.node_id
        @test v[:,id] ≈ f[:,id] # invariant of logically valid circuits
    end

    cpu_gpu_agree_approx(input) do d
        satisfies_flows(r, d)[1] # same value
    end

    cpu_gpu_agree_approx(input) do d
        satisfies_flows(r, d)[2][3:end,:] # same flows (except constants)
    end

    l_a = LogicCircuit(Lit(1))
    l_na = LogicCircuit(Lit(-1))
    l_b = LogicCircuit(Lit(2))
    l_nb = LogicCircuit(Lit(-2))
    l_c = LogicCircuit(Lit(3))
    l_nc = LogicCircuit(Lit(-3))

    o_c = (l_a & l_b & l_c) | (l_a & l_nb & l_c) | (l_na & l_b & l_c) | (l_na & l_nb & l_c)
    o_nc = (l_a & l_b & l_nc) | (l_a & l_nb & l_nc) | (l_na & l_b & l_nc) | (l_na & l_nb & l_nc)
    r = (o_c | o_nc)

    input = DataFrame(rand(Float64, (10,3)))

    @test all(r(input) .≈ 1.0)

    v, f = satisfies_flows(r, input)
    @test all(v[:,end] .≈ 1.0)
    @test all(f[:,end].≈ 1.0)
    @test all(v[:,o_c.data.node_id] .≈ input[:, 3])

    @test all(f[:,o_c.data.node_id] .≈ input[:, 3])
    @test all(v[:,o_nc.data.node_id] .+ input[:, 3] .≈ 1.0)
    @test all(f[:,o_nc.data.node_id] .+ input[:, 3] .≈ 1.0)

    foreach(literal_nodes(r)) do n
        id = n.data.node_id
        @test v[:,id] ≈ f[:,id] # invariant of logically valid circuits
    end

    cpu_gpu_agree_approx(input) do d
        satisfies_flows(r, d)[1] # same value
    end

    cpu_gpu_agree_approx(input) do d
        satisfies_flows(r, d)[2][3:end,:] # same flows (except constants)
    end

    o_c = (l_a & l_b & l_c) | (l_na & l_b & l_c)
    o_nc = (l_a & l_nb & l_nc) | (l_na & l_nb & l_nc)
    r = (o_c | o_nc)

    input = DataFrame(rand(Float64, (10,3)))
    input[1:5, 3] .= 0.0

    v, f = satisfies_flows(r, input)
    @test all(v[:,o_c.data.node_id] .≈ input[:, 2] .* input[:, 3])
    @test all(f[:,o_c.data.node_id] .≈ v[:,o_c.data.node_id])
    @test all(v[:,o_nc.data.node_id] .≈ (1.0 .- input[:, 2]) .* (1.0 .- input[:, 3]))
    @test all(f[:,o_nc.data.node_id] .≈ v[:,o_nc.data.node_id])
    @test all(f[:,l_a.data.node_id] .≈ v[:,l_a.data.node_id] .* f[:,r.data.node_id])
    @test all(f[:,l_na.data.node_id] .≈ v[:,l_na.data.node_id] .* f[:,r.data.node_id])
    @test all(f[:,l_b.data.node_id] .≈ input[:, 2] .* input[:, 3])
    @test all(f[:,l_b.data.node_id] .≈ f[:,l_c.data.node_id])
    @test all(f[:,l_nb.data.node_id] .≈ (1.0 .- input[:, 2]) .* (1.0 .- input[:, 3]))
    @test all(f[:,l_nb.data.node_id] .≈ f[:,l_nc.data.node_id])

    cpu_gpu_agree_approx(input) do d
        satisfies_flows(r, d)[1] # same value
    end

    cpu_gpu_agree_approx(input) do d
        satisfies_flows(r, d)[2][3:end,:] # same flows (except constants)
    end

end


text:
module GridWorlds

import Requires
import Crayons
import MacroTools:@forward
import Random
import DataStructures
const DS = DataStructures
import ReinforcementLearningBase
import ReinforcementLearningBase:RLBase
import StaticArrays
const SA = StaticArrays

const GW = GridWorlds

include("directions.jl")
include("actions.jl")
include("objects.jl")
include("grid_world_base.jl")
include("abstract_grid_world.jl")
include("envs/envs.jl")
include("textual_rendering.jl")

function __init__()
    Requires.@require Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("graphical_rendering.jl")
end

end


text:
# MIT License
# Copyright (c) 2017: Xavier Gandibleux, Anthony Przybylski, Gauthier Soleilhac, and contributors.
mutable struct mono_problem
    pb::problem
    psize::Int
    p::Vector{Int} #combined profit vector
    p1::Vector{Int}
    p2::Vector{Int}
    w::Vector{Int}
    c::Int
    C0::Vector{Int} #Variables that have been fixed to 0 for the mono_problem only
    C1::Vector{Int} #Variables that have been fixed to 1 for the mono_problem only
    ω::Int #Minimum weight
    min_profit_λ::Int
    min_profit_1::Int
    min_profit_2::Int
    eff::Vector{Float64}
    variables::Vector{Int} #Order of variables with decreasing utiliy wrt to the combined objective
    ub_z1_i0::Vector{Int} #Upper bound on z1 for each variable i fixed to 0
    ub_z1_i1::Vector{Int} #Upper bound on z1 for each variable i fixed to 1
    ub_z2_i0::Vector{Int} #Upper bound on z2 for each variable i fixed to 0
    ub_z2_i1::Vector{Int} #Upper bound on z2 for each variable i fixed to 1
end

#Printer
Base.show(io::IO, pb::mono_problem) =
    print(io, "Mono-objective Knapsack Problem with $(size(pb))($(pb.psize)) variables :
\tω=$(pb.ω)\n\tc=$(pb.c)\n\tvars=$(pb.variables)\n\tp=$(pb.p[pb.variables])\n\tw=$(pb.w[pb.variables])")

size(pb::mono_problem)::Int = length(pb.variables)
variables(pb::mono_problem)::Vector{Int} = pb.variables

function fix!(mp::mono_problem, C0::Vector{Int}, C1::Vector{Int})
    for i in C0
        v = mp.variables[i]
        push!(mp.C0, v)
    end
    for i in C1
        v = mp.variables[i]
        push!(mp.C1, v)
        mp.min_profit_λ += mp.p[v]
        mp.min_profit_1 += mp.p1[v]
        mp.min_profit_2 += mp.p2[v]
        mp.ω += mp.w[v]
    end
    deleteat!(mp.variables, sort(vcat(C0,C1)))
end

#Constructor from a bi_problem
function mono_problem(p::problem, λ1::Int, λ2::Int, calculate_efficiency::Bool = true)

    profits = λ1*p.p1 + λ2*p.p2
    eff = calculate_efficiency ? profits ./ p.w : Float64[]
    vars = variables(p)[:]
    calculate_efficiency && sort!(vars, by = x -> eff[x], rev=true)

    return mono_problem(
        p,
        size(p), #Number of variables
        profits, #Combined profit vector
        p.p1,
        p.p2,
        p.w, #Weight vector
        p.c, #Capacity
        Int[], #C0
        Int[], #C1
        p.ω, #ω
        λ1*p.min_profit_1 + λ2*p.min_profit_2,#min_profit
        p.min_profit_1,
        p.min_profit_2,
        eff, #Utility vector
        vars, #variables
        p.ub_z1_i0,
        p.ub_z1_i1,
        p.ub_z2_i0,
        p.ub_z2_i1
        )
end

function reduce!(mp::mono_problem, Δ, output::Bool)

    #Here, C0 and C1 will contain not the variables to be fixed to 0/1 but their indices in mp.variables
    #It makes their removal faster : we don't have to use findfirst(...) to find the values we have to remove from mp.variables
    lb_zλ, lb_z1, lb_z2 = Δ.lb, obj_1(Δ.xr), obj_2(Δ.xs)
    vars = variables(mp)

    ### Reduction with lb_z1
    C0,C1 = Vector{Int}(), Vector{Int}()
    for i = 1:size(mp)
        if mp.ub_z1_i0[vars[i]] <= lb_z1
            push!(C1, i)
        elseif mp.ub_z1_i1[vars[i]] <= lb_z1
            push!(C0, i)
        end
    end
    ##@assert isempty(intersect(C1,C0))
    # @show C0
    # @show C1
    fix!(mp, C0, C1)
    empty!(C0) ; empty!(C1)

    #Reduction with lb_z2
    for i = 1:size(mp)
        if mp.ub_z2_i0[vars[i]] <= lb_z2
            push!(C1, i)
        end
        if mp.ub_z2_i1[vars[i]] <= lb_z2
            push!(C0, i)
        end
    end
    ##@assert isempty(intersect(C1,C0))
    fix!(mp, C0, C1)
    empty!(C0) ; empty!(C1)

    #TODO : find the break_item and only calculate C0 for variables before the break_item and C1 for those after.

    #Reduction with lb_zλ
    for i=1:size(mp)
        #Calculate relaxation with variable i set to 0
        w = mp.ω
        ub = mp.min_profit_λ
        j = 1
        while j <= size(mp) && (i==j || w + mp.w[vars[j]] <= mp.c)
            v = vars[j]
            if j != i
                w += mp.w[v]
                ub += mp.p[v]
            end
            j += 1
        end

        if j <= size(mp)
            v = vars[j]
            cleft = mp.c - w
            if j+1 == i || j == size(mp)
                if j+2 <= size(mp)
                    U0 = ub + floor(Int, cleft * mp.p[vars[j+2]]/mp.w[vars[j+2]])
                else
                    U0 = ub
                end
            else
                U0 = ub + floor(Int, cleft * mp.p[vars[j+1]]/mp.w[vars[j+1]])
            end
            if j-1 == i || j == 1
                if j-2 >= 1
                    U1 = ub + floor(Int, mp.p[v] - (mp.w[v] - cleft)*(mp.p[vars[j-2]]/mp.w[vars[j-2]]))
                else
                    U1 = 0
                end
            else
                U1 = ub + floor(Int, mp.p[v] - (mp.w[v] - cleft)*(mp.p[vars[j-1]]/mp.w[vars[j-1]]))
            end
        ub = max(U0,U1)
        end
        if ub < lb_zλ
            push!(C1, i)
        end

    end

    for i=1:size(mp)
        #Calculate relaxation with variable i set to 1
        w = mp.ω + mp.w[vars[i]]
        ub = mp.min_profit_λ + mp.p[vars[i]]
        j = 1
        while j <= size(mp) && (i==j || w + mp.w[vars[j]] <= mp.c)
            v = vars[j]
            if j != i
                w += mp.w[v]
                ub += mp.p[v]
            end
            j += 1
        end

        if j <= size(mp)
            v = vars[j]
            cleft = mp.c - w
            if j+1 == i || j == size(mp)
                if j+2 <= size(mp)
                    U0 = ub + floor(Int, cleft * mp.p[vars[j+2]]/mp.w[vars[j+2]])
                else
                    U0 = ub
                end
            else
                U0 = ub + floor(Int, cleft * mp.p[vars[j+1]]/mp.w[vars[j+1]])
            end
            if j-1==i || j == 1
                if j-2 >= 1
                    U1 = ub + floor(Int, mp.p[v] - (mp.w[v] - cleft)*(mp.p[vars[j-2]]/mp.w[vars[j-2]]))
                else
                    U1 = 0
                end
            else
                U1 = ub + floor(Int, mp.p[v] - (mp.w[v] - cleft)*(mp.p[vars[j-1]]/mp.w[vars[j-1]]))
            end
        ub = max(U0,U1)
        end
        if ub < lb_zλ
            push!(C0, i)
        end
    end

    fix!(mp, C0, C1)
    reduce_again!(mp, Δ)
    output && mp.psize != size(mp) && println("Problem reduced from $(mp.psize) to $(size(mp)) variables")

    nothing
end

function reduce_again!(mp::mono_problem, Δ)

    size0 = size(mp)

    lb_zλ = Δ.lb
    vars = variables(mp)
    C0,C1 = Vector{Int}(), Vector{Int}()
    M = 500000

    for i = 1:length(vars)
        v = vars[i]
        p_old = mp.p[v] ; mp.p[v] = 0 #fix_0!(p, v)
        ms = solve_mono(mp)
        if obj(ms) < lb_zλ
            push!(C1, i)
        end
        mp.p[v] = p_old
    end

    for i = 1:length(vars)
        v = vars[i]
        if mp.ω + mp.w[v] > mp.c
            push!(C0, i)
            continue
        end

        mp.p[v] += M ; mp.min_profit_λ -= M # fix_1!(mp, v)
        ms = solve_mono(mp)
        if obj(ms) < lb_zλ
            push!(C0, i)
        end
        mp.p[v] -= M ; mp.min_profit_λ += M #unfix_0
    end

    fix!(mp, C0, C1)

    return nothing
end

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
struct LogBook{T<:AbstractFloat}              # provide the number format as T
    logbook::Array{UInt64,1}
end

# in case no type T is provided use Float16
LogBook(lb::Array{UInt64,1}) = LogBook{Float16}(lb)

const n16 = 2^16
const nlogbooks = 32
const logbooks = [LogBook(zeros(UInt64,n16)) for _ in 1:nlogbooks]

"""Change the histogram array of logbook #id with UInt64 array x."""
function set_logbook(x::Array{UInt64,1},id::Int)
    if length(x) == n16
        logbooks[id].logbook .= x
    else
        error("Array of length $(length(x)) instead of $n16.")
    end
    return nothing
end

"""Change the type T of logbook #id. Must be a 16bit type."""
function set_logbook(T::DataType,id::Int=1)
    if sizeof(T) == 2   # must be a 16bit type
        # copy over the logbook but change type to T
        logbooks[id] = LogBook{T}(logbooks[id].logbook)
    else
        error("$(sizeof(T)*8)bit type $T was provided but must be a 16bit type.")
    end
    return nothing
end

get_logbooks() = logbooks
get_logbook(id::Int=1) = logbooks[id]
reset_logbook(id::Int=1) = set_logbook(zeros(UInt64,n16),id)
reset_logbooks() = for id in 1:nlogbooks reset_logbook(id) end
set_logbooks(T::DataType) = for id in 1:nlogbooks set_logbook(T,id) end

"""Count arithmetic result in logbook."""
function log_it(::Type{T},x::AbstractFloat,id::Int) where {T<:AbstractFloat}
    @inbounds logbooks[id].logbook[reinterpret(UInt16,T(x)) + 1] += 1
end

Base.sum(lb::LogBook)::Int64 = sum(lb.logbook)
entropy(lb::LogBook,b::Real) = entropy(lb.logbook/sum(lb),b)
entropy(lb::LogBook) = entropy(lb.logbook/sum(lb))
Base.length(lb::LogBook) = length(lb.logbook)
Base.getindex(lb::LogBook,i::Integer) = Int(getindex(lb.logbook,i))
Base.getindex(lb::LogBook,u::UnitRange) = Int.(getindex(lb.logbook,u))
Base.getindex(lb::LogBook,s::StepRange) = Int.(getindex(lb.logbook,s))
Base.lastindex(lb::LogBook) = Int(lb.logbook[end])

function Base.show(io::IO,lb::LogBook{T}) where T
    n = length(lb)
    if n > 10
        print(io,"$n-element LogBook{$T}(")
        [print(io,string(Int(i)),", ") for i in lb.logbook[1:5]]
        print(io,"… , ")
        [print(io,string(Int(i)),", ") for i in lb.logbook[end-5:end-1]]
        print(io,string(Int(lb.logbook[end])),")")
    else
        print(io,"$n-element LogBook{$T}(")
        [print(io,string(Int(i)),", ") for i in lb.logbook[1:end-1]]
        print(io,string(Int(lb.logbook[end])),")")
    end
end

"""Find the largest number that was logged into the logbook `lb`."""
function Base.maximum(lb::LogBook{T}) where T
    n_half = length(lb)÷2
    i = n_half - findfirst(x->x>0,lb[n_half:-1:1])
    return reinterpret(T,UInt16(i))
end

"""Find the smallest number that was logged into the logbook `lb`."""
function Base.minimum(lb::LogBook{T}) where T
    n = length(lb)
    i = n - findfirst(x->x>0,lb[n:-1:1])
    return reinterpret(T,UInt16(i))
end

"""Sort a logbook as [-NaNs,-∞,...,-0,0,...,∞,NaNs]"""
function Base.sort(lb::LogBook)
    @views positives = lb.logbook[1:2^15]
    @views negatives = lb.logbook[2^15+1:end]
    lb_sorted = vcat(reverse(negatives),positives)
    return LogBook(lb_sorted)
end

"""Convert a Float16/BFloat16 value to the corresponding UInt16 in a sorted logbook"""
function uint16_sort(x::AbstractFloat)
    ui = reinterpret(UInt16,x)      # will trigger an error for non-16bit floats
    return signbit(x) ? 0x8000 - (ui - 0x7fff) : ui + 0x8000
end


text:
# Autogenerated wrapper script for opencv_jll for x86_64-linux-gnu
export OpenCV_jl, libopencv_calib3d, libopencv_core, libopencv_dnn, libopencv_features2d, libopencv_flann, libopencv_highgui, libopencv_imgcodecs, libopencv_imgproc, libopencv_julia, libopencv_objdetect, libopencv_videoio

using Qt_jll
using Libglvnd_jll
using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("opencv")
JLLWrappers.@declare_file_product(OpenCV_jl)
JLLWrappers.@declare_library_product(libopencv_calib3d, "libopencv_calib3d.so.4.5")
JLLWrappers.@declare_library_product(libopencv_core, "libopencv_core.so.4.5")
JLLWrappers.@declare_library_product(libopencv_dnn, "libopencv_dnn.so.4.5")
JLLWrappers.@declare_library_product(libopencv_features2d, "libopencv_features2d.so.4.5")
JLLWrappers.@declare_library_product(libopencv_flann, "libopencv_flann.so.4.5")
JLLWrappers.@declare_library_product(libopencv_highgui, "libopencv_highgui.so.4.5")
JLLWrappers.@declare_library_product(libopencv_imgcodecs, "libopencv_imgcodecs.so.4.5")
JLLWrappers.@declare_library_product(libopencv_imgproc, "libopencv_imgproc.so.4.5")
JLLWrappers.@declare_library_product(libopencv_julia, "libopencv_julia.so")
JLLWrappers.@declare_library_product(libopencv_objdetect, "libopencv_objdetect.so.4.5")
JLLWrappers.@declare_library_product(libopencv_videoio, "libopencv_videoio.so.4.5")
function __init__()
    JLLWrappers.@generate_init_header(Qt_jll, Libglvnd_jll, libcxxwrap_julia_jll)
    JLLWrappers.@init_file_product(
        OpenCV_jl,
        "OpenCV.jl.tar",
    )

    JLLWrappers.@init_library_product(
        libopencv_calib3d,
        "lib/libopencv_calib3d.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_core,
        "lib/libopencv_core.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_dnn,
        "lib/libopencv_dnn.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_features2d,
        "lib/libopencv_features2d.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_flann,
        "lib/libopencv_flann.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_highgui,
        "lib/libopencv_highgui.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_imgcodecs,
        "lib/libopencv_imgcodecs.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_imgproc,
        "lib/libopencv_imgproc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_julia,
        "lib/libopencv_julia.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_objdetect,
        "lib/libopencv_objdetect.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopencv_videoio,
        "lib/libopencv_videoio.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
for N in [Float64, Rational{Int}, Float32]
    # approximate rectification by rectifying the vertices
    # - exact approximation
    X = Ball1(N[1, 1], N(1))
    Y = approximate(Rectification(X))
    @test isequivalent(Y, X)
    # - underapproximation
    X = Ball1(N[2, 2], N(3))
    Y = approximate(Rectification(X))
    Z = VPolygon([N[5, 2], N[2, 5], N[2, 0], N[0, 2]])
    @test Y ⊆ Z && Z ⊆ Y
    # - overapproximation
    X = Ball1(N[-1, -1], N(2))
    Y = approximate(Rectification(X))
    Z = VPolygon([N[0, 0], N[0, 1], N[1, 0]])
    @test Y ⊆ Z && Z ⊆ Y
    # - neither over- nor underapproximation
    X = VPolygon([N[3, 2], N[1, -1], N[-1, 2], N[-1, 4]])
    Y = approximate(Rectification(X))
    Z = VPolygon([N[3, 2], N[1, 0], N[0, 2], N[0, 4]])
    @test Y ⊆ Z && Z ⊆ Y
end


text:
module SkeletonCellFieldsTests

using Test
using Gridap
import Gridap: ∇

model = CartesianDiscreteModel(partition=(3,4,2))

order = 1
diritag = "boundary"
fespace = H1ConformingFESpace(Float64,model,order,diritag)

trian = Triangulation(model)

tags = [27,]
strian = SkeletonTriangulation(model,tags)
squad = CellQuadrature(strian,degree=2)
s = coordinates(squad)

ufun(x) = 2.0*x[1]^2 * x[2]^2
grad_ufun(x) = VectorValue(4.0*x[1]*x[2]^2, 4.0*x[1]^2*x[2])
∇(::typeof(ufun)) = grad_ufun

u = CellField(trian,ufun)

uh = interpolate(fespace,ufun)

su = restrict(u,strian)

suh = restrict(uh,strian)

jsu = jump(su)
jsu_s = evaluate(jsu,s)

msu = mean(su)
msu_s = evaluate(msu,s)

jsuh = jump(suh)
jsuh_s = evaluate(jsuh,s)

msuh = mean(suh)
msuh_s = evaluate(msuh,s)

jsu_grad = jump(∇(su))
jsu_grad_s = evaluate(jsu_grad,s)

msu_grad = mean(∇(su))
msu_grad_s = evaluate(msu_grad,s)

jsuh_grad = jump(∇(suh))
jsuh_grad_s = evaluate(jsuh_grad,s)

msuh_grad = mean(∇(suh))
msuh_grad_s = evaluate(msuh_grad,s)

#writevtk(strian,"strian",cellfields=[
# "jumpu"=>jsu,"meanu"=>msu,"jumpu_grad"=>jsu_grad,"meanu_grad"=>msu_grad,
# "jumph"=>jsuh,"meanh"=>msuh,"jumpuh_grad"=>jsuh_grad,"meanh_grad"=>msuh_grad ])

v = CellBasis(trian)
du = v

u_Γ = restrict(u,strian)
du_Γ = restrict(v,strian)
v_Γ = restrict(v,strian)

cm = inner(jump(v_Γ),jump(u_Γ))

cv = integrate(cm,strian,squad)
_ = collect(cv.cellvector1)
_ = collect(cv.cellvector2)

cm = inner(mean(v_Γ),jump(u_Γ))

cv = integrate(cm,strian,squad)
_ = collect(cv.cellvector1)
_ = collect(cv.cellvector2)

cm = inner(jump(v_Γ),jump(du_Γ))

cv = integrate(cm,strian,squad)
_ = collect(cv.cellmatrix11)
_ = collect(cv.cellmatrix12)
_ = collect(cv.cellmatrix21)
_ = collect(cv.cellmatrix22)

cm = inner(mean(v_Γ),jump(du_Γ))

cv = integrate(cm,strian,squad)
_ = collect(cv.cellmatrix11)
_ = collect(cv.cellmatrix12)
_ = collect(cv.cellmatrix21)
_ = collect(cv.cellmatrix22)

cm = inner(jump(u_Γ),jump(u_Γ))
cv = integrate(cm,strian,squad)
_ = collect(cv)

phi = CellGeomap(trian)
phi_Γ = restrict(phi,strian)

x = evaluate(jump(phi_Γ),s)
z = Point(0.0,0.0,0.0)
zs = fill(z,4)
@test all([ xi .+ 1 ≈ zs .+ 1 for xi in x ])

x = evaluate(mean(phi_Γ),s)
x2 = evaluate(phi_Γ.cellfield1,s)
@test x ≈ x2

end # module


text:
# Unit Test for `SolveQuadraticProgram()`

using LinearAlgebra
using SparseArrays

@enum ProblemClass randomQp = 1 inequalityConstrainedQp equalityConstrainedQp optimalControl portfolioOptimization lassoOptimization huberFitting supportVectorMachine isotonicRegression

function GenerateRandomQP(problemClass::ProblemClass, numElements::Int64 = 1000; numConstraints::Int64 = 0)

    if (problemClass ∈ (randomQp, inequalityConstrainedQp, equalityConstrainedQp, optimalControl))
        densityFctr = 0.15;
        α           = 1e-2;

        mM = sprandn(numElements, numElements, densityFctr);
        mP = (mM' * mM) + sparse(α * I, numElements, numElements);
        vQ = randn(numElements);
        if (problemClass == inequalityConstrainedQp)
            numConstraints ≠ 0 || (numConstraints = 10 * numElements);
            mA = sprandn(numConstraints, numElements, densityFctr);
            vL = -rand(numConstraints);
            vU = rand(numConstraints);
        elseif (problemClass == equalityConstrainedQp)
            numConstraints ≠ 0 || (numConstraints = numElements ÷ 2);
            mA = sprandn(numConstraints, numElements, densityFctr);
            vL = randn(numConstraints);
            vU = copy(vL);
        else
            numConstraints ≠ 0 || (numConstraints = numElements ÷ 2);
            mA = sprandn(numConstraints, numElements, densityFctr);
            vL = -rand(numConstraints);
            vU = rand(numConstraints);
            vI = rand(numConstraints) .<= 0.15;
            vL[vI] .= vU[vI];
            vI = rand(numConstraints) .<= 0.15;
            vU[vI] .= vI[vI];
        end
    elseif (problemClass == portfolioOptimization)
        densityFctr = 0.5;

        numConstraints ≠ 0 || (numConstraints = max(5, numElements ÷ 100));
        mD = spdiagm(rand(numElements) * sqrt(numConstraints));
        mP = [mD spzeros(numElements, numConstraints); spzeros(numConstraints, numElements) sparse(I, numConstraints, numConstraints)];
        vQ = [randn(numElements); zeros(numConstraints)];
        mF = sprandn(numElements, numConstraints, densityFctr);
        mA = [mF' sparse(-I, numConstraints, numConstraints); ones(1, numElements) spzeros(1, numConstraints); sparse(I, numElements, numElements) spzeros(numElements, numConstraints)];
        vL = [zeros(numConstraints); 1; zeros(numElements)];
        vU = [zeros(numConstraints); 1; ones(numElements)];
    elseif (problemClass == lassoOptimization)
        densityFctr = 0.15;

        numConstraints ≠ 0 || (numConstraints = numElements * 100);
        mAd = sprandn(numConstraints, numElements, densityFctr);
        vXX = (randn(numElements) ./ sqrt(numElements)) .* (rand(numElements) .> 0.5);
        vB  = mAd * vXX + randn(numConstraints);
        λ   = norm(mAd' * vB, Inf) / 5.0;

        mP = blockdiag(spzeros(numElements, numElements), sparse(2I, numConstraints, numConstraints), spzeros(numElements, numElements));
        vQ = [zeros(numElements + numConstraints); λ .* ones(numElements)];
        mA = [mAd sparse(-I, numConstraints, numConstraints) spzeros(numConstraints, numElements); sparse(I, numElements, numElements) spzeros(numElements, numConstraints) sparse(-I, numElements, numElements); sparse(I, numElements, numElements) spzeros(numElements, numConstraints) sparse(I, numElements, numElements)];
        vL = [vB; -Inf .* ones(numElements); zeros(numElements)];
        vU = [vB; zeros(numElements); Inf .* ones(numElements)];
    elseif (problemClass == huberFitting)
        densityFctr = 0.15;

        numConstraints ≠ 0 || (numConstraints = numElements * 100);
        mAd = sprandn(numConstraints, numElements, densityFctr);
        vXX = randn(numElements) ./ sqrt(numElements);
        vI = rand(numConstraints) .< 0.95;
        vB = (mAd * vXX) + (0.5 .* vI .* randn(numConstraints)) + (10 .* .!vI .* rand(numConstraints));

        mP = blockdiag(spzeros(numElements, numElements), sparse(2I, numConstraints, numConstraints), spzeros(2 * numConstraints, 2 * numConstraints));
        vQ = [zeros(numElements + numConstraints); 2 .* ones(2 * numConstraints)];
        mIm = sparse(I, numConstraints, numConstraints);
        mA = [mAd -mIm -mIm mIm; spzeros(numConstraints, numElements + numConstraints) mIm spzeros(numConstraints, numConstraints); spzeros(numConstraints, numElements + numConstraints + numConstraints) mIm];
        vL = [vB; zeros(2 * numConstraints)];
        vU = [vB; Inf .* ones(2 * numConstraints)];
    elseif (problemClass == supportVectorMachine)
        densityFctr = 0.15;

        (numConstraints ≠ 0) || (numConstraints = numElements * 100);
        numClassA       = numConstraints ÷ 2;
        λ               = 1;
        vB = [ones(numClassA); -ones(numClassA)];
        mAu = sprandn(numClassA, numElements, densityFctr);
        mAl = sprandn(numClassA, numElements, densityFctr);
        mAd = [(mAu ./ sqrt(numConstraints)) + ((mAu .≠ 0)./ numConstraints); (mAl ./ sqrt(numConstraints)) - ((mAl .≠ 0) ./ numConstraints)];

        mP = blockdiag(sparse(2I, numElements, numElements), spzeros(numConstraints, numConstraints));
        vQ = λ .* [zeros(numElements); ones(numConstraints)];
        mA = [(spdiagm(0 => vB) * mAd) sparse(-I, numConstraints, numConstraints); spzeros(numConstraints, numElements) sparse(I, numConstraints, numConstraints)];
        vL = [-Inf .* ones(numConstraints); zeros(numConstraints)];
        vU = [-ones(numConstraints); Inf .* ones(numConstraints)];
    elseif (problemClass == isotonicRegression)
        densityFctr = 0.25;
        α           = 1e-2;

        mM = sprandn(numElements, numElements, densityFctr);
        mP = (mM' * mM) + sparse(α * I, numElements, numElements);
        vQ = randn(numElements);

        if (rand() >= 0.5)
            # Monotonic Non Increasing
            mA = spdiagm(numElements - 1, numElements, 0 => ones(numElements - 1), 1 => -ones(numElements - 1));
        else
            # Monotonic Non Decreasing
            mA = spdiagm(numElements - 1, numElements, 0 => -ones(numElements - 1), 1 => ones(numElements - 1));
        end
        vL = zeros(numElements - 1);
        vU = 10 * ones(numElements - 1);
    end

    return mP, vQ, mA, vL, vU;

end

max_stars_repo_name: RoyiAvital/QuadraticProgramSolver

max_issues_repo_name: RoyiAvital/QuadraticProgramSolver

max_forks_repo_name: RoyiAvital/QuadraticProgramSolver


text:
export solveDDSerial,setupDDSerial#,getSubMeshOfCell;#,solveDDJac,solveDDRB

function computeResidualAtIdx(AT::SparseMatrixCSC,b::Array,x::Array,Idxs::Array{DDIndType})
r = b[Idxs];
for t = 1:length(Idxs)
	for tt = AT.colptr[Idxs[t]]:AT.colptr[Idxs[t]+1]-1
		r[t] = r[t] - conj(AT.nzval[tt])*x[AT.rowval[tt]];
	end
end

#void computeResidualAtIdx(long long *rowptr , double complex *valA ,long long *colA,double complex *b,double complex *x,long long* Idxs, double complex *local_r, int blockSize){
# local_r = r*0.0;
# ccall((:computeResidualAtIdx,lib),Void,(Ptr{Int64},Ptr{Complex128},Ptr{Int64},Ptr{Complex128},Ptr{Complex128},Ptr{Int64},Ptr{Complex128},Int16,)
								 # ,AT.colptr,AT.nzval,AT.rowval,b,x,Idxs,local_r,convert(Int16,length(Idxs)));
# if norm(local_r - r) > 1e-14
	# error("Fix computeResidualAtIdx in C: ",norm(local_r - r));
# end
return r;
end

function performSetup(i::Array{Int64}, AI::SparseMatrixCSC,  DDparam::DomainDecompositionParam)
	# println("setupping sub domain: ",i," by ",myid());
	if DDparam.getSubDomainMass!=identity
		AI = AI + DDparam.getSubDomainMass(DDparam,i);
	end
	DDPrec = DomainDecompositionPreconditionerParam([],i,AI,[],copySolver(DDparam.Ainv))
	DDPrec.Ainv = setupSolver(AI,DDPrec.Ainv);

	if isa(DDPrec.Ainv,ParallelJuliaSolver.parallelJuliaSolver)
		DDPrec.A_i = spzeros(0,0);
	end
	# println("setup for sub domain: ",i," all-done by ",myid());
	return DDPrec;
end

function performSetup(i::Array{Int64}, AI::DomainDecompositionOperatorConstructor,  DDparam::DomainDecompositionParam)
	AII = AI.getOperator(AI.problem_param);
	if AI.getDirichletMass !=identity
		DirichletMass = AI.getDirichletMass(DDparam,AI.problem_param,i)[:];
		AII = AII + sparse(Diagonal(DirichletMass));
	end
	DDPrec = DomainDecompositionPreconditionerParam(AI.problem_param,i,AII,DirichletMass,copySolver(DDparam.Ainv))
	DDPrec.Ainv = setupSolver(AII,DDPrec.Ainv);
	# println("setup for sub domain: ",i," by ",myid());
	return DDPrec;
end

function computeResidual(precParam::DomainDecompositionPreconditionerParam, x, r)
	# r is initialized by b and overrun here
	r .-= precParam.A_i'*x;

	# AT = precParam.A_i
	# for t = 1:length(r)
		# for tt = AT.colptr[t]:AT.colptr[t+1]-1
			# r[t] = r[t] - conj(AT.nzval[tt])*x[AT.rowval[tt]];
		# end
	# end
	return r;
end

function setupDDSerial(AT::Union{SparseMatrixCSC,DomainDecompositionOperatorConstructor},DDparam::DomainDecompositionParam)
M 			= DDparam.Mesh;
numDomains 	= DDparam.numDomains;
overlap 	= DDparam.overlap;
DDPreconditioners 	= Array{DomainDecompositionPreconditionerParam}(undef,prod(numDomains));
DDparam.GlobalIndices = Array{Array{DDIndType}}(undef,prod(numDomains));
n = M.n;
for ii = 1:prod(numDomains)
	i = cs2loc(ii,numDomains);
	IIp = DDparam.getIndicesOfCell(numDomains,overlap, i,n);
	if isa(AT,SparseMatrixCSC)==true
		AI = sparse(AT[IIp,IIp]');
		DDPreconditioners[ii] = performSetup(i, AI,  DDparam);
	else
		subparams = AT.getSubParams(AT.problem_param, M,i,numDomains,overlap);
		AI = DomainDecompositionOperatorConstructor(subparams,AT.getSubParams,AT.getOperator,AT.getDirichletMass);
		DDPreconditioners[ii] = performSetup(i, AI,  DDparam);
	end
	IIp = convert(Array{DDIndType},IIp);
	DDparam.GlobalIndices[ii] = IIp;
end
DDparam.PrecParams = DDPreconditioners;
return DDparam;
end

function solveDDSerial(AT::SparseMatrixCSC,b::Array,x::Array,DDparam::DomainDecompositionParam,niter=1,doTranspose::Int64=0)

ncells = DDparam.Mesh.n;
dim = length(ncells);
numDomains = DDparam.numDomains;
overlap = DDparam.overlap;

for k=1:niter
	for color = 1:2^dim
		for ic = 1:prod(numDomains);
			icloc = cs2loc(ic,numDomains);
			if cellColor(icloc)==color
				#IIp = DDparam.getIndicesOfCell(numDomains,overlap,icloc,ncells);
				IIp = DDparam.GlobalIndices[ic];
				if isa(AT,SparseMatrixCSC)==true
					r = computeResidualAtIdx(AT,b,x,IIp);
				else
					r = computeResidualAtIdx(AT,b,x,IIp);
					# r = computeResidual(DDparam.PrecParams[ic],x[IIp],b[IIp]);
				end
				(t,DDparam.PrecParams[ic]) = solveSubDomain(r,DDparam.PrecParams[ic],doTranspose)
				x[IIp] = x[IIp] + t;
			end
		end
	end
	if niter > 1
		println(k,": ",norm(b - AT'*x)/norm(b))
	end
end

return x,DDparam;
end

function solveGSDDSerial(AT::SparseMatrixCSC,b::Array,x::Array,DDparam::DomainDecompositionParam,niter=1,doTranspose::Int64=0)

ncells = DDparam.Mesh.n;
dim = length(ncells);
numDomains = DDparam.numDomains;
overlap = DDparam.overlap;

for k=1:niter
	for ic = 1:prod(numDomains);
		# icloc = cs2loc(ic,numDomains);
		#IIp = DDparam.getIndicesOfCell(numDomains,overlap,icloc,ncells);
		IIp = DDparam.GlobalIndices[ic];
		if isa(AT,SparseMatrixCSC)==true
			r = computeResidualAtIdx(AT,b,x,IIp);
		else
			r = computeResidualAtIdx(AT,b,x,IIp);
			# r = computeResidual(DDparam.PrecParams[ic],x[IIp],b[IIp]);
		end
		(t,DDparam.PrecParams[ic]) = solveSubDomain(r,DDparam.PrecParams[ic],doTranspose)
		x[IIp] = x[IIp] + t;
	end
	for ic = prod(numDomains):-1:1;
		# icloc = cs2loc(ic,numDomains);
		#IIp = DDparam.getIndicesOfCell(numDomains,overlap,icloc,ncells);
		IIp = DDparam.GlobalIndices[ic];
		if isa(AT,SparseMatrixCSC)==true
			r = computeResidualAtIdx(AT,b,x,IIp);
		else
			r = computeResidualAtIdx(AT,b,x,IIp);
			# r = computeResidual(DDparam.PrecParams[ic],x[IIp],b[IIp]);
		end
		(t,DDparam.PrecParams[ic]) = solveSubDomain(r,DDparam.PrecParams[ic],doTranspose)
		x[IIp] = x[IIp] + t;
	end
end
if niter > 1
	println(k,": ",norm(b - AT'*x)/norm(b))
end
return x,DDparam;
end

function solveSubDomain(r,prec::DomainDecompositionPreconditionerParam,doTranspose::Int64)
	t = zeros(eltype(r),size(r));
	t,prec.Ainv = solveLinearSystem!(prec.A_i,r,t,prec.Ainv,doTranspose);
	return t,prec;
end


text:
function FieldOrekitFixedStepHandlerMultiplexer()
    return FieldOrekitFixedStepHandlerMultiplexer(())
end

function add(obj::FieldOrekitFixedStepHandlerMultiplexer, arg0::FieldOrekitFixedStepHandler)
    return jcall(obj, "add", void, (FieldOrekitFixedStepHandler,), arg0)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function handle_step(obj::FieldOrekitFixedStepHandlerMultiplexer, arg0::FieldSpacecraftState, arg1::jboolean)
    return jcall(obj, "handleStep", void, (FieldSpacecraftState, jboolean), arg0, arg1)
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function init(obj::FieldOrekitFixedStepHandlerMultiplexer, arg0::FieldSpacecraftState, arg1::FieldAbsoluteDate, arg2::RealFieldElement)
    return jcall(obj, "init", void, (FieldSpacecraftState, FieldAbsoluteDate, RealFieldElement), arg0, arg1, arg2)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
export Shamanskii_MA57

function Shamanskii_MA57(nlp_stop, PData::PDataFact, dₙ::Vector, gt::Vector)
    nlp_at_x = nlp_stop.current_state
    x = nlp_at_x.x
    xdemi = x + dₙ
    gtemp = grad(nlp_stop.pb, xdemi)
    ϵ2 = sqrt(eps(eltype(gt)))
    Γ = max.(abs.(PData.Δ), ϵ2)

    sg = PData.s .* gtemp
    ĝ = PData.L \ (sg[PData.pp])
    g̃ = PData.Q' * ĝ
    d̃ = -(PData.Δ) .\ g̃
    d = AInv(PData, d̃)

    dₕₒ = dₙ .+ d
    return dₕₒ
end


text:
# Julia 1.0 compat
if VERSION < v"1.1"
    @eval isnothing(::Any) = false
    @eval isnothing(::Nothing) = true
end

const possessive_pronouns = Dict(
    "they/them" => "their",
    "she/her" => "her",
    "he/him" => "his"
)

"""
    generate_card(line)

Return the generated moderation card from the given `line`.

# Examples
```jldoctest
julia> println(generate_card("»How I learned to say Farewell« – Madeline (she/her). Start: 13:40, Q&A: 14:00, End: 14:05"))
- Our next speaker is Madeline
- Her talk is called »How I learned to say Farewell«
- She will answer your questions in the Q&A session at the end of the talk, starting at 14:00

13:40 - 14:00 - 14:05
```
"""
function generate_card(line)
    # You can also write this as a multi-line regex string:
    # re = r"""»(.*)«                     # title
    #          \ –\ (.*?)                 # name
    #          (?:\ \((.*?)\)|)           # optional pronouns
    #          (?:,\ from\ (.*?)|).       # optional organisation
    #          \ Start:\ (.*?),\          # start time
    #          (?:Q&A:\ (.*?),\ |)        # optional Q&A time
    #          End:\ (.*)                 # end time
    #       """x
    re = r"»(.*?)« – (.*?)(?: \((.*?)\)|)(?:, from (.*?)|). Start: (.*?), (?:Q&A: (.*?), |)End: (.*)"

    title, name, pronouns, org, start_time, qa_time, end_time = match(re, line).captures

    they = isnothing(pronouns) ? first(split(name)) : titlecase(first(split(pronouns, '/')))
    their = isnothing(pronouns) ? first(split(name)) * "'s" : titlecase(possessive_pronouns[pronouns])
    qa_line = (isnothing(qa_time) ?
               "There will not be a Q&A session" :
               "$they will answer your questions in the Q&A session at the end of the talk, starting at $qa_time")

    """
    - Our next speaker is $(name * (isnothing(org) ? "" : ", from $org"))
    - $their talk is called »$title«
    - $qa_line

    $start_time - $(isnothing(qa_time) ? "NO Q&A -" : "$qa_time -") $end_time
    """
end


text:
#md # !!! note
#
#     This tutorial is under construction, but the code below is already functional.
#
#
# ## Problem statement

using Gridap
using LinearAlgebra: inv, det
using LineSearches: BackTracking

# Material parameters
const λ = 100.0
const μ = 1.0

# Deformation Gradient
F(∇u) = one(∇u) + ∇u'

J(F) = sqrt(det(C(F)))

#Green strain

#E(F) = 0.5*( F'*F - one(F) )

@law dE(∇du,∇u) = 0.5*( ∇du⋅F(∇u) + (∇du⋅F(∇u))' )

# Right Cauchy-green deformation tensor

C(F) = (F')⋅F

# Constitutive law (Neo hookean)

@law function S(∇u)
  Cinv = inv(C(F(∇u)))
  μ*(one(∇u)-Cinv) + λ*log(J(F(∇u)))*Cinv
end

@law function dS(∇du,∇u)
  Cinv = inv(C(F(∇u)))
  _dE = dE(∇du,∇u)
	λ*(Cinv⊙_dE)*Cinv + 2*(μ-λ*log(J(F(∇u))))*Cinv⋅_dE⋅(Cinv')
end

# Cauchy stress tensor

@law σ(∇u) = (1.0/J(F(∇u)))*F(∇u)⋅S(∇u)⋅(F(∇u))'

# Weak form

res(u,v) = dE(∇(v),∇(u)) ⊙ S(∇(u))

jac_mat(u,du,v) =  dE(∇(v),∇(u)) ⊙ dS(∇(du),∇(u))

jac_geo(u,du,v) = ∇(v) ⊙ ( S(∇(u))⋅∇(du) )

jac(u,du,v) = jac_mat(u,v,du) + jac_geo(u,v,du)

# Model
domain = (0,1,0,1)
partition = (20,20)
model = CartesianDiscreteModel(domain,partition)

# Define new boundaries
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"diri_0",[1,3,7])
add_tag_from_tags!(labels,"diri_1",[2,4,8])

# Construct the FEspace
V = TestFESpace(
  model=model,valuetype=VectorValue{2,Float64},
  reffe=:Lagrangian,conformity=:H1,order=1,
  dirichlet_tags = ["diri_0", "diri_1"])

# Setup integration
trian = Triangulation(model)
degree = 2
quad = CellQuadrature(trian,degree)

# Setup weak form terms
t_Ω = FETerm(res,jac,trian,quad)

# Setup non-linear solver
nls = NLSolver(
  show_trace=true,
  method=:newton,
  linesearch=BackTracking())

solver = FESolver(nls)

function run(x0,disp_x,step,nsteps)

  g0 = VectorValue(0.0,0.0)
  g1 = VectorValue(disp_x,0.0)
  U = TrialFESpace(V,[g0,g1])

  #FE problem
  op = FEOperator(U,V,t_Ω)

  println("\n+++ Solving for disp_x $disp_x in step $step of $nsteps +++\n")

  uh = FEFunction(U,x0)

  uh, = solve!(uh,solver,op)

  writevtk(trian,"results_$(lpad(step,3,'0'))",cellfields=["uh"=>uh,"sigma"=>σ(∇(uh))])

  return get_free_values(uh)

end

function runs()

 disp_max = 0.75
 disp_inc = 0.02
 nsteps = ceil(Int,abs(disp_max)/disp_inc)

 x0 = zeros(Float64,num_free_dofs(V))

 for step in 1:nsteps
   disp_x = step * disp_max / nsteps
   x0 = run(x0,disp_x,step,nsteps)
 end

end

#Do the work!
runs()

# Picture of the last load step
# ![](../assets/hyperelasticity/neo_hook_2d.png)
#
# ##  Extension to 3D
#
# Extending this tutorial to the 3D case is straightforward. It is leaved as an exercise.
#
# ![](../assets/hyperelasticity/neo_hook_3d.png)

max_stars_repo_name: andre-massing/Tutorials

max_issues_repo_name: andre-massing/Tutorials

max_forks_repo_name: andre-massing/Tutorials


text:
"""
    function addCleanInteraction_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

Generate the functions related with Cleaning Interaction Updates.
"""
function addCleanInteraction_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

    if "UpdateLocalInteraction" in keys(abm.declaredUpdates) || "UpdateInteraction" in keys(abm.declaredUpdates)

        #Construct cleanup function
        s = []
        for i in ["UpdateInteraction"]
            if i in keys(abm.declaredUpdates)
                syms = symbols_(abm,abm.declaredUpdates[i])
                syms = syms[Bool.((syms[:,"placeDeclaration"] .== :Model) .* syms[:,"updated"]),"Symbol"]
                append!(s,syms)
            end
        end
        unique!(s)
        up = quote end
        for i in s
            pos = findfirst(abm.declaredSymbols["Local"] .== i)
            push!(up.args,:(localV[ic1_,$pos]=0))
        end
        fclean = simpleFirstLoopWrapInFunction_(platform,:cleanInteraction_!,up)
        push!(p.declareF.args,fclean)
    end

    return nothing
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
using Gtk.ShortNames, GtkReactive, Graphics, Cairo, Colors

mutable struct Hexagon
    center::Point
    radius::Int
    letter::String
    color::Colorant
end

const offset = 50
const hgt = 450
const wid = 400
const hcombdim = (rows = 5, cols = 4)
const randletters = reshape(string.(Char.(shuffle(UInt8('A'):UInt8('Z'))))[1:20], Tuple(hcombdim))
const win = Window("Honeycombs", wid, hgt)
const can = Canvas()
const honeycomb = Dict{Point, Hexagon}()
const chosen = Vector{String}()

function hexmat(p, rad)
    shor = rad * 0.5
    long = rad * sqrt(3.0) / 2.0
    mat = reshape([shor, long, -shor, long, Float64(-rad), 0.0, -shor, -long, shor, -long, Float64(rad), 0.0], 2, 6)
    [Point(mat[1, n] + p.x, mat[2, n] + p.y) for n in 1:6]
end

function whichclicked(clickpos)
    centers = [c for c in keys(honeycomb)]
    (maybeclicked, idx) = findmin(map(c -> sqrt((clickpos.x - c.x)^2 + (clickpos.y - c.y)^2), centers))
    return maybeclicked < offset * sqrt(3) / 2.0 ? centers[idx] : nothing
end

whichtyped(ch) = (for (k, v) in honeycomb if v.letter == ch return k end end; nothing)

function hexagon(ctx, pos, rad, ltr, colr = colorant"yellow")
    set_source(ctx, colr)
    points = hexmat(pos, rad)
    set_line_width(ctx, 4)
    polygon(ctx, points)
    close_path(ctx)
    fill(ctx)
    set_source(ctx, colorant"black")
    polygon(ctx, points)
    close_path(ctx)
    stroke(ctx)
    move_to(ctx, pos.x - (ltr == "I" ? 7 : 18), pos.y + 15)
    set_source(ctx, colr == colorant"yellow" ? colorant"red" : colorant"black")
    set_font_size(ctx, offset)
    show_text(ctx, ltr)
    Hexagon(pos, rad, ltr, colr)
end

hexagon(ctx, h::Hexagon) = hexagon(ctx, h.center, h.radius, h.letter, h.color)

function makehoneycomb(ctx)
    centers = fill(Point(0, 0), hcombdim.rows, hcombdim.cols)
    xdelta = 75.0
    ydelta = 90.0
    for i in 1:hcombdim.rows, j in 1:hcombdim.cols
        center = Point((i - 1) * xdelta + offset, (j - 1) * ydelta + ((i - 1 ) % 2 + 1) * offset)
        centers[i, j] = center
        honeycomb[center] = hexagon(ctx, center, offset, randletters[i, j])
    end
    centers
end

@guarded draw(can) do widget
    ctx = getgc(can)
    if length(honeycomb) == 0
        makehoneycomb(ctx)
    else
        map(c -> hexagon(ctx, honeycomb[c]), collect(keys(honeycomb)))
    end
end

""" At entry to this function we have just found out what letter was chosen."""
function changecolor(colr)
    h = honeycomb[colr]
    h.color = colorant"violet"
    hexagon(getgc(can), h)
    reveal(win, true)
    push!(chosen, h.letter)
    if all(map(k -> honeycomb[k].color == colorant"violet", collect(keys(honeycomb))))
        println("All hexagons ($chosen, and the last letter was $(chosen[end])) have been chosen. Exiting.")
        exit(0)
    end
end

signal_connect(win, "key-press-event") do widget, event
    if (whichhexgon = whichtyped(string(uppercase(Char(event.keyval))))) != nothing
        changecolor(whichhexgon)
    end
end

can.mouse.button1press = @guarded (widget, event) -> begin
    if (whichhexgon = whichclicked(Point(event.x, event.y))) != nothing
        changecolor(whichhexgon)
    end
end

push!(win, can)
show(can)
condition = Condition()
endit(w) = notify(condition)
signal_connect(endit, win, :destroy)
show(win)
wait(condition)
exit()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_index

"""
    ''_index(;kwargs...)
    ''_index(children::Any;kwargs...)
    ''_index(children_maker::Function;kwargs...)

An index component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional)
- `id` (String; optional)
- `width` (Real; optional)
"""
function ''_index(; kwargs...)
        available_props = Symbol[:children, :id, :width]
        wild_props = Symbol[]
        return Component("''_index", "index", "python_version", available_props, wild_props; kwargs...)
end

''_index(children::Any; kwargs...) = ''_index(;kwargs..., children = children)
''_index(children_maker::Function; kwargs...) = ''_index(children_maker(); kwargs...)

max_stars_repo_licenses: CNRI-Python Xnet X11

max_issues_repo_licenses: CNRI-Python Xnet X11

max_forks_repo_licenses: CNRI-Python Xnet X11


text:
__precompile__()

module Bokeh

include("bokehjs.jl")
include("glyphs.jl")
include("display.jl")
include("generate.jl")
include("plot.jl")

# if DEBUG == true, displays indented JSON, uses unminified js and saves the raw
# JSON to "bokeh_models.json" if true
DEBUG = false
debug(b::Bool) = (global DEBUG = b)
debug() = DEBUG

# whether or not to show the plot immediately after `plot`
AUTOOPEN = isdefined(Main, :IJulia)
autoopen(b::Bool) = (global AUTOOPEN = b)
autoopen() = AUTOOPEN

# default width of plot
WIDTH = 800
width(w::Int) = (global WIDTH = w)
width() = WIDTH

# default height of plot
HEIGHT = 600
height(h::Int) = (global HEIGHT = h)
height() = HEIGHT

# default axis types
X_AXIS_TYPE = :auto
Y_AXIS_TYPE = :auto

# default glyph type
DEFAULT_GLYPHS_STR = "b|r|g|k|y|c|m|b--|r--|g--|k--|y--|c--|m--|--"
DEFAULT_GLYPHS = convert(Vector{Glyph}, DEFAULT_GLYPHS_STR)
glyphs(gs::Vector{Glyph}) = (global DEFAULT_GLYPHS = gs)
glyphs(s::AbstractString) = (global DEFAULT_GLYPHS = convert(Vector{Glyph},s))
glyphs() = DEFAULT_GLYPHS

# default glyph size
DEFAULT_SIZE = 6
glyphsize(gs::Int) = (global DEFAULT_SIZE = gs)
glyphsize() = DEFAULT_SIZE

# default alpha value for filled glyphs
DEFAULT_FILL_ALPHA = 0.7

# default filename
PLOTFILE = "bokeh_plot.html"
function warn_overwrite()
    if ispath(PLOTFILE) && isinteractive() && !isdefined(Main, :IJulia)
        warn("""
        $PLOTFILE already exists, it will be overwritten when a plot is generated
        Change the output file with plotfile(<new file name>)
        """)
    end
end

warn_overwrite()

function plotfile(fn::AbstractString)
    global PLOTFILE = fn
    warn_overwrite()
    nothing
end
plotfile() = PLOTFILE

# default plot title
TITLE = "Bokeh Plot"
title(t::AbstractString) = (global TITLE = t)
title() = TITLE

# number of points at which to evaluate functions
COUNTEVAL = 500
counteval(count::Int) = (global COUNTEVAL = count)
counteval() = COUNTEVAL

# hold on to plots
HOLD = false

function hold(h::Bool, clear::Union{Bool, Void}=nothing)
    if (!h && clear == nothing) || clear == true
        global CURPLOT = nothing
    end
    global HOLD = h
end

hold() = hold(!HOLD)  # method to flip hold status

# current plot object
CURPLOT = nothing
curplot(cp::Plot) = (global CURPLOT = cp)
curplot() = CURPLOT

# tools to add to display
TOOLS = [:pan, :wheelzoom, :boxzoom, :resize, :reset]
tools(t::Vector{Symbol}) = (global TOOLS = t)
tools() = TOOLS

# to avoid giving the same warning lots of times remember the file we've
# just warned about overwriting
WARN_FILE = nothing

# this overrides autoopen and disables opening html files used for travis CI,
# shouldn't be necessary elsewhere
NOSHOW = false
noshow(b::Bool) = (global NOSHOW = b)

# generally js and css is not copied into HTML pages but rather the files
# referenced, the only common exception is IJulia plots, this allows that to be
# ovwritten so js and css are always included in pages
INCLUDE_JS = false
includejs(ijs::Bool) = (global INCLUDE_JS = ijs)

# filesystem warning switch
FILE_WARNINGS = true
filewarnings(warn::Bool) = (global FILE_WARNINGS = warn)

# if we're in IJulia call setupnotebook to load js and css
if isdefined(Main, :IJulia) && Main.IJulia.inited
    setupnotebook()
end

export plot,
       setupnotebook,
       Glyph,
       Plot,
       BokehDataSet,
       renderplot,
       genplot,
       showplot,
       autoopen,
       width,
       height,
       glyphs,
       glyphsize,
       plotfile,
       title,
       counteval,
       hold,
       curplot,
       tools

end  # module


text:
##########
## table
##########

d = DataVector[1,1,2,4,5,9,NA,NA]
t = table(d)
@assert t[1] == 2
@assert t[NA] == 2

x = LETTERS[[1,1,9,5,9,9,9]]
t = table(x)
@assert t["A"] == 2
@assert t["I"] == 4
@assert length(t) == 3

##########
## paste
##########

@assert paste(["a", "b"], "X", [1:2]) == ["aX1", "bX2"]
@assert paste(["a", "b"], "X", [1:4]) == ["aX1", "bX2", "aX3", "bX4"]

##########
## cut
##########

@assert isequal(cut([2, 3, 5], [1, 3, 6]), PooledDataArray(["(1,3]", "(1,3]", "(3,6]"]))
@assert isequal(cut([2, 3, 5], [3, 6]), PooledDataArray(["[2,3]", "[2,3]", "(3,6]"]))
@assert isequal(cut([2, 3, 5, 6], [3, 6]), PooledDataArray(["[2,3]", "[2,3]", "(3,6]", "(3,6]"]))
@assert isequal(cut([1, 2, 4], [1, 3, 6]), PooledDataArray(["[1,3]", "[1,3]", "(3,6]"]))
@assert isequal(cut([1, 2, 4], [3, 6]), PooledDataArray(["[1,3]", "[1,3]", "(3,6]"]))
@assert isequal(cut([1, 2, 4], [3]), PooledDataArray(["[1,3]", "[1,3]", "(3,4]"]))
@assert isequal(cut([1, 5, 7], [3, 6]), PooledDataArray(["[1,3]", "(3,6]", "(6,7]"]))

ages = [20, 22, 25, 27, 21, 23, 37, 31, 61, 45, 41, 32]
bins = [18, 25, 35, 60, 100]
cats = cut(ages, bins)
pdv = PooledDataArray(["(18,25]", "(18,25]", "(18,25]",
                       "(25,35]", "(18,25]", "(18,25]",
                       "(35,60]", "(25,35]", "(60,100]",
                       "(35,60]", "(35,60]", "(25,35]"])
@assert isequal(cats, pdv)

##########
## rep
##########

@assert rep(3, 2) == [3,3]
@assert rep([3,4], 2) == [3,4,3,4]
@assert rep([3,4], 1, 2) == [3,3,4,4]
@assert rep([3,4], each = 2) == [3,3,4,4]
@assert rep([3,4], times = 2) == [3,4,3,4]
@assert rep([3,4], times = [2,3]) == [3,3,4,4,4]
@assert rep([3,4], [2,3]) == [3,3,4,4,4]
@assert isequal(rep(DataVector[NA,3,4], 2), DataVector[NA,3,4,NA,3,4])
@assert isequal(rep(DataVector[NA,3,4], [2,1,2]), DataVector[NA,NA,3,4,4])
@assert isequal(rep(DataVector[NA,3,4], [2,1,0]), DataVector[NA,NA,3])

##########
## findat
##########

@assert findat([4,1,9,10,1], [1:6]) == [4,1,0,0,1]
@assert findat([1:6], [1,5,9,3,1]) == [1,0,4,0,2,0]


text:
# SPDX-License-Identifier: MIT

# Half-vectorization of a matrix.
function vech(A::AbstractMatrix, k=0)
    @assert size(A, 1) == size(A, 2)

    D = size(A, 1)
    out = similar(A, Int((D-k) + ((D-k)^2 - (D-k)) / 2))
    i = 0
    for c = 1:(D-k)
        for r = (c+k):D
            @inbounds out[i += 1] = A[r,c]
        end
    end
    out
end

struct CompressedLowerTriangular{T,S<:AbstractVector{T}} <: AbstractMatrix{T}
    D::Int64
    k::Int64
    vech::S
end

function CompressedLowerTriangular(M::AbstractMatrix, k = 0)
    T = eltype(M)
    v = vech(M, k)
    CompressedLowerTriangular{T, typeof(v)}(size(M, 1), k, v)
end

@Zygote.adjoint CompressedLowerTriangular(D, k, v) =
    CompressedLowerTriangular(D, k, v), c -> (0, 0, vech(c, k))

@Zygote.adjoint vech(M, k) =
    vech(M, k), c -> (CompressedLowerTriangular(size(M, 1), k, c), 0)

Base.size(M::CompressedLowerTriangular) = (M.D, M.D)

function Base.getindex(M::CompressedLowerTriangular, r::Int, c::Int)
    (r > M.k && c <= (M.D - M.k) && c <= (r - M.k)) || return zero(eltype(M))
    s = ((c-1) * c) ÷ 2
    i = (c-1) * (M.D - M.k) + (r - M.k) - s
    M.vech[i]
end

function Base.replace_in_print_matrix(M::CompressedLowerTriangular, r::Integer,
                                      c::Integer, s::AbstractString)
    if (r > M.k && c <= (M.D - M.k) && c <= (r - M.k))
        return s
    else
        Base.replace_with_centered_mark(s)
    end
end

CompressedSymmetric(D, k, v) = Symmetric(CompressedLowerTriangular(D, k, v), :L)

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
function uppersum(M)
    n = size(M, 1)
    s = zero(eltype(M))
    for i in 2:n
        @simd for j in 1:(i-1)
            @inbounds s += M[j, i]
        end
    end
    s
end

@time m = randn(10_000, 10_000);
m
@time sum(m)
@time uppersum(M)

max_stars_repo_name: bkamins/WarszawskieForumJulia

max_issues_repo_name: bkamins/WarszawskieForumJulia

max_forks_repo_name: bkamins/WarszawskieForumJulia


text:
# the following code assumes you have a Kafka broker running 127.0.0.1:9092
# with topic "test" and one partition

using Kafka

# create KafkaClient using single bootstrap broker
# KafkaClient will then fetch information of all other brokers
# and initialize connections to all of them
kc = KafkaClient("127.0.0.1", 9092)

# get metadata about a topic(s)
# all requests in Kafka.jl return channels that hold actual result
# this way requests may be chained asynchronously and results taken later
md_channel = metadata(kc, ["test"])
md = take!(md_channel)
# or, if you prefer synchronous logic, use one-linear
take!(metadata(kc, ["test"]))
# if you don't pass any topics, metadata for all of them will be returned
take!(metadata(kc))

# get earliest and latest available offsets for topic "test" and partition 0
take!(earliest_offset(kc, "test", 0))
take!(latest_offset(kc, "test", 0))
# or list offsets before particular time
# WARNING: this method follows the protocol, but I couldn't make it working
#          please, report if you can or see an error in usage!
now = floor(Int64, time())
max_number_of_offsets = 100
take!(list_offsets(kc, "test", 0, now, max_number_of_offsets))

# produce new messages
# each message is a key-value pair where both key and value are byte arrays
keys = [unsafe_wrap(Vector{UInt8}, key) for key in ["1", "2", "3"]]
values = [unsafe_wrap(Vector{UInt8}, value) for value in ["feel", "good", "inc."]]
messages = collect(zip(keys, values))
# messages are produced to a specific topic (e.g. "test") and partition (e.g. 0)
# take!(produce(...)) returns an offset of the first message
offset = take!(produce(kc, "test", 0, messages))

# fetch messages
# in addition to topic name and partition id, fetch() accepts
# offset to start reading from
# return value is an array of triples (offset, key, value)
# note that due to storage internals, brokers may actually return messages
# with offsets less than start_offset; client is responsible to filter them out
start_offset = 0
offset_messages = take!(fetch(kc, "test", 0, start_offset))
# fetch() supports several options to control troughput/latency tradeoff:
# * max_wait_time - max time to wait for new messages (in milliseconds);
#   default is 100
# * min_bytes - minimum number of bytes that broker should have before returning
#   response to a client; default is 1KB
# * max_bytes - maximum number of bytes to include into response;
#   default is 1M
offset_messages = take!(fetch(kc, "test", 0, start_offset, max_wait_time=5000))

# Kafka.jl tries to make convenient wrapper around Kafka protocol,
# but you can always get access to raw responses using methods with "_" prefix
# for example:
produce_resp = take!(_produce(kc, "test", 0, messages))
fetch_resp = take!(_fetch(kc, "test", 0, start_offset))


text:
const cube_edges = ( 0, 1, 0, 2, 0, 4, 1, 3, 1, 5, 2, 3,
                     2, 6, 3, 7, 4, 5, 4, 6, 5, 7, 6, 7 )
const sn_edge_table = (7, 25, 30, 98, 101, 123, 124, 168, 175, 177, 182, 202,
                       205, 211, 212, 772, 771, 797, 794, 870, 865, 895, 888,
                       940, 939, 949, 946, 974, 969, 983, 976, 1296, 1303, 1289,
                       1294, 1394, 1397, 1387, 1388, 1464, 1471, 1441, 1446,
                       1498, 1501, 1475, 1476, 1556, 1555, 1549, 1546, 1654,
                       1649, 1647, 1640, 1724, 1723, 1701, 1698, 1758, 1753,
                       1735, 1728, 2624, 2631, 2649, 2654, 2594, 2597, 2619,
                       2620, 2792, 2799, 2801, 2806, 2698, 2701, 2707, 2708,
                       2372, 2371, 2397, 2394, 2342, 2337, 2367, 2360, 2540,
                       2539, 2549, 2546, 2446, 2441, 2455, 2448, 3920, 3927,
                       3913, 3918, 3890, 3893, 3883, 3884, 4088, 4095, 4065,
                       4070, 3994, 3997, 3971, 3972, 3156, 3155, 3149, 3146,
                       3126, 3121, 3119, 3112, 3324, 3323, 3301, 3298, 3230,
                       3225, 3207, 3200, 3200, 3207, 3225, 3230, 3298, 3301,
                       3323, 3324, 3112, 3119, 3121, 3126, 3146, 3149, 3155,
                       3156, 3972, 3971, 3997, 3994, 4070, 4065, 4095, 4088,
                       3884, 3883, 3893, 3890, 3918, 3913, 3927, 3920, 2448,
                       2455, 2441, 2446, 2546, 2549, 2539, 2540, 2360, 2367,
                       2337, 2342, 2394, 2397, 2371, 2372, 2708, 2707, 2701,
                       2698, 2806, 2801, 2799, 2792, 2620, 2619, 2597, 2594,
                       2654, 2649, 2631, 2624, 1728, 1735, 1753, 1758, 1698,
                       1701, 1723, 1724, 1640, 1647, 1649, 1654, 1546, 1549,
                       1555, 1556, 1476, 1475, 1501, 1498, 1446, 1441, 1471,
                       1464, 1388, 1387, 1397, 1394, 1294, 1289, 1303, 1296,
                       976, 983, 969, 974, 946, 949, 939, 940, 888, 895, 865,
                       870, 794, 797, 771, 772, 212, 211, 205, 202, 182, 177,
                       175, 168, 124, 123, 101, 98, 30, 25, 7)