text:
struct BlockMap{T,As,Rs<:Tuple{Vararg{Int}},Rranges<:Tuple{Vararg{UnitRange{Int}}},Cranges<:Tuple{Vararg{UnitRange{Int}}}} <: LinearMap{T}
    maps::As
    rows::Rs
    rowranges::Rranges
    colranges::Cranges
    # Store BlockMap size explicitly, to allow south-east corner to be
    # empty.
    M::Int
    N::Int

    function BlockMap{T,R,S}(maps::R, rows::S) where {T, R<:Tuple{Vararg{LinearMap}}, S<:Tuple{Vararg{Int}}}
        for A in maps
            promote_type(T, eltype(A)) == T || throw(InexactError())
        end
        rowranges, colranges = rowcolranges(maps, rows)
        M,N = last(last(rowranges)), last(last(colranges))
        return new{T,R,S,typeof(rowranges),typeof(colranges)}(maps, rows, rowranges, colranges, M, N)
    end

    BlockMap{T}(maps::As, rows::Rs, rowranges::Rranges, colranges::Cranges,
                M::Integer, N::Integer) where {T,As,Rs,Rranges,Cranges} =
        new{T,As,Rs,Rranges,Cranges}(maps, rows, rowranges, colranges, M, N)
end

BlockMap{T}(maps::As, rows::S) where {T,As<:Tuple{Vararg{LinearMap}},S} = BlockMap{T,As,S}(maps, rows)

"""
    BlockMap{T}(M, N, (m, n), maps, is, js)

Constructor for a `BlockMap` of `eltype(T)` for the case when all rows
have the same column dimensions; these are specified using the vectors
`m` and `n`. Only those blocks `maps` that are actually occupied need
to be specified, their coordinates are given by the vectors `is` and `js`.
"""
function BlockMap{T}(M::Int, N::Int, (m,n), maps, is, js) where T
    all(m .> 0) && all(n .> 0) ||
        throw(ArgumentError("Block sizes must be positive integers"))

    sum(m) == M && sum(n) == N ||
        throw(DimensionMismatch("Cumulative block dimensions $(sum(m))×$(sum(n)) do not agree with overall size $(M)×$(N)"))

    length(maps) == length(is) == length(js) ||
        throw(ArgumentError("Must provide block coordinates for $(length(maps)) blocks"))

    allunique(zip(is, js)) ||
        throw(ArgumentError("Not all block indices unique"))

    for (A,i,j) in zip(maps,is,js)
        promote_type(T, eltype(A)) == T || throw(InexactError())
        0 < i ≤ length(m) || throw(ArgumentError("Invalid block row $(i) ∉ 1:$(length(m))"))
        0 < j ≤ length(n) || throw(ArgumentError("Invalid block column $(j) ∉ 1:$(length(n))"))
        size(A) == (m[i],n[j]) ||
            throw(DimensionMismatch("Block of size $(size(A)) does not match size at $((i,j)): $((m[i],n[j]))"))
    end

    rows = zeros(Int, length(m))
    colranges = ()
    colstarts = 1 .+ vcat(0,cumsum(n))
    for p = sortperm(collect(zip(is,js)))
        A,i,j = maps[p],is[p],js[p]
        rows[i] += 1
        colranges = (colranges..., colstarts[j]:colstarts[j+1]-1)
    end
    rows = (rows...,)

    rowranges = ()
    i = 1
    for (mm,r) in zip(m,rows)
        iprev = i
        i += mm
        rowranges = (rowranges...,iprev:i-1)
    end

    maps = map(convert_to_lmaps_, maps)

    BlockMap{T}(maps, rows, rowranges, colranges, M, N)
end

function Base.show(io::IO, B::BlockMap{T}) where T
    nrows = length(B.rows)
    # One block in one row may actually span multiple block columns
    # for other rows, but it's nice to get a sense of the structure
    # anyway.
    ncols = maximum(B.rows)
    M,N = size(B)
    write(io, "$(nrows)×$(ncols)-blocked $(M)×$(N) BlockMap{$T}")
end

MulStyle(A::BlockMap) = MulStyle(A.maps...)

function check_dim(A::LinearMap, dim, n)
    n == size(A, dim) || throw(DimensionMismatch("Expected $n, got $(size(A, dim))"))
    return nothing
end

"""
    rowcolranges(maps, rows)

Determines the range of rows for each block row and the range of columns for each
map in `maps`, according to its position in a virtual matrix representation of the
block linear map obtained from `hvcat(rows, maps...)`.
"""
function rowcolranges(maps, rows)
    rowranges = ()
    colranges = ()
    mapind = 0
    rowstart = 1
    for row in rows
        xinds = vcat(1, map(a -> size(a, 2), maps[mapind+1:mapind+row])...)
        cumsum!(xinds, xinds)
        mapind += 1
        rowend = rowstart + size(maps[mapind], 1) - 1
        rowranges = (rowranges..., rowstart:rowend)
        colranges = (colranges..., xinds[1]:xinds[2]-1)
        for colind in 2:row
            mapind +=1
            colranges = (colranges..., xinds[colind]:xinds[colind+1]-1)
        end
        rowstart = rowend + 1
    end
    return rowranges::NTuple{length(rows), UnitRange{Int}}, colranges::NTuple{length(maps), UnitRange{Int}}
end

Base.size(A::BlockMap) = (A.M, A.N)

############
# concatenation
############

for k in 1:8 # is 8 sufficient?
    Is = ntuple(n->:($(Symbol(:A,n))::UniformScaling), Val(k-1))
    L = :($(Symbol(:A,k))::LinearMap)
    args = ntuple(n->Symbol(:A,n), Val(k))

    @eval Base.hcat($(Is...), $L, As::Union{LinearMap,UniformScaling}...) = _hcat($(args...), As...)
    @eval Base.vcat($(Is...), $L, As::Union{LinearMap,UniformScaling}...) = _vcat($(args...), As...)
    @eval Base.hvcat(rows::Tuple{Vararg{Int}}, $(Is...), $L, As::Union{LinearMap,UniformScaling}...) = _hvcat(rows, $(args...), As...)
end

############
# hcat
############
"""
    hcat(As::Union{LinearMap,UniformScaling}...)

Construct a `BlockMap <: LinearMap` object, a (lazy) representation of the
horizontal concatenation of the arguments. `UniformScaling` objects are promoted
to `LinearMap` automatically. To avoid fallback to the generic [`Base.hcat`](@ref),
there must be a `LinearMap` object among the first 8 arguments.

# Examples
```jldoctest; setup=(using LinearMaps)
julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;

julia> L = [CS LinearMap(ones(Int, 3, 3))]::LinearMaps.BlockMap;

julia> L * ones(Int, 6)
3-element Array{Int64,1}:
 4
 5
 6
```
"""
Base.hcat

function _hcat(As::Union{LinearMap,UniformScaling}...)
    T = promote_type(map(eltype, As)...)
    nbc = length(As)

    nrows = -1
    # find first non-UniformScaling to detect number of rows
    for A in As
        if !(A isa UniformScaling)
            nrows = size(A, 1)
            break
        end
    end
    nrows == -1 && throw(ArgumentError("hcat of only UniformScaling objects cannot determine the linear map size"))
    return BlockMap{T}(promote_to_lmaps(fill(nrows, nbc), 1, 1, As...), (nbc,))
end

############
# vcat
############
"""
    vcat(As::Union{LinearMap,UniformScaling}...)

Construct a `BlockMap <: LinearMap` object, a (lazy) representation of the
vertical concatenation of the arguments. `UniformScaling` objects are promoted
to `LinearMap` automatically. To avoid fallback to the generic [`Base.vcat`](@ref),
there must be a `LinearMap` object among the first 8 arguments.

# Examples
```jldoctest; setup=(using LinearMaps)
julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;

julia> L = [CS; LinearMap(ones(Int, 3, 3))]::LinearMaps.BlockMap;

julia> L * ones(Int, 3)
6-element Array{Int64,1}:
 1
 2
 3
 3
 3
 3
```
"""
Base.vcat

function _vcat(As::Union{LinearMap,UniformScaling}...)
    T = promote_type(map(eltype, As)...)
    nbr = length(As)

    ncols = -1
    # find first non-UniformScaling to detect number of columns
    for A in As
        if !(A isa UniformScaling)
            ncols = size(A, 2)
            break
        end
    end
    ncols == -1 && throw(ArgumentError("vcat of only UniformScaling objects cannot determine the linear map size"))

    return BlockMap{T}(promote_to_lmaps(fill(ncols, nbr), 1, 2, As...), ntuple(i->1, nbr))
end

############
# hvcat
############
"""
    hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling}...)

Construct a `BlockMap <: LinearMap` object, a (lazy) representation of the
horizontal-vertical concatenation of the arguments. The first argument specifies
the number of arguments to concatenate in each block row. `UniformScaling` objects
are promoted to `LinearMap` automatically. To avoid fallback to the generic
[`Base.hvcat`](@ref), there must be a `LinearMap` object among the first 8 arguments.

# Examples
```jldoctest; setup=(using LinearMaps)
julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;

julia> L = [CS CS; CS CS]::LinearMaps.BlockMap;

julia> L.rows
(2, 2)

julia> L * ones(Int, 6)
6-element Array{Int64,1}:
 2
 4
 6
 2
 4
 6
```
"""
Base.hvcat

function _hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling}...)
    nr = length(rows)
    T = promote_type(map(eltype, As)...)
    sum(rows) == length(As) || throw(ArgumentError("mismatch between row sizes and number of arguments"))
    n = fill(-1, length(As))
    j = 0
    for i in 1:nr # infer UniformScaling sizes from row counts, if possible:
        ni = -1 # number of rows in this block-row, -1 indicates unknown
        for k in 1:rows[i]
            if !isa(As[j+k], UniformScaling)
                na = size(As[j+k], 1)
                ni >= 0 && ni != na &&
                    throw(DimensionMismatch("mismatch in number of rows"))
                ni = na
            end
        end
        if ni >= 0
            for k = 1:rows[i]
                n[j+k] = ni
            end
        end
        j += rows[i]
    end
    # check for consistent total column number
    nc = -1
    j = 0
    for i in 1:nr
        nci = 0
        rows[i] > 0 && n[j+1] == -1 && (j += rows[i]; continue)
        for k = 1:rows[i]
            nci += isa(As[j+k], UniformScaling) ? n[j+k] : size(As[j+k], 2)
        end
        nc >= 0 && nc != nci && throw(DimensionMismatch("mismatch in number of columns"))
        nc = nci
        j += rows[i]
    end
    nc == -1 && throw(ArgumentError("sizes of UniformScalings could not be inferred"))
    j = 0
    for i in 1:nr
        if rows[i] > 0 && n[j+1] == -1 # this row consists entirely of UniformScalings
            nci, r = divrem(nc, rows[i])
            r != 0 && throw(DimensionMismatch("indivisible UniformScaling sizes"))
            for k = 1:rows[i]
                n[j+k] = nci
            end
        end
        j += rows[i]
    end

    return BlockMap{T}(promote_to_lmaps(n, 1, 1, As...), rows)
end

promote_to_lmaps_(n::Int, dim, J::UniformScaling) = UniformScalingMap(J.λ, n)
promote_to_lmaps_(n::Int, dim, A::LinearMap) = (check_dim(A, dim, n); A)
promote_to_lmaps(n, k, dim) = ()
promote_to_lmaps(n, k, dim, A) = (promote_to_lmaps_(n[k], dim, A),)
@inline promote_to_lmaps(n, k, dim, A, B, Cs...) =
    (promote_to_lmaps_(n[k], dim, A), promote_to_lmaps_(n[k+1], dim, B), promote_to_lmaps(n, k+2, dim, Cs...)...)

############
# basic methods
############

function isblocksquare(A::BlockMap)
    rows = A.rows
    N = length(rows)
    return all(==(N), rows)
end

# the following rules are sufficient but not necessary
function LinearAlgebra.issymmetric(A::BlockMap)
    isblocksquare(A) || return false
    N = length(A.rows)
    maps = A.maps
    symindex = vec(permutedims(reshape(collect(1:N*N), N, N)))
    for i in 1:N*N
        if (i == symindex[i] && !issymmetric(maps[i]))
            return false
        elseif (maps[i] != transpose(maps[symindex[i]]))
            return false
        end
    end
    return true
end

LinearAlgebra.ishermitian(A::BlockMap{<:Real}) = issymmetric(A)
function LinearAlgebra.ishermitian(A::BlockMap)
    isblocksquare(A) || return false
    N = length(A.rows)
    maps = A.maps
    symindex = vec(permutedims(reshape(collect(1:N*N), N, N)))
    for i in 1:N*N
        if (i == symindex[i] && !ishermitian(maps[i]))
            return false
        elseif (maps[i] != adjoint(maps[symindex[i]]))
            return false
        end
    end
    return true
end

############
# comparison of BlockMap objects, sufficient but not necessary
############

Base.:(==)(A::BlockMap, B::BlockMap) = (eltype(A) == eltype(B) && A.maps == B.maps && A.rows == B.rows)

# special transposition behavior

LinearAlgebra.transpose(A::BlockMap) = TransposeMap(A)
LinearAlgebra.adjoint(A::BlockMap)  = AdjointMap(A)

############
# multiplication helper functions
############

@inline function _blockmul!(y, A::BlockMap, x, α, β)
    maps, rows, yinds, xinds = A.maps, A.rows, A.rowranges, A.colranges
    mapind = 0
    @views @inbounds for (row, yi) in zip(rows, yinds)
        iszero(row) && continue
        yrow = selectdim(y, 1, yi)
        mapind += 1
        mul!(yrow, maps[mapind], selectdim(x, 1, xinds[mapind]), α, β)
        for _ in 2:row
            mapind +=1
            mul!(yrow, maps[mapind], selectdim(x, 1, xinds[mapind]), α, true)
        end
    end
    return y
end

@inline function _transblockmul!(y, A::BlockMap, x, α, β, transform)
    maps, rows, xinds, yinds = A.maps, A.rows, A.rowranges, A.colranges
    @views @inbounds begin
        # first block row (rowind = 1) of A, meaning first block column of A', fill all of y
        xcol = selectdim(x, 1, first(xinds))
        for rowind in 1:first(rows)
            mul!(selectdim(y, 1, yinds[rowind]), transform(maps[rowind]), xcol, α, β)
        end
        mapind = first(rows)
        # subsequent block rows of A (block columns of A'),
        # add results to corresponding parts of y
        # TODO: think about multithreading
        for (row, xi) in zip(Base.tail(rows), Base.tail(xinds))
            xcol = selectdim(x, 1, xi)
            for _ in 1:row
                mapind +=1
                mul!(selectdim(y, 1, yinds[mapind]), transform(maps[mapind]), xcol, α, true)
            end
        end
    end
    return y
end

############
# multiplication with vectors & matrices
############

Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) =
    mul!(y, A, x)

Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::TransposeMap{<:Any,<:BlockMap}, x::AbstractVector) =
    mul!(y, A, x)

Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) =
    mul!(y, transpose(A), x)

Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::AdjointMap{<:Any,<:BlockMap}, x::AbstractVector) =
    mul!(y, A, x)

Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) =
    mul!(y, adjoint(A), x)

for Atype in (AbstractVector, AbstractMatrix)
    @eval Base.@propagate_inbounds function LinearAlgebra.mul!(y::$Atype, A::BlockMap, x::$Atype,
                        α::Number=true, β::Number=false)
        require_one_based_indexing(y, x)
        @boundscheck check_dim_mul(y, A, x)
        return _blockmul!(y, A, x, α, β)
    end

    for (maptype, transform) in ((:(TransposeMap{<:Any,<:BlockMap}), :transpose), (:(AdjointMap{<:Any,<:BlockMap}), :adjoint))
        @eval Base.@propagate_inbounds function LinearAlgebra.mul!(y::$Atype, wrapA::$maptype, x::$Atype,
                        α::Number=true, β::Number=false)
            require_one_based_indexing(y, x)
            @boundscheck check_dim_mul(y, wrapA, x)
            return _transblockmul!(y, wrapA.lmap, x, α, β, $transform)
        end
    end
end

############
# materialization into matrices
############

function materialize!(M::MT, A::BlockMap) where {MT<:AbstractArray}
    axes(M) == axes(A) ||
        throw(DimensionMismatch("Cannot materialize BlockMap of size $(size(A)) into matrix of size $(size(M))"))

    M .= false

    maps, rows, yinds, xinds = A.maps, A.rows, A.rowranges, A.colranges
    mapind = 0
    @views @inbounds for (row, yi) in zip(rows, yinds)
        Mrow = selectdim(M, 1, yi)
        for _ in 1:row
            mapind +=1
            copyto!(selectdim(Mrow, 2, xinds[mapind]), convert(Matrix, maps[mapind]))
        end
    end

    M
end

############
# show methods
############

# block2string(b, s) = string(join(map(string, b), '×'), "-blocked ", Base.dims2string(s))
# Base.summary(a::BlockMap) = string(block2string(nblocks(a), size(a)), " ", typeof(a))
# # _show_typeof(io, a) = show(io, typeof(a))
# function Base.summary(io::IO, a::AbstractBlockMap)
#     print(io, block2string(nblocks(a), size(a)))
#     print(io, ' ')
#     _show_typeof(io, a)
# end
# function _show_typeof(io::IO, a::AbstractBlockMap{T}) where {T}
#     Base.show_type_name(io, typeof(a).name)
#     print(io, '{')
#     show(io, T)
#     print(io, '}')
# end

############
# BlockDiagonalMap
############

struct BlockDiagonalMap{T,As<:Tuple{Vararg{LinearMap}},Ranges<:Tuple{Vararg{UnitRange{Int}}}} <: LinearMap{T}
    maps::As
    rowranges::Ranges
    colranges::Ranges
    function BlockDiagonalMap{T,As}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}}
        for A in maps
            promote_type(T, eltype(A)) == T || throw(InexactError())
        end
        # row ranges
        inds = vcat(1, size.(maps, 1)...)
        cumsum!(inds, inds)
        rowranges = ntuple(i -> inds[i]:inds[i+1]-1, Val(length(maps)))
        # column ranges
        inds[2:end] .= size.(maps, 2)
        cumsum!(inds, inds)
        colranges = ntuple(i -> inds[i]:inds[i+1]-1, Val(length(maps)))
        return new{T,As,typeof(rowranges)}(maps, rowranges, colranges)
    end
end

BlockDiagonalMap{T}(maps::As) where {T,As<:Tuple{Vararg{LinearMap}}} =
    BlockDiagonalMap{T,As}(maps)
BlockDiagonalMap(maps::LinearMap...) =
    BlockDiagonalMap{promote_type(map(eltype, maps)...)}(maps)

for k in 1:8 # is 8 sufficient?
    Is = ntuple(n->:($(Symbol(:A,n))::AbstractMatrix), Val(k-1))
    # yields (:A1, :A2, :A3, ..., :A(k-1))
    L = :($(Symbol(:A,k))::LinearMap)
    # yields :Ak
    mapargs = ntuple(n -> :(LinearMap($(Symbol(:A,n)))), Val(k-1))
    # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1)))

    @eval begin
        SparseArrays.blockdiag($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...) =
            BlockDiagonalMap($(mapargs...), $(Symbol(:A,k)), convert_to_lmaps(As...)...)
        function Base.cat($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...; dims::Dims{2})
            if dims == (1,2)
                return BlockDiagonalMap($(mapargs...), $(Symbol(:A,k)), convert_to_lmaps(As...)...)
            else
                throw(ArgumentError("dims keyword in cat of LinearMaps must be (1,2)"))
            end
        end
    end
end

Base.size(A::BlockDiagonalMap) = (last(A.rowranges[end]), last(A.colranges[end]))

LinearAlgebra.issymmetric(A::BlockDiagonalMap) = all(issymmetric, A.maps)
LinearAlgebra.ishermitian(A::BlockDiagonalMap{<:Real}) = all(issymmetric, A.maps)
LinearAlgebra.ishermitian(A::BlockDiagonalMap) = all(ishermitian, A.maps)

LinearAlgebra.adjoint(A::BlockDiagonalMap{T}) where {T} = BlockDiagonalMap{T}(map(adjoint, A.maps))
LinearAlgebra.transpose(A::BlockDiagonalMap{T}) where {T} = BlockDiagonalMap{T}(map(transpose, A.maps))

Base.:(==)(A::BlockDiagonalMap, B::BlockDiagonalMap) = (eltype(A) == eltype(B) && A.maps == B.maps)

Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) =
    mul!(y, A, x, true, false)

Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) =
    mul!(y, transpose(A), x, true, false)

Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) =
    mul!(y, adjoint(A), x, true, false)

for Atype in (AbstractVector, AbstractMatrix)
    @eval Base.@propagate_inbounds function LinearAlgebra.mul!(y::$Atype, A::BlockDiagonalMap, x::$Atype,
                        α::Number=true, β::Number=false)
        require_one_based_indexing(y, x)
        @boundscheck check_dim_mul(y, A, x)
        return _blockscaling!(y, A, x, α, β)
    end
end

@inline function _blockscaling!(y, A::BlockDiagonalMap, x, α, β)
    maps, yinds, xinds = A.maps, A.rowranges, A.colranges
    # TODO: think about multi-threading here
    @views @inbounds for i in eachindex(yinds, maps, xinds)
        mul!(selectdim(y, 1, yinds[i]), maps[i], selectdim(x, 1, xinds[i]), α, β)
    end
    return y
end


text:
module TESTOCCA
using OCCA
using Base.Test

include("testtypes.jl");

#Test vector add kernel.
include("test_okl_ofl_vectoradd.jl");
@test test_okl_ofl_vectoradd(testserial) == true;
@test test_okl_ofl_vectoradd(testopenmp) == true;
@test test_okl_ofl_vectoradd(testopencl) == true;
@test test_okl_ofl_vectoradd(testpthreads) == true;
@test test_okl_ofl_vectoradd(testcuda) == true;

@test test_okl_ofl_vectoradd(testserial) == true;
@test test_okl_ofl_vectoradd(testopencl) == true;
@test test_okl_ofl_vectoradd(testpthreads) == true;
@test test_okl_ofl_vectoradd(testcuda) == true;

#Test reduction kernel.
include("test_okl_reduction.jl");
@test test_okl_reduction(testserial) == true;
@test test_okl_reduction(testopenmp) == true;
@test test_okl_reduction(testopencl) == true;
@test test_okl_reduction(testpthreads) == true;
@test test_okl_reduction(testcuda) == true;

end

ext: jl


text:
using Test
using MathLink

@testset "README" begin

    ###The operations mentioned for MathLink the README
    sin1 = W"Sin"(1.0)
    @test sin1 == W"Sin"(1.0)

    sinx = W"Sin"(W"x")
    @test sinx == W"Sin"(W"x")
    @test W`Sin[1]` == W"Sin"(1)

    @test weval(sin1) == 0.8414709848078965

    @test weval(sinx) == W"Sin"(W"x")

    @test weval(W"Integrate"(sinx, (W"x", 0, 1))) == W"Plus"(1, W"Times"(-1, W"Cos"(1)))
end

using MathLinkExtras

@testset "GreedyEval" begin
    ###Testing turning on and turning of the greedy evaluation
    ###The default is "false"
    @test W"a"+W"b" == W"Plus"(W"a",W"b")
    @test W"a"+W"a" == W"Plus"(W"a",W"a")
    @test W"a"-W"a" == W"Plus"(W"a",W"Minus"(W"a"))
    set_GreedyEval(true)
    @test W"a"+W"b" == W"Plus"(W"a",W"b")
    @test W"a"+W"a" == W"Times"(2,W"a")
    @test W"a"-W"a" == 0
    set_GreedyEval(false)
    @test W"b"+W"b" == W"Plus"(W"b",W"b")
    set_GreedyEval(true)
end

@testset "Rationals" begin
    #### Test Rationals parts
    @test (4//5)*W"a" == weval(W`4 a/5`)
    @test W"a"*(4//5) == weval(W`4 a/5`)
    @test (4//5)/W"a" == weval(W`4/(a 5)`)
    @test W"a"/(4//5) == weval(W`5 a/4`)
    @test weval(1//2) == weval(W`1/2`)
    @test weval([1//2,W`a`]) == W"List"(weval(W`1/2`),W`a`)
end

@testset "Complex Numbers" begin
    #### Test imaginary parts
    @test im*W"a" == weval(W`I * a`)
    @test (2*im)*W"a" == weval(W`2 I a`)
    @test im/W"a" == weval(W`I / a`)
    @test W"a"/(2* im) == weval(W`- I a/2`)
    @test im*(im*W"c") == weval(W`-c`)

    #####Testing that complex nubmers can be put in weval
    @test weval(im+2) == weval(W`I+2`)
    @test weval(im*2) == weval(W`I*2`)
    @test weval(im) == weval(W`I`)

    @test (3*im)*(2*im)*W"a" == weval(W`-6 a`)
    @test (3*im) + (2*im)*W"a" == weval(W`3 I + 2 I a`)
    @test (3*im) - (2*im)*W"a" == weval(W`3 I - 2 I a`)
end

@testset "Unary operators" begin
    @test +W"b" == W"b"
    @test +W`a+b` == W`a+b`
    @test -W"b" == W`-b`
    @test -W`a-b` == W`-a+b`
end

@testset "W2Mstr" begin
    ###Test of a naive MathLink to Mathematica converter function (to resuts can be copied into mathematica directly"

    @testset "Basic Algebra" begin
        @test W2Mstr(W"a") == "a"
        @test W2Mstr(W"x") == "x"
        @test W2Mstr(W"x"+W"y") == "(x + y)"
        @test W2Mstr(W`Sqrt[a + b]`) == "Sqrt[(a + b)]"
        @test W2Mstr(W`Pow[x,2]`) == "Pow[x,2]"
        @test W2Mstr(W`x^2`) == "(x^2)"
        @test W2Mstr(W`a+b`) == "(a + b)"
        @test W2Mstr(weval(W`a + c + v`)) == "(a + c + v)"
        @test W2Mstr(2) == "2"
        @test W2Mstr(W`x`) == "x"
        @test W2Mstr(W"Sin"(W"x")) == "Sin[x]"
        @test W2Mstr(W`Sin[x]`) == "Sin[x]"

    end
    @testset "Nested functions" begin
        @test W2Mstr(weval(W`a + c*b + v`)) == "(a + (b*c) + v)"
        @test W2Mstr(weval(W`(a + c)*(b + v)`)) == "((a + c)*(b + v))"
        @test W2Mstr(weval(W`a^(b+c)`)) == "(a^(b + c))"
        @test W2Mstr(weval(W`a^2`)) == "(a^2)"
        @test W2Mstr(weval(W`e+a^(b+c)`)) == "((a^(b + c)) + e)"
        @test W2Mstr(weval(W`a + c + v + Sin[2 + x + Cos[q]]`)) == "(a + c + v + Sin[(2 + x + Cos[q])])"
        @test W2Mstr(W"a"+W"c"+W"v"+W"Sin"(2 +W"x" + W"Cos"(W"q"))) == "(a + c + v + Sin[(2 + x + Cos[q])])"
        @test W2Mstr(W`Sqrt[x+Sin[y]+z^(3/2)]`) == "Sqrt[(x + Sin[y] + (z^(3*(2^-1))))]"

    end

    @testset "Complex values" begin
        @test W2Mstr(weval(W`2*I`)) == "(2*I)"
        @test W2Mstr(weval(W`2/I`)) == "(-2*I)"
        @test W2Mstr(W`2 + 0*I`) == "(2 + (0*I))"
        @test W2Mstr(W"Complex"(W"c",0)) == "c"
        @test W2Mstr(weval(W"Complex"(W"c",0))) == "c"
        @test W2Mstr(weval(W"Complex"(W"c",W"b"))) == "(c+b*I)"
        @test W2Mstr(im) == "(1*I)"
        @test W2Mstr(2*im) == "(2*I)"
    end

    @testset "Factions" begin
        @test W2Mstr(W`3/4`) == "(3*(4^-1))"
        @test W2Mstr(3//4) == "(3/4)"
        @test W2Mstr(W`b/c`) == "(b*(c^-1))"
        @test W2Mstr(W`b/(c^(a+c))`) == "(b*((c^(a + c))^-1))"
        @test W2Mstr(W`(b^2)/(c^3)`) == "((b^2)*((c^3)^-1))"
        @test W2Mstr(weval(W`(b^2)/(c^3)`)) == "((b^2)*(c^-3))"
    end

    @testset "Lists & Arrays" begin
        @test W2Mstr([W`x`,W`a`]) == "{x,a}"
        @test W2Mstr([W`x`]) == "{x}"
        @test W2Mstr([W`x` W`y`; W`z` W`x`]) == "{{x,y},{z,x}}"
    end

end

@testset "Matrix Multiplication" begin

    P12 = [ 0 1 ; 1 0 ]
    @test P12 * [W"a" W"b" ; W`a+b` 2] == [ W`a+b` 2 ; W"a" W"b"]
    @test [W"a" W"b" ; W`a+b` 2] * P12  == [ W"b" W"a" ; 2 W`a+b`]

    #### test larger matrix
    @test P12 * [W"a" W"b" ; W`a+b` W`v+2`] == [ W`a+b` W`2+v` ; W"a" W"b"]
    @test [W"a" W"b" ; W`a+b` W`v+2`] * P12  == [ W"b" W"a" ; W`2+v` W`a+b`]

    #### test larger matrix
    P13 = fill(0,(3,3))
    P13[1,3]=1
    P13[3,1]=1
    P13[2,2]=1
    Mat = fill(W`a+d`,3,3)
    Mat[:,:] = [W`a+d` W`a+d` W`f*g`; W`a+b` W`v+2` W`f*g` ;  W`d+b`  W`a+b`  W`a+b`]
    P13 * Mat * P13
    #HM2 = P13*Mat*P13

    ###A real live eample
    P14 = fill(0,(4,4))
    P14[1,4]=1
    P14[4,1]=1
    P14[2,2]=1
    P14[3,3]=1
    Mat = MathLink.WExpr[W"Plus"(W"J1245", W"J1346", W"J2356") W"Plus"(W"Times"(W"Complex"(0, 1), W"J1356"), W"Times"(W"Complex"(0, -1), W"J2346")) W"Plus"(W"Times"(W"Complex"(0, -1), W"J1256"), W"Times"(W"Complex"(0, 1), W"J2345")) W"Plus"(W"J1246", W"Times"(-1, W"J1345")); W"Plus"(W"Times"(W"Complex"(0, -1), W"J1356"), W"Times"(W"Complex"(0, 1), W"J2346")) W"Plus"(W"J1245", W"Times"(-1, W"J1346"), W"Times"(-1, W"J2356")) W"Plus"(W"J1246", W"J1345") W"Plus"(W"Times"(W"Complex"(0, -1), W"J1256"), W"Times"(W"Complex"(0, -1), W"J2345")); W"Plus"(W"Times"(W"Complex"(0, 1), W"J1256"), W"Times"(W"Complex"(0, -1), W"J2345")) W"Plus"(W"J1246", W"J1345") W"Plus"(W"Times"(-1, W"J1245"), W"J1346", W"Times"(-1, W"J2356")) W"Plus"(W"Times"(W"Complex"(0, -1), W"J1356"), W"Times"(W"Complex"(0, -1), W"J2346")); W"Plus"(W"J1246", W"Times"(-1, W"J1345")) W"Plus"(W"Times"(W"Complex"(0, 1), W"J1256"), W"Times"(W"Complex"(0, 1), W"J2345")) W"Plus"(W"Times"(W"Complex"(0, 1), W"J1356"), W"Times"(W"Complex"(0, 1), W"J2346")) W"Plus"(W"Times"(-1, W"J1245"), W"Times"(-1, W"J1346"), W"J2356")]
    ####WE just want to see that the numbers can be computed
    Mat * P14
    P14 * Mat
    P14 * Mat* P14
end

@testset "Int128" begin
    set_GreedyEval(false)
    ###These tets check that things work when the integer is larger than int64.
    A0=3^3 * 199 * 3797 * Int128(372928745438657)
    A=7608224128671509719617
    @test A0 == A
    B=1216448257343019439234
    AB = A+B
    @test typeof(A) == Int128
    @test typeof(B) == Int128
    AB1=weval(W"Plus"(W"a", 7608224128671509719617))
    AB2=weval(W"a"+A)
    @test AB1 == AB2
    ###For numbers of this size there is a special WInteger used
    @test weval(AB1 - W"a").value == "$A"
    @test weval(AB1 - W"a"+B).value == "$AB"
end

@testset "W2Tex - LaTex conversion" begin
    @test W2Tex(W`(a+b)^(b+x)`) == "(a+b)^{b+x}"
    @test W2Tex(W`a`) == "a"
    @test W2Tex(W`ab`) == "\\text{ab}"
    @test W2Tex(W`ab*cd`) == "\\text{ab} \\text{cd}"

    ###Testing that MIME form exists for the text/latex option of show.
    @test show(stdout,MIME"text/latex"(),W"a")  == nothing
    @test show(stdout,MIME"text/latex"(),W`a+b`)  == nothing
end


text:
#####
##### apply flux boundary conditions to tracer tendency `Gc`
#####

# a positive west flux is associated with an *increase* in `Gc` near the west boundary.
# same for south and bottom.
@inline function apply_west_bc!(Gc, west_flux, j, k, iter, grid, ΔT)
    @inbounds Gc[1+grid.Hx, j, k] += getbc(west_flux, j, k, iter) * ΔT * Ax(1+grid.Hx, j, k, grid) / volume(1+grid.Hx, j, k, grid)
    return nothing
end

@inline function apply_south_bc!(Gc, south_flux, i, k, iter, grid, ΔT)
    @inbounds Gc[i, 1+grid.Hy, k] += getbc(south_flux, i, k, iter) * ΔT * Ay(i, 1+grid.Hy, k, grid) / volume(i, 1+grid.Hy, k, grid)
    return nothing
end

@inline function apply_bottom_bc!(Gc, bottom_flux, i, j, iter, grid, ΔT)
    @inbounds Gc[i, j, grid.Nz+grid.Hz] += getbc(bottom_flux, i, j, iter) * ΔT * Az(i, j, grid.Nz+grid.Hz, grid) / volume(i, j, grid.Nz+grid.Hz, grid)
    return nothing
end

# a positive east flux is associated with an *decrease* in `Gc` near the east boundary.
# same for north and top.
@inline function apply_east_bc!(Gc, east_flux, j, k, iter, grid, ΔT)
    @inbounds Gc[grid.Nx+grid.Hx, j, k] += getbc(east_flux, j, k, iter) * ΔT * Ax(grid.Nx+grid.Hx, j, k, grid) / volume(grid.Nx+grid.Hx, j, k, grid)
    return nothing
end

@inline function apply_north_bc!(Gc, north_flux, i, k, iter, grid, ΔT)
    @inbounds Gc[i, grid.Ny+grid.Hy, k] += getbc(north_flux, i, k, iter) * ΔT * Ay(i, grid.Ny+grid.Hy, k, grid) / volume(i, grid.Ny+grid.Hy, k, grid)
    return nothing
end

@inline function apply_top_bc!(Gc, top_flux, i, j, iter, grid, ΔT)
    @inbounds Gc[i, j, 1+grid.Hz] += getbc(top_flux, i, j, iter) * ΔT * Az(i, j, 1+grid.Hz, grid) / volume(i, j, 1+grid.Hz, grid)
    return nothing
end

##### apply flux boundary conditions in x direction (west and east)
@kernel function apply_x_bcs_kernel!(Gc, grid, west_bc, east_bc, iter, ΔT)
    j, k = @index(Global, NTuple)
    jj = j + grid.Hy
    kk = k + grid.Hz
    apply_west_bc!(Gc, west_bc, jj, kk, iter, grid, ΔT)
    apply_east_bc!(Gc, east_bc, jj, kk, iter, grid, ΔT)
end
function apply_x_bcs!(Gc, grid, west_bc, east_bc, iter, ΔT, arch, dep)
    kernel! = apply_x_bcs_kernel!(device(arch), (16,16), (grid.Ny, grid.Nz))
    event = kernel!(Gc, grid, west_bc, east_bc, iter, ΔT, dependencies=dep )
    wait(device(arch), event)
    return nothing
end

##### apply flux boundary conditions in y direction (south and north)
@kernel function apply_y_bcs_kernel!(Gc, grid, south_bc, north_bc, iter, ΔT)
    i, k = @index(Global, NTuple)
    ii = i + grid.Hx
    kk = k + grid.Hz
    apply_south_bc!(Gc, south_bc, ii, kk, iter, grid, ΔT)
    apply_north_bc!(Gc, north_bc, ii, kk, iter, grid, ΔT)
end
function apply_y_bcs!(Gc, grid, south_bc, north_bc, iter, ΔT, arch, dep)
    kernel! = apply_y_bcs_kernel!(device(arch), (16,16), (grid.Nx, grid.Nz))
    event = kernel!(Gc, grid, south_bc, north_bc, iter, ΔT, dependencies=dep )
    wait(device(arch), event)
    return nothing
end

##### apply flux boundary conditions in z direction (bottom and top)
@kernel function apply_z_bcs_kernel!(Gc, grid, bottom_bc, top_bc, iter, ΔT)
    i, j = @index(Global, NTuple)
    ii = i + grid.Hx
    jj = j + grid.Hx
    apply_bottom_bc!(Gc, bottom_bc, ii, jj, iter, grid, ΔT)
    apply_top_bc!(Gc, top_bc, ii, jj, iter, grid, ΔT)
end
function apply_z_bcs!(Gc, grid, bottom_bc, top_bc, iter, ΔT, arch, dep)
    kernel! = apply_z_bcs_kernel!(device(arch), (16,16), (grid.Nx, grid.Ny))
    event = kernel!(Gc, grid, bottom_bc, top_bc, iter, ΔT, dependencies=dep )
    wait(device(arch), event)
    return nothing
end

##### apply flux boundary conditions in all three directions, six faces.
function apply_bcs!(Gcs, nuts, grid, iter, ΔT, arch)
    barrier = Event(device(arch))
    events = []
    for name in nut_names
        x_event = apply_x_bcs!(Gcs[name].data, grid, nuts[name].bc.west,   nuts[name].bc.east,  iter, ΔT, arch, barrier)
        y_event = apply_y_bcs!(Gcs[name].data, grid, nuts[name].bc.south,  nuts[name].bc.north, iter, ΔT, arch, barrier)
        z_event = apply_z_bcs!(Gcs[name].data, grid, nuts[name].bc.bottom, nuts[name].bc.top,   iter, ΔT, arch, barrier)
        push!(events, x_event, y_event, z_event)
    end

    # filter out ::Nothing
    events = filter(e -> typeof(e) <: Event, events)

    wait(device(arch), MultiEvent(Tuple(events)))
    return nothing
end

# Avoid some computation / memory accesses when no flux boundary conditions will be applied
apply_x_bcs!(Gc, grid, ::Nothing, ::Nothing, iter, ΔT, arch, dep) = nothing
apply_y_bcs!(Gc, grid, ::Nothing, ::Nothing, iter, ΔT, arch, dep) = nothing
apply_z_bcs!(Gc, grid, ::Nothing, ::Nothing, iter, ΔT, arch, dep) = nothing


text:
using projecttwo, Test

@testset "Single variable optimizers" begin

	@testset "Tolerance tests" begin
		for tolerance in [1, 0.5, 1e-1]
			@test isapprox(line_optimize(x -> x^2 - 1, -5, 5), 0.0, atol = tolerance)
			@test isapprox(line_optimize(x -> x^2 - 2x, -5, 5), 1.0, atol = tolerance)
		end
	end
end

max_stars_repo_name: madzias/Julia

max_issues_repo_name: madzias/Julia

max_forks_repo_name: madzias/Julia


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

julia_version = v"1.6.0"

name = "Libtask"
version = v"0.4.3"

# Collection of sources required to build Libtask
sources = [
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
# create output directory
mkdir -p "${libdir}"

# compile library
cd "${WORKSPACE}/srcdir/"

CFLAGS="-I${includedir} -I${includedir}/julia -O2 -shared -std=gnu99 -fPIC"
if [[ "${target}" == i686-* ]]; then
  CFLAGS="${CFLAGS} -march=pentium4"
fi
if [[ "${target}" == *-mingw* ]]; then
  CFLAGS="${CFLAGS} -Wl,--export-all-symbols"
elif [[ "${target}" == *-linux-* ]]; then
  CFLAGS="${CFLAGS} -Wl,--export-dynamic"
fi

LDFLAGS="-L${libdir} -ljulia"
if [[ "${target}" == *-mingw* ]]; then
  LDFLAGS="${LDFLAGS} -lopenlibm"
fi

$CC $CFLAGS $LDFLAGS task.c -o "${libdir}/libtask_julia.${dlext}"

# install license
install_license LICENSE.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")
platforms = libjulia_platforms(julia_version)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version)),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat = "~$(julia_version.major).$(julia_version.minor)",
               lock_microarchitecture = false)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
Base.@kwdef struct DiscretizedFunction{S<:AbstractRealLineDiscretizer, D<:AbstractDictionary, E<:AbstractDictionary}
    discretizer::S
    dictionary::D
    index_dictionary::E = Dictionary(keys(dictionary), Base.OneTo(length(dictionary)))
end

function DiscretizedFunction(discr::AbstractRealLineDiscretizer, vals::AbstractVector)
    DiscretizedFunction(;discretizer=discr, dictionary = Dictionary(discr, vals))
end

function DiscretizedFunction(discr::AbstractRealLineDiscretizer, f)
    DiscretizedFunction(discr, f.(discr))
end

function DiscretizedFunction(discr::AbstractRealLineDiscretizer, flift, f)
    dictionary = Dictionary( flift.(discr),  f.(discr))
    DiscretizedFunction(;discretizer=discr, dictionary = dictionary)
end

function (f::DiscretizedFunction)(x)
    f.dictionary[f.discretizer(x)]
end

function (f::DiscretizedFunction)(x, default_value)
    get(f.dictionary, f.discretizer(x), default_value)
end


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Test
using Behavior
using Behavior.Gherkin
using Behavior: transformoutline

@testset "Scenario Outline     " begin
    @testset "Transform; Outline description is \"Some description\"; Result description is same" begin
        outline = ScenarioOutline("Some description", String[],
            ScenarioStep[Given("placeholder <foo>")],
            ["foo"],
            ["bar"])

        scenarios = transformoutline(outline)

        @test scenarios[1].description == "Some description"
    end

    @testset "Transform; Outline tags are @foo @bar; Result tags are @foo @bar" begin
        outline = ScenarioOutline("", ["@foo", "@bar"],
            ScenarioStep[Given("some <foo>")],
            ["foo"],
            ["bar"])

        scenarios = transformoutline(outline)

        @test scenarios[1].tags == ["@foo", "@bar"]
    end

    @testset "Transform; Scenario Outline has one example; One Scenario" begin
        outline = ScenarioOutline("", String[],
            ScenarioStep[Given("placeholder <foo>")],
            ["foo"],
            ["bar"])

        scenarios = transformoutline(outline)

        @test length(scenarios) == 1
        scenario = scenarios[1]
        @test scenario.steps[1] == Given("placeholder bar")
    end

    @testset "Transform; Placeholder is quux; quux is replaced by example" begin
        outline = ScenarioOutline("", String[],
            ScenarioStep[Given("placeholder <quux>")],
            ["quux"],
            ["baz"])

        scenarios = transformoutline(outline)

        scenario = scenarios[1]
        @test scenario.steps[1] == Given("placeholder baz")
    end

    @testset "Transform; Two placeholders foo, quux; foo and quux are replaced" begin
        outline = ScenarioOutline("", String[],
            ScenarioStep[Given("placeholders <foo> <quux>")],
            ["foo", "quux"],
            ["bar"; "baz"])

        scenarios = transformoutline(outline)

        scenario = scenarios[1]
        @test scenario.steps[1] == Given("placeholders bar baz")
    end

    @testset "Transform; Steps Given and When; Both steps are transformed" begin
        steps = ScenarioStep[Given("placeholder <quux>"), When("other <quux>")]
        outline = ScenarioOutline("", String[],
            steps,
            ["quux"],
            ["baz"])

        scenarios = transformoutline(outline)

        scenario = scenarios[1]
        @test scenario.steps[1] == Given("placeholder baz")
        @test scenario.steps[2] == When("other baz")
    end

    @testset "Transform; Steps Given and When again; Both steps are transformed" begin
        steps = ScenarioStep[Given("placeholder <quux>"), When("another step <quux>")]
        outline = ScenarioOutline("", String[],
            steps,
            ["quux"],
            ["baz"])

        scenarios = transformoutline(outline)

        scenario = scenarios[1]
        @test scenario.steps[1] == Given("placeholder baz")
        @test scenario.steps[2] == When("another step baz")
    end

    @testset "Transform; Step Then; Step is transformed" begin
        steps = ScenarioStep[Then("step <quux>")]
        outline = ScenarioOutline("", String[],
            steps,
            ["quux"],
            ["baz"])

        scenarios = transformoutline(outline)

        scenario = scenarios[1]
        @test scenario.steps[1] == Then("step baz")
    end

    @testset "Transform; Two examples; Two scenarios in the result" begin
        outline = ScenarioOutline("", String[],
            ScenarioStep[Given("step <quux>")],
            ["quux"],
            ["bar" "baz"])

        scenarios = transformoutline(outline)

        scenario = scenarios[1]
        @test scenarios[1].steps[1] == Given("step bar")
        @test scenarios[2].steps[1] == Given("step baz")
    end

    @testset "Transform; Placeholders in the block_text; Placeholders are replaced with examples" begin
        outline = ScenarioOutline("", String[],
            [Given(""; block_text="given <quux>"),
             When(""; block_text="when <quux>"),
             Then(""; block_text="then <quux>")],
            ["quux"],
            ["bar"])

        scenarios = transformoutline(outline)

        scenario = scenarios[1]
        @test scenarios[1].steps[1] == Given(""; block_text="given bar")
        @test scenarios[1].steps[2] == When(""; block_text="when bar")
        @test scenarios[1].steps[3] == Then(""; block_text="then bar")
    end

    @testset "Transform; Outline examples are AbstractStrings; Interpolation works" begin
        outline = ScenarioOutline("Some description", String[],
            ScenarioStep[Given("placeholder <foo>")],
            ["foo"],
            AbstractString["bar"])

        scenarios = transformoutline(outline)

        @test scenarios[1].description == "Some description"
    end

    # TODO: Mismatching placeholders
end


text:
using Knet, Base.Test

# Uncomment these if you want lots of messages:
# import Base.Test: default_handler, Success, Failure, Error
# default_handler(r::Success) = info("$(r.expr)")
# default_handler(r::Failure) = warn("FAIL: $(r.expr)")
# default_handler(r::Error)   = warn("$(r.err): $(r.expr)")

load_only = true
isapprox3(a,b,c)=all(map((x,y,z)->isapprox(x,y;rtol=z), a,b,c))

include(Knet.dir("examples","linreg.jl"))
#@test LinReg.main("") == (0.0005497846637255735,32.77257400591496, 0.11265426200775067) #gpu
#@test isapprox3(LinReg.main(""), (0.0005439055920768844,32.772149551563935,0.11210599283551381), (1e-4,1e-4,1e-4)) #failed on 2016-05-16
@test isapprox3(LinReg.main(""),  (0.0005438227392702326,32.772139750826454,0.11215015965120426), (1e-4,1e-4,1e-4))

# (1,0.30519268f0,33.105003f0,5.6772943f0,0.16634920246297397,0.9490333333333334,0.1688508650672793,0.9490000000000001)
# ERROR: LoadError: LoadError: test failed: isapprox3(MNIST2D.main("--epochs 1"),(0.30533937f0,33.110886f0,5.677294f0),(0.0001,0.0001,0.0001))

include(Knet.dir("examples","mnist2d.jl"))
#@test MNIST2D.main("--epochs 1") == (0.3204898f0, 32.93997f0, 4.614684f0) #gpu
#@test isapprox3(MNIST2D.main("--epochs 1"), (0.30533937f0,33.110886f0,5.677294f0), (1e-4,1e-4,1e-4))
@test isapprox3(MNIST2D.main("--epochs 1"), (0.30519268f0,33.105003f0,5.6772943f0), (1e-4,1e-4,1e-4))

#@test MNIST2D.main("--epochs 1 --ysparse") == (0.3204898f0, 32.93997f0, 4.614684f0) #gpu
#@test isapprox3(MNIST2D.main("--epochs 1 --ysparse"), (0.30533937f0,33.110886f0,5.677294f0), (1e-4,1e-4,1e-4))
@test isapprox3(MNIST2D.main("--epochs 1 --ysparse"), (0.30519268f0,33.105003f0,5.6772943f0), (1e-4,1e-4,1e-4))

warn("Need to implement: A_mul_B!(::Array{Float32,2}, ::Array{Float32,2}, ::SparseMatrixCSC{Float32,Int32}); the other direction is implemented in sparse/linalg.jl")
# @test isapprox3(MNIST2D.main("--epochs 1 --xsparse"), (0.32048982f0,32.93997f0,4.6146836f0), (0.01,0.001,0.1))
# @test isapprox3(MNIST2D.main("--epochs 1 --xsparse --ysparse"), (0.3204898f0,32.93997f0,4.614684f0), (0.01,0.001,0.1))

warn("cpu conv not implemented yet")
# include(Knet.dir("examples","mnist4d.jl"))
# @test isapprox3(MNIST4D.main("--epochs 1"), (0.20520441331629022,65.88371276855469,114.9615707397461), (0.02,0.001,0.2))

# Differences below due mostly to random initialization

include(Knet.dir("examples","adding.jl"))
#@test Adding.main("--epochs 1") == (0.22713459f0,3.3565507f0,5.3267756f0) #gpu
@test isapprox3(Adding.main("--epochs 1"), (0.21368636f0,3.4062932f0,6.621738f0), (1e-4,1e-4,1e-4))

#@test Adding.main("--epochs 1 --nettype lstm") == (0.24768005f0,3.601481f0,1.2290705f0) #gpu
@test isapprox3(Adding.main("--epochs 1 --nettype lstm"), (0.24432828f0,3.6435652f0,1.1989625f0), (1e-4,1e-4,1e-4))

seqdata = Knet.dir("data","seqdata.txt")

include(Knet.dir("examples","rnnlm.jl"))
#@test RNNLM.main("--max_max_epoch 1 --dense seqdata.txt") == (30.65637197763184, 110.81809997558594,29.36899185180664) #gpu
@test isapprox3(RNNLM.main("--max_max_epoch 1 --dense $seqdata"), (30.349935446042096,110.69515228271484,26.056880950927734), (1e-4,1e-4,1e-4))

include(Knet.dir("examples","copyseq.jl"))
#@test      CopySeq.main("--epochs 1 --dense seqdata.txt") == (40.00286169097269, 30.352935791015625,1.646486520767212)
#@test isapprox3(CopySeq.main("--epochs 1 --dense $seqdata"), (40.00211618458305,29.903160095214844,1.4534363746643066), (1e-4,1e-4,1e-4)) # failed 2016-05-16
@test  isapprox3(CopySeq.main("--epochs 1 --dense $seqdata"), (40.000922734021565,22.803560256958008,0.6356764435768127), (1e-4,1e-4,1e-4))

### DEAD CODE:
# using Compat
# using Knet
# @date include("testdense.jl")
# @date include("testsparse.jl")
# @date include("testconvert.jl")
# @date include("testlinalg.jl")
# @date include("testcolops.jl")
# @date include("testmnist.jl")
# @date include("testlayers.jl")
# @date include("testperceptron.jl")
# @date include("testkperceptron.jl")
# # @date include("tutorial.jl")


text:
using OhMyREPL
using Revise

enable_autocomplete_brackets(false)

max_stars_repo_name: mgarbiso/markus-config

max_issues_repo_name: mgarbiso/markus-config

max_forks_repo_name: mgarbiso/markus-config


text:
#: -------------------------------
IntRange = UnitRange{Int}

TCode = Int
TPattern = Vector{TCode}
global const __DEFAULT_PATT__ = TCode[]
empty_TPattern() = TCode[]
global const __M0_PATT__  = TCode[]
one_elem_TPattern(x) = TCode[x,]
global const __M1_PATT__  = one_elem_TPattern(-1)
# reserved token 999999999 for all number line
# used by tokenize
global const __PATT__all_number_line__ = one_elem_TPattern(999999999)

global const __DEFAULT_HASH_CODE__ = UInt64(0x0)
global const __DEFAULT__R__ = (__DEFAULT_HASH_CODE__, TPattern[])  #+ modify
global const RTYPE = typeof(__DEFAULT__R__)

global const SIMILARITY_LEVEL = 0.7

import Base.copy

copy(x::RTYPE) = (x[1],copy(x[2]))

TokenValueLine = Vector{Pair{String,Any}}
DATATYPE = Vector{Pair{IntRange,Vector{TokenValueLine}}}

#: -------------------------------

max_stars_repo_name: algorithmx/SmartParser

max_issues_repo_name: algorithmx/SmartParser

max_forks_repo_name: algorithmx/SmartParser


text:
"""
# module LUFactorization

# Examples

```jldoctest
julia>
```
"""
module LUFactorization

using Reexport

include("Substitution.jl")

include("Doolittle.jl")
@reexport using .Doolittle

include("LUP.jl")
@reexport using .LUP

end


text:
# FIELDS:
#==============================================================================#
"""
    abstract type Field end

abstract type at the top of the type hierarchy for denoting fields over which
vector spaces can be defined. Two common fields are `ℝ` and `ℂ`, representing the
field of real or complex numbers respectively.
"""
abstract type Field end

struct RealNumbers <: Field end
struct ComplexNumbers <: Field end

const ℝ = RealNumbers()
const ℂ = ComplexNumbers()

Base.show(io::IO, ::RealNumbers) = print(io,"ℝ")
Base.show(io::IO, ::ComplexNumbers) = print(io,"ℂ")

Base.in(::Any, ::Field) = false
Base.in(::Real, ::RealNumbers) = true
Base.in(::Number, ::ComplexNumbers) = true

Base.@pure Base.issubset(::Type, ::Field) = false
Base.@pure Base.issubset(::Type{<:Real}, ::RealNumbers) = true
Base.@pure Base.issubset(::Type{<:Number}, ::ComplexNumbers) = true
Base.@pure Base.issubset(::RealNumbers, ::ComplexNumbers) = true

# VECTOR SPACES:
#==============================================================================#
"""
    abstract type VectorSpace end

abstract type at the top of the type hierarchy for denoting vector spaces
"""
abstract type VectorSpace end

"""
    function fieldtype(V::VectorSpace) -> Field

Returns the field type over which a vector space is defined.
"""
function fieldtype end
fieldtype(V::VectorSpace) = fieldtype(typeof(V))

# Basic vector space methods
#----------------------------
"""
    space(a) -> VectorSpace

Returns the vector space associated to object `a`.
"""
function space end

"""
    dim(V::VectorSpace) -> Int
Returns the total dimension of the vector space `V` as an Int.
"""
function dim end

"""
    dual(V::VectorSpace) -> VectorSpace
Returns the dual space of `V`; also obtained via `V'`. It is assumed that
`typeof(V) == typeof(V')`.
"""
function dual end

# convenience definitions:
Base.adjoint(V::VectorSpace) = dual(V)
Base.:*(V1::VectorSpace, V2::VectorSpace) = ⊗(V1, V2)

# Hierarchy of elementary vector spaces
#---------------------------------------
"""
    abstract type ElementarySpace{k} <: VectorSpace end

Elementary finite-dimensional vector space over a field `k` that can be used as the
index space corresponding to the indices of a tensor. Every elementary vector space
should respond to the methods `conj` and `dual`, returning the complex conjugate space
and the dual space respectively. The complex conjugate of the dual space is obtained
as `dual(conj(V)) === conj(dual(V))`. These different spaces should be of the same type,
so that a tensor can be defined as an element of a homogeneous tensor product
of these spaces.
"""
abstract type ElementarySpace{k} <: VectorSpace end
const IndexSpace = ElementarySpace

fieldtype(::Type{<:ElementarySpace{k}}) where {k} = k

"""
    conj(V::ElementarySpace) -> ElementarySpace
Returns the conjugate space of `V`. For `fieldtype(V)==ℝ`, `conj(V) == V`
It is assumed that `typeof(V) == typeof(conj(V))`.
"""
Base.conj(V::ElementarySpace{ℝ}) = V

"""
    abstract type InnerProductSpace{k} <: ElementarySpace{k} end

abstract type for denoting vector with an inner product and a corresponding
metric, which can be used to raise or lower indices of tensors
"""
abstract type InnerProductSpace{k} <: ElementarySpace{k} end

"""
    abstract type EuclideanSpace{k<:Union{ℝ,ℂ}} <: InnerProductSpace{k} end

abstract type for denoting real or complex spaces with a standard (Euclidean)
inner product (i.e. orthonormal basis), such that the dual space is naturally
isomorphic to the conjugate space (in the complex case) or even to the
space itself (in the real case)
"""
abstract type EuclideanSpace{k} <: InnerProductSpace{k} end # k should be ℝ or ℂ

dual(V::EuclideanSpace) = conj(V)
# dual space is naturally isomorphic to conjugate space for inner product spaces

# representation spaces: we restrict to complex Euclidean space supporting unitary representations
"""
    abstract type AbstractRepresentationSpace{G<:Sector} <: EuclideanSpace{ℂ} end

Complex Euclidean space with a direct sum structure corresponding to different
different superselection sectors of type `G<:Sector`, e.g. the elements or
irreps of a compact or finite group, or the labels of a unitary fusion category.
"""
abstract type AbstractRepresentationSpace{G<:Sector} <: EuclideanSpace{ℂ} end

"""
    function sectortype(a) -> Sector

Returns the type of sector over which object `a` (e.g. a representation space or
an invariant tensor) is defined. Also works in type domain.
"""
sectortype(V::VectorSpace) = sectortype(typeof(V))
sectortype(::Type{<:ElementarySpace}) = Trivial
sectortype(::Type{<:AbstractRepresentationSpace{G}}) where {G} = G

"""
    function sectors(a)

    Returns the different sectors of object `a`( e.g. a representation space or
    an invariant tensor).
"""
sectors(::ElementarySpace) = (Trivial(),)
dim(V::ElementarySpace, ::Trivial) = sectortype(V) == Type{Trivial} ? dim(V) : throw(SectorError())

# Composite vector spaces
#-------------------------
"""
    abstract type CompositeSpace{S<:ElementarySpace} <: VectorSpace end

Abstract type for composite spaces that are defined in terms of a number of
elementary vector spaces of a homogeneous type `S<:ElementarySpace{k}`.
"""
abstract type CompositeSpace{S<:ElementarySpace} <: VectorSpace end

fieldtype(::Type{<:CompositeSpace{S}}) where {S<:ElementarySpace} = fieldtype(S)
sectortype(::Type{<:CompositeSpace{S}}) where {S<:ElementarySpace} = sectortype(S)

# Specific realizations of ElementarySpace types
#------------------------------------------------
# spaces without internal structure
include("cartesianspace.jl")
include("complexspace.jl")
include("generalspace.jl")
include("representationspace.jl")

# # Specific realizations of CompositeSpace types
# #-----------------------------------------------
include("productspace.jl")
# include("superspace.jl")
# include("invariantspace.jl")

# Other examples might include:
# braidedspace and fermionspace
# symmetric and antisymmetric subspace of a tensor product of identical vector spaces
# ...


text:
using Base.Test
using NullableArrays
using Lifting

_a_ = [Int, Nullable{Float64}]
@test lift.(_a_) == [ Int, Float64 ]
@test unlift.(_a_) == [ Nullable{Int}, Nullable{Float64} ]

@test lift(Nullable(1)) == 1
@test lift(Nullable("hello")) == "hello"

@test get(unlift(1) == Nullable(1))
@test get(unlift("hello") == Nullable("hello"))

@test lift(NullableVector{Int}) == Array{Int,1}
@test unlift(Vector{String}) == NullableArray{String, 1}

@test lift(NullableArray([1,2,3])) == [1, 2, 3]

a = unlift([1, 2, 3])
b = NullableArray([1,2,3])

for i in 1:3
	@test get(a[i] == b[i])
end


text:
using LotMassing
using Test

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


text:
struct PlanarTrivial <: TensorKit.Sector
end

struct 𝔹 end
Base.:^(::Type{𝔹}, d::Int) = GradedSpace(PlanarTrivial()=>d);

Base.one(::Type{PlanarTrivial})  = PlanarTrivial();
Base.conj(c::PlanarTrivial) = PlanarTrivial();
TensorKit.:⊗(c1::PlanarTrivial, c2::PlanarTrivial) = (PlanarTrivial(),)

Base.IteratorSize(::Type{TensorKit.SectorValues{PlanarTrivial}}) = TensorKit.HasLength()
Base.length(::TensorKit.SectorValues{PlanarTrivial}) = 1
Base.iterate(::TensorKit.SectorValues{PlanarTrivial}, i = 0) = i==0 ? (PlanarTrivial(),1) : nothing;
Base.getindex(::TensorKit.SectorValues{PlanarTrivial}, i::Int) = i == 1 ? PlanarTrivial() : ArgumentError("lol");
TensorKit.findindex(::TensorKit.SectorValues{PlanarTrivial}, c::PlanarTrivial)  = 1
Base.isless(::PlanarTrivial, ::PlanarTrivial) = false;
TensorKit.BraidingStyle(::Type{PlanarTrivial}) = TensorKit.NoBraiding();
TensorKit.FusionStyle(::Type{PlanarTrivial}) = TensorKit.UniqueFusion();
TensorKit.Fsymbol(args::Vararg{PlanarTrivial,6}) = 1
TensorKit.Nsymbol(args::Vararg{PlanarTrivial,3}) = 1

#take a normal mpo hamiltonian and change its spaces to be \bbB, therefore disabling non planar operations
function force_planar(x::AbstractTensorMap)
    t = TensorMap(zeros,eltype(x),reduce(*,map(i->𝔹^dim(space(x,i)),codomainind(x))),reduce(*,map(i->𝔹^dim(space(x,i)),domainind(x))))
    copyto!(blocks(t)[PlanarTrivial()],convert(Array,x));
    t
end

function force_planar(mpo::MPOHamiltonian)
    MPOHamiltonian(map(Iterators.product(1:mpo.period,1:mpo.odim,1:mpo.odim)) do (i,j,k)
        force_planar(mpo[i,j,k])
    end)
end
force_planar(mpo::InfiniteMPO) = InfiniteMPO(force_planar.(mpo.opp))


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "FMUComplianceChecker"
version = v"2.0.4"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/modelica-tools/FMUComplianceChecker/archive/f0fd0b2ca78c415a725b028deec46c1962799f1d.zip",
                  "db4607e7b6230d9b7f2e394e6ec132e90cee19ae6625997c09d946ab61dda156"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/FMUComplianceChecker*/
atomic_patch -p1 ../patches/forward-cmake-toolchain.patch
atomic_patch -p1 ../patches/windows-lowercase-header-file.patch
export CFLAGS="-I${includedir}"
if [[ "${target}" == *-apple-* ]]; then
    # Replace `/usr/bin/libtool` with `libtool` from our toolchain to fix error
    #     libtool:   error: unrecognised option: '-static'
    ln -sf $(which libtool) /usr/bin/libtool
fi
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ..
make -j${nproc}
install -Dm755 fmuCheck.* "${bindir}/fmuCheck${exeext}"
# Delete expat stuff
rm -rf "${includedir}" "${prefix}/lib"
install_license "../LICENSE"
"""

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

# The products that we will ensure are always built
products = [
    ExecutableProduct("fmuCheck", :fmucheck),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

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

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
Build = Command(
   "build",
   "pangraph build <options> [arguments]",
   "align genomes into a multiple sequence alignment graph",
   """zero, one or more fasta files.
      files can be optionally gzipped.
      multiple records within one file are treated as seperate genomes.
      if no file is given, reads fasta record from stdin""",
   [
    Arg(
        Int,
        "minimum length",
        (short="-l", long="--len"),
        "minimum block size for alignment graph (in nucleotides)",
        100,
    ),
    Arg(
        Float64,
        "block junction cost",
        (short="-a", long="--alpha"),
        "energy cost for introducing junction due to alignment merger",
        100,
    ),
    Arg(
        Float64,
        "block diversity cost",
        (short="-b", long="--beta"),
        "energy cost for interblock diversity due to alignment merger",
        20,
    ),
    Arg(
        Bool,
        "circular genomes",
        (short="-c", long="--circular"),
        "toggle if input genomes are circular",
        false,
    ),
    Arg(
        Bool,
        "enforce uppercase",
        (short="-u", long="--upper-case"),
        "transforms all sequence to upper case",
        false,
    ),
    Arg(
        String,
        "pairwise sensitivity",
        (short="-s", long="--sensitivity"),
        "used to set pairwise alignment sensitivity\n\trecognized options: [5, 10, 20]",
        "10",
    ),
    Arg(
        Int,
        "maximum self maps",
        (short="-x", long="--max-self-map"),
        "maximum number of self mappings to consider per pairwise graph merger",
        100,
    ),
    Arg(
        String,
        "distance calculator",
        (short="-d", long="--distance-backend"),
        "backend to use to estimate pairwise distance for guide tree\n\trecognized options: [native, mash]",
        "native",
    ),
   ],

   (args) -> let
       files = parse(Build, args)
       files = if files === nothing || length(files) == 0
           ["/dev/stdin"]
       else
           files
       end

       minblock  = arg(Build, "-l")
       circular  = arg(Build, "-c")
       uppercase = arg(Build, "-u")

       α = arg(Build, "-a")
       β = arg(Build, "-b")

       energy = (aln) -> let
           len = aln.length
           len < minblock && return Inf

           cuts(hit) = (hit.start > minblock) + ((hit.length-hit.stop) > minblock)

           ncuts = cuts(aln.qry)+cuts(aln.ref)
           nmuts = aln.divergence*aln.length

           return -len + α*ncuts + β*nmuts
       end

       maxiter = arg(Build, "-x")
       sensitivity = @match arg(Build, "-s") begin
           "5"  => "asm5"
           "10" => "asm10"
           "20" => "asm20"
            _   => begin
                usage(Build)
                exit(1)
            end
       end

       graph(io) = graphs(io; circular=circular, upper=uppercase)
       isolates  = (G for file in files for G ∈ open(graph,file))

       compare = @match arg(Build, "-d") begin
           "native" => Graphs.Mash.distance
           "mash"   => begin
               if !Graphs.havecommand("mash")
                   panic("external command mash not found. either install or use native backend\n")
               end
               Graphs.mash
           end
           _        => begin
                # XXX: hacky...
                Build.arg[8].value = "native"

                usage(Build)
                exit(1)
           end
       end

       graph = Graphs.align(isolates...;
            compare     = compare,
            energy      = energy,
            minblock    = minblock,
            maxiter     = maxiter,
            sensitivity = sensitivity,
       )
       finalize!(graph)

       marshal(stdout, graph; fmt=:json)
   end
)

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
@testset "FFTGS" begin
  𝒟 = CartesianGrid(100,100)
  problem = SimulationProblem(𝒟, :z=>Float64, 3)

  Random.seed!(2019)
  solver = FFTGS(:z => (variogram=GaussianVariogram(range=10.),))
  sol = solve(problem, solver)

  if visualtests
    @test_reference "data/FFT-iso.png" plot(sol,size=(900,300))
  end

  Random.seed!(2019)
  solver = FFTGS(:z => (variogram=GaussianVariogram(distance=aniso2distance([20.,5.],[0.])),))
  sol = solve(problem, solver)

  if visualtests
    @test_reference "data/FFT-aniso.png" plot(sol,size=(900,300))
  end
end


text:
"""
Molecular structure manipulation and analysis
"""
module Dorothy

using LinearAlgebra
using Statistics
using Formats
using FormatStreams

export
		MolecularModelHeader, MolecularModel, MolecularModelView,
		ParticleCollection, Particle, MolecularTrajectory, pbcstrategy,
		pbcstrategy!, SelectionCache, SelectionMode, Selector, Selectors,
		DorothyIO, @DorothyAll

include("utils.jl")
include("properties.jl")
include("headers.jl")
include("graphs.jl")
include("multicolls.jl")
include("geometry.jl")
include("pbc.jl")
include("neighbors.jl")

using .Utils
using .Properties
using .Headers
using .Graphs
using .Multicollections
using .Geometry
using .PBC
using .Neighbors

include("models.jl")
include("atoms.jl")
include("hierarchies.jl")

using .Atoms
using .Hierarchies

include("selections.jl")

using .Selectors

struct DorothyIO <: Formats.FormatHandler end

include("formats/gromos87.jl")
include("formats/pdb.jl")
include("formats/xdr.jl")
include("formats/trr.jl")
include("formats/xtc.jl")
include("formats/ndx.jl")

using .Gromos87
using .PDB
using .TRR
using .XTC
using .NDX

include("ss.jl")
include("topology.jl")

using .SecondaryStructure
using .Topology

macro DorothyAll()
	return :( using LinearAlgebra,
			Statistics,
			Formats,
			FormatCodecs,
			FormatStreams,
			Dorothy,
			Dorothy.Utils,
			Dorothy.Properties,
			Dorothy.Graphs,
			Dorothy.Multicollections,
			Dorothy.Geometry,
			Dorothy.PBC,
			Dorothy.Neighbors,
			Dorothy.Atoms,
			Dorothy.Hierarchies,
			Dorothy.Selectors,
			Dorothy.Gromos87,
			Dorothy.PDB,
			Dorothy.TRR,
			Dorothy.XTC,
			Dorothy.NDX,
			Dorothy.SecondaryStructure,
			Dorothy.Topology )
end

datapath = joinpath(@__DIR__, "..", "data")

function __init__()
	register_gromos87()
	register_pdb()
	register_trr()
	register_xtc()
	register_ndx()
end

end # module


text:
module day4

function readInput()
    result = [Dict{String, String}()]
    for line in readlines(joinpath(@__DIR__, "input.txt"))
        if line == ""
            push!(result, Dict{String, String}())
        else
            push!(result[end], map(x -> Pair(split(x, ":")...), split(line, " "))...) # ¯\_(ツ)_/¯
        end
    end
    result
end

checkerP1(d::AbstractDict) = all(key -> haskey(d, key), ("byr","iyr","eyr","hgt","hcl","ecl","pid"))

part1(input) = count(checkerP1, input)

function checkerP2(d::AbstractDict)
    !checkerP1(d) && return false

    byr = tryparse(Int, d["byr"]) ; (byr === nothing || byr < 1920 || byr > 2002) && return false

    iyr = tryparse(Int, d["iyr"]) ; (iyr === nothing || iyr < 2010 || iyr > 2020) && return false

    eyr = tryparse(Int, d["eyr"]) ; (eyr === nothing || eyr < 2020 || eyr > 2030) && return false

    hgt = match(r"^(\d+)(in|cm)$", d["hgt"]) ; hgt === nothing && return false
    h, unit = parse(Int, hgt.captures[1]), hgt.captures[2]
    unit == "cm" && (h < 150 || h > 193) && return false
    unit == "in" && (h < 59 || h > 76) && return false

    hcl = match(r"^#[a-f0-9]{6}$", d["hcl"]) ; hcl === nothing && return false

    d["ecl"] in ("amb", "blu", "brn", "gry", "grn", "hzl", "oth") || return false

    pid = match(r"^\d{9}$", d["pid"]) ; pid === nothing && return false

    true
end

part2(input) = count(checkerP2, input)

end


text:
module Scrabble

include("setup.jl")
include("iterators.jl")
include("anagrams.jl")
include("jumble.jl")

end # module


text:
module Memoize
using MacroTools: isexpr, combinedef, namify, splitarg, splitdef
export @memoize, memoize_cache

cache_name(f) = Symbol("##", f, "_memoized_cache")

macro memoize(args...)
    if length(args) == 1
        dicttype = :(IdDict)
        ex = args[1]
    elseif length(args) == 2
        (dicttype, ex) = args
    else
        error("Memoize accepts at most two arguments")
    end

    cache_dict = isexpr(dicttype, :call) ? dicttype : :(($dicttype)())

    def_dict = try
        splitdef(ex)
    catch
        error("@memoize must be applied to a method definition")
    end

    # a return type declaration of Any is a No-op because everything is <: Any
    rettype = get(def_dict, :rtype, Any)
    f = def_dict[:name]
    def_dict_unmemoized = copy(def_dict)
    def_dict_unmemoized[:name] = u = Symbol("##", f, "_unmemoized")

    args = def_dict[:args]
    kws = def_dict[:kwargs]
    # Set up arguments for tuple
    tup = [splitarg(arg)[1] for arg in vcat(args, kws)]

    # Set up identity arguments to pass to unmemoized function
    identargs = map(args) do arg
        arg_name, typ, slurp, default = splitarg(arg)
        if slurp || namify(typ) === :Vararg
            Expr(:..., arg_name)
        else
            arg_name
        end
    end
    identkws = map(kws) do kw
        arg_name, typ, slurp, default = splitarg(kw)
        if slurp
            Expr(:..., arg_name)
        else
            Expr(:kw, arg_name, arg_name)
        end
    end

    fcachename = cache_name(f)
    mod = __module__
    fcache = isdefined(mod, fcachename) ?
             getfield(mod, fcachename) :
             Core.eval(mod, :(const $fcachename = $cache_dict))

    body = quote
        get!($fcache, ($(tup...),)) do
            $u($(identargs...); $(identkws...))
        end
    end

    if length(kws) == 0
        def_dict[:body] = quote
            $(body)::Core.Compiler.return_type($u, typeof(($(identargs...),)))
        end
    else
        def_dict[:body] = body
    end

    esc(quote
        $(combinedef(def_dict_unmemoized))
        empty!($fcache)
        Base.@__doc__ $(combinedef(def_dict))
    end)

end

function memoize_cache(f::Function)
    # This will fail in certain circumstances (eg. @memoize Base.sin(::MyNumberType) = ...) but I don't think there's
    # a clean answer here, because we can already have multiple caches for certain functions, if the methods are
    # defined in different modules.
    getproperty(parentmodule(f), cache_name(f))
end

end


text:
module disambig
export adagramDisambig

using AdaGram

GC.enable(false)
vm, dict = load_model("model/model_output_full.out");

function adagramDisambig(target_word, context)

	counter = 0
	maxval = 0.0
	maxcnt = 0

	arr = disambiguate(vm, dict, target_word, split(context))

	for w in arr
		counter += 1
		if (w > maxval)
			maxval = w
			maxcnt = counter
		end
	end
	return maxcnt
end

end

max_stars_repo_name: benandow/UiRef

max_issues_repo_name: benandow/UiRef

max_forks_repo_name: benandow/UiRef


text:
using NumericalMethodsforEngineers

#=
@sym begin
  LagrangePolynomial(xi_, yi_) := Module([N, sum, num, den],
    begin
      N = Length(xi)
      sum = 0
      For( i=1, i <= N, Increment(i),
        begin
          num = 1
          den = 1
          For( j=1, j <= N, Increment(j),
            begin
              If(j != i, num *= (x-xi[j]))
              If(j != i, den *= (xi[i]-xi[j]))
            end
          )
          sum += yi[i] * num/den
        end
      ),
      Return(Simplify(sum))
    end
  )
end
=#

println("\n\n--------1------------\n")

@sym begin
  ClearAll(xi, yi, Y, R, sol, ytilde)
  xi = [0, 1//2, 1]
  yi = [0, a, 1]
  Y(x_) := LagrangePolynomial(xi, yi)
  Println(Y(x))
end

println("\n---------2-----------\n")

@sym begin
  ClearAll(xi, yi, Y1, R, sol, ytilde)
  Y1(x_) := x^3 -3x^2 +6x - 9
  R(x_) = Simplify(D(Y1(x), x, 2))
  Println(R(x))
end

println("\n----------3----------\n")

@sym begin
  ClearAll(xi, yi, Y, R, sol, ytilde)
  xi = [0, 1//2, 1]
  yi = [0, a, 1]
  Y(x_) := LagrangePolynomial(xi, yi)
  R(x_) = Simplify(D(Y(x), x, 2) - 3*x - 4*Y(x))
  R(x_) = Simplify(R(x) ./ (a => -C1/4))
  Println(R(x))
  SetJ(r, ToString(Simplify(R(x))))
end

println("\n----------4----------\n")

println(r)

println("\n----------5----------\n")

println("\n\nExample 7.14: y'' = 3x + 4y, y(0)=0, y(1)=1")
println("Residual for Weighted Residual Method using 2 point Lagragian Polynomial")
@sym Println("\nY(x): ", Expand(Y(x)), "\n")
println("( Example 7.14 gives: ytilde = 1/2*x^3*(27a-27b+9)-x^2*(45a-36b+9)+x*(18a-9b+2) )")
@sym Println("\nR(x) = ", R(x), "\n")


text:
module CellQuadratures

# Dependencies of this module

using Gridap
using Gridap.Helpers

# Functionality provided by this module

export CellQuadrature
export ConstantCellQuadrature
import Gridap: coordinates
import Gridap: weights

# Abstract types and interfaces

"""
Abstract type representing a collection of quadratures, one for each cell
"""
const CellQuadrature{D} = CellValue{<:Quadrature{D}}

coordinates(::CellQuadrature{D} where D )::CellPoints{D} = @abstractmethod

weights(::CellQuadrature)::CellValues{Float64} = @abstractmethod

function CellQuadrature end

# Pretty printing

import Base: show

function show(io::IO,self::CellQuadrature{D}) where {D}
  print(io,"CellQuadrature object")
end

function show(io::IO,::MIME"text/plain",quad::CellQuadrature{D}) where {D}
  show(io,quad)
  print(io,":")
  print(io,"\n dim: $D")
  print(io,"\n ncells: $(length(quad))")
end

# Concrete structs

"""
A concrete implementation of CellQuadrature for the particular case
that all cells have the same quadrature
"""
const ConstantCellQuadrature{D} = ConstantCellValue{Quadrature{D}}

function ConstantCellQuadrature(quad::Quadrature{D}, l::Int) where D
  ConstantCellQuadrature{D}(quad, l)
end

function coordinates(self::ConstantCellQuadrature)
  ConstantCellValue(coordinates(self.value), self.length)
end

function weights(self::ConstantCellQuadrature)
  ConstantCellValue(weights(self.value), self.length)
end

end # module CellQuadratures


text:
## sigmoidal functions

export logistic, logistic_prime

logistic(x) = 1/(1 + exp(-x))
logistic_prime(x) = exp(-x)/(1 + exp(-x))^2


text:
# SQL identifier (possibly qualified).

mutable struct IdentifierClause <: AbstractSQLClause
    over::Union{SQLClause, Nothing}
    name::Symbol

    IdentifierClause(;
                     over = nothing,
                     name::Union{Symbol, AbstractString}) =
        new(over, Symbol(name))
end

IdentifierClause(name; over = nothing) =
    IdentifierClause(over = over, name = name)

"""
    ID(; over = nothing, name)
    ID(name; over = nothing)

A SQL identifier.  Specify `over` or use the `|>` operator to make a qualified
identifier.

# Examples

```jldoctest
julia> c = ID(:person);

julia> print(render(c))
"person"
```

```jldoctest
julia> c = ID(:p) |> ID(:birth_datetime);

julia> print(render(c))
"p"."birth_datetime"
```
"""
ID(args...; kws...) =
    IdentifierClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(ID), pats::Vector{Any}) =
    dissect(scr, IdentifierClause, pats)

Base.convert(::Type{AbstractSQLClause}, name::Symbol) =
    IdentifierClause(name)

Base.convert(::Type{AbstractSQLClause}, qname::Tuple{Symbol, Symbol}) =
    IdentifierClause(qname[2], over = IdentifierClause(qname[1]))

function PrettyPrinting.quoteof(c::IdentifierClause, ctx::QuoteContext)
    ex = Expr(:call, nameof(ID), quoteof(c.name))
    if c.over !== nothing
        ex = Expr(:call, :|>, quoteof(c.over, ctx), ex)
    end
    ex
end

rebase(c::IdentifierClause, c′) =
    IdentifierClause(over = rebase(c.over, c′), name = c.name)


text:
if !@isdefined(TabLevel)
    TabLevel = ""
end
println(TabLevel*"Open OpEqual.jl")
TabLevel=TabLevel*"    "

import Base.==

==(x::Union{MajoranaOp,MajoranaSp,MajElem},
   y::Union{MajoranaOp,MajoranaSp,MajElem}) = do_OpEqual(x,y)

==(x::Number,y::Union{MajoranaOp,MajoranaSp,MajElem}) = do_OpEqual(x,y)
==(y::Union{MajoranaOp,MajoranaSp,MajElem},x::Number) = do_OpEqual(y,x)

function OpEqual(a,b;strict=true)
    same = do_OpEqual(a,b)
    if strict && !same
        throw(ErrorException("a="*OpToStr(a)*"\nis not equal to\nb="*OpToStr(b)))
    end
    return same
end

isapprox(a::MajoranaOp,b::MajoranaOp;atol=sqrt(eps())) = do_OpEqual(a,b,atol=atol)
isapprox(a::MajoranaSp,b::MajoranaOp;atol=sqrt(eps())) = do_OpEqual(a,b,atol=atol)
isapprox(a::MajoranaOp,b::MajoranaSp;atol=sqrt(eps())) = do_OpEqual(a,b,atol=atol)
isapprox(a::MajoranaSp,b::MajoranaSp;atol=sqrt(eps())) = do_OpEqual(a,b,atol=atol)

do_OpEqual(b::MajoranaOp,a::Number)=do_OpEqual(a,b)
function do_OpEqual(a::Number,b::MajoranaOp)
    ###There two can only be the same if the majorna operator is actually a number (of if both are zero)
    if a == 0
        ###If the scale is zero they are the same
        return b.scale == 0
    else
        ##They can still be zero if there the operators list is empty and there are ket numbers
        if b.KetState
            return false
        else
            return a == b.scale && length(b.OpList)==0
        end
    end
end

do_OpEqual(b::MajoranaSp,a::Number) =  do_OpEqual(a,b)
do_OpEqual(a::Number,b::MajoranaSp) = (MajScalar(a) == b)

function do_OpEqual(a::MajoranaOp,b::MajoranaOp;atol=0)
    #println("1: Check a="*OpToStr(a)*" vs b="*OpToStr(b))
    same=OpProp(a,b)
    if same
        if atol!=0
            return isapprox(a.scale,b.scale,atol=atol)
        else
            return (a.scale==b.scale)
        end
    else
        return false
    end
end

do_OpEqual(a::MajoranaSp,b::MajoranaOp;atol=0) = do_OpEqual(b,a,atol=atol)
function do_OpEqual(a::MajoranaOp,b::MajoranaSp;atol=0)
    #println("2: Check:\na=$a vs \nb=$b")
    if length(b.SpList)!=1
        if atol==0
            return false
        else ###It could be that some elements are (approximately) equal to zero
            ###That should be ok
            EqOp=false
            for j in 1:length(b.SpList)
                #println("Tets j=$j")
                TestEqOp=do_OpEqual(b.SpList[j],a,atol=atol)
                if TestEqOp ###This was equal to the one term
                    EqOp=true
                else ###if now it zould be approximately zero
                    TestZeroOp=isapprox(b.SpList[j].scale,0.0,atol=atol)
                    if !TestZeroOp
                        return false
                    end
                end
            end
            return EqOp
        end
    else
        return do_OpEqual(b.SpList[1],a,atol=atol)
    end
end
function do_OpEqual(a::MajoranaSp,b::MajoranaSp;atol=0)
    #println("3: Check a=\n"*OpToStr(a)*" vs b=\n"*OpToStr(b))
    if length(b.SpList)!=length(a.SpList)
        return false
    else
        for indx in 1:length(a.SpList)
            if !do_OpEqual(a.SpList[indx],b.SpList[indx],atol=atol)
                return false
            end
        end
    end
    return true
end

function do_OpEqual(a::MajElem,b::MajElem)
    return a.Type == b.Type && a.Index == b.Index
end

do_OpEqual(a::MajElem,b::MajoranaOp) = MajoranaOp(1,[a])==b
do_OpEqual(b::MajoranaOp,a::MajElem) = MajoranaOp(1,[a])==b

TabLevel=TabLevel[1:end-4]
println(TabLevel*"Close OpEqual.jl")


text:
using ImagineHardware
using Test

include("piezos.jl")
#include("amplifiers.jl")
#include("pid.jl")
#include("piezosystems.jl")


text:
@testset "Errors" begin
  @testset "Error message for Humans(TM)" begin
    let errmsg = "Couldn't resolve host name"
      server = "smtp://nonexists"
      body = IOBuffer("test")

      try
        send(server, ["nobody@earth"], "nobody@earth", body)
        @assert false, "send should fail"
      catch e
        @test occursin(string(errmsg), string(e))
      end
    end
  end

  @testset "Non-blocking send" begin
    let errmsg = "Couldn't resolve host name"
      server = "smtp://nonexists"
      body = IOBuffer("test")

      t = @async send(server, ["nobody@earth"], "nobody@earth", body)
      try
        wait(t)
      catch e
        @test e isa TaskFailedException
        @test occursin(errmsg, e.task.exception.msg)
      end
    end
  end
end  # @testset "Errors"


text:
# function annotate(minefield::Array{<:AbstractString,1}) # Sensible IMO but fails the first test
function annotate(minefield::Vector)
    isempty(minefield) && return minefield

    # Julia gets upset / slow here as matrix can be Array{,2} or Array{,1}
    # TODO: work out how to ensure it is always Array{,2} for hopeful speedup
    matrix = map(x->x=='*',mapreduce(collect,hcat,minefield))  # |> permutedims |> permutedims (this is a crazy way of doing this)
    flagged = flag_mines(matrix)
    [[d == -1 ? '*' : d == 0 ? ' ' : codepoint('0') + Char(d) for d in s] |> String for s in eachcol(flagged)]
end

# This function works for hypercube minefields too, which is pretty cool.
# If you decide to construct your own hypercube minefield, bear in mind that
# the curse of dimensionality means that mines become useless as the number of
# dimensions increases to even moderate numbers.
#
# (A more useful metric for 'danger from mines' is the percentage of neighbouring
# cells which contain mines).
function flag_mines(matrix::Array)
    flagged = zeros(Int,size(matrix))
    @inbounds for inds in Tuple.(CartesianIndices(matrix))
        flagged[inds...] = matrix[inds...] == 1 ? -1 : sum(window(matrix,inds,ones(Int,length(size(matrix)))))
    end
    flagged
end

function window(mat,coords,δs)
    dims = size(mat)
    view(mat, (max(1,c-δ):min(dims[dim],c+δ) for (dim,(δ,c)) in enumerate(zip(δs,coords)))...)
end

# Testing
using BenchmarkTools
mines(n,m) = [rand([repeat(' ',4)...; '*'],m) |> String for _ in 1:n]
@benchmark annotate(m) setup=(m=mines(20,20))


text:
## Test naive implementation of entropy calculations towards efficient calculation of entropy on a manifold

using LinearAlgebra
# using DocStringExtensions
using KernelDensityEstimate

using Cairo
using Gadfly, Colors
using Random

using ApproxManifoldProducts

using TransformUtils
const TU = TransformUtils

setForceEvalDirect!(true)

## single Bimodal density

pts = [0.1*randn(50); TU.wrapRad.(0.1*randn(50).+pi)]
shuffle!(pts)

pp = kde!(pts, [0.05], (addtheta,), (difftheta,))

plotKDECircular(pp)

##

pp = kde!(pts, [0.05], (addtheta,), (difftheta,))

x = 2pi*rand(1,1000).-pi
y = pp(x)

Gadfly.plot(x=x, y=y, Geom.line)

## create two densities

pts1 = TU.wrapRad.(0.1*randn(100).-pi.+0.5)
pts2 = TU.wrapRad.(0.1*randn(100).+pi.-0.5)

# pc1 = kde!_CircularNaiveCV(pts1)
# pc2 = kde!_CircularNaiveCV(pts2)
pc1 = kde!(pts1, [0.1;], (addtheta,) , (difftheta,))
pc2 = kde!(pts2, [0.1;], (addtheta,) , (difftheta,))

pl = plotKDECircular([pc1;pc2])

## Calculate the approximate product between the densities

# TODO: make circular KDE
dummy = kde!(rand(100),[1.0;], (addtheta,), (difftheta,));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2], nothing, nothing, Niter=2,
                          addop=(addtheta,), diffop=(difftheta,), getMu=(getCircMu,));

pc12 = kde!_CircularNaiveCV(pGM[:])
getBW(pc12)[1,1]
# pc12 = kde!(pGM, [0.05;], (addtheta,), (difftheta,))

pl = plotKDECircular([pc1;pc2; pc12], scale=0.07, c=["cyan";"blue";"red"])

# x = reshape(collect(range(-pi, stop=pi, length=1000)),1,:)
# y = pc12(x, false, 1e-3, (addtheta,), (difftheta,))
# Gadfly.plot(x=x, y=y, Geom.line)

# Gadfly.push_theme(:default)
pl |> PDF("/tmp/test.pdf",10cm,10cm)
@async run(`evince /tmp/test.pdf`)

pl.coord = Coord.Cartesian(xmin = -1.8,
xmax = 1.8,
ymin = -1.8,
ymax = 1.8,
aspect_ratio=1.0)

pl.theme
pl.theme.line_width = 2.0mm

pl
## create two densities at +Y

N = 100

pts1 = TU.wrapRad.(0.1*randn(N).+pi/2.0.+0.5)
pts2 = TU.wrapRad.(0.1*randn(N).+pi/2.0.-0.5)

# pc1 = kde!_CircularNaiveCV(pts1)
# pc2 = kde!_CircularNaiveCV(pts2)
pc1 = kde!(pts1, [0.1;], (addtheta,), (difftheta,))
pc2 = kde!(pts2, [0.1;], (addtheta,), (difftheta,))

# TODO: make circular KDE
dummy = kde!(rand(N),[1.0;], (addtheta,), (difftheta,));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2], nothing, nothing, Niter=2,
                          addop=(addtheta,), diffop=(difftheta,), getMu=(getCircMu,));

pc12 = kde!_CircularNaiveCV(pGM[:])
# getBW(pc12)[1,1]
# pc12 = kde!(pGM, [0.1;], addtheta, difftheta)

pl = plotKDECircular([pc1;pc2; pc12])

## Different two

N=100

pts1 = TU.wrapRad.(0.1*randn(N).-0.5)
pts2 = TU.wrapRad.(0.1*randn(N).+0.5)

# pc1 = kde!_CircularNaiveCV(pts1)
# pc2 = kde!_CircularNaiveCV(pts2)
pc1 = kde!(pts1, [0.1;], (addtheta,), (difftheta,))
pc2 = kde!(pts2, [0.1;], (addtheta,), (difftheta,))

# TODO: make circular KDE
dummy = kde!(rand(N),[1.0;], (addtheta,), (difftheta,));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2], nothing, nothing, Niter=1, addop=(addtheta,), diffop=(difftheta,), getMu=(getCircMu,));

# pc12 = kde!_CircularNaiveCV(pGM[:])
# getBW(pc12)[:,1][1]

pc12 = kde!(pGM,[0.05], (addtheta,), (difftheta,));

pl = plotKDECircular([pc1;pc2; pc12])

## another two equally spaced

pts1 = TU.wrapRad.(0.1*randn(100).-pi/2.0)
pts2 = TU.wrapRad.(0.1*randn(100).+pi/2.0)

# pc1 = kde!_CircularNaiveCV(pts1)
# pc2 = kde!_CircularNaiveCV(pts2)
pc1 = kde!(pts1, [0.1;], (addtheta,), (difftheta,))
pc2 = kde!(pts2, [0.1;], (addtheta,), (difftheta,))

# TODO: make circular KDE
dummy = kde!(rand(100),[1.0;], (addtheta,), (difftheta,));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2], nothing, nothing, Niter=2, addop=(addtheta,), diffop=(difftheta,), getMu=(getCircMu,));

pc12 = kde!_CircularNaiveCV(pGM[:])
pl = plotKDECircular([pc1;pc2; pc12])

## handrolically trying product of 3

pts1 = TU.wrapRad.(0.1*randn(100))
pts2 = TU.wrapRad.(0.1*randn(100).-2pi/3.0.+1.0)
pts3 = TU.wrapRad.(0.1*randn(100).+2pi/3.0)

# pc1 = kde!_CircularNaiveCV(pts1)
# pc2 = kde!_CircularNaiveCV(pts2)
pc1 = kde!(pts1, [0.1;], (addtheta,), (difftheta,))
pc2 = kde!(pts2, [0.1;], (addtheta,), (difftheta,))
pc3 = kde!(pts3, [0.1;], (addtheta,), (difftheta,))

# pl = plotKDECircular([pc1;pc2;pc3])

dummy = kde!(rand(100),[1.0;], (addtheta,), (difftheta,));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2; pc3], nothing, nothing, Niter=1, addop=(addtheta,), diffop=(difftheta,), getMu=(getCircMu,));

pc123 = kde!_CircularNaiveCV(pGM[:])

pl = plotKDECircular([pc1;pc2;pc3; pc123])

## product of 3

pts1 = TU.wrapRad.(0.1*randn(100))
pts2 = TU.wrapRad.(0.1*randn(100).-2pi/3.0.+1.0)
pts3 = TU.wrapRad.(0.1*randn(100).+2pi/3.0.-1.0)

# pc1 = kde!_CircularNaiveCV(pts1)
# pc2 = kde!_CircularNaiveCV(pts2)
pc1 = kde!(pts1, [0.1;], (addtheta,), (difftheta,))
pc2 = kde!(pts2, [0.1;], (addtheta,), (difftheta,))
pc3 = kde!(pts3, [0.1;], (addtheta,), (difftheta,))

# pl = plotKDECircular([pc1;pc2;pc3])

dummy = kde!(rand(100),[1.0;], (addtheta,), (difftheta,));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2; pc3], nothing, nothing, Niter=1, addop=(addtheta,), diffop=(difftheta,), getMu=(getCircMu,));

pc123 = kde!_CircularNaiveCV(pGM[:])

pl = plotKDECircular([pc1;pc2;pc3; pc123])

## product of 3

pts1 = TU.wrapRad.(0.1*randn(100))
pts2 = TU.wrapRad.(0.1*randn(100).-2pi/3.0)
pts3 = TU.wrapRad.(0.1*randn(100).+2pi/3.0)

# pc1 = kde!_CircularNaiveCV(pts1)
# pc2 = kde!_CircularNaiveCV(pts2)
pc1 = kde!(pts1, [0.1;], (addtheta,), (difftheta,))
pc2 = kde!(pts2, [0.1;], (addtheta,), (difftheta,))
pc3 = kde!(pts3, [0.1;], (addtheta,), (difftheta,))

# pl = plotKDECircular([pc1;pc2;pc3])

dummy = kde!(rand(100),[1.0;], (addtheta,), (difftheta,));

pGM, = prodAppxMSGibbsS(dummy, [pc1; pc2; pc3], nothing, nothing, Niter=1, addop=(addtheta,), diffop=(difftheta,), getMu=(getCircMu,));

pc123 = kde!_CircularNaiveCV(pGM[:])

pl = plotKDECircular([pc1;pc2;pc3; pc123])

##

Gadfly.push_theme(:default)

pl |> PDF("/tmp/test.pdf", 20cm,15cm)

#


text:
@testset "julian_versions" begin

    # test filter
    sc = SparkContext(master="local")

    rdd = parallelize(sc, 1:30)
    pairs = map_pair(rdd, x -> (x,1))

    f = (prt, idx) -> sum(prt)
    @test (map_partitions_with_index(rdd, f) |> collect ==
           map_partitions_with_index(f, rdd) |> collect)

    f = prt -> sum(prt)
    @test (map_partitions(rdd, f) |> collect ==
           map_partitions(f, rdd) |> collect)

    f = prt -> ((x, 1) for x in prt)
    @test (map_partitions_pair(rdd, f) |> collect ==
           map_partitions_pair(f, rdd) |> collect)

    f = x -> x + 1
    @test map(rdd, f) |> collect == map(f, rdd) |> collect

    f = x -> (x,1)
    @test map_pair(rdd, f) |> collect == map_pair(f, rdd) |> collect

    f = x -> [x, x]
    @test flat_map(rdd, f) |> collect == flat_map(f, rdd) |> collect

    # currently broken
    # f -> x -> [x, x]
    # @test flat_map_pair(rdd, f) |> collect == flat_map_pair(f, rdd) |> collect

    f = x -> mod(x, 2) == 1
    @test filter(rdd, f) |> collect == filter(f, rdd) |> collect

    # currently broken
    # f = (x, y) -> mod(x, 2) == 1
    # @test filter(pairs, f) |> collect == filter(f, pairs) |> collect

    f = (+)
    @test reduce(rdd, f) |> collect == reduce(f, rdd) |> collect

    pairs2 = map_pair(parallelize(sc, ["a", "b", "a", "c"]), x -> (x, 1))
    f = (+)
    @test reduce_by_key(pairs2, f) |> collect == reduce_by_key(f, pairs2) |> collect

    close(sc)

end


text:
function lolrng(x    :: Integer,
                y    :: Integer,
                z    :: Integer,
                seed :: Integer)
    r1 = UInt32(x) * UInt32(0x1B873593)
    r2 = UInt32(y) * UInt32(0x19088711)
    r3 = UInt32(z) * UInt32(0xB2D05E13)

    r  = UInt32(seed) + r1 + r2 + r3
    r ⊻= r >> UInt32(5)
    r *= UInt32(0xCC9E2D51)
    return r / 0xffffffff
end

function octave(x    :: Real,
                y    :: Real,
                z    :: Real,
                oct  :: Integer,
                seed :: Integer)
    divisor = 2.0^(-oct)

    xi = (x / divisor) |> floor |> Int
    yi = (y / divisor) |> floor |> Int
    zi = (z / divisor) |> floor |> Int

    δx = rem(x, divisor) / divisor
    δy = rem(y, divisor) / divisor
    δz = rem(z, divisor) / divisor

    v000 = lolrng(xi,     yi,     zi, seed)
    v001 = lolrng(xi + 1, yi,     zi, seed)
    v010 = lolrng(xi,     yi + 1, zi, seed)
    v011 = lolrng(xi + 1, yi + 1, zi, seed)

    v100 = lolrng(xi,     yi,     zi + 1, seed)
    v101 = lolrng(xi + 1, yi,     zi + 1, seed)
    v110 = lolrng(xi,     yi + 1, zi + 1, seed)
    v111 = lolrng(xi + 1, yi + 1, zi + 1, seed)

    inter(v1, v2, x) = v1 + (v2 - v1)*x
    v00 = inter(v000, v001, δx)
    v01 = inter(v010, v011, δx)
    v10 = inter(v100, v101, δx)
    v11 = inter(v110, v111, δx)

    v0 = inter(v00, v01, δy)
    v1 = inter(v10, v11, δy)

    v = inter(v0, v1, δz)
    return v
end

"""
    value_noise(x :: Real, y :: Real, z :: Real, octaves :: Integer, seed :: Integer)

Calculate value noise at the point (x,y,z).

Calaculated noise consists of `octaves` octaves (higher
harmonics). `seed` is a seed for random number generator.

# Examples
```jddoctest
julia> [value_noise(i/10, 0, 0, 10, 1343) for i in 1:20]
20-element Vector{Float64}:
 0.6782011856880998
 0.6811537717432827
 0.6749341884205909
 0.6837415201436138
 0.7459242823636113
 0.6684127212046574
 0.6642874704554677
 0.6833569738117538
 0.7214221918747463
 0.8051134343235813
 0.6986164658038979
 0.647098867616866
 0.6158556023606773
 0.6356918681071072
 0.6617448759515002
 0.6710733979362564
 0.7094478163638382
 0.7446686007256124
 0.7987984329263587
 0.9086695925335629
```
"""
function value_noise(x       :: Real,
                     y       :: Real,
                     z       :: Real,
                     octaves :: Integer,
                     seed    :: Integer)
    mapreduce((x, octave) -> x ./ 2.0^octave, +,
              (octave(x, y, z, o, seed) for o in 0:octaves-1),
              countfrom(0)) ./ 2*(1 - 2.0^(-octaves))
end


text:
multiple_of(x, n) = isinteger(x/n)
multiple_of(n) = MultipleOf(n)


text:
abstract type AbstractChebyshevBasis{P} <: AbstractGegenbauerBasis{P} end

polynomial_type(::Type{<:AbstractChebyshevBasis}, V::Type) = MP.polynomialtype(V, Float64)

reccurence_first_coef(::Type{<:AbstractChebyshevBasis}, degree) = 2
reccurence_third_coef(::Type{<:AbstractChebyshevBasis}, degree) = -1
reccurence_deno_coef(::Type{<:AbstractChebyshevBasis}, degree) = 1

"""
    struct ChebyshevBasisFirstKind{P} <: AbstractChebyshevBasis{P}
        polynomials::Vector{P}
    end

Orthogonal polynomial with respect to the univariate weight function ``w(x) = \\frac{1}{\\sqrt{1 - x^2}}`` over the interval ``[-1, 1]``.
"""
struct ChebyshevBasisFirstKind{P} <: AbstractChebyshevBasis{P}
    polynomials::Vector{P}
end

const ChebyshevBasis{P} = ChebyshevBasisFirstKind{P}

degree_one_univariate_polynomial(::Type{<:ChebyshevBasisFirstKind}, variable::MP.AbstractVariable) = MA.@rewrite(variable + 0)

"""
    struct ChebyshevBasisSecondKind{P} <: AbstractChebyshevBasis{P}
        polynomials::Vector{P}
    end

Orthogonal polynomial with respect to the univariate weight function ``w(x) = \\sqrt{1 - x^2}`` over the interval ``[-1, 1]``.
"""
struct ChebyshevBasisSecondKind{P} <: AbstractChebyshevBasis{P}
    polynomials::Vector{P}
end

degree_one_univariate_polynomial(::Type{<:ChebyshevBasisSecondKind}, variable::MP.AbstractVariable) = MA.@rewrite(2variable + 0)

max_stars_repo_name: JuliaAlgebra/PolynomialBases

max_issues_repo_name: JuliaAlgebra/PolynomialBases

max_forks_repo_name: JuliaAlgebra/PolynomialBases


text:
module DataFormatter

using DataFrames: DataFrame, groupby, Not, insertcols!

include("types.jl")
include("simulation.jl")
include("vizualize.jl")

end

max_stars_repo_name: Leooop/DataFormatter

max_issues_repo_name: Leooop/DataFormatter

max_forks_repo_name: Leooop/DataFormatter


text:
module test_julia_stdlib_unicode

using Test
using Unicode

@test Unicode.normalize("거북이", :NFD) == "거북이"
@test Unicode.normalize("거북이", :NFC) == "거북이"

end # module test_julia_stdlib_unicode


text:
module DNC

using Flux
using LinearAlgebra
using Zygote
using TensorCast

export Dnc

include("utils.jl")
include("access.jl")
include("addressing.jl")
include("computer.jl")

end # module


text:
struct ModelParam{T, S <: NTuple{N, Int} where N}
  m::T
  log_s::T
  support::String # real, unit, simplex, positive
  size::S
end

# scalar param
function ModelParam(support::String; m=nothing, s=nothing)
  @assert(support in ["real", "unit", "positive"])
  m = param(m==nothing ? rand(1) : m)
  log_s = param(s==nothing ?  randn(1) : log.(s))
  return ModelParam(m, log_s, support, (1, ))
end

# Vector param
function ModelParam(K::Integer, support::String; m=nothing, s=nothing)
  @assert(support in ["real", "unit", "simplex", "positive"])
  m = param(m==nothing ? rand(K) : m)
  log_s = param(s==nothing ? randn(K) : log.(s))
  return ModelParam(m, log_s, support, (K, ))
end

# ND-Array param
function ModelParam(D::Tuple, support::String; m=nothing, s=nothing)
  @assert(support in ["real", "unit", "simplex", "positive"])
  m = param(m==nothing ? rand(D...) : m)
  log_s = param(s==nothing ? randn(D...) : log.(s))
  return ModelParam(m, log_s, support, D)
end

function log_q(mp::ModelParam, r::R; args...) where R
  m, s = vp(mp; args...)
  return sum(lpdf_normal.(r, m, s))
end

"""
Get variational parameters
"""
function vp(mp::ModelParam; m_min::Float64=-10.0, m_max::Float64=10.0, s_max::Float64=10.0)
  if mp.support in ("unit", "simplex")
    param_range = m_max - m_min
    m = sigmoid.(mp.m) * param_range .+ m_min
    s = sigmoid.(mp.log_s) * s_max
  else
    m = mp.m
    s = exp.(mp.log_s)
  end

  return (m, s)
end

function logabsdetJ(mp::ModelParam, real::R, tran::T) where {R, T}
  if mp.support == "simplex"
    return SB_logabsdetJ(real, tran)
  elseif mp.support == "unit"
    return log.(tran) + log1p.(-tran)
  elseif mp.support == "positive"
    return real
  elseif mp.support == "real"
    return zero(mp.m)
  else
    ErrorException("ADVI.ModelParam.logabsdet is not implemented for support=$(mp.support)")
  end
end

function transform(mp::ModelParam, real::T) where T
  if mp.support == "simplex"
    return SB_transform(real)
  elseif mp.support == "unit"
    return sigmoid.(real)
  elseif mp.support == "positive"
    return exp.(real)
  elseif mp.support == "real"
    return real
  else
    ErrorException("ADVI.ModelParam.transform is not implemented for support=$(mp.support)")
  end
end

"""
Reparameterized sampling from variational distribution
"""
function rsample(mp::ModelParam)
  # TODO: Optimize this
  m, s = vp(mp)
  return randn(mp.size) .* s + m
end

#= TEST
mp = ModelParam(3, "unit")
vp(mp)
rsample(mp)
=#

"""
Get variational parameters
"""
vparams(mp::ModelParam) = Flux.params(mp.m, mp.log_s)

function vparams(s::S) where S
  ps = []
  for key in fieldnames(S)
    if isdefined(s, key)
      f = getfield(s, key)
      if typeof(f) <: ModelParam
        append!(ps, [f.m])
        append!(ps, [f.log_s])
      elseif Tracker.istracked(f)
        append!(ps, [f])
      end
    else
      ErrorException("Field $key is undefined in $(S.name)")
    end
  end
  return Flux.params(ps...)
end


text:
using EffectiveWaves, Test

@testset "match wave low volfrac" begin

using LinearAlgebra
## Low volume fraction
    medium = Medium(1.0,1.0+0.0im)
    specie = Specie(ρ=0.5, r=0.4, c=0.5, volfrac=0.001)
    species = [specie]
    θin = 0.2
    tol = 1e-9
    hankel_order = 2

    ωs = [0.2,1.2]
    radius_multiplier = 1.005

    wave_effs_arr = [
        effective_waves(ω, medium, [specie];
            hankel_order=hankel_order, tol = tol,  θin = θin,
            radius_multiplier = radius_multiplier,
            num_wavenumbers = 1
            #, mesh_points = 10, mesh_size = 2.0 #, max_Rek = 20.0, max_Imk = 20.0
            , extinction_rescale = false)
    for ω in ωs]

    k_eff_φs = wavenumber_low_volfrac(ωs, medium, species; tol=tol,
        hankel_order=hankel_order,
        radius_multiplier = radius_multiplier)

    wave_eff_φs = [
        EffectiveWave(ωs[i], k_eff_φs[i], medium, species; tol = tol,
            hankel_order=hankel_order, θin = θin,
            extinction_rescale = true,
            radius_multiplier = radius_multiplier)
    for i in eachindex(ωs)]

    @test maximum(abs(wave_effs_arr[i][1].k_eff/k_eff_φs[i] - 1.0) for i in eachindex(ωs)) < 5e-8

    ds = [abs(dot(wave_eff_φs[i].amplitudes[:], wave_effs_arr[i][1].amplitudes[:])) for i in eachindex(ωs)]
    ds2 = [norm(wave_eff_φs[i].amplitudes[:])*norm(wave_effs_arr[i][1].amplitudes[:]) for i in eachindex(ωs)]
    @test maximum(ds./ds2 .- 1.0) < tol

    match_ws = [
        MatchWave(ωs[i], medium, specie; θin = θin,
            radius_multiplier = radius_multiplier,
            max_size=100,
            tol = tol, wave_effs = wave_effs_arr[i])
    for i in eachindex(ωs)]

    @test maximum(match_error.(match_ws)) < 1e-6

    avg_wave_φs = [
        AverageWave(match_ws[i].average_wave.x, wave_eff_φs[i])
    for i in eachindex(ωs)]

    Rs_near = [reflection_coefficient(ωs[i], match_ws[i].average_wave, medium, specie) for i in eachindex(ωs)]
    Rs_near_φ = [reflection_coefficient(ωs[i], avg_wave_φs[i], medium, specie) for i in eachindex(ωs)]
    @test maximum(abs.(Rs_near_φ - Rs_near)) < 5e-6

    Rs = [reflection_coefficient(ωs[i], match_ws[i], medium, specie) for i in eachindex(ωs)]
    Rs_φ = [reflection_coefficient(ωs[i], wave_eff_φs[i], medium, species) for i in eachindex(ωs)]
    @test maximum(abs.(Rs_near_φ - Rs_near)) < 5e-6

    @test maximum(
        norm(match_ws[i].average_wave.amplitudes[:] .- avg_wave_φs[i].amplitudes[:])/norm(avg_wave_φs[i].amplitudes[:])
    for i in eachindex(ωs)) < 1e-3
end


text:
using PlasmoPlots
using Test

@testset "PlasmoPlots.jl" begin
    @testset "$(file)" for file in filter(f -> endswith(f, ".jl"), readdir(@__DIR__))
        if file == "runtests.jl"
            continue
        end
        include(file)
    end
end


text:
struct BanditSim
    b::Bandit
    G::ObjectiveFunc
    metadata::Dict
end
function BanditSim(b::Bandit, G::ObjectiveFunc)
    BanditSim(b, G, Dict(:sim_algorithm=>string(b), :sim_objective=>string(G)))
end
POMDPs.simulate(sim::BanditSim) = POMDPs.solve(sim.b, sim.G)