text:
module AlgebraInterfacesTests

using Gridap.Algebra
using Test

a = allocate_vector(Vector{Int},10)
@test isa(a,Vector{Int})
@test length(a) == 10

a = allocate_vector(Vector{Int},1:10)
@test isa(a,Vector{Int})
@test length(a) == 10

a = zeros(4,6)

b = allocate_in_range(Vector{Int},a)
@test length(b) == 4

b = allocate_in_domain(Vector{Int},a)
@test length(b) == 6

fill_entries!(b,4)
@test all( b .== 4 )

a = rand(6)
b = rand(6)

copy_entries!(a,b)
@test all( a .== b )

a = rand(6)
c = copy(a)
b = rand(6)

a .+= b
@test all( a .== ( c .+ b) )

a = rand(6)
c = copy(a)
b = rand(6)

a .-= b
@test all( a .== ( c .- b) )

a = rand(6)
c = copy(a)
scale_entries!(a,10)
@test all( a .== 10*c)

a = rand(4,6)
c = rand(4)
d = copy(c)
b = rand(6)

muladd!(c,a,b)

@test all( c .≈ (d .+ a*b ) )

A = Vector{Float64}
n = 10
rows = Base.OneTo(n)
a = nz_counter(ArrayBuilder(A),(rows,))
@test LoopStyle(a) == DoNotLoop()
add_entry!(a,1.0,1)
add_entries!(a,[1.0,-1.0],[1,1])
add_entries!(a,nothing,[1,1])
b = nz_allocation(a)
@test LoopStyle(b) == DoNotLoop()
@test isa(b,Vector{Float64})
@test length(b) == n
add_entry!(b,1.0,1)
add_entry!(b,1.0,1)
add_entry!(b,1.0,4)
add_entries!(b,[1.0,-1.0],[1,1])
add_entries!(b,nothing,[1,1])
r = zeros(n)
r[1] = 2
r[4] = 1
@test b == r
c = create_from_nz(b)
@test c === b
add_entries!(c,[1.0,-1.0],[1,1])
add_entries!(c,nothing,[1,1])

using SparseArrays
using SparseMatricesCSR

m = 6
n = 9
rows = Base.OneTo(m)
cols = Base.OneTo(n)
builder = SparseMatrixBuilder(SparseMatrixCSC{Float64,Int},MinMemory())
counter = nz_counter(builder,(rows,cols))
@test LoopStyle(counter) == Loop()

builder = SparseMatrixBuilder(SparseMatrixCSC{Float64,Int},MinMemory(4))
counter = nz_counter(builder,(rows,cols))
@test counter.colnnzmax == fill(4,n)
@test LoopStyle(counter) == DoNotLoop()

builder = SparseMatrixBuilder(SparseMatrixCSC{Float64,Int},MinCPU())
counter = nz_counter(builder,(rows,cols))
@test LoopStyle(counter) == Loop()

for A in (SparseMatrixCSC{Float64,Int},SparseMatrixCSC{Float64,Int32})
  builder = SparseMatrixBuilder(A,MinMemory())
  m = 6
  n = 9
  rows = Base.OneTo(m)
  cols = Base.OneTo(n)
  a = nz_counter(builder,(rows,cols))
  @test LoopStyle(a) == Loop()
  add_entry!(a,1.0,1,1)
  add_entry!(a,nothing,1,1)
  add_entry!(a,nothing,3,1)
  add_entry!(a,nothing,2,1)
  add_entry!(a,3.0,4,9)
  @test a.colnnzmax == [4, 0, 0, 0, 0, 0, 0, 0, 1]
  b = nz_allocation(a)
  @test LoopStyle(b) == Loop()
  @test b.colnnz == [0, 0, 0, 0, 0, 0, 0, 0, 0]
  @test b.colptr == [1, 5, 5, 5, 5, 5, 5, 5, 5, 6]
  add_entry!(b,1.0,1,1)
  add_entry!(b,nothing,1,1)
  add_entry!(b,4.0,3,1)
  add_entry!(b,2.0,3,1)
  add_entry!(b,8.0,2,1)
  add_entry!(b,3.0,4,9)
  @test b.colnnz == [3, 0, 0, 0, 0, 0, 0, 0, 1]
  c = create_from_nz(b)
  @test LoopStyle(c) == DoNotLoop()
  @test isa(c,get_array_type(builder))
  @test isa(c,A)
  I,J,V = findnz(c)
  @test I == [1,2,3,4]
  @test J == [1,1,1,9]
  @test V == Float64[1,8,6,3]
  add_entry!(c,1.0,1,1)
  add_entry!(c,nothing,1,1)
  add_entry!(c,nothing,3,1)
  add_entry!(c,3.0,4,9)

  a = nz_counter(builder,(rows,cols))
  add_entries!(a,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(a,nothing,[1,1],[1,-1])
  b = nz_allocation(a)
  add_entries!(b,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(b,nothing,[1,1],[1,-1])
  c = create_from_nz(b)
  add_entries!(c,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(c,nothing,[1,1],[1,-1])
end

for A in (SparseMatrixCSC{Float64,Int},SparseMatrixCSC{Float64,Int32})
  builder = SparseMatrixBuilder(A,MinCPU())
  m = 6
  n = 9
  rows = Base.OneTo(m)
  cols = Base.OneTo(n)
  a = nz_counter(builder,(rows,cols))
  @test LoopStyle(a) == Loop()
  add_entry!(a,1.0,1,1)
  add_entry!(a,nothing,1,1)
  add_entry!(a,4.0,3,1)
  add_entry!(a,2.0,3,1)
  add_entry!(a,8.0,2,1)
  add_entry!(a,3.0,4,9)
  @test a.rowptrs == [0, 2, 1, 2, 1, 0, 0]
  b = nz_allocation(a)
  @test LoopStyle(b) == Loop()
  add_entry!(b,1.0,1,1)
  add_entry!(b,nothing,1,1)
  add_entry!(b,4.0,3,1)
  add_entry!(b,2.0,3,1)
  add_entry!(b,8.0,2,1)
  add_entry!(b,3.0,4,9)
  @test b.rowptrs == [3, 4, 6, 7, 7, 7, 7]
  @test b.colvals == [1, 1, 1, 1, 1, 9]
  c = create_from_nz(b)
  @test LoopStyle(c) == DoNotLoop()
  @test isa(c,get_array_type(builder))
  @test isa(c,A)
  I,J,V = findnz(c)
  @test I == [1,2,3,4]
  @test J == [1,1,1,9]
  @test V == Float64[1,8,6,3]
  add_entry!(c,1.0,1,1)
  add_entry!(c,nothing,1,1)
  add_entry!(c,nothing,3,1)
  add_entry!(c,3.0,4,9)

  a = nz_counter(builder,(rows,cols))
  add_entries!(a,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(a,nothing,[1,1],[1,-1])
  b = nz_allocation(a)
  add_entries!(b,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(b,nothing,[1,1],[1,-1])
  c = create_from_nz(b)
  add_entries!(c,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(c,nothing,[1,1],[1,-1])
end

m = 6
n = 9
rows = Base.OneTo(m)
cols = Base.OneTo(n)
builder = SparseMatrixBuilder(SparseMatrixCSR{0,Float64,Int},MinMemory())
counter = nz_counter(builder,(rows,cols))
@test counter.bi == Val(0)
@test LoopStyle(counter) == Loop()

builder = SparseMatrixBuilder(SparseMatrixCSR{1,Float64,Int},MinMemory(4))
counter = nz_counter(builder,(rows,cols))
@test counter.csc.colnnzmax == fill(4,m)
@test LoopStyle(counter) == DoNotLoop()
@test counter.bi == Val(1)

builder = SparseMatrixBuilder(SparseMatrixCSR{0,Float64,Int},MinCPU())
counter = nz_counter(builder,(rows,cols))
@test LoopStyle(counter) == Loop()

for A in (SparseMatrixCSR{1,Float64,Int},SparseMatrixCSR{0,Float64,Int32})
  builder = SparseMatrixBuilder(A,MinMemory())
  m = 9
  n = 6
  rows = Base.OneTo(m)
  cols = Base.OneTo(n)
  a = nz_counter(builder,(rows,cols))
  @test LoopStyle(a) == Loop()
  add_entry!(a,1.0,1,1)
  add_entry!(a,nothing,1,1)
  add_entry!(a,nothing,1,3)
  add_entry!(a,nothing,1,2)
  add_entry!(a,3.0,9,4)
  @test a.csc.colnnzmax == [4, 0, 0, 0, 0, 0, 0, 0, 1]
  b = nz_allocation(a)
  @test LoopStyle(b) == Loop()
  @test b.csc.colnnz == [0, 0, 0, 0, 0, 0, 0, 0, 0]
  @test b.csc.colptr == [1, 5, 5, 5, 5, 5, 5, 5, 5, 6]
  add_entry!(b,1.0,1,1)
  add_entry!(b,nothing,1,1)
  add_entry!(b,4.0,1,3)
  add_entry!(b,2.0,1,3)
  add_entry!(b,8.0,1,2)
  add_entry!(b,3.0,9,4)
  @test b.csc.colnnz == [3, 0, 0, 0, 0, 0, 0, 0, 1]
  c = create_from_nz(b)
  @test LoopStyle(c) == DoNotLoop()
  @test isa(c,get_array_type(builder))
  @test isa(c,A)
  I,J,V = findnz(c)
  @test I == [1,1,1,9]
  @test J == [1,2,3,4]
  @test V == Float64[1,8,6,3]
  add_entry!(c,1.0,1,1)
  add_entry!(c,nothing,1,1)
  add_entry!(c,nothing,1,3)
  add_entry!(c,3.0,9,4)

  a = nz_counter(builder,(rows,cols))
  add_entries!(a,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(a,nothing,[1,1],[1,-1])
  b = nz_allocation(a)
  add_entries!(b,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(b,nothing,[1,1],[1,-1])
  c = create_from_nz(b)
  add_entries!(c,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(c,nothing,[1,1],[1,-1])
end

for A in (SparseMatrixCSR{0,Float64,Int},SparseMatrixCSR{1,Float64,Int32})
  builder = SparseMatrixBuilder(A,MinCPU())
  m = 9
  n = 6
  rows = Base.OneTo(m)
  cols = Base.OneTo(n)
  a = nz_counter(builder,(rows,cols))
  @test LoopStyle(a) == Loop()
  add_entry!(a,1.0,1,1)
  add_entry!(a,nothing,1,1)
  add_entry!(a,4.0,1,3)
  add_entry!(a,2.0,1,3)
  add_entry!(a,8.0,1,2)
  add_entry!(a,3.0,9,4)
  @test a.csc.rowptrs == [0, 2, 1, 2, 1, 0, 0]
  b = nz_allocation(a)
  @test LoopStyle(b) == Loop()
  add_entry!(b,1.0,1,1)
  add_entry!(b,nothing,1,1)
  add_entry!(b,4.0,1,3)
  add_entry!(b,2.0,1,3)
  add_entry!(b,8.0,1,2)
  add_entry!(b,3.0,9,4)
  @test b.csc.rowptrs == [3, 4, 6, 7, 7, 7, 7]
  @test b.csc.colvals == [1, 1, 1, 1, 1, 9]
  c = create_from_nz(b)
  @test LoopStyle(c) == DoNotLoop()
  @test isa(c,get_array_type(builder))
  @test isa(c,A)
  I,J,V = findnz(c)
  @test I == [1,1,1,9]
  @test J == [1,2,3,4]
  @test V == Float64[1,8,6,3]
  add_entry!(c,1.0,1,1)
  add_entry!(c,nothing,1,1)
  add_entry!(c,nothing,3,1)
  add_entry!(c,3.0,9,4)

  a = nz_counter(builder,(rows,cols))
  add_entries!(a,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(a,nothing,[1,1],[1,-1])
  b = nz_allocation(a)
  add_entries!(b,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(b,nothing,[1,1],[1,-1])
  c = create_from_nz(b)
  add_entries!(c,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(c,nothing,[1,1],[1,-1])
end

for A in (
  SymSparseMatrixCSR{1,Float64,Int},
  SymSparseMatrixCSR{0,Float64,Int})

  builder = SparseMatrixBuilder(A)
  m = 6
  n = 6
  rows = Base.OneTo(m)
  cols = Base.OneTo(n)
  a = nz_counter(builder,(rows,cols))
  @test LoopStyle(a) == Loop()
  add_entry!(a,1.0,1,1)
  add_entry!(a,nothing,1,1)
  add_entry!(a,nothing,3,1)
  add_entry!(a,nothing,1,3)
  add_entry!(a,3.0,2,2)
  add_entry!(a,3.0,2,6)
  add_entry!(a,3.0,6,2)
  @test a.nnz == 5
  b = nz_allocation(a)
  add_entry!(b,1.0,1,1)
  add_entry!(b,nothing,1,1)
  add_entry!(b,nothing,3,1)
  add_entry!(b,nothing,1,3)
  add_entry!(b,3.0,2,2)
  add_entry!(b,3.0,2,6)
  add_entry!(b,3.0,6,2)
  @test LoopStyle(b) == Loop()
  @test length(b.I) == a.nnz
  @test length(b.J) == a.nnz
  @test length(b.V) == a.nnz
  c = create_from_nz(b)
  @test LoopStyle(c) == DoNotLoop()
  @test isa(c,A)
  I,J,V = findnz(c)
  @test I == [1,1,2,2,3,4,5,6]
  @test J == [1,3,2,6,3,4,5,6]
  @test V == Float64[1,0,3,3,0,0,0,0]
  add_entry!(c,1.0,1,1)
  add_entry!(c,nothing,1,1)
  add_entry!(c,nothing,3,1)
  add_entry!(c,nothing,1,3)
  add_entry!(c,3.0,2,2)
  add_entry!(c,3.0,2,6)
  add_entry!(c,3.0,6,2)

  a = nz_counter(builder,(rows,cols))
  add_entries!(a,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(a,nothing,[1,1],[1,-1])
  @test a.nnz == 3
  b = nz_allocation(a)
  add_entries!(b,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(b,nothing,[1,1],[1,-1])
  c = create_from_nz(b)
  add_entries!(c,[1.0 -1.0; -1.0 1.0],[1,-1],[-1,1])
  add_entries!(c,nothing,[1,1],[1,-1])
end

end # module


text:
bg = BenchmarkGroup()
SUITE["core"] = bg

function bench_iteredges(g::AbstractGraph)
  i = 0
  for e in edges(g)
    i += 1
  end
  return i
end

function bench_has_edge(g::AbstractGraph)
    srand(1)
    nvg = nv(g)
    srcs = rand([1:nvg;], cld(nvg, 4))
    dsts = rand([1:nvg;], cld(nvg, 4))
    i = 0
    for (s, d) in zip(srcs, dsts)
        if has_edge(g, s, d)
            i += 1
        end
    end
    return i
end

EDGEFNS = [
    bench_iteredges,
    bench_has_edge
]

for fun in EDGEFNS
    for (name, g) in GRAPHS
        bg["edges","$fun","graph","$name"] = @benchmarkable $fun($g)
    end

    for (name, g) in DIGRAPHS
        bg["edges","$fun","digraph","$name"] = @benchmarkable $fun($g)
    end

end

ext: jl

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
using Primes

isalmostprime(n::Integer, k::Integer) = sum(values(factor(n))) == k

function almostprimes(N::Integer, k::Integer) # return first N almost-k primes
    P = Vector{typeof(k)}(undef,N)
    i = 0; n = 2
    while i < N
        if isalmostprime(n, k) P[i += 1] = n end
        n += 1
    end
    return P
end

for k in 1:5
    println("$k-Almost-primes: ", join(almostprimes(10, k), ", "), "...")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Copyright (c) 2015, Intel Corporation
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#     * Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#     * Neither the name of Intel Corporation nor the names of its contributors
#       may be used to endorse or promote products derived from this software
#       without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

export BatchNormalizationLayer

type BatchNormalizationNeuron <: Neuron
    value :: Float32
    ∇     :: Float32
end

type BatchNormalizationEnsemble <: NormalizationEnsemble
    name         :: Symbol
    neurons      :: Array{BatchNormalizationNeuron}
    connections  :: Vector{Connection}
    num_out      :: Int
    eps          :: Float32
    momentum     :: Float32
    phase        :: Phase
    net_subgroup :: Cint
end

function BatchNormalizationEnsemble(name::Symbol, shape, num_out, eps::Float32, momentum::Float32)
    neurons = Array(BatchNormalizationNeuron, shape...)
    for i = 1:length(neurons)
        neurons[i] = BatchNormalizationNeuron(0.0f0, 0.0f0)
    end

    BatchNormalizationEnsemble(name, neurons, [], num_out, eps, momentum, Train, 1)
end

function init(ensemble::BatchNormalizationEnsemble, net::Net)
    shape = size(ensemble)
    set_buffer(net, symbol(ensemble.name,:value), zeros(Float32, shape..., net.batch_size))
    set_buffer(net, symbol(ensemble.name,:∇), zeros(Float32, shape..., net.batch_size))
    set_buffer(net, symbol(ensemble.name,:running_mean), zeros(Float32, net.batch_size))
    set_buffer(net, symbol(ensemble.name,:running_var), ones(Float32, net.batch_size))
    set_buffer(net, symbol(ensemble.name,:weight), rand(Float32, net.batch_size))
    set_buffer(net, symbol(ensemble.name,:∇weight), zeros(Float32, net.batch_size))
    set_buffer(net, symbol(ensemble.name,:bias), zeros(Float32, net.batch_size))
    set_buffer(net, symbol(ensemble.name,:∇bias), zeros(Float32, net.batch_size))
    set_buffer(net, symbol(ensemble.name,:momentum), Float32[ensemble.momentum])
    set_buffer(net, symbol(ensemble.name,:eps), Float32[ensemble.eps])
end

function get_forward_args(ens::BatchNormalizationEnsemble)
    return [symbol(ens.name,:value),
            symbol(ens.name,:weight),
            symbol(ens.name,:bias),
            symbol(ens.name,:running_mean),
            symbol(ens.name,:running_var),
            symbol(ens.name,:momentum),
            symbol(ens.name,:eps)]
end

function get_backward_args(ens::BatchNormalizationEnsemble)
    return [
        symbol(ens.name, :∇),
        symbol(ens.name, :weight),
        symbol(ens.name, :∇weight),
        symbol(ens.name, :∇bias),
        symbol(ens.name, :momentum),
        symbol(ens.name, :eps)]
end

function BatchNormalizationLayer(name::Symbol, net::Net, input_ens::AbstractEnsemble, eps::Float32, momentum::Float32)
    @assert(ndims(input_ens) == 3)
    ens = BatchNormalizationEnsemble(name, size(input_ens), 3, eps, momentum)
    add_ensemble(net, ens)
    add_connections(net, input_ens, ens, (i, j, k) -> (i:i, j:j, k:k))
    ens
end

@output BatchNormalizationEnsemble function forward(output::Array,
    weight::Array, bias::Array, running_mean::Array, running_var::Array,
    momentum::Array, epsilon::Array, x::Array)

    _momentum = momentum[1]
    _eps = epsilon[1]

    N = size(x, 4)
    inner_size = size(x, 1) * size(x, 2) * size(x, 3)
    for n in 1:N
        x_sum = 0.0f0
        for i_3 in 1:size(x, 3)
            for i_2 in 1:size(x, 2)
                for i_1 in 1:size(x, 1)
                    x_sum += x[i_1, i_2, i_3, n]
                end
            end
        end
        x_mean = x_sum / inner_size
        x_var = 0.0f0
        for i_3 in 1:size(x, 3)
            for i_2 in 1:size(x, 2)
                for i_1 in 1:size(x, 1)
                    dev = x[i_1, i_2, i_3, n] - x_mean
                    x_var += dev * dev
                end
            end
        end
        x_var /= inner_size

        invstd = 1.0f0 / sqrt(x_var + _eps)
        running_mean[n] = _momentum * x_mean + (1 - _momentum) * running_mean[n]
        running_var[n] = _momentum * x_var + (1 - _momentum) * running_var[n]
        for i_3 in 1:size(x, 3)
            for i_2 in 1:size(x, 2)
                for i_1 in 1:size(x, 1)
                    output[i_1, i_2, i_3, n] = ((x[i_1, i_2, i_3, n] - x_mean) * invstd) * weight[n] + bias[n]
                end
            end
        end
    end
end

@output BatchNormalizationEnsemble function backward(output_diff::Array,
    weight::Array, weight_diff::Array, bias_diff::Array,
    momentum::Array, eps::Array, x::Array, x_diff::Array)

    _momentum = momentum[1]
    _eps = eps[1]

    #= x_sum = sum(x, 4) =#
    #= x_var = var(x, 4) =#
    #= xhat = (x - x_sum) ./ sqrt(x_var .+ eps) =#
    #= grad_weight[:] = sum(xhat .* grad_output, 4) =#
    #= grad_bias[:] = sum(grad_output, ndims(x)) =#
    #= xtmp = (xhat .* grad_weight .+ grad_bias) / size(x, 4) =#
    #= grad_input[:] = weight .* (grad_output - xtmp) / sqrt(x_var + eps) =#
end


text:
module GLDebugControl

export printError
export checkError
export hasError
export hasNoError
export getInfoLog
export validate

using ModernGL

using ..GLLists

function printError(title=nothing)
 result=false
	while (err = checkError())[1]
		if typeof(title) == String
			err=err[2]
			errstr = haskey(LIST_ERROR,err) ? LIST_ERROR[err] : ""
			warn("OpenGL Error! $(title): $(errstr)")
		end
		if !result result=true end
	end
	result
end

function checkError()
	err = glGetError()
	(err != GL_NO_ERROR, err)
end

hasError() = checkError()[1]
hasNoError() = !hasError()

function getInfoLog(s::Symbol, id::GLuint)
	len = zeros(GLint,1)
	LIST_STATUS[s][:STATE](id, GL_INFO_LOG_LENGTH, len)
	len = len[]
	#log := strings.Repeat("\x00", int(logLength+1))
	#d[::INFO](id, logLength, nil, glStr(log))
	#glGetShaderInfoLog(id, logLength, nil, gl.Str(log))
	#glGetProgramInfoLog(id, logLength, nil, gl.Str(log))

	if len > 0
		buffer = zeros(GLchar, len)
		sizei = zeros(GLsizei,1)
		LIST_STATUS[s][:INFO](id, len, sizei, buffer)
		len = sizei[]
		unsafe_string(pointer(buffer),len) # bytestring(pointer(buffer), len)
	else
		""
	end
end

validate(s::Symbol, id::GLuint, status::GLenum) = (result = zeros(GLint,1); LIST_STATUS[s][:STATE](id, status, result); result[] == GL_TRUE)

end #GLDebugControl

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
module Day4

mutable struct BingoBoard
    x::Array{Int64, 2}
    checked::Array{Bool, 2}
end

function Base.parse(::Type{BingoBoard}, s)
    s = strip(s)
    b = Array{Int64, 2}(undef, 5, 5)
    checked = zeros(Bool, 5, 5)
    for (i, line) in enumerate(split(s, "\n"))
        row = split(line, " ") |>
            x -> filter(!isempty, x) |>
            x -> parse.(Int64, x)
        b[i, :] = row
    end
    BingoBoard(b, checked)
end

function read(f::String)
    s = Base.read(f, String) |> s -> split(s, "\n\n")
    numbers = popfirst!(s) |>
        a -> split(a, ",") |>
        a -> parse.(Int64, a)
    boards = parse.(BingoBoard, s)
    (boards, numbers)
end

function mark!(board::BingoBoard, value::Int64)
    board.checked = board.checked .| (board.x .== value)
    board
end

check(board::BingoBoard) = (sum(board.checked, dims=1) |> a -> any(x -> x == 5, a)) |
    (sum(board.checked, dims=2) |> a -> any(x -> x == 5, a))

score(board::BingoBoard, last::Int64) = sum(.!(board.checked) .* board.x) * last

function first_winning_board!(boards::Vector{BingoBoard}, numbers::Vector{Int64})
    for num in numbers
        mark!.(boards, num)
        idx = findfirst(check, boards)
        if idx != nothing
            return (boards[idx], num)
        end
    end
    nothing
end

run_part1(data_dir::String) = joinpath(data_dir, "day4.txt") |>
    read |>
    a -> first_winning_board!(a...) |>
    a -> score(a...)

function last_winning_board!(boards::Vector{BingoBoard}, numbers::Vector{Int64})
    for num in numbers
        mark!.(boards, num)
        if (length(boards) == 1) && check(first(boards))
            return (first(boards), num)
        end
        findall(check, boards) |> idxs -> deleteat!(boards, idxs)
    end
    nothing
end

run_part2(data_dir::String) = joinpath(data_dir, "day4.txt") |>
    read |>
    a -> last_winning_board!(a...) |>
    a -> score(a...)

end # module


text:
using Test
using LinearMaps
using SparseArrays
using LinearAlgebra

import Base: *
import LinearAlgebra: issymmetric, mul!

function myft(v::AbstractVector)
    # not so fast fourier transform
    N = length(v)
    w = zeros(complex(eltype(v)), N)
    for k = 1:N
        kappa = (2*(k-1)/N)*pi
        for n = 1:N
            w[k] += v[n]*exp(kappa*(n-1)*im)
        end
    end
    return w
end

A = 2 * rand(ComplexF64, (20, 10)) .- 1
v = rand(ComplexF64, 10)
w = rand(ComplexF64, 20)
V = rand(ComplexF64, 10, 3)
W = rand(ComplexF64, 20, 3)
α = rand()
β = rand()

# test wrapped map for matrix
M = LinearMap(A)
N = LinearMap(M)
Av = A * v
AV = A * V
@test M * v == Av
@test N * v == Av
@test mul!(copy(w), M, v) == Av
@test mul!(copy(w), N, v) == Av
@test mul!(copy(W), M, V) ≈ AV
@test typeof(M * V) <: LinearMap

# test of mul!
@test mul!(copy(w), M, v, 0, 0) == zero(w)
@test mul!(copy(w), M, v, 0, 1) == w
@test mul!(copy(w), M, v, 0, β) == β * w
@test mul!(copy(w), M, v, 1, 1) ≈ Av + w
@test mul!(copy(w), M, v, 1, β) ≈ Av + β * w
@test mul!(copy(w), M, v, α, 1) ≈ α * Av + w
@test mul!(copy(w), M, v, α, β) ≈ α * Av + β * w
@test mul!(copy(w), M, v, α)    ≈ α * Av

# test matrix-mul!
@test mul!(copy(W), M, V, α, β) ≈ α * AV + β * W
@test mul!(copy(W), M, V, α) ≈ α * AV

# test transposition and Matrix
@test M' * w == A' * w
@test mul!(copy(V), adjoint(M), W) ≈ A' * W
@test transpose(M) * w == transpose(A) * w
@test transpose(M') * v ≈ transpose(A') * v
@test mul!(copy(V), transpose(M), W) ≈ transpose(A) * W

@test Matrix(M) == A
@test Array(M) == A
@test convert(Matrix, M) == A
@test convert(Array, M) == A

@test Matrix(M') == A'
@test Matrix(transpose(M)) == copy(transpose(A))

B = LinearMap(Symmetric(rand(10, 10)))
@test transpose(B) == B
@test B == transpose(B)

B = LinearMap(Hermitian(rand(ComplexF64, 10, 10)))
@test adjoint(B) == B
@test B == B'

# test sparse conversions
@test sparse(M) == sparse(Array(M))
@test convert(SparseMatrixCSC, M) == sparse(Array(M))

B = copy(A)
B[rand(1:length(A), 30)] .= 0
MS = LinearMap(B)
@test sparse(MS) == sparse(Array(MS))

# test function map
F = LinearMap(cumsum, 2)
@test Matrix(F) == [1. 0.; 1. 1.]
@test Array(F) == [1. 0.; 1. 1.]

N = 100
F = LinearMap{ComplexF64}(myft, N) / sqrt(N)
U = Matrix(F) # will be a unitary matrix
@test U'U ≈ Matrix{eltype(U)}(I, N, N)

F = LinearMap(cumsum, 10; ismutating=false)
@test F * v == cumsum(v)
@test *(F, v) == cumsum(v)
@test_throws ErrorException F' * v

F = LinearMap((y,x) -> y .= cumsum(x), 10; ismutating=true)
@test F * v == cumsum(v)
@test *(F, v) == cumsum(v)
@test_throws ErrorException F'v

# Test fallback methods:
L = LinearMap(x -> x, x-> x, 10)
v = randn(10);
@test (2*L)' * v ≈ 2 * v

# test linear combinations
A = 2 * rand(ComplexF64, (10, 10)) .- 1
B = rand(size(A)...)
M = LinearMap(A)
N = LinearMap(B)
v = rand(ComplexF64, 10)

@test Matrix(3 * M) == 3 * A
@test Array(M + A) == 2 * A
@test Matrix(-M) == -A
@test Array(3 * M' - F) == 3 * A' - Array(F)
@test (3 * M - 1im * F)' == 3 * M' + 1im * F'

@test (2 * M' + 3 * I) * v ≈ (2 * A' + 3 * I) * v
@test transpose(LinearMap(2 * M' + 3 * I)) * v ≈ transpose(2 * A' + 3 * I) * v
@test LinearMap(2 * M' + 3 * I)' * v ≈ (2 * A' + 3 * I)' * v

# test composition
@test (F * F) * v == F * (F * v)
@test (F * A) * v == F * (A * v)
@test Matrix(M * transpose(M)) ≈ A * transpose(A)
@test !isposdef(M * transpose(M))
@test isposdef(M * M')
@test issymmetric(N * N')
@test ishermitian(N * N')
@test !issymmetric(M' * M)
@test ishermitian(M' * M)
@test isposdef(transpose(F) * F)
@test isposdef((M * F)' * M * F)
@test transpose(M * F) == transpose(F) * transpose(M)
L = 3 * F + 1im * A + F * M' * F
LF = 3 * Matrix(F) + 1im * A + Matrix(F) * Matrix(M)' * Matrix(F)
@test Array(L) ≈ LF
R1 = rand(ComplexF64, 10, 10)
R2 = rand(ComplexF64, 10, 10)
R3 = rand(ComplexF64, 10, 10)
CompositeR = prod(R -> LinearMap(R), [R1, R2, R3])
Lt = transpose(LinearMap(CompositeR))
@test Lt * v ≈ transpose(R3) * transpose(R2) * transpose(R1) * v
Lc = adjoint(LinearMap(CompositeR))
@test Lc * v ≈ R3' * R2' * R1' * v

# test inplace operations
w = similar(v)
mul!(w, L, v)
@test w ≈ LF * v

# test new type
struct SimpleFunctionMap <: LinearMap{Float64}
    f::Function
    N::Int
end

Base.size(A::SimpleFunctionMap) = (A.N, A.N)
LinearAlgebra.issymmetric(A::SimpleFunctionMap) = false
*(A::SimpleFunctionMap, v::Vector) = A.f(v)
mul!(y::Vector, A::SimpleFunctionMap, x::Vector) = copyto!(y, *(A, x))

F = SimpleFunctionMap(cumsum, 10)
@test ndims(F) == 2
@test size(F, 1) == 10
@test length(F) == 100
w = similar(v)
mul!(w, F, v)
@test w == F * v
@test_throws MethodError F' * v
@test_throws MethodError transpose(F) * v
@test_throws MethodError mul!(w, adjoint(F), v)
@test_throws MethodError mul!(w, transpose(F), v)

# test composition of several maps with shared data #31
sizes = ( (5, 2), (3, 3), (3, 2), (2, 2), (9, 2), (7, 1) )
N = length(sizes) - 1
Lf = []
Lt = []
Lc = []

# build list of operators [LN, ..., L2, L1] for each mode
for (fi, i) in [ (Symbol("f$i"), i) for i in 1:N]
    @eval begin
        function ($fi)(source)
            dest = ones(prod(sizes[$i + 1]))
            tmp = reshape(source, sizes[$i])
            return conj.($i * dest)
        end
        insert!(Lf, 1, LinearMap($fi, prod(sizes[$i + 1]), prod(sizes[$i])))
        insert!(Lt, 1, LinearMap(x -> x, $fi, prod(sizes[$i]), prod(sizes[$i + 1])))
        insert!(Lc, 1, LinearMap{ComplexF64}(x -> x, $fi, prod(sizes[$i]), prod(sizes[$i + 1])))
    end
end

# multiply as composition and as recursion
v1 = ones(prod(sizes[1]))
u1 = ones(prod(sizes[1]))
w1 = im.*ones(ComplexF64, prod(sizes[1]))
for i = N:-1:1
    v2 = prod(Lf[i:N]) * ones(prod(sizes[1]))
    u2 = transpose(LinearMap(prod(Lt[N:-1:i]))) * ones(prod(sizes[1]))
    w2 = adjoint(LinearMap(prod(Lc[N:-1:i]))) * ones(prod(sizes[1]))

    global v1 = Lf[i] * v1
    global u1 = transpose(Lt[i]) * u1
    global w1 = adjoint(Lc[i]) * w1

    @test v1 == v2
    @test u1 == u2
    @test w1 == w2
end


text:
"""
		spikeCorr(arr1, binWidth, T)
		spikeCorr(arr1, arr2, binWidth, T)

Spike timing correlation function. If one array is given calculates autocorrelation if two arrays are given calculates cross-correlation.

# Arguments
arr1=1-dimensional array of spike timings\\
\\
binWidth=size of bin interval in seconds\\
\\
T=Total time overwhich to calulate correlation\\

# Output
corrArr=A 1-dimensional array of correlation coefficients

# Example
```julia-repl
julia> spikeCorr(arr1,0.002, 0.1)
Float64[0.089335, 0.01514, 0.0151317, 0.02228, ...]

```
```julia-repl
julia> spikeCorr(arr1, arr2, 0.002, 0.1)
Float64[0.0893333, 0.0207067, 0.00866, 0.00900833, ...]
```
"""
function spikeCorr(arr1, binWidth, T)
	corrArr = zeros(Int(T/binWidth))

	for i in 1:Int(T/binWidth)
		corrArr[i] = sum(fit(Histogram, arr1, (0:0.002:1200)).weights .*
						 fit(Histogram, arr1 .+ (i-1)*binWidth, 											(0:0.002:1200)).weights)/(1200/0.002)
	end
	return corrArr
end

function spikeCorr(arr1, arr2, binWidth, T)
	corrArr = zeros(Int(T/binWidth))

	for i in 1:Int(T/binWidth)
		corrArr[i] = sum(fit(Histogram, arr1, (0:0.002:1200)).weights .*
						 fit(Histogram, arr2 .+ (i-1)*binWidth, 											(0:0.002:1200)).weights)/(1200/0.002)
	end
	return corrArr
end


text:
using StatsBase
using MXNet

import POMDPs: MDP, Policy

#
# Define a policy that will act as a deep crossentropy policy
#
mutable struct DQNPolicy <: Policy
    action_space::Vector{Int64}
    net::MXNet.mx.FeedForward
    optimizer
    epsilon::Float64
    epsilon_discount::Float64
end

#
# Default uniform intializer based on a games input
#
function DQNPolicy(pomdp::MDP; nnet_spec::Vector{Int64} = [48, 24], nnet_optimizer = mx.ADAM(), epsilon = 0.5, epsilon_discount = 0.999)

    action_count = n_actions(pomdp)
    state_length = n_states(pomdp)

    return DQNPolicy(state_length, action_count, nnet_spec, nnet_optimizer, epsilon, epsilon_discount)
end

#
# Model initializer by specifying number of states, actions and neural network parameters
#
function DQNPolicy(state_length::Int64, action_count::Int64, nnet_spec::Vector{Int64}, nnet_optimizer, epsilon::Float64, epsilon_discount::Float64)

    action_space = collect(0:(action_count-1))

    mlp = @mx.chain mx.Variable(:data)
    for layer_size in nnet_spec
        mlp = mx.FullyConnected(mlp, num_hidden=layer_size)
        mlp = mx.Activation(mlp, act_type=:relu)
    end

    mlp = mx.FullyConnected(mlp, num_hidden=action_count)
    mlp = mx.LinearRegressionOutput(mlp, mx.Variable(:label))

    model = mx.FeedForward(mlp, context = mx.gpu())

    # initialize data provider with identical dataset for all possible outcomes
    # this is done to have initial weight initialization
    random_examples = action_count
    data_provider = mx.ArrayDataProvider(:data => zeros(state_length, random_examples), :label => zeros(action_count, random_examples))
    mx.fit(model, nnet_optimizer, data_provider, initializer = mx.UniformInitializer(0.1), n_epoch = 1, callbacks = [mx.speedometer()])

    DQNPolicy(action_space, model, nnet_optimizer, epsilon, epsilon_discount)
end

#
# Select random action from a policy based on a state and action weight
#
function action(policy::DQNPolicy, s::Vector{Float64})

    # we need to reshape our single state into a state vector of size 1
    # so it is supported by mxnet predict function
    current_values = values(policy, s)
    value, index = findmax(current_values)

    selected_action = if policy.epsilon < rand() index - 1 else rand(policy.action_space) end
    selected_action, current_values
end

function values(policy::DQNPolicy, s::Vector{Float64})
    data_provider = mx.ArrayDataProvider(:data => reshape(s, (length(s), 1)))
    return mx.predict(policy.net, data_provider; verbosity = 0)[:, 1]
end

function update(policy::DQNPolicy, s::Vector{Float64}, v::Vector{Float32})
    data_provider = mx.ArrayDataProvider(:data => repmat(s, 1, length(policy.action_space)), :label => repmat(v, 1, length(policy.action_space)))
    return mx.fit(policy.net, policy.optimizer, data_provider; n_epoch = 3, verbosity = 0)
end

max_stars_repo_name: dmitrijsc/practical-rl

max_issues_repo_name: dmitrijsc/practical-rl

max_forks_repo_name: dmitrijsc/practical-rl


text:
@testset "ch06" begin
    Random.seed!(1234)
    rng = MersenneTwister(1234)
    function create_toy_data(func, n_samples::Int64, std = 1.0, domain = [0.0, 1.0])
        x = collect(range(domain[1], domain[2], length = n_samples))
        shuffle!(rng, x)
        noise = rand(Uniform(-std, std), n_samples)
        return x, func.(x) + noise
    end

    function sinusoidal(x)
        return sin.(2 * pi .* x)
    end

    @testset "ch06_1" begin
        x_train, y_train = create_toy_data(sinusoidal, 10, 0.2)

        x_train = reshape(x_train, 1, size(x_train)[1])
        model = GaussianProcessRegressor(PolynomialKernel(3, 1.0), 1e10)
        fitting(model, x_train, y_train)

        x = reshape(collect(range(0.0, 1.0, length = 100)), 1, 100)
        y, y_sigmas = predict(model, x, true)
    end

    @testset "ch06_2" begin
        x_train, y_train = create_toy_data(sinusoidal, 7, 0.1, [0.0, 0.7])
        x_train = reshape(x_train, 1, size(x_train)[1])
        model = GaussianProcessRegressor(RBFKernel([1.0, 15.0]), 100)
        fitting(model, x_train, y_train)

        x = reshape(collect(range(0.0, 1.0, length = 100)), 1, 100)
        y, y_sigmas = predict(model, x, true)

        y_std = sqrt.(reshape(y_sigmas, length(y_sigmas)))
    end
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
push!(LOAD_PATH, "../src");
push!(LOAD_PATH, ".");

import SocialTest: sdcompare, plotsdcompare, plotRiskvsConc, plotHomovsHetero

#plotsdcompare(trials = 10, simtime = 2000)
plotRiskvsConc()
#plotHomovsHetero()

max_stars_repo_name: josephbakarji/socio-hydro-abm

max_issues_repo_name: josephbakarji/socio-hydro-abm

max_forks_repo_name: josephbakarji/socio-hydro-abm


text:
@testset "Empirical" begin
  # homogeneous field has zero variogram
  sdata = georef((z=ones(3),), Matrix(1.0I, 3, 3))
  γ = EmpiricalVariogram(sdata, :z, nlags=2, maxlag=2.)
  x, y, n = values(γ)
  @test x ≈ [1/2, √2]
  @test y[2] == 0.
  @test n == [0, 3]

  # basic test on number of lags
  sdata = georef((z=[1.,0.,1.],), [25. 50. 75.; 25. 75. 50.])
  γ = EmpiricalVariogram(sdata, :z, nlags=20, maxlag=1.)
  x, y, n = values(γ)
  @test length(x) == 20
  @test length(y) == 20
  @test length(n) == 20

  # empirical variogram on integer coordinates
  sdata = georef((z=ones(3),), Matrix(1I, 3, 3))
  γ = EmpiricalVariogram(sdata, :z, nlags=2, maxlag=2, algo=:full)
  x, y, n = values(γ)
  @test x ≈ [1/2, √2]
  @test y[2] == 0.
  @test n == [0, 3]

  # empirical variogram with only missing data
  X = rand(2,3)
  z = Union{Float64,Missing}[missing,missing,missing]
  𝒟 = georef((z=z,), X)
  γ = EmpiricalVariogram(𝒟, :z, maxlag=1., nlags=5)
  x, y, n = values(γ)
  @test x == [.1, .3, .5, .7, .9]
  @test all(iszero.(n))

  # accumulation algorithms give the same result
  Random.seed!(2021)
  sdata = georef((z=rand(1000),), rand(3,1000))
  γ₁ = EmpiricalVariogram(sdata, :z, maxlag=0.01, algo=:full)
  γ₂ = EmpiricalVariogram(sdata, :z, maxlag=0.01, algo=:ball)
  @test isequal(values(γ₁), values(γ₂))

  # custom distance is recorded
  sdata = georef((z=rand(1000),), rand(2,1000))
  γ = EmpiricalVariogram(sdata, :z, distance=Haversine(6371.), algo=:full)
  @test distance(γ) == Haversine(6371.)

  # print methods
  rng = MersenneTwister(123)
  d = georef((z=rand(rng,100,100),))
  γ = EmpiricalVariogram(d, :z)
  @test sprint(show, γ) == "EmpiricalVariogram"
  @test sprint(show, MIME"text/plain"(), γ) == "EmpiricalVariogram\n  abscissa: (0.3535533905932738, 13.84261778461484)\n  ordinate: (0.0, 0.08323850196902784)\n  N° pairs: 2790126"

  if visualtests
    wl = geostatsimage("WalkerLake")
    TI = asarray(wl, :Z)[1:20,1:20]
    d = georef((z=TI,))
    γ = EmpiricalVariogram(d, :z, maxlag=15.)
    @test_reference "data/empirical.png" plot(γ)
  end

  # test variography with compositional data
  data = georef((z=rand(Composition{3}, 100),), rand(2, 100))
  γ = EmpiricalVariogram(data, :z, maxlag = 1.0, algo=:full)
  x, y, n = values(γ)
  @test all(≥(0), x)
  @test all(≥(0), y)
  @test all(>(0), n)

  # test variography with unitful data
  data = georef((z=[1*u"K" for i in 1:100],), rand(2, 100))
  γ = EmpiricalVariogram(data, :z, nlags=20)
  x, y, n = values(γ)
  @test all(≥(0), x)
  @test y == fill(0.0*u"K^2", 20)
end


text:
using Libdl

const build_loopback = true

# Build SoapyLoopback
# TODO: Yggdrasil once stable
if build_loopback
    include("build_tarballs.jl")
end
cd("products")
loopback_tar = readdir()[1]
loopback = splitext(loopback_tar)[1]
run(`tar -xzf $(loopback_tar)`) # BB without tar output?
dlopen("lib/SoapySDR/modules0.8/libsoapyloopback.so")


text:
using Zlib
using Libz

abstract Compressor

# Min length of input str for compression to happen
minlen(c::Compressor) = c.minlen

compress{S <: AbstractString}(c::Compressor, str::S) = throw("Missing compress function")

clen{S <: AbstractString}(c::Compressor, str::S) = length(compress(c, str))
clen{S <: AbstractString}(c::Compressor, strs::Vector{S}) = clen(c, join(strs))
clen{S <: AbstractString}(c::Compressor, x::S, y::S) = clen(c, x * y)

type ZlibCompressor <: Compressor
  minlen::Int
  level::Int
  # No compression if less than 5 chars.
  ZlibCompressor(; level = 9) = new(5, level)
end

compress{S <: AbstractString}(zc::ZlibCompressor, str::S) = Zlib.compress(str, zc.level, false, true)

# Some compressors have streaming API's and then we can use that to speed things
# up
abstract StreamCompressor <: Compressor

# Create and return one IO stream for output and one compression stream (for adding input).
# Note that the IO stream is for internal use only so can be reused between multiple calls.
newstreams(sc::StreamCompressor) = raise("Subtypes have not defined a newstreams function")
addtostream(cs, str) = raise("Subtypes have not defined a addtostream function")

function compress{S <: AbstractString}(sc::StreamCompressor, str::S)
	ios, cs = newstreams(sc)
	addtostream(cs, str)
	close(cs)
	takebuf_string(ios)
end

function clen{S <: AbstractString}(sc::StreamCompressor, str::S)
	ios, cs = newstreams(sc)
	addtostream(cs, str)
	close(cs)
	position(ios)
end

function clen{S <: AbstractString}(sc::StreamCompressor, x::S, y::S)
	ios, cs = newstreams(sc)
	addtostream(cs, x)
	addtostream(cs, y)
	close(cs)
	position(ios)
end

function clen{S <: AbstractString}(sc::StreamCompressor, strs::Vector{S})
	ios, cs = newstreams(sc)
	for str in strs
		addtostream(cs, str)
	end
	close(cs)
	position(ios)
end

# Libz is a faster (than Zlib) julia wrapper of Zlib
type LibzCompressor <: StreamCompressor
  minlen::Int
  level::Int
  iobf::IOBuffer
  LibzCompressor(; level = 9) = new(0, level, IOBuffer())
end

function newstreams(lzc::LibzCompressor)
	# Reset the io buffer to its first position
	seek(lzc.iobf, 0)
	zs = ZlibDeflateOutputStream(lzc.iobf; level = lzc.level, gzip = false) # Don't write gzip header
	return lzc.iobf, zs
end

addtostream{S <: AbstractString}(cs, str::S) = write(cs, str)

# Convenience function for finding the minlen for a compressor
function find_min_len_for_compression(compressor, seedStr = "a", maxReps = 1000)
  prevclen = currclen = clen(compressor, seedStr)
  for i in 1:maxReps
    prevclen = currclen
    str = seedStr ^ i
    currclen = clen(compressor, str)
    if currclen < prevclen
      return i, length(str), currclen, length(str)-length(seedStr), clen(compressor, seedStr^(i-1))
    end
  end
  return -1, str
end

###############################################################################
# Only load Blosc-related functions if the Blosc.jl package is installed
###############################################################################
if false # isinstalled("Blosc") # exclude Blosc for now since it is not currently used because of its long min length before it compresses

using Blosc

blosc_compress(s, level = 9, shuffle = false, compressor = "blosclz") = begin
  Blosc.set_compressor(compressor)
  Blosc.compress(s; level = level, shuffle = shuffle)
end

type BloscCompressor <: Compressor
  minlen::Int
  level::Int
  compressorName
  # Limit for Blosc compression to start seems to be 128 chars.
  BloscCompressor(name; minlen = 128, level = 9) = new(minlen, level, name)
end

compress{S <: AbstractString}(bc::BloscCompressor, str::S) = blosc_compress(str, bc.level,
	false, bc.compressorName)

compress_str(c::BloscCompressor, str) = blosc_compress(str, c.level, false, c.compressorName)

# Define specific compressors which can be used directly
ZlibC = ZlibCompressor()
BloscLz = BloscCompressor("blosclz")
Lz4 = BloscCompressor("lz4")
Lz4hc = BloscCompressor("lz4hc")
Snappy = BloscCompressor("snappy")
BloscZlib = BloscCompressor("zlib")

QLz4hc = BloscCompressor("lz4hc"; level = 1)
QSnappy = BloscCompressor("snappy"; level = 1)

zliblength(x) = clen(ZlibC, x)

end # end if isinstalled("Blosc")


text:
__precompile__(true)

module AOPE

using Dates, BusinessDays, Roots
# using LinearAlgebra
using Reexport
@reexport using Dates
@reexport using Roots

export
# Types
  Following,           # bd-conventions.jl
  ModifiedFollowing,   # bd-conventions.jl
  Preceding,           # bd-conventions.jl
  ModifiedPreceding,   # bd-conventions.jl
  EndOfMonth,          # bd-conventions.jl
  AllBD,               # bd-conventions.jl
  ActualActualISDA,    # dc-conventions.jl
  ActualActualICMA,    # dc-conventions.jl
  ActualActualISMA,    # bd-conventions.jl
  ActualActualExcel,   # dc-conventions.jl
  Actual365Fixed,      # dc-conventions.jl
  Actual360,           # dc-conventions.jl
  Thirty360,           # dc-conventions.jl
  ThirtyE360,          # dc-conventions.jl
  ThirtyE360ISDA,      # dc-conventions.jl
  Thirty360Excel,      # dc-conventions.jl
  Discrete,            # AOPE.jl
  Continuous,          # AOPE.jl

# Functions
  #transpose,   # LinearAlgebra.jl
  #det,         # LinearAlgebra.jl
  #diag,        # LinearAlgebra.jl
  #pinv,        # LinearAlgebra.jl
  liqd,         # bd-conventions.jl
  frac,         # dc-conventions.jl
  utild,        # bd-conventions.jl
  bdarray,      # bd-conventions.jl
  nfcaixa,      # valuation.jl
  fcaixa,       # valuation.jl
  spreço,       # valuation.jl
  yca,          # yields.jl
  ytm,          # yields.jl
  ytmisda,      # yields.jl
  zspread       # yields.jl

# Operators
  #\,           # LinearAlgebra.jl
  #/,           # LinearAlgebra.jl

  # Constants
  #I            # LinearAlgebra.jl

abstract type Composition end

struct Discrete <: Composition end
struct Continuous <: Composition end

const D = Discrete
const C = Continuous
#const uteis = bd
#const spreço = bondprice
#const ysimples = simpleyield
#const nfcaixa = nfc
#const fcaixa = fc

include("dc-conventions.jl")
include("excel-conventions.jl")
include("bd-conventions.jl")
include("valuation.jl")
include("yields.jl")

println("\nMestrado Análise Financeira: AOPE|2019 - versão 0.1.0 <2019-03-05>\n")

end # module AOPE


text:
# Block-coordinate and incremental aggregated proximal
# gradient methods for nonsmooth nonconvex problems
# Algorithm 2

struct FINITO_basic_iterable{R<:Real,C<:Union{R,Complex{R}},Tx<:AbstractArray{C},Tf,Tg} <: CIAO_iterable
    F::Array{Tf}            # smooth term
    g::Tg                   # nonsmooth term
    x0::Tx                  # initial point
    N::Int                  # number of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}}  # Lipschitz moduli of nabla f_i
    γ::Maybe{Union{Array{R},R}}  # stepsizes
    sweeping::Int8          # to only use one stepsize γ
    batch::Int              # batch size
    α::R                    # in (0, 1), e.g.: 0.99
end

mutable struct FINITO_basic_state{R<:Real,Tx}
    s::Array{Tx}            # table of x_j- γ_j/N nabla f_j(x_j)
    γ::Array{R}             # stepsize parameters
    hat_γ::R                # average γ
    av::Tx                  # the running average
    z::Tx
    ind::Array{Array{Int}} # running index set
    # some extra placeholders
    d::Int                  # number of batches
    ∇f_temp::Tx             # placeholder for gradients
    idxr::Int               # running idx in the iterate
    idx::Int                # location of idxr in 1:N
    inds::Array{Int}        # needed for shuffled only
end

function FINITO_basic_state(s, γ, hat_γ::R, av::Tx, z::Tx, ind, d) where {R,Tx} #? why it is not inside?
    return FINITO_basic_state{R,Tx}(
        s,
        γ,
        hat_γ,
        av,
        z,
        ind,
        d,
        zero(av),
        Int(1),
        Int(0),
        collect(1:d),
    )
end

function Base.iterate(iter::FINITO_basic_iterable{R,C,Tx}) where {R,C,Tx}
    N = iter.N
    # define the batches
    r = iter.batch # batch size
    # create index sets
    if iter.sweeping == 1
        ind = [collect(1:r)] # placeholder
    else
        ind = Vector{Vector{Int}}(undef, 0)
        d = Int(floor(N / r))
        for i = 1:d
            push!(ind, collect(r*(i-1)+1:i*r))
        end
        r * d < N && push!(ind, collect(r*d+1:N))
    end
    d = cld(N, r) # number of batches
    # updating the stepsize
    if iter.γ === nothing
        if iter.L === nothing
            @warn "--> smoothness parameter absent"
            return nothing
        else
            γ = zeros(R, N)
            for i = 1:N
                isa(iter.L, R) ? (γ = fill(iter.α * R(iter.N) / iter.L, (N,))) :
                (γ[i] = iter.α * R(N) / iter.L[i])
            end
        end
    else
        isa(iter.γ, R) ? (γ = fill(iter.γ, (N,))) : (γ = iter.γ) #provided γ
    end
    # computing the gradients and updating the table s
    s = Vector{Tx}(undef, 0)
    for i = 1:N
        ∇f, ~ = gradient(iter.F[i], iter.x0)
        push!(s, iter.x0 - γ[i] / N * ∇f) # table of x_i
    end
    #initializing the vectors
    hat_γ = 1 / sum(1 ./ γ)
    av = hat_γ * (sum(s ./ γ)) # the running average
    z, ~ = prox(iter.g, av, hat_γ)

    state = FINITO_basic_state(s, γ, hat_γ, av, z, ind, d)

    return state, state
end

# Block-coordinate and incremental aggregated proximal
# gradient methods for nonsmooth nonconvex problems
# Algorithm 2

function Base.iterate(
    iter::FINITO_basic_iterable{R},
    state::FINITO_basic_state{R},
) where {R}
    # manipulating indices
    # if iter.sweeping == 1 # uniformly random
    #     state.ind = [sample(1:iter.N, iter.batch, replace = false)]
    # elseif iter.sweeping == 2  # cyclic
    #     state.idxr = mod(state.idxr, state.d) + 1
    # elseif iter.sweeping == 3  # shuffled cyclic
    #     if state.idx == state.d
    #         state.inds = randperm(state.d)
    #         state.idx = 1
    #     else
    #         state.idx += 1
    #     end
    #     state.idxr = state.inds[state.idx]
    # end

    for j = 1:iter.N
        state.idxr = mod(state.idxr,iter.N)+1
        i = state.idxr
        gradient!(state.∇f_temp, iter.F[i], state.z)
        state.∇f_temp .*= -(state.γ[i] / iter.N)
        state.∇f_temp .+= state.z
        @. state.av += (state.∇f_temp - state.s[i]) * (state.hat_γ / state.γ[i])
        state.s[i] .= state.∇f_temp  #update x_i

        # the iterate
        # for i in collect(1:iter.N)
        #     # perform the main steps
        #     gradient!(state.∇f_temp, iter.F[i], state.z)
        #     state.∇f_temp .*= -(state.γ[i] / iter.N)
        #     state.∇f_temp .+= state.z
        #     @. state.av += (state.∇f_temp - state.s[i]) * (state.hat_γ / state.γ[i])
        #     state.s[i] .= state.∇f_temp  #update x_i
        # end

        prox!(state.z, iter.g, state.av, state.hat_γ)
    end

    return state, state
end

solution(state::FINITO_basic_state) = state.z

# count(state::FINITO_basic_state) = []

#TODO list
## in cyclic/shuffled minibatchs are static


text:
module IO
module Client
    using HTTP
    using DataStructures

    function ws_msg(msg="Hello!", wss="ws://164.90.144.33:8080")
        HTTP.WebSockets.open(wss) do ws
            write(ws, msg)
            x = readavailable(ws)
            println(String(x))
        end;
    end
end # module

module Server
    using HTTP

    @async HTTP.WebSockets.listen("164.90.144.33", UInt16(8080)) do ws
        while !eof(ws)
            data = readavailable(ws)
            write(ws, data)
        end
    end
end # module

end # module

max_stars_repo_name: dustyirwin/PlaymatSimulator

max_issues_repo_name: dustyirwin/PlaymatSimulator

max_forks_repo_name: dustyirwin/PlaymatSimulator


text:
"""
    A = JopDerivative(R::JetSpace[; dim=1, accuracy=4, delta=1.0)

`A*m` is the centered finite different approximation to the derivative of `m`
along dimension `dim`.  The `accuracy` of the approximation is either `4` for
a 4th order accurate estimate of the derivative, or `8` for an 8th order accurate
derivative.
"""
function JopDerivative(sp::JetSpace{T}; dim=1, accuracy=4, delta=1.0) where {T}
    local a
    if accuracy == 4
        a = T[1/12, -2/3, 0, 2/3, -1/12]
    elseif accuracy == 8
        a = T[1/280, -4/105, 1/5, -4/5, 0, 4/5, -1/5, 4/105, -1/280]
    else
        error("order of accuracy must be either 4 or 8, got accuracy=$(accuracy)")
    end
    a ./= delta
    JopLn(df! = JopDerivative_df!, df′! = JopDerivative_df′!, dom = sp, rng = sp, s = (a=a,dim=dim))
end
export JopDerivative

JopDerivative_df!(d::AbstractArray, m::AbstractArray; a, dim, kwargs...) = derivative_helper(d, m, a, dim)
JopDerivative_df′!(m::AbstractArray, d::AbstractArray; a, dim, kwargs...) = derivative_helper(m, d, -a, dim)

function derivative_helper_dim(d, m, a, n, Rpre, Rpost)
    d .= 0
    M = div(length(a)-1,2)
    for Ipost in Rpost
        for i = (M+1):(n-M)
            for j = 1:length(a)
                k = i-M+j-1
                if M < k <= n-M
                    for Ipre in Rpre
                        d[Ipre, i, Ipost] += a[j] * m[Ipre, k, Ipost]
                    end
                end
            end
        end
    end
    d
end
function derivative_helper(d::AbstractArray, m::AbstractArray, a::Vector, dim::Integer)
    Rpre = CartesianIndices(size(m)[1:dim-1])
    Rpost = CartesianIndices(size(m)[dim+1:end])
    n = size(m, dim)
    derivative_helper_dim(d, m, a, n, Rpre, Rpost)
end


text:
"""
	SeisDiff(d ; <keyword arguments>)

Apply differentiation to seismic traces in freq. It can also be
used to apply a phase rotation

# Arguments
- `d`: Input 2D data array in tx domain. First dimension is time.

# Keyword arguments
- `delay=0.1`: desired time delay in seconds.
- `pow=-2`: order of derivative
- `rot=0`: constant phase shift or rotation

# Example
```
julia> d = SeisLinearEvents(); SeisPlotTX(d);
julia> d2 = SeisDiff(d;dt=0.004,pow=1); SeisPlotTX(d);
```
"""
function SeisDiff(d;dt=0.004,pow=-2,rot=0)

	nt = size(d,1)
	D = fft(d,1)
	dw = 2. * pi/nt/dt
	nw = convert(Int,floor(nt/2)) + 1
	eps = pow < 0 ? 10.0 : 0.0
	for iw=1:nw
		D[iw,:] *= exp(rot*1im)*((iw*dw) + eps)^pow
	end
	# symmetries
	for iw=nw+1:nt
		D[iw,:] = conj(D[nt-iw+2,:])
	end
	d = real(ifft(D,1))
	return d

end


text:
module Day1

export parse_input
export part_1
export part_2

function parse_input(input::String)
	nums = split(input)
	return map(x -> parse(Int, x), nums)
end

function part_1(numbers::Array{Int})
	return sum(numbers)
end

function part_2(numbers::Array{Int})
	current_frequency = 0
	visited_frequencies = Set()
	idx = 1
	while !in(current_frequency, visited_frequencies)
		push!(visited_frequencies, current_frequency)
		current_frequency += numbers[idx]
		idx = mod1(idx + 1, length(numbers))
	end
	return current_frequency
end

end


text:
## -------------  Ranking students for sequential college admissions

"""
	$(SIGNATURES)

Abstract type for switches governing student ranking.
"""
abstract type AbstractRankingSwitches{F1 <: Real} end

"""
	$(SIGNATURES)

Abstract student ranking type.
"""
abstract type AbstractRanking{F1 <: Real} <: ModelObject end

# StructLH.describe(x :: AbstractRanking) = nothing;
StructLH.describe(switches :: AbstractRankingSwitches) =
    ["Generic student ranking"];

Lazy.@forward AbstractRanking.switches (
    StructLH.describe, high_draws_first, calibrate_weights
)

"""
	$(SIGNATURES)

Initialize an `AbstractRanking` object from its switches.
"""
function make_student_ranking end

"""
	$(SIGNATURES)

Retrieve draws for endowment `eName`. Needs to be defined for objects passed into ranking functions.
"""
retrieve_draws(draws, eName) =
    error("Caller must define `retrieve_draws` for input $(typeof(draws))");

"""
	$(SIGNATURES)

Number of individuals in endowment draws. Needs to be defined for objects passed into ranking functions.
"""
n_draws(draws) =
    error("Caller must define `n_draws` for input $(typeof(draws))");

"""
	$(SIGNATURES)

Rank students. Returns indices of students in rank order from best to worst.

# Arguments
- `draws`: must support `retrieve_draws`.
"""
function rank_students(e :: AbstractRanking{F1}, draws) where F1
    return sortperm(score_students(e, draws), rev = high_draws_first(e));
end

"""
	$(SIGNATURES)

Score students. Higher scores are better.
Scores are simply weighted sums of endowment draws.
Not scaled. Scores can be positive or negative.
"""
function score_students(e :: AbstractRanking{F1}, draws) where F1
    scoreV = zeros(F1, n_draws(draws));
    wtV = weights(e);
    nameV = endow_names(e);
    for (j, eName) in enumerate(nameV)
        scoreV .+= retrieve_draws(draws, eName) .* wtV[j];
    end
    # scaled  &&  scale_scores(e, scoreV);
    return scoreV
end

"""
	$(SIGNATURES)

Range of scores, given ranges for the scoring variables.
"""
function range_of_scores(e :: AbstractRanking{F1};
    lbV = lower_bounds(e), ubV = upper_bounds(e)) where F1
    wtV = weights(e);
    lb = sum(min.(wtV .* lbV, wtV .* ubV));
    ub = sum(max.(wtV .* lbV, wtV .* ubV));
    return lb, ub
end

"""
	$(SIGNATURES)

Scale scores to lie in [0, 1].
"""
function scale_scores(e :: AbstractRanking{F1}, scoreV;
    lbV = lower_bounds(e),
    ubV = upper_bounds(e)) where F1

    lb, ub = range_of_scores(e; lbV, ubV);
    scoreV = (scoreV .- lb) ./ (ub .- lb);
    return scoreV
end

"""
	$(SIGNATURES)

Lower bounds of scoring variables.
"""
function lower_bounds(e :: AbstractRanking{F1}) where F1 end

"""
	$(SIGNATURES)

Upper bounds of scoring variables.
"""
function upper_bounds(e :: AbstractRanking{F1}) where F1 end

"""
	$(SIGNATURES)

Return endowment names.
"""
endow_names(e :: AbstractRankingSwitches{F1}) where F1 = e.eNameV;

"""
	$(SIGNATURES)

Return endowment names.
"""
endow_names(e :: AbstractRanking{F1}) where F1 = endow_names(e.switches);

"""
	$(SIGNATURES)

Validate an `AbstractRanking`.
"""
function validate_ranking end

"""
	$(SIGNATURES)

Validate switches for an `AbstractRanking`.
"""
function validate_ranking_switches end

# -----------------


text:
type Moments <: ContinuousMultivariateStreamStat
    m1::Float64
    m2::Float64
    m3::Float64
    m4::Float64
    n::Int
end

Moments() = Moments(0.0, 0.0, 0.0, 0.0, 0)

function update!(stat::Moments, x::Real)
    stat.n += 1

    n = stat.n

    δ = x - stat.m1
    δ_n = δ / stat.n
    δ_n_sq = δ_n * δ_n
    term1 = δ * δ_n * (n - 1)

    stat.m1 += δ_n
    stat.m4 += term1 * δ_n_sq * (n * n - 3 * n + 3) +
               6 * δ_n_sq * stat.m2 - 4 * δ_n * stat.m3
    stat.m3 += term1 * δ_n * (n - 2) - 3 * δ_n * stat.m2
    stat.m2 += term1

    return
end

function Base.mean(stat::Moments)
    return stat.m1
end

function Base.var(stat::Moments)
    return stat.m2 / (stat.n - 1)
end

function StatsBase.skewness(stat::Moments)
    return sqrt(stat.n) * stat.m3 / stat.m2^1.5
end

function StatsBase.kurtosis(stat::Moments)
    return stat.n * stat.m4 / (stat.m2 * stat.m2) - 3.0
end

nobs(stat::Moments) = stat.n

function state(stat::Moments)
    return (mean(stat), var(stat), StatsBase.skewness(stat), StatsBase.kurtosis(stat))
end

Base.copy(stat::Moments) = Moments(stat.m1, stat.m2, stat.m3, stat.m4, stat.n)

function Base.merge(a::Moments, b::Moments)
    merged = Moments()

    merged.n = a.n + b.n

    δ = b.m1 - a.m1
    δ2 = δ * δ
    δ3 = δ * δ2
    δ4 = δ2 * δ2

    merged.m1 = (a.n * a.m1 + b.n * b.m1) / merged.n

    merged.m2 = a.m2 + b.m2 + δ2 * a.n * b.n / merged.n

    merged.m3 = a.m3 + b.m3 +
                  δ3 * a.n * b.n * (a.n - b.n) / (merged.n * merged.n)
    merged.m3 += 3.0 * δ * (a.n * b.m2 - b.n * a.m2) / merged.n

    merged.m4 = a.m4 + b.m4 + δ4 * a.n * b.n *
                  (a.n * a.n - a.n * b.n + b.n * b.n) /
                  (merged.n * merged.n * merged.n)
    merged.m4 += 6.0 * δ2 * (a.n * a.n * b.m2 + b.n * b.n * a.m2) /
                   (merged.n * merged.n) +
                   4.0 * δ * (a.n * b.m3 - b.n * a.m3) / merged.n

    return merged
end

function Base.empty!(stat::Moments)
    stat.n = 0
    stat.m1 = 0.0
    stat.m2 = 0.0
    stat.m3 = 0.0
    stat.m4 = 0.0
    return
end

function Base.show(io::IO, stat::Moments)
    m, v, s, k = state(stat)
    @printf(io, "Online Moments\n")
    @printf(io, " * Mean:     %f\n", m)
    @printf(io, " * Variance: %f\n", v)
    @printf(io, " * Skewness: %f\n", s)
    @printf(io, " * Kurtosis: %f\n", k)
    @printf(io, " * N:        %d\n", nobs(stat))
    return
end

Base.mean(stat::Moments) = stat.m1

Base.var(stat::Moments) = stat.m2 / (stat.n - 1)

Base.std(stat::Moments) = sqrt(var(stat))

StatsBase.skewness(stat::Moments) = sqrt(stat.n) * stat.m3 / stat.m2^1.5

StatsBase.kurtosis(stat::Moments) = stat.n * stat.m4 / (stat.m2 * stat.m2) - 3.0


text:
include("util.jl")
include("math.jl")
include("rand.jl")
include("lalg.jl")
include("blockinv.jl")


text:
function visualize(states::AbstractArray, states2::AbstractArray; statenames = string.(1:length(states)), statenames2 = string.(1:length(states2)), aspect = false, resolution = (2880, 1080), statistics = false, title = "Field = ", title2 = "Field = ", units1 = ["" for i in eachindex(states)], units2 = ["" for i in eachindex(states)])
    # Create scene
    scene, layout = layoutscene(resolution = resolution)

    lscene = layout[2:4, 2:4] = LScene(scene)
    lscene2 = layout[2:4, 5:7] = LScene(scene)
    width = round(Int, resolution[1] / 6) # make menu 1/4 of preliminary resolution

    # Create choices and nodes
    stateindex = collect(1:length(states))
    statenode = Node(stateindex[1])

    colorchoices = [:balance, :thermal, :dense, :deep, :curl, :thermometer]
    colornode = Node(colorchoices[1])

    stateindex2 = collect(1:length(states2))
    statenode2 = Node(stateindex2[1])

    colorchoices2 = [:balance, :thermal, :dense, :deep, :curl, :thermometer]
    colornode2 = Node(colorchoices2[1])

    if statistics
        llscene = layout[4,1] = Axis(scene, xlabel = @lift(statenames[$statenode] * units1[$statenode]),
                         xlabelcolor = :black, ylabel = "pdf",
                         ylabelcolor = :black, xlabelsize = 40, ylabelsize = 40,
                         xticklabelsize = 25, yticklabelsize = 25,
                         xtickcolor = :black, ytickcolor = :black,
                         xticklabelcolor  = :black, yticklabelcolor = :black)
        layout[3, 1] = Label(scene, "Statistics", width = width, textsize = 50)

        llscene2 = layout[4,8] = Axis(scene, xlabel = @lift(statenames2[$statenode2] * units2[$statenode2]),
                         xlabelcolor = :black, ylabel = "pdf",
                         ylabelcolor = :black, xlabelsize = 40, ylabelsize = 40,
                         xticklabelsize = 25, yticklabelsize = 25,
                         xtickcolor = :black, ytickcolor = :black,
                         xticklabelcolor  = :black, yticklabelcolor = :black)
        layout[3, 8] = Label(scene, "Statistics", width = width, textsize = 50)
    end

    # x,y,z are for determining the aspect ratio of the box
    if (typeof(aspect) <: Tuple) & (length(aspect) == 3)
        x, y, z = aspect
    else
        x, y, z = size(states[1])
    end

    # Clim sliders
    upperclim_slider = Slider(scene, range = range(0, 1, length = 101), startvalue = 0.99)
    upperclim_node = upperclim_slider.value
    lowerclim_slider = Slider(scene, range = range(0, 1, length = 101), startvalue = 0.01)
    lowerclim_node = lowerclim_slider.value

    upperclim_slider2 = Slider(scene, range = range(0, 1, length = 101), startvalue = 0.99)
    upperclim_node2 = upperclim_slider2.value
    lowerclim_slider2 = Slider(scene, range = range(0, 1, length = 101), startvalue = 0.01)
    lowerclim_node2 = lowerclim_slider2.value

    # Lift Nodes
    state = @lift(states[$statenode])
    statename = @lift(statenames[$statenode])
    clims = @lift((quantile($state[:], $lowerclim_node) , quantile($state[:], $upperclim_node)))
    cmap_rgb = @lift(to_colormap($colornode))
    titlename = @lift(title * $statename) # use padding and appropriate centering

    state2 = @lift(states2[$statenode2])
    statename2 = @lift(statenames2[$statenode2])
    clims2 = @lift((quantile($state2[:], $lowerclim_node2) , quantile($state2[:], $upperclim_node2)))
    cmap_rgb2 = @lift(to_colormap($colornode2))
    titlename2 = @lift(title2 * $statename2) # use padding and appropriate centering

    # Statistics
    if statistics
        histogram_node = @lift(histogram($state, bins = 300))
        xs = @lift($histogram_node[1])
        ys = @lift($histogram_node[2])
        pdf = GLMakie.AbstractPlotting.barplot!(llscene, xs, ys, color = :red,
                        strokecolor = :red,
                        strokewidth = 1)
        @lift(GLMakie.AbstractPlotting.xlims!(llscene, extrema($state)))
        @lift(GLMakie.AbstractPlotting.ylims!(llscene, extrema($histogram_node[2])))
        vlines!(llscene, @lift($clims[1]), color = :black, linewidth = width / 100)
        vlines!(llscene, @lift($clims[2]), color = :black, linewidth = width / 100)
        # 2
        histogram_node2 = @lift(histogram($state2, bins = 300))
        xs2 = @lift($histogram_node2[1])
        ys2 = @lift($histogram_node2[2])
        pdf2 = GLMakie.AbstractPlotting.barplot!(llscene2, xs2, ys2, color = :red,
                        strokecolor = :red,
                        strokewidth = 1)
        @lift(GLMakie.AbstractPlotting.xlims!(llscene2, extrema($state2)))
        @lift(GLMakie.AbstractPlotting.ylims!(llscene2, extrema($histogram_node2[2])))
        vlines!(llscene2, @lift($clims2[1]), color = :black, linewidth = width / 100)
        vlines!(llscene2, @lift($clims2[2]), color = :black, linewidth = width / 100)
    end

    # Volume Plot
    volume!(lscene, 0..x, 0..y, 0..z, state,
            camera = cam3d!,
            colormap = cmap_rgb,
            colorrange = clims)
    volume!(lscene2, 0..x, 0..y, 0..z, state2,
        camera = cam3d!,
        colormap = cmap_rgb2,
        colorrange = clims2)
    # Title
    supertitle = layout[1, 2:4] = Label(scene, titlename , textsize = 50, color = :black)
    supertitle2 = layout[1, 5:7] = Label(scene, titlename2 , textsize = 50, color = :black)

    # Menus 1
    statemenu = Menu(scene, options = zip(statenames, stateindex))
    on(statemenu.selection) do s
        statenode[] = s
    end

    colormenu = Menu(scene, options = zip(colorchoices, colorchoices))
    on(colormenu.selection) do s
        colornode[] = s
    end
    lowerclim_string = @lift("quantile = " *  @sprintf("%0.2f", $lowerclim_node) * ", value = " * @sprintf("%0.1e", $clims[1]))
    upperclim_string = @lift("quantile = " *  @sprintf("%0.2f", $upperclim_node) * ", value = " * @sprintf("%0.1e", $clims[2]))
    # depends on makie version, vbox for old, vgrid for new
    layout[2, 1] = vgrid!(
        Label(scene, "State", width = nothing),
        statemenu,
        Label(scene, "Color", width = nothing),
        colormenu,
        Label(scene, lowerclim_string, width = nothing),
        lowerclim_slider,
        Label(scene, upperclim_string, width = nothing),
        upperclim_slider,
    )
    layout[1,1] = Label(scene, "Menu", width = width, textsize = 50)

    # Menus 2
    statemenu2 = Menu(scene, options = zip(statenames2, stateindex2))
    on(statemenu2.selection) do s
        statenode2[] = s
    end

    colormenu2 = Menu(scene, options = zip(colorchoices2, colorchoices2))
    on(colormenu2.selection) do s
        colornode2[] = s
    end
    lowerclim_string2 = @lift("quantile = " *  @sprintf("%0.2f", $lowerclim_node2) * ", value = " * @sprintf("%0.1e", $clims2[1]))
    upperclim_string2 = @lift("quantile = " *  @sprintf("%0.2f", $upperclim_node2) * ", value = " * @sprintf("%0.1e", $clims2[2]))
    # depends on makie version, vbox for old, vgrid for new
    layout[2, 8] = vgrid!(
        Label(scene, "State", width = nothing),
        statemenu2,
        Label(scene, "Color", width = nothing),
        colormenu2,
        Label(scene, lowerclim_string2, width = nothing),
        lowerclim_slider2,
        Label(scene, upperclim_string2, width = nothing),
        upperclim_slider2,
    )
    layout[1,8] = Label(scene, "Menu", width = width, textsize = 50)

    ## Axis
        # Modify Axis
    axis = scene.children[1][OldAxis]
    axis[:names][:axisnames] = ("↓ Zonal [m] ", "Meriodonal [m]↓ ", "Depth [m]↓ ")
    axis[:names][:align] = ((:left, :center), (:right, :center), (:right, :center))
    # need to adjust size of ticks first and then size of axis names
    axis[:names][:textsize] = (50.0, 50.0, 50.0)
    axis[:ticks][:textsize] = (00.0, 00.0, 00.0)
    # axis[:ticks][:ranges_labels].val # current axis labels
    xticks = collect(range(-0, aspect[1], length = 2))
    yticks = collect(range(-0, aspect[2], length = 6))
    zticks = collect(range(-0, aspect[3], length = 2))
    ticks = (xticks, yticks, zticks)
    axis[:ticks][:ranges] = ticks
    xtickslabels = [@sprintf("%0.1f", (xtick)) for xtick in xticks]
    xtickslabels[end] = "1e6"
    ytickslabels = ["", "south","", "", "north", ""]
    ztickslabels = [@sprintf("%0.1f", (xtick)) for xtick in xticks]
    labels = (xtickslabels, ytickslabels, ztickslabels)
    axis[:ticks][:labels] = labels

        # Modify Axis
    axis = scene.children[2][OldAxis]
    axis[:names][:axisnames] = ("↓ Zonal [m] ", "Meriodonal [m]↓ ", "Depth [m]↓ ")
    axis[:names][:align] = ((:left, :center), (:right, :center), (:right, :center))
    # need to adjust size of ticks first and then size of axis names
    axis[:names][:textsize] = (50.0, 50.0, 50.0)
    axis[:ticks][:textsize] = (00.0, 00.0, 00.0)
    # axis[:ticks][:ranges_labels].val # current axis labels
    xticks = collect(range(-0, aspect[1], length = 2))
    yticks = collect(range(-0, aspect[2], length = 6))
    zticks = collect(range(-0, aspect[3], length = 2))
    ticks = (xticks, yticks, zticks)
    axis[:ticks][:ranges] = ticks
    xtickslabels = [@sprintf("%0.1f", (xtick)) for xtick in xticks]
    xtickslabels[end] = "1e6"
    ytickslabels = ["", "south","", "", "north", ""]
    ztickslabels = [@sprintf("%0.1f", (xtick)) for xtick in xticks]
    labels = (xtickslabels, ytickslabels, ztickslabels)
    axis[:ticks][:labels] = labels

    display(scene)
    return scene
end

function grabtitle(filename)
    resolutionnumber = parse(Int64, split(filename, "_")[end-2])
    return string(round(Int, 1000 * 16 / 192 / resolutionnumber)) * " km resolution, "
end

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
@testset "McDonald Kreitman" begin

    @testset "mkt" begin
        testseqs = [[dna"AAAGGGCCC", dna"AGAGGGCGC"],
                    [dna"AACGGGCCG", dna"ACAGGGCCC"]]

        @test mkt(NaturalSelection.aligned_codons(testseqs[1]),
                  NaturalSelection.aligned_codons(testseqs[2])) == (2, 3, 0, 1, 1.0)
    end

end


text:
"""
Redundent code, should replace this in runtests when I get the chance
"""

small_square = QNetwork()
A = BasicNode("A")
B = BasicNode("B")
C = BasicNode("C")
D = BasicNode("D")
AB = BasicChannel(A, B)
BC = BasicChannel(B, C)
AD = BasicChannel(A, D)
DC = BasicChannel(D, C)

AB.costs = unit_costvector()
BC.costs = unit_costvector()
AD.costs = unit_costvector()
DC.costs = unit_costvector()

for i in [A, B, C, D, AB, BC, AD, DC]
    add(small_square, i)
end

refresh_graph!(small_square)

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
#=

# 301: 3D Laplace equation
([source code](SOURCE_URL))

=#

module Example301_Laplace3D

using VoronoiFVM,ExtendableGrids
using GridVisualize

## Flux function which describes the flux
## between neigboring control volumes
function g!(f,u,edge)
    f[1]=u[1,1]-u[1,2]
end

function s(f,node)
    n=view(node.coord,:,node.index)
    f[1]=n[1]*sin(5.0*n[2])*exp(n[3])
end

function main(;Plotter=nothing,n=5)
    nspecies=1
    ispec=1
    X=collect(0:1/n:1)
    grid=VoronoiFVM.Grid(X,X,X)
    physics=VoronoiFVM.Physics(flux=g!,source=s)
    sys=VoronoiFVM.System(grid,physics)
    enable_species!(sys,ispec,[1])
    boundary_dirichlet!(sys,ispec,5,0.0)
    boundary_dirichlet!(sys,ispec,6,0.0)
    inival=unknowns(sys,inival=0)
    solution=unknowns(sys)
    solve!(solution,inival,sys)
    scalarplot(grid,solution[1,:],Plotter=Plotter,zplane=0.5, flevel=0.5)
    return solution[43]
end

## Called by unit test

function test()
    main() ≈ 0.012234524449380824
end

end


text:
import Random

"""
    permute_params(params; filter)

Creates a list of dictionaries with all possible combinates of the provided
parameters `params`, where each element in `params` represents one parameter and
all its possible values. Optionally, a validator function can be supplied
to filter invalid or unwanted parameter combinations.

# Arguments
- `params`:  Dictionary where each element holds a list with all possible values.
- `filter`:  Function for filtering parameter combinations. It is called with one argument of type `Dict{Any, Any}` and expects a `Bool` return value.
             Example: `filter=x -> x[:key1] > 0.5 || x[:key2] < 1.0`
- `shuffle`: Indicates whether to randomly shuffle combinations.

# Returns
List of dictionaries with all possible parameter combinations.

# Examples
```jldoctest
julia> permute_params(Dict{Any, Vector{Any}}("wnd" => [1,2,3], :mode => ["A"], "coef" => [0.1, 0.5, 1.0]))

9-element Array{Dict{Any,Any},1}:
    Dict("wnd" => 1,:mode => "A","coef" => 0.1)
    Dict("wnd" => 1,:mode => "A","coef" => 0.5)
    Dict("wnd" => 1,:mode => "A","coef" => 1.0)
    Dict("wnd" => 2,:mode => "A","coef" => 0.1)
    Dict("wnd" => 2,:mode => "A","coef" => 0.5)
    Dict("wnd" => 2,:mode => "A","coef" => 1.0)
    Dict("wnd" => 3,:mode => "A","coef" => 0.1)
    Dict("wnd" => 3,:mode => "A","coef" => 0.5)
    Dict("wnd" => 3,:mode => "A","coef" => 1.0)
julia>
julia> params = Dict{Any, Vector{Any}}(:wnd => [1,2], :mode => ["A", "B"], :coef => [0.1, 0.5]);
julia> filter = x -> x[:mode] != "A" || x[:wnd] > 1;
julia> permute_params(params; filter=filter)
6-element Array{Dict{Any,Any},1}:
    Dict(:mode => "A",:wnd => 2,:coef => 0.1)
    Dict(:mode => "A",:wnd => 2,:coef => 0.5)
    Dict(:mode => "B",:wnd => 1,:coef => 0.1)
    Dict(:mode => "B",:wnd => 1,:coef => 0.5)
    Dict(:mode => "B",:wnd => 2,:coef => 0.1)
    Dict(:mode => "B",:wnd => 2,:coef => 0.5)
```
"""
function permute_params(params::Dict{Any, Vector{Any}}; filter::Function=x -> true, shuffle::Bool=false)::Vector{Dict{Any, Any}}
    result = Vector{Dict{Any, Any}}()
    permute_params_internal(params, Dict{Any, Any}(), filter, result)
    if shuffle
        Random.shuffle!(result)
    end
    result
end

function permute_params_internal(
    params          ::Dict{Any, Vector{Any}},
    fixed_params    ::Dict{Any, Any},
    filter          ::Function,
    result          ::Vector{Dict{Any, Any}}
)::Nothing
    if !isempty(params)
        params2 = copy(params)
        key = collect(keys(params))[1]
        all_key_values = params2[key]
        delete!(params2, key)
        for val in all_key_values
            fixed_params2 = copy(fixed_params)
            fixed_params2[key] = val
            permute_params_internal(params2, fixed_params2, filter, result)
        end
    else
        if filter(fixed_params)
            push!(result, fixed_params)
        end
    end
    nothing
end


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

mutable struct OrchestrationServiceStateInput <: SwaggerModel
    serviceName::Any # spec type: Union{ Nothing, String } # spec name: serviceName
    action::Any # spec type: Union{ Nothing, String } # spec name: action

    function OrchestrationServiceStateInput(;serviceName=nothing, action=nothing)
        o = new()
        validate_property(OrchestrationServiceStateInput, Symbol("serviceName"), serviceName)
        setfield!(o, Symbol("serviceName"), serviceName)
        validate_property(OrchestrationServiceStateInput, Symbol("action"), action)
        setfield!(o, Symbol("action"), action)
        o
    end
end # type OrchestrationServiceStateInput

const _property_map_OrchestrationServiceStateInput = Dict{Symbol,Symbol}(Symbol("serviceName")=>Symbol("serviceName"), Symbol("action")=>Symbol("action"))
const _property_types_OrchestrationServiceStateInput = Dict{Symbol,String}(Symbol("serviceName")=>"String", Symbol("action")=>"String")
Base.propertynames(::Type{ OrchestrationServiceStateInput }) = collect(keys(_property_map_OrchestrationServiceStateInput))
Swagger.property_type(::Type{ OrchestrationServiceStateInput }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_OrchestrationServiceStateInput[name]))}
Swagger.field_name(::Type{ OrchestrationServiceStateInput }, property_name::Symbol) =  _property_map_OrchestrationServiceStateInput[property_name]

const _allowed_OrchestrationServiceStateInput_serviceName = ["AutomaticRepairs"]

const _allowed_OrchestrationServiceStateInput_action = ["Resume", "Suspend"]

function check_required(o::OrchestrationServiceStateInput)
    true
end

function validate_property(::Type{ OrchestrationServiceStateInput }, name::Symbol, val)
    if name === Symbol("serviceName")
        Swagger.validate_param(name, "OrchestrationServiceStateInput", :enum, val, _allowed_OrchestrationServiceStateInput_serviceName)
    end
    if name === Symbol("action")
        Swagger.validate_param(name, "OrchestrationServiceStateInput", :enum, val, _allowed_OrchestrationServiceStateInput_action)
    end
end


text:
using MemPool
using Serialization

export domain, UnitDomain, project, alignfirst, ArrayDomain

import Base: isempty, getindex, intersect, ==, size, length, ndims

"""
    domain(x::T)

Returns metadata about `x`. This metadata will be in the `domain`
field of a Chunk object when an object of type `T` is created as
the result of evaluating a Thunk.
"""
function domain end

"""
Default domain -- has no information about the value
"""
struct UnitDomain end

"""
If no `domain` method is defined on an object, then
we use the `UnitDomain` on it. A `UnitDomain` is indivisible.
"""
domain(x::Any) = UnitDomain()

###### Chunk ######

"""
A reference to a piece of data located on a remote worker. `Chunk`s are typically created with `Dagger.tochunk(data)`, and the data can then be accessed from any worker with `collect(::Chunk)`. `Chunk`s are serialization-safe, and use distributed refcounting (provided by `MemPool.DRef`) to ensure that the data referenced by a `Chunk` won't be GC'd, as long as a reference exists on some worker.

Each `Chunk` is associated with a given `Dagger.Processor`, which is (in a sense) the processor that "owns" or contains the data. Calling `collect(::Chunk)` will perform data movement and conversions defined by that processor to safely serialize the data to the calling worker.
"""
mutable struct Chunk{T, H, P<:Processor}
    chunktype::Type{T}
    domain
    handle::H
    processor::P
    persist::Bool
    function (::Type{Chunk{T,H,P}})(::Type{T}, domain, handle, processor, persist) where {T,H,P}
        c = new{T,H,P}(T, domain, handle, processor, persist)
        finalizer(x -> @async(myid() == 1 && free!(x)), c)
        c
    end
end

domain(c::Chunk) = c.domain
chunktype(c::Chunk) = c.chunktype
persist!(t::Chunk) = (t.persist=true; t)
shouldpersist(p::Chunk) = t.persist
affinity(c::Chunk) = affinity(c.handle)

function unrelease(c::Chunk{<:Any,DRef})
    # set spilltodisk = true if data is still around
    try
        destroyonevict(c.handle, false)
        return c
    catch err
        if isa(err, KeyError)
            return nothing
        else
            rethrow(err)
        end
    end
end
unrelease(c::Chunk) = c

Base.:(==)(c1::Chunk, c2::Chunk) = c1.handle == c2.handle
Base.hash(c::Chunk, x::UInt64) = hash(c.handle, x)

collect_remote(chunk::Chunk) =
    move(chunk.processor, OSProc(), poolget(chunk.handle))

function collect(ctx::Context, chunk::Chunk; options=nothing)
    # delegate fetching to handle by default.
    if chunk.handle isa DRef && !(chunk.processor isa OSProc)
        return remotecall_fetch(collect_remote, chunk.handle.owner, chunk)
    elseif chunk.handle isa FileRef
        return poolget(chunk.handle)
    else
        return move(chunk.processor, OSProc(), chunk.handle)
    end
end
collect(ctx::Context, ref::DRef; options=nothing) =
    move(OSProc(ref.owner), OSProc(), ref)
collect(ctx::Context, ref::FileRef; options=nothing) =
    poolget(ref) # FIXME: Do move call

# Unwrap Chunk, DRef, and FileRef by default
move(from_proc::Processor, to_proc::Processor, x::Chunk) =
    move(from_proc, to_proc, x.handle)
move(from_proc::Processor, to_proc::Processor, x::Union{DRef,FileRef}) =
    move(from_proc, to_proc, poolget(x))

# Determine from_proc when unspecified
move(to_proc::Processor, chunk::Chunk) =
    move(chunk.processor, to_proc, chunk)
move(to_proc::Processor, d::DRef) =
    move(OSProc(d.owner), to_proc, d)
move(to_proc::Processor, x) =
    move(OSProc(), to_proc, x)

### ChunkIO
affinity(r::DRef) = Pair{OSProc, UInt64}[OSProc(r.owner) => r.size]
function affinity(r::FileRef)
    if haskey(MemPool.who_has_read, r.file)
        Pair{OSProc, UInt64}[OSProc(dref.owner) => r.size for dref in MemPool.who_has_read[r.file]]
    elseif r.force_pid[] !== nothing
        Pair{OSProc, UInt64}[OSProc(r.force_pid[]) => 1]
    else
        Pair{OSProc, UInt64}[OSProc(w) => r.size for w in MemPool.get_workers_at(r.host)]
    end
end

"""
    tochunk(x, proc; persist=false, cache=false) -> Chunk

Create a chunk from sequential object `x` which resides on `proc`.
"""
function tochunk(x::X, proc::P=OSProc(); persist=false, cache=false) where {X,P}
    ref = poolset(x, destroyonevict=persist ? false : cache)
    Chunk{X, typeof(ref), P}(X, domain(x), ref, proc, persist)
end
tochunk(x::Union{Chunk, Thunk}, proc=nothing) = x

# Check to see if the node is set to persist
# if it is foce can override it
function free!(s::Chunk{X, DRef, P}; force=true, cache=false) where {X,P}
    if force || !s.persist
        if cache
            try
                destroyonevict(s.handle, true) # keep around, but remove when evicted
            catch err
                isa(err, KeyError) || rethrow(err)
            end
        end
    end
end
free!(x; force=true,cache=false) = x # catch-all for non-chunks

function savechunk(data, dir, f)
    sz = open(joinpath(dir, f), "w") do io
        serialize(io, MemPool.MMWrap(data))
        return position(io)
    end
    fr = FileRef(f, sz)
    proc = OSProc()
    Chunk{typeof(data), typeof(fr), typeof(proc)}(typeof(data), domain(data), fr, proc, true)
end

Base.@deprecate_binding AbstractPart Union{Chunk, Thunk}
Base.@deprecate_binding Part Chunk
Base.@deprecate parts(args...) chunks(args...)
Base.@deprecate part(args...) tochunk(args...)
Base.@deprecate parttype(args...) chunktype(args...)


text:
__precompile__()

module AxisArrays

using Base: tail
using RangeArrays, IntervalSets
using Compat

export AxisArray, Axis, axisnames, axisvalues, axisdim, axes, atindex

# From IntervalSets:
export ClosedInterval, ..
Base.@deprecate_binding Interval ClosedInterval

include("core.jl")
include("intervals.jl")
include("search.jl")
include("indexing.jl")
include("sortedvector.jl")
include("combine.jl")

end


text:
export getLaplacianMatrix, getConvectionMatrix

"""
function ddxc(n)

centered one-dimensional finite difference operator.

Input:

    n  - number of cells

Output:

    dx - sparse matrix for 1D derivative size(dx) = (n,n)
"""
function ddxc(n::Int)
    dx = spdiagm(-1=>fill(-1.0,n-1), 1=>fill(1.0,n-1))
    return dx
end

"""
function ddxfwd(n)

forward one-dimensional finite difference operator.

Input:

    n  - number of cells

Output:

    dx - sparse matrix for 1D derivative size(dx) = (n,n)
"""
function ddxfwd(n::Int)
    dx = spdiagm(0=>fill(-1.0,n),1=>fill(1.0,n-1))
    return dx
end

"""
function ddxbwd(n)

backward one-dimensional finite difference operator.

Input:

    n  - number of cells

Output:

    dx - sparse matrix for 1D derivative size(dx) = (n,n)
"""
function ddxbwd(n::Int)
    dx = spdiagm(-1=>fill(-1.0,n-1),0=>fill(1.0,n))
    return dx
end

"""
function getLaplacian(M::RegularMesh)

builds 2D/3D Laplacian for padded cell-centered vector including
ghost points and corners.

Input:

	M  - computational mesh

Output:

	Lap - 2D/3D Laplacian

"""
function getLaplacianMatrix(M::RegularMesh)
    if M.dim==2
		dx = Mesh.ddx(M.n[1]+1)/M.h[1]
		dy = Mesh.ddx(M.n[2]+1)/M.h[2]

		Dxx = kron(sparse(1.0I,M.n[2]+2, M.n[2]+2),dx'*dx)
		Dyy = kron(dy'*dy, sparse(I, M.n[1]+2, M.n[1]+2))
		return Dxx+Dyy
	elseif M.dim==3
		dx = Mesh.ddx(M.n[1]+1)/M.h[1]
		dy = Mesh.ddx(M.n[2]+1)/M.h[2]
		dz = Mesh.ddx(M.n[3]+1)/M.h[3]

		Dxx = kron(sparse(1.0I, M.n[3]+2, M.n[3]+2),kron(sparse(1.0I, M.n[2]+2, M.n[2]+2),dx'*dx))
		Dyy = kron(sparse(1.0I, M.n[3]+2, M.n[3]+2),kron(dy'*dy, sparse(1.0I, M.n[1]+2, M.n[1]+2)))
		Dzz = kron(dz'*dz, kron(sparse(1.0I, M.n[2]+2, M.n[2]+2), sparse(1.0I, M.n[1]+2, M.n[1]+2)))

		return Dxx + Dyy + Dzz
	end
end

"""
function getConvectionMatrix(M,v)

builds stationary convection operator on extended mesh

2D: C(v_1,v_2)     = Diagonal(v_1)*Dx + Diagonal(v_2)*Dy

3D: C(v_1,v_2,v_3) = Diagonal(v_1)*Dx + Diagonal(v_2)*Dy + Diagonal(v_3)*Dz

Input:

	M :: RegularMesh - original domain (not padded)
	v                - velocities discretized on cell-centers of extended mesh
	 				   or function that returns values (X -> vFun) where
					   X has x,y(,z) coordinates column wise.

Output:

	Conv - 2D Convection operator

"""
function getConvectionMatrix(M::RegularMesh,v::Array{Float64})

	n   = M.n
	h   = M.h
	dim = M.dim
	v  = reshape(v,prod(n.+2),dim)

	if dim==2
    	dxf = ddxfwd(n[1]+2)./(h[1])
    	dxb = ddxbwd(n[1]+2)./(h[1])
    	dyf = ddxfwd(n[2]+2)./(h[2])
    	dyb = ddxbwd(n[2]+2)./(h[2])

    	Dxf = kron(sparse(1.0I, n[2]+2, n[2]+2), dxf)
    	Dxb = kron(sparse(1.0I, n[2]+2, n[2]+2), dxb)
    	Dyf = kron(dyf, sparse(1.0I, n[1]+2, n[1]+2))
    	Dyb = kron(dyb, sparse(1.0I, n[1]+2, n[1]+2))
    	return Diagonal(max.(v[:,1],0))*Dxb + Diagonal(min.(v[:,1],0))*Dxf +
		       Diagonal(max.(v[:,2],0))*Dyb + Diagonal(min.(v[:,2],0))*Dyf
	elseif dim==3
		dxf = ddxfwd(n[1]+2)./(M.h[1])
		dxb = ddxbwd(n[1]+2)./(M.h[1])
		dyf = ddxfwd(n[2]+2)./(M.h[2])
		dyb = ddxbwd(n[2]+2)./(M.h[2])
		dzf = ddxfwd(n[3]+2)./(M.h[3])
		dzb = ddxbwd(n[3]+2)./(M.h[3])

		Dxf = kron(sparse(1.0I, n[3]+2, n[3]+2),kron(sparse(1.0I, n[2]+2, n[2]+2), dxf));
		Dxb = kron(sparse(1.0I, n[3]+2, n[3]+2),kron(sparse(1.0I, n[2]+2, n[2]+2), dxb));
		Dyf = kron(sparse(1.0I, n[3]+2, n[3]+2),kron(dyf, sparse(1.0I, n[1]+2, n[1]+2)))
		Dyb = kron(sparse(1.0I, n[3]+2, n[3]+2),kron(dyb, sparse(1.0I, n[1]+2, n[1]+2)))
		Dzf = kron(dzf, sparse(1.0I, prod(n[1:2].+2), prod(n[1:2].+2)))
		Dzb = kron(dzb, sparse(1.0I, prod(n[1:2].+2), prod(n[1:2].+2)))
    	return Diagonal(max.(v[:,1],0))*Dxb + Diagonal(min.(v[:,1],0))*Dxf +
               Diagonal(max.(v[:,2],0))*Dyb + Diagonal(min.(v[:,2],0))*Dyf +
               Diagonal(max.(v[:,3],0))*Dzb + Diagonal(min.(v[:,3],0))*Dzf
	end
end

function getConvectionMatrix(M::RegularMesh,v::Function)
	xc = getNodalGrid(getPaddedMesh(M))
	return getConvectionMatrix(M,v(xc))
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
import Base.*
# struct NGramIterator{T} where {T<:Union{Base.CodeUnits{UInt8,S} where S,Vector{I} where I<:Integer}}
"""
struct NGramIterator{T}
s::T
n::Int
b::Int
end

Iterates and enumerates ngrams of collection of integers `s::T` with zero padding. Enumeration is computed as in positional number systems, where items of `s` are digits and `b` is the base.

In order to reduce collisions when mixing ngrams of different order one should avoid zeros and negative integers in `s` and should set base `b` to be equal to the expected number of unique tokkens in `s`.

# Examples
```jldoctest
julia> it = Mill.NGramIterator(collect(1:10), 3, 10)
Mill.NGramIterator{Array{Int64,1}}([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 3, 10)

julia> [gram for gram in it]
12-element Array{Int64,1}:
1
12
123
234
345
456
567
678
789
900
100
10

julia> sit = Mill.NGramIterator(codeunits("deadbeef"), 3, 256)    # creates collisions as codeunits returns tokens from 0x00:0xff
Mill.NGramIterator{Base.CodeUnits{UInt8,String}}(UInt8[0x64, 0x65, 0x61, 0x64, 0x62, 0x65, 0x65, 0x66], 3, 256)

julia> [gram for gram in sit]
10-element Array{Int64,1}:
100
25701
6579553
6644068
6382690
6578789
6448485
6645094
25958
102
```
"""
struct NGramIterator{T}
    s::T
    n::Int
    b::Int
    # function NGramIterator{T}(x, n, b) where {T}
    # new(s, n, b)
    # end
end

# NGramIterator(s::AbstractString) = NGramIterator(codeunits(s), 3, 256)
# NGramIterator(s::AbstractString, n, b) = NGramIterator(codeunits(s), n, b)

Base.length(it::NGramIterator) = length(it.s) + it.n - 1

function Base.iterate(it::NGramIterator, s = (0, 1))
    idx, i = s
    b, n = it.b, it.n
    if i <= length(it.s)
        idx = idx * b + it.s[i]
        idx = (i>n) ? idx - it.s[i - n]*b^n : idx
        return(idx, (idx, i + 1))
    elseif i < length(it.s) + n
        idx = (i>n) ? idx - it.s[i - n]*b^(n - (i - length(it.s))) : idx
        return(idx, (idx, i + 1))
    else
        return(nothing)
    end
end

"""
ngrams!(o,x,n::Int,b::Int)

store indexes of `n` grams of `x` with base `b` to `o`

"""
function ngrams!(o,x::T,n::Int,b::Int) where {T<:Union{Base.CodeUnits{UInt8,S} where S,Vector{I} where I<:Integer}}
    @assert length(o) >= length(x) + n - 1
    for (i, idx) in enumerate(NGramIterator(x, n, b))
        o[i] = idx
    end
    o
end

"""
ngrams(x,n::Int,b::Int)

indexes of `n` grams of `x` with base `b`

"""
ngrams(x::T,n::Int,b::Int) where {T<:Union{Base.CodeUnits{UInt8,S} where S,Vector{I} where I<:Integer}} =   collect(NGramIterator(x, n, b))
ngrams(x::T,n::Int,b::Int) where {T<:AbstractString} = ngrams(codeunits(x),n,b)

"""
function countngrams!(o,x,n::Int,b::Int)

counts number of of `n` grams of `x` with base `b` to `o` and store it to o

"""
function countngrams!(o,x::T,n::Int,b::Int) where {T<:Union{Base.CodeUnits{UInt8,S} where S,Vector{I} where I<:Integer}}
    for idx in NGramIterator(x, n, b)
        o[mod(idx, length(o))+1] += 1
    end
    o
end

countngrams!(o,x::T,n::Int,b::Int) where {T<:AbstractString} = countngrams!(o,codeunits(x),n,b)

"""
function countngrams(x,n::Int,b::Int)

counts number of of `n` grams of `x` with base `b` to `o`

"""
countngrams(x,n::Int,b::Int,m) = countngrams!(zeros(Int,m),x,n,b)
function countngrams(x::Vector{T},n::Int,b::Int,m) where {T<:AbstractString}
    o = zeros(Int,m,length(x))
    for (i,s) in enumerate(x)
        countngrams!(view(o,:,i),x[i],n,b)
    end
    o
end

string2ngrams(x::T, n, m) where {T <: AbstractArray{I} where I<: AbstractString} = countngrams(Vector(x[:]),n,256,m)
string2ngrams(x::T, n, m) where {T<: AbstractString} = countngrams(x, n, 256, m)
string2ngrams(x, n, m) = x

"""
struct NGramMatrix{T}
s::Vector{T}
n::Int
b::Int
m::Int
end

Represents strings stored in array `s` as ngrams of cardinality `n`. Strings are internally stored as strings and the multiplication with
dense matrix is overloaded and `b` is a base for calculation of trigrams. Finally `m` is the modulo applied on indexes of ngrams.

The structure essentially represents module one-hot representation of strings, where each columns contains one observation (string).
Therefore the structure can be viewed as a matrix with `m` rows and `length(s)` columns
"""
struct NGramMatrix{T} <: AbstractMatrix{T}
    s::Vector{T}
    n::Int
    b::Int
    m::Int
end

Base.show(io::IO, n::NGramMatrix) = (print(io, "NGramMatrix($(n.b), $(n.m))"); show(io, n.s))
Base.show(io::IO, ::MIME{Symbol("text/plain")}, n::NGramMatrix) = Base.show(io, n)

NGramIterator(a::NGramMatrix{T}, i) where {T<:AbstractString} = NGramIterator(codeunits(a.s[i]), a.n, a.b)
NGramIterator(a::NGramMatrix{T}, i) where {T<:Vector{U}} where {U<:Integer} = NGramIterator(a.s[i], a.n, a.b)

Base.length(a::NGramMatrix) = length(a.s)
Base.size(a::NGramMatrix) = (a.m, length(a.s))
Base.size(a::NGramMatrix, d) = (d == 1) ? a.m : length(a.s)
subset(a::NGramMatrix{T}, i) where {T<:AbstractString} = NGramMatrix(a.s[i], a.n, a.b, a.m)
Base.getindex(a::NGramMatrix{T}, i) where {T<:AbstractString} = NGramMatrix(a.s[i], a.n, a.b, a.m)
Base.reduce(::typeof(catobs), a::Vector{S}) where {S<:NGramMatrix} = _catobs(a[:])
Base.reduce(::typeof(catobs), a::Matrix{S}) where {S<:NGramMatrix} = _catobs(a[:])
Base.reduce(::typeof(hcat), a::Vector{S}) where {S<:NGramMatrix} = _catobs(a[:])
Base.reduce(::typeof(hcat), a::Matrix{S}) where {S<:NGramMatrix} = _catobs(a[:])
Base.hcat(a::NGramMatrix...) = reduce(catobs, collect(a))
catobs(a::NGramMatrix...) = _catobs(collect(a))
_lastcat(a::Array{S}) where {S<:NGramMatrix} = _catobs(a)
_catobs(a::AbstractVecOrTuple{NGramMatrix}) = NGramMatrix(reduce(vcat, [i.s for i in a]), a[1].n, a[1].b, a[1].m)

Base.Matrix(x::NGramMatrix) = Matrix(SparseMatrixCSC(x))

SparseArrays.SparseMatrixCSC(x::NGramMatrix) = SparseArrays.SparseMatrixCSC{Float32, UInt}(x)

function SparseArrays.SparseMatrixCSC{Tv, Ti}(x::NGramMatrix) where {Tv, Ti <: Integer}
    l = sum(map(i -> length(NGramIterator(x, i)), 1:length(x)))
    I = zeros(Ti, l)
    J = zeros(Ti, l)
    V = ones(Tv, l)
    vid = 1
    for j in 1:length(x)
        for i in NGramIterator(x, j)
            I[vid] = mod(i, x.m) + 1
            J[vid] = j
            vid += 1
        end
    end
    sparse(I, J, V, size(x,1), size(x,2))
end

function mulkernel!(C, A, jB, mA, nA, idxs)
    for iB in idxs
        miB = mod(iB, nA) + 1
        @inbounds for iA in 1:mA
            C[iA, jB] += A[iA, miB]
        end
    end
end

*(A::Matrix, B::NGramMatrix) = mul(A, B)
function mul(A::AbstractMatrix, B::NGramMatrix)
    mA, nA = size(A)
    @assert nA == size(B,1)
    nB = length(B)
    C = zeros(eltype(A), mA, nB)
    for jB in 1:length(B)
        mulkernel!(C, A, jB, mA, nA, NGramIterator(B, jB))
    end
    return C
end

function multkernel!(C, A, jB, mA, bm, idxs)
    for iB in idxs
        miB = mod(iB, bm) + 1
        @inbounds for iA in 1:mA
            C[iA, miB] += A[iA, jB]
        end
    end
end

function multrans(A::AbstractMatrix, B::NGramMatrix)
    mA, nA = size(A)
    mB = length(B)
    C = zeros(eltype(A), mA, B.m)
    for jB in 1:length(B)
        multkernel!(C, A, jB, mA, B.m, NGramIterator(B, jB))
    end
    return C
end

Zygote.@adjoint function mul(a::AbstractMatrix, b::NGramMatrix)
    return mul(a,b) , Δ -> (multrans(Δ, b),nothing)
end

Zygote.@adjoint function *(a::AbstractMatrix, b::NGramMatrix)
    return mul(a,b) , Δ -> (multrans(Δ, b),nothing)
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/concave_increasing.jl
# Contains definitions of commonly used activation functions:
# relu, parametric relu, leaky relu, maxsig, elu, selu, gelu, sigmoid,
# bisigmoid, swish, maxtanh, pentanh, softsign, softplus.
#############################################################################

# RELU DEFINITION
relu_kernel(t::ANYRELAX, x::MCNoGrad, z::Interval{Float64}) = max_kernel(t, x, 0.0, z)
function relu(t::ANYRELAX, x::MCNoGrad)
    relu_Intv = max(x.Intv, 0.0)
    z, cvi, cci, dcv, dcc = relu_kernel(t, x, relu_Intv)
    return z
end

@inline leaky_relu_kernel(t::ANYRELAX, x::MCNoGrad, z::Interval{Float64}) = param_relu_kernel(t, x, 0.01, z)
@inline leaky_relu(t::ANYRELAX, x::MCNoGrad) = leaky_relu_kernel(t, x, param_relu(x.Intv, 0.01))

function param_relu_kernel(t::ANYRELAX, x::MCNoGrad, α::Float64, z::Interval{Float64})
    @assert α >= 0.0
    xLc = z.lo
    xUc = z.hi
    xL = x.Intv.lo
    xU = x.Intv.hi
    midcv, cvi = mid3(x.cc, x.cv, xL)
    midcc, cci = mid3(x.cc, x.cv, xU)
    dcv = param_relu_deriv(midcv, α)
    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
    u = param_relu(midcv, α)
    o = dcc*(midcc - xL) + xLc
    return MCNoGrad(u, o, z, x.cnst), cvi, cci, dcv, dcc
end
function param_relu(t::ANYRELAX, x::MCNoGrad, α::Float64)
    param_relu_kernel_Intv = param_relu(x.Intv, α)
    z, cvi, cci, dcv, dcc = param_relu_kernel(t, x, α, param_relu_kernel_Intv)
    return z
end

function maxtanh_kernel(t::ANYRELAX, x::MCNoGrad, z::Interval{Float64})
    xLc = z.lo
    xUc = z.hi
    xL = x.Intv.lo
    xU = x.Intv.hi
    midcv, cvi = mid3(x.cc, x.cv, xL)
    midcc, cci = mid3(x.cc, x.cv, xU)
    dcv = maxtanh_deriv(midcv)
    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
    u = maxtanh(midcv)
    o = dcc*(midcc - xL) + xLc
    return MCNoGrad(u, o, z, x.cnst), cvi, cci, dcv, dcc
end
function maxtanh(t::ANYRELAX, x::MCNoGrad)
    maxtanh_Intv = maxtanh(x.Intv)
    z, cvi, cci, dcv, dcc = maxtanh_kernel(t, x, maxtanh_Intv)
    return z
end

function elu_kernel(t::ANYRELAX, x::MCNoGrad, α::Float64, z::Interval{Float64})
    xLc = z.lo
    xUc = z.hi
    xL = x.Intv.lo
    xU = x.Intv.hi
    midcv, cvi = mid3(x.cc, x.cv, xL)
    midcc, cci = mid3(x.cc, x.cv, xU)
    dcv = elu_deriv(midcv, α)
    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
    u = elu(midcv, α)
    o = dcc*(midcc - xL) + xLc
    return MCNoGrad(u, o, z, x.cnst), cvi, cci, dcv, dcc
end
function elu(t::ANYRELAX, x::MCNoGrad, α::Float64)
    elu_Intv = elu(x.Intv, α)
    z, cvi, cci, dcv, dcc = elu_kernel(t, x, α, elu_Intv)
    return z
end

function softplus_kernel(t::ANYRELAX, x::MCNoGrad, z::Interval{Float64})
    xLc = z.lo
    xUc = z.hi
    xL = x.Intv.lo
    xU = x.Intv.hi
    midcv, cv_id = mid3(x.cc, x.cv, xL)
    midcc, cc_id = mid3(x.cc, x.cv, xU)
    dcv = softplus_deriv(midcv)
    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
    u = softplus(midcv)
    o = dcc*(midcc - xL) + xLc
    return MCNoGrad(u, o, z, x.cnst), cvi, cci, dcv, dcc
end
function softplus(t::ANYRELAX, x::MCNoGrad)
    maxtanh_Intv = softplus(x.Intv)
    z, cvi, cci, dcv, dcc = softplus_kernel(t, x, maxtanh_Intv)
    return z
end

function maxsig_kernel(t::ANYRELAX, x::MCNoGrad, z::Interval{Float64})
    xLc = z.lo
    xUc = z.hi
    xL = x.Intv.lo
    xU = x.Intv.hi
    midcv, cv_id = mid3(x.cc, x.cv, xL)
    midcc, cc_id = mid3(x.cc, x.cv, xU)
    dcv = maxsig_deriv(midcv)
    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
    u = maxsig(midcv)
    o = dcc*(midcc - xL) + xLc
    return MCNoGrad(u, o, z, x.cnst), cvi, cci, dcv, dcc
end
function maxsig(t::ANYRELAX, x::MCNoGrad)
    maxsig_Intv = maxsig(x.Intv)
    z, cvi, cci, dcv, dcc = maxsig_kernel(t, x, maxsig_Intv)
    return z
end

# define kernel and operator for sigmoid, bisigmoid, softsign, gelu
for expri in (:pentanh, :sigmoid, :bisigmoid, :softsign)
    expri_cv = Symbol("cv_"*String(expri))
    expri_cc = Symbol("cc_"*String(expri))
    expri_kernel = Symbol(String(expri)*"_kernel")
    eps_min = :xL
    eps_max = :xU
    @eval @inline function ($expri_kernel)(t::ANYRELAX, x::MCNoGrad, y::Interval{Float64}, cv_p::Float64, cc_p::Float64)
        xL = x.Intv.lo
        xU = x.Intv.hi
        midcv, cv_id = mid3(x.cc, x.cv, $eps_min)
        midcc, cc_id = mid3(x.cc, x.cv, $eps_max)
        cv, dcv, cv_p = $(expri_cv)(midcv, xL, xU, cv_p)
        cc, dcc, cc_p = $(expri_cc)(midcc, xL, xU, cc_p)
        return MCNoGrad(cv, cc, y, x.cnst), cv_id, cc_id, cv_p, cc_p
    end
    @eval @inline function ($expri)(t::ANYRELAX, x::MCNoGrad)
        z, tp1, tp2 = ($expri_kernel)(t, x, ($expri)(x.Intv), Inf, Inf)
        return z
    end
end

for expri in (:swish, :gelu)
    expri_cv = Symbol("cv_"*String(expri))
    expri_cc = Symbol("cc_"*String(expri))
    expri_kernel = Symbol(String(expri)*"_kernel")
    if expri == swish
        eps_min = :(SWISH1_MIN < xL ? xL : (SWISH1_MIN > xU ? xU : SWISH1_MIN))
        eps_max = :(swish(xL) < swish(xU) ? xU : xL)
    else
        eps_min = :(GELU_MIN > xU ? xU : (GELU_MIN < xL ? xL : GELU_MIN))
        eps_max = :(gelu(xL) < gelu(xU) ? xU : xL)
    end
    @eval @inline function ($expri_kernel)(t::ANYRELAX, x::MCNoGrad, y::Interval{Float64},
                            cv_p::Float64, cc_p::Float64, cv_p2::Float64, cc_p2::Float64)
        xL = x.Intv.lo
        xU = x.Intv.hi
        midcv, cv_id = mid3(x.cc, x.cv, $eps_min)
        midcc, cc_id = mid3(x.cc, x.cv, $eps_max)
        cv, dcv, cv_p, cv_p2 = $(expri_cv)(midcv, xL, xU, cv_p, cv_p2)
        cc, dcc, cc_p, cc_p2 = $(expri_cc)(midcc, xL, xU, cc_p, cc_p2)
        return MCNoGrad(cv, cc, y, x.cnst), cv_id, cc_id, cv_p, cc_p, cv_p, cc_p, cv_p2, cc_p2
    end
    @eval @inline function ($expri)(t::ANYRELAX, x::MCNoGrad)
        z, cv_id, cc_id, dcv, dcc, tp1, tp2, tp3, tp4 = ($expri_kernel)(t, x, ($expri)(x.Intv), Inf, Inf, Inf, Inf)
        return z
    end
end

@inline function logcosh_kernel(t::ANYRELAX, x::MCNoGrad, y::Interval{Float64})
    xL = x.Intv.lo
    xU = x.Intv.hi
    eps_min = in(0.0, x) ? 0.0 : (xU <= 0.0 ? xU : xL)
    eps_max = (abs(xL) < abs(xU)) ? xU : xL
    midcv, cv_id = mid3(x.cc, x.cv, eps_min)
    midcc, cc_id = mid3(x.cc, x.cv, eps_max)
    cv, dcv = cv_logcosh(midcv, xL, xU)
    cc, dcc = cc_logcosh(midcc, xL, xU)
    return MCNoGrad(cv, cc, y, x.cnst), cv_id, cc_id, dcv, dcc
end
function logcosh(t::ANYRELAX, x::MCNoGrad)
    logcosh_Intv = logcosh(x.Intv)
	z, cv_id, cc_id, dcv, dcc = logcosh_kernel(t, x, logcosh_Intv)
    return z
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions