text:
struct HermitianPlot{A,B}
    x::A
    y::B
    function HermitianPlot(x::AbstractVector{T},y::AbstractVector{R}) where {R<:AbstractMatrix{S}} where {S<:Number,T<:Number}
        new{typeof(x),typeof(y)}(x,y)
    end
end
@recipe function f(h::HermitianPlot)
    x,y = h.x,h.y
    y_new = Vector{Matrix{Float64}}(undef, length(y))
    for i = 1:length(y)
        y_temp = zeros(size(y[i]))
        for j = 1:size(y[i],1), k = 1:size(y[i],2)
            y_temp[j,k] = j >= k ? real(y[i][j,k]) : imag(y[i][j,k])
        end
        y_new[i] = y_temp
    end
    MatrixPlot(x,y_new)
end

struct MatrixPlot{A,B}
    x::A
    y::B
    function MatrixPlot(x::AbstractVector{T},y::AbstractVector{R}) where {R<:AbstractMatrix{S}} where {S<:Number,T<:Number}
        all(size(yi)==size(y[1]) for yi in y) || error("Matricies given to MatrixPlot are not all the same size.")
        new{typeof(x),typeof(y)}(x,y)
    end
end
@recipe function f(m::MatrixPlot)
    x,y = m.x,m.y
    layout --> size(y[1])
    link --> :all
    label --> false
    count = 1
    for i = 1:size(y[1],1), j = 1:size(y[1],2)
        @series begin
            subplot --> count
            count += 1
            seriestype := :line
            x, getindex.(y,i,j)
        end
    end
end

@userplot PlotSdf
@recipe function f(p::PlotSdf)
    if length(p.args) == 1
        p.args[1] isa TimeSeriesModel || throw(ArgumentError("If only one argument is provided to plotsdf, it must be a TimeSeriesModel."))
        Ω = range(-π,π,length=200)
    elseif length(p.args) == 2
        p.args[1] isa TimeSeriesModel || throw(ArgumentError("First argument to plotsdf must be a TimeSeriesModel."))
        p.args[2] isa AbstractVector{T} where {T<:Number} || throw(ArgumentError("Second argument to plotasdf should be the frequencies, a vector of numbers."))
        Ω = p.args[2]
    else
        throw(ArgumentError("A maximum of two numbers should be passed to plotsdf."))
    end
    model = p.args[1]
    y = sdf.(model, Ω)
    if ndims(model) > 1
        return HermitianPlot(Ω, y)
    else
        return Ω, real.(y)
    end
end

"""
    plotsdf(model)
    plotsdf!(model)
    plotsdf(model, Ω)
    plotsdf!(model, Ω)

Plot the spectral density function of a model at the frequencies `Ω`.
By default, `Ω = range(-π,π,length=200)` and `Δ = 1`.
"""
plotsdf

"""
    checkposreal(x)

Checks if x is a positive real number."""
checkposreal(x) = x isa Real && x > 0

@userplot PlotAsdf
@recipe function f(p::PlotAsdf)
    if length(p.args) == 1
        p.args[1] isa TimeSeriesModel || throw(ArgumentError("If only one argument is provided to plotasdf, it must be a TimeSeriesModel."))
        Ω = range(-π,π,length=200)
        Δ = 1
    elseif length(p.args) == 3
        p.args[1] isa TimeSeriesModel || throw(ArgumentError("First argument to plotasdf must be a TimeSeriesModel."))
        p.args[2] isa AbstractVector{T} where {T<:Number} || throw(ArgumentError("Second argument to plotasdf should be the frequencies, a vector of numbers."))
        checkposreal(p.args[3]) || throw(ArgumentError("Third argument to plotasdf should be the sampling rate, a positive real number."))
        Ω = p.args[2]
        Δ = p.args[3]
    else
        throw(ArgumentError("Either only a model should be provided to plotasdf or a model, vector of frequencies and sampling rate Δ."))
    end
    model = p.args[1]
    y = asdf.(model, Ω, Δ)
    if ndims(model) > 1
        return HermitianPlot(Ω, y)
    else
        return Ω, real.(y)
    end
end

"""
    plotasdf(model)
    plotasdf!(model)
    plotasdf(model, Ω, Δ)
    plotasdf!(model, Ω, Δ)

Plot the aliased spectral density function of a model at the frequencies Ω.
By default, `Ω = range(-π,π,length=200)` and `Δ = 1`.
"""
plotasdf

@userplot PlotAcv
@recipe function f(p::PlotAcv)
    if length(p.args) == 1
        p.args[1] isa TimeSeriesModel || throw(ArgumentError("If only one argument is provided to plotacv, it must be a TimeSeriesModel."))
        n = 30
        Δ = 1
        τ = (-n+1:n-1).*Δ
    elseif length(p.args) == 2
        p.args[1] isa TimeSeriesModel || throw(ArgumentError("First argument to plotsdf must be a TimeSeriesModel."))
        p.args[2] isa AbstractVector{T} where {T<:Number} || throw(ArgumentError("If two arguments provided to the second must be a vector of time lags."))
        !(p.args[1] isa UnknownAcvTimeSeriesModel) || throw(ArgumentError("Custom lag vector only possible if acv of model is known, pass the number of observations and sampling rate instead."))
        τ = p.args[2]
    elseif length(p.args) == 3
        p.args[1] isa TimeSeriesModel || throw(ArgumentError("First argument to plotsdf must be a TimeSeriesModel."))
        checkposreal(p.args[2]) || throw(ArgumentError("Second argument to plotasdf should be the number of lags, a positive real number."))
        checkposreal(p.args[3]) || throw(ArgumentError("Third argument to plotasdf should be the sampling rate, a positive real number."))
        n = p.args[2]
        Δ = p.args[3]
        τ = (-n+1:n-1).*Δ
    else
        throw(ArgumentError("Either only a model should be provided to plotacv or a model, number of observations and sampling rate Δ (or vector of lags if appropriate)."))
    end
    model = p.args[1]
    y = length(p.args)==2 ? acv(model, τ) : acv(model, n, Δ)
    if ndims(model) > 1
        return MatrixPlot(τ, y)
    else
        return τ, real.(y)
    end
end

"""
    plotacv(model)
    plotacv!(model)
    plotacv(model, τ)
    plotacv!(model, τ)
    plotacv(model, n, Δ)
    plotacv!(model, n, Δ)

Plot the aliased spectral density function of a model at lags `τ`.
If the model does not have known acv, then the number of lags `n` and sampling rate `Δ` should be provided.
In this case, the lags are `-(n-1)*Δ:Δ:(n-1)*Δ`.
If unprovided, `n = 30` and `Δ = 1`.
"""
plotacv

@userplot PlotEI
@recipe function f(p::PlotEI)
    if length(p.args) == 1
        p.args[1] isa TimeSeriesModel || throw(ArgumentError("If only one argument is provided to plotei, it must be a TimeSeriesModel."))
        n = 1024
        Δ = 1
    elseif length(p.args) == 3
        p.args[1] isa TimeSeriesModel || throw(ArgumentError("First argument to plotei must be a TimeSeriesModel."))
        checkposreal(p.args[2]) || throw(ArgumentError("Second argument to plotei should be the number of observations, a positive real number."))
        checkposreal(p.args[3]) || throw(ArgumentError("Third argument to plotei should be the sampling rate, a positive real number."))
        n = p.args[2]
        Δ = p.args[3]
    else
        throw(ArgumentError("Either only a model should be provided to plotei or a model, number of observations and sampling rate Δ."))
    end
    Ω = fftshift(fftfreq(n, 2π/Δ))
    model = p.args[1]
    y = EI(model, n, Δ)
    if ndims(model) > 1
        return HermitianPlot(Ω, y)
    else
        return Ω, real.(y)
    end
end

"""
    plotei(model)
    plotei!(model)
    plotei(model, n, Δ)
    plotei!(model, n, Δ)

Plot the aliased spectral density function of a model at the angular Fourier frequencies `fftfreq(n,2π/Δ)`.
By default, `n = 1024` and `Δ = 1`.
"""
plotei

@userplot PlotCoh
@recipe function f(p::PlotCoh)
    if length(p.args) == 1
        p.args[1] isa TimeSeriesModel || throw(ArgumentError("If only one argument is provided to plotsdf, it must be a TimeSeriesModel."))
        Ω = range(-π,π,length=200)
    elseif length(p.args) == 2
        p.args[1] isa TimeSeriesModel || throw(ArgumentError("First argument to plotsdf must be a TimeSeriesModel."))
        p.args[2] isa AbstractVector{T} where {T<:Number} || throw(ArgumentError("Second argument to plotasdf should be the frequencies, a vector of numbers."))
        Ω = p.args[2]
    else
        throw(ArgumentError("A maximum of two numbers should be passed to plotsdf."))
    end
    model = p.args[1]
    if ndims(model) > 1
        S = sdf.(model, Ω)
        C = coherance.(model, Ω)
        G = groupdelay.(model, Ω)
        y = [
            [i==j ? real(s[i,j]) : (i<j ? c[i,j] : g[i,j]) for i in 1:ndims(model), j in 1:ndims(model)]
        for (s,c,g) in zip(S,C,G)]
        layout --> size(y[1])
        link --> :none
        label --> false
        count = 1
        for i = 1:size(y[1],1), j = 1:size(y[1],2)
            @series begin
                subplot --> count
                count += 1
                seriestype := :line
                if i<j
                    yscale := :identity
                    ylims --> (0,1)
                elseif i>j
                    yscale := :identity
                    ylims --> (-π,π)
                end
                Ω, getindex.(y,i,j)
            end
        end
    else
        return Ω, real.(sdf.(model, Ω))
    end
end

"""
    plotcoh(model)
    plotcoh!(model)
    plotcoh(model, Ω)
    plotcoh!(model, Ω)

Plot the coherance and group delay of a model at the frequencies `Ω`.
If multivariate, above the diagonal is the coherance, and below is the group delay
and on the diagonal is the spectral density function.
By default, `Ω = range(-π,π,length=200)` and `Δ = 1`.
"""
plotcoh


text:
function coefs(R, xroot, lroot, ϵroot, κ, norm)

    @info "Calculating coefficients..."

    energ = @. xroot*ħ/R
    j = @. lroot+ϵroot/2

    norm² = norm * norm'
    uu̅(r) = @. sphericalbessel(xroot/R*r,lroot)^2 + sphericalbessel(xroot/R*r,lroot+ϵroot)^2
    uu̅uu̅(r) = uu̅(r) * uu̅(r)' * r^2
    ∫uu̅uu̅ = quadgk.(uu̅uu̅, 0, R)[1]

    coef = @. norm²*∫uu̅uu̅/(4*π)/2
    avgcoef = sum(coef)/κ^2

    @info "Done!"
    @debug "Corresponding j: " j
    @debug "Corresponding E: " energ
    @debug "Average value of coefficients: " avgcoef

    return coef, avgcoef, energ, j

end

ext: jl

max_stars_repo_name: avramescudana/QuarkDroplets

max_issues_repo_name: avramescudana/QuarkDroplets

max_forks_repo_name: avramescudana/QuarkDroplets


text:
# meant to be a faster version of column reader
path = "C:/git/parquet-data-collection/dsd50p.parquet"

filemetadata = metadata(path)

io = open(path)

seek(io, 4)

read_thrift(io, PAR2.PageHeader)

using Snappy

compressed_page = read(io, 2461)

uncompressed_page = Snappy.uncompress(compressed_page)

dict = reinterpret(Float64, uncompressed_page)

next_page = read_thrift(io, PAR2.PageHeader)

dict

a = read_parquet("C:/git/parquet-data-collection/dsd50p.parquet")

compressed_page = read(io, 36328)

uncompressed_page = Snappy.uncompress(compressed_page)

Int.(uncompressed_page)

dict[140]

Int(uncompressed_page[5])

# every number here represents a encoding
bit_width = uncompressed_page[1]

iobuffer = IOBuffer(uncompressed_page)

bit_width = Int(read(iobuffer, UInt8))

data_length = Int(ntoh(read(iobuffer, UInt32)))

using LittleEndianBase128

LittleEndianBase128.decode(uncompressed_page[6:9])

using LittleEndianBase128: encode

encode(UInt(1))

uncompressed_page[6:9]

Int(28998//2) # bytes are using to read from dictionary

indices = (uncompressed_page[2:Int(28998//2)+1] .<< 4) .>> 4

indices = uncompressed_page[2:Int(28998//2)+1] .& UInt8(15)

indices2 = uncompressed_page[2:Int(28998//2)+1] .>> 4

dict[indices .+ 1]

dict[indices2 .+ 1]

dict[(uncompressed_page[2:Int(28998//2)+1] .>> 4) .+ 1]


text:
"""
    Board

A `DataType` representing the state of a chess board.
"""
mutable struct Board
    squares::MVector{64,Piece}
    pieces::MVector{6,Bitboard}
    colors::MVector{2,Bitboard}
    checkers::Bitboard
    pinned::Bitboard
    turn::Color
    castling::UInt8
    enpass::UInt8
    halfmovecount::UInt16
    movecount::UInt16
    hash::ZobristHash
    history::Vector{ZobristHash}
    psqteval::Int32
    pkhash::ZobristHash
end

function Board(squares::AbstractArray{Piece}, pieces::AbstractArray{Bitboard}, colors::AbstractArray{Bitboard},
    checkers::Bitboard, pinned::Bitboard, turn::Color, castling::UInt8, enpass::UInt8, halfmovecount::UInt16, movecount::UInt16,
    hash::ZobristHash, history::Vector{ZobristHash}, psqteval::Int32, pkhash::ZobristHash)

    return Board(MVector(squares...), MVector(pieces...), MVector(colors...),
        checkers, pinned, turn, castling, enpass, halfmovecount, movecount,
        hash, history, psqteval, pkhash)
end
Board() = Board(repeat([BLANK], 64), repeat([EMPTY], 6), repeat([EMPTY], 2),
    EMPTY, EMPTY, WHITE, zero(UInt8), zero(UInt8), zero(UInt16), zero(UInt16),
    ZobristHash(zero(UInt64)), ZobristHash.(zeros(UInt64, 512)), zero(Int32), ZobristHash(zero(UInt64)))

"""
    copy!(board_1::Board, board_2::Board)

Copy the contents of `board_2` to `board_1`.
"""
function copy!(board_1::Board, board_2::Board)
    copyto!(board_1.squares, board_2.squares)
    copyto!(board_1.pieces, board_2.pieces)
    copyto!(board_1.colors, board_2.colors)
    board_1.checkers = board_2.checkers
    board_1.pinned = board_2.pinned
    board_1.turn = board_2.turn
    board_1.castling = board_2.castling
    board_1.enpass = board_2.enpass
    board_1.halfmovecount = board_2.halfmovecount
    board_1.movecount = board_2.movecount
    board_1.hash = board_2.hash
    board_1.psqteval = board_2.psqteval
    board_1.pkhash = board_2.pkhash
    copy!(board_1.history, board_2.history)
end

getindex(board::Board, color::Color) = @inbounds board.colors[color.val]
getindex(board::Board, type::PieceType) = @inbounds board.pieces[type.val]
getindex(board::Board, piece::Piece) = @inbounds board[type(piece)] & board[color(piece)]
getindex(board::Board, idx::Integer) = @inbounds board.squares[idx]

setindex!(board::Board, bb::Bitboard, type::PieceType) = setindex!(board.pieces, bb, type.val)
setindex!(board::Board, bb::Bitboard, color::Color) = setindex!(board.colors, bb, color.val)
setindex!(board::Board, piece::Piece, idx::Integer) = setindex!(board.squares, piece, idx)

"""
    add!(board::Board, piece::Piece, bb::Bitboard)
    add!(board::Board, piece::Piece, sqr::Integer)

Add a `piece` to the `board` square given by an `Integer` `sqr`, or a `Bitboard` `bb`.
"""
function add!(board::Board, piece::Piece, bb::Bitboard, sqr::Integer)
    @inbounds board[type(piece)] |= bb
    @inbounds board[color(piece)] |= bb
    @inbounds board[sqr] = piece
    board.hash ⊻= zobkey(piece, sqr)
    board.psqteval += psqt(piece, sqr)
    if type(piece) == PAWN || type(piece) == KING
        board.pkhash ⊻= zobkey(piece, sqr)
    end
end
add!(board::Board, piece::Piece, bb::Bitboard) = add!(board, piece, bb, square(bb))
add!(board::Board, piece::Piece, sqr::Integer) = add!(board, piece, Bitboard(sqr), sqr)

"""
    remove!(board::Board, bb::Bitboard)
    remove!(board::Board, sqr::Integer)

Remove a `piece` from the `board` at the square given by an `Integer` `sqr`, or a `Bitboard` `bb`.
"""
function remove!(board::Board, bb::Bitboard, sqr::Integer)
    piece = board.squares[sqr]
    board[type(piece)] &= ~bb
    board[color(piece)] &= ~bb
    board.squares[sqr] = BLANK
    board.hash ⊻= zobkey(piece, sqr)
    board.psqteval -= psqt(piece, sqr)
    if type(piece) == PAWN
        board.pkhash ⊻= zobkey(piece, sqr)
    end
end
remove!(board::Board, bb::Bitboard) = remove!(board, bb, square(bb))
remove!(board::Board, sqr::Integer) = remove!(board, Bitboard(sqr), sqr)

"""
    addremove!(board::Board, piece::Piece, bb::Bitboard)
    addremove!(board::Board, piece::Piece, sqr::Integer)

Remove a `piece` from the `board` at the square given by an `Integer` `sqr`, or a `Bitboard` `bb`, and then add the given `piece` in its place.
"""
function addremove!(board::Board, piece::Piece, bb::Bitboard, sqr::Integer)
    captured = board.squares[sqr]
    board[type(captured)] &= ~bb
    board[color(captured)] &= ~bb
    board[type(piece)] |= bb
    board[color(piece)] |= bb
    board[sqr] = piece
    board.hash ⊻= zobkey(captured, sqr)
    board.hash ⊻= zobkey(piece, sqr)
    board.psqteval += psqt(piece, sqr)
    board.psqteval -= psqt(captured, sqr)
    if type(piece) == PAWN
        board.pkhash ⊻= zobkey(piece, sqr)
    end
    if type(captured) == PAWN
        board.pkhash ⊻= zobkey(captured, sqr)
    end
end
addremove!(board::Board, piece::Piece, bb::Bitboard) = addremove!(board, piece, bb, square(bb))
addremove!(board::Board, piece::Piece, sqr::Integer) = addremove!(board, piece, Bitboard(sqr), sqr)

"""
    ischeck(board::Board)

Returns `true` if the king is in check, `false` otherwise.
"""
ischeck(board::Board) = isempty(board.checkers) === false

"""
    isdoublecheck(board::Board)

Returns `true` if the king is in check, `false` otherwise.
"""
isdoublecheck(board::Board) = ismany(board.checkers)

"""
    switchturn(board::Board)

Switches the `Color` of the side to move on a given `Board` object.
"""
switchturn!(board::Board) = board.turn = !board.turn

"""
    enemies(board::Board)

Return the positions, as a `Bitboard`, of all the enemies on the board.
"""
enemy(board::Board) = @inbounds board[!board.turn]

"""
    friendly(board::Board)

Return the positions, as a `Bitboard`, of all the friendly units on the board.
"""
friendly(board::Board) = @inbounds board[board.turn]

"""
    occupied(board::Board)

Return the positions, as a `Bitboard`, of all the occupied squares on the board.
"""
occupied(board::Board) = white(board) | black(board)

"""
    empty(board::Board)

Return the positions, as a `Bitboard`, of all the empty squares on the board.
"""
empty(board::Board) = ~occupied(board)

"""
    SYMBOLS

A constant holding the glyphs for each chess piece.
"""
const SYMBOLS = ['♙','♘','♗','♖','♕','♔'] .+ 6

# custom show for board type
function Base.show(io::IO, board::Board)
    for row in 1:8
        for col in 1:8
            if col == 1
                print(9 - row, " ")
            end
            sqr = square(FILE[col] & RANK[9 - row])
            pcolor = color(board[sqr])
            ptype = type(board[sqr])
            sym = (ptype.val > zero(UInt8)) ? SYMBOLS[ptype.val] : ' '
            foreground = (pcolor == WHITE) ? :white : :black
            background = isodd(row + col) ? :blue : :light_blue
            print(Crayon(foreground = foreground, background = background), sym, " ")
            if (row == 1) && (col == 8) && (board.turn == WHITE)
                print(Crayon(reset = true), " White to move...")
            elseif (row == 1) && (col == 8) && (board.turn == BLACK)
                print(Crayon(reset = true), " Black to move...")
            end
        end
        print(Crayon(reset = true), "\n")
    end
    println("  A B C D E F G H")
end

"""
    white(board::Board)

Get the location of all the white pieces on the board, as a `Bitboard`.
"""
white(board::Board) = @inbounds board.colors[1]

"""
    black(board::Board)

Get the location of all the black pieces on the board, as a `Bitboard`.
"""
black(board::Board) = @inbounds board.colors[2]

"""
    pawns(board::Board)

Get the location of all the pawns on the `board`, as a `Bitboard`.
"""
pawns(board::Board) = @inbounds board.pieces[1]

"""
    knights(board::Board)

Get the location of all the knights on the `board`, as a `Bitboard`.
"""
knights(board::Board) = @inbounds board.pieces[2]

"""
    bishops(board::Board)

Get the location of all the bishops on the `board`, as a `Bitboard`.
"""
bishops(board::Board) = @inbounds board.pieces[3]

"""
    rooks(board::Board)

Get the location of all the rooks on the `board`, as a `Bitboard`.
"""
rooks(board::Board) = @inbounds board.pieces[4]

"""
    queens(board::Board)

Get the location of all the queens on the `board`, as a `Bitboard`.
"""
queens(board::Board) = @inbounds board.pieces[5]

"""
    kings(board::Board)

Get the location of all the kings on the `board`, as a `Bitboard`.
"""
kings(board::Board) = @inbounds board.pieces[6]

"""
    rooklike(board::Board)

Get the location of all the rooks and queens on the `board`, as a `Bitboard`.
"""
rooklike(board::Board) = rooks(board) | queens(board)

"""
    bishoplike(board::Board)

Get the location of all the bishops and queens on the `board`, as a `Bitboard`.
"""
bishoplike(board::Board) = bishops(board) | queens(board)

"""
    piece(board::Board, sqr::Integer)

Get the `Piece` located at a square, `sqr`.
"""
piece(board::Board, sqr::Integer) = @inbounds board[sqr]

"""
    checkers(board::Board)

Return the `Bitboard` of all the pieces giving check.
"""
checkers(board::Board) = board.checkers

"""
    pinned(board::Board)

Return the `Bitboard` of all the pieces that are pinned.
"""
pinned(board::Board) = board.pinned

# Castling as follows
# white kingside UInt8(1) << 0 (color val - 1)
# black kingside UInt8(1) << 1 (color val - 1)
# white queenside UInt8(1) << 2 (color val + 1)
# black queenside UInt8(1) << 3(color val + 1)
"""
    cancastlekingside(board::Board, color::Color)
    cancastlekingside(board::Board)

Return a `Bool` which denotes if the player of `color` can castle kingside.
"""
cancastlekingside(board::Board, color::Color) = isone((board.castling >> (color.val - 1)) & 1)
cancastlekingside(board::Board) = cancastlekingside(board, board.turn)

"""
    cancastlequeenside(board::Board, color::Color)
    cancastlequeenside(board::Board)

Return a `Bool` which denotes if the player of `color` can castle queenside. If no color is given, it assumes the `color` of the current turn.
"""
cancastlequeenside(board::Board, color::Color) = isone((board.castling >> (color.val + 1)) & 1)
cancastlequeenside(board::Board) = cancastlequeenside(board, board.turn)

"""
    has_non_pawn_material(board::Board)

Check if a position has non-pawn material left.
"""
has_non_pawn_material(board::Board) = !isempty(knights(board)) || !isempty(bishops(board)) || !isempty(rooks(board)) || !isempty(queens(board))

"""
    islegal(board::Board)

A safety check to see if the position on the board is legal. Returns `true` if legal.
"""
function islegal(board::Board)
    switchturn!(board)
    bool = isempty(kingAttackers(board))
    switchturn!(board)
    return bool
end

"""
    isdraw(board::Board)

Detects checkmate. Returns `true` if checkmated.
"""
function ischeckmate(board::Board)
    if ischeck(board)
        ml = MoveStack(50)
        gen_moves!(ml, board)
        if length(ml) == 0
            return true
        else
            return false
        end
    else
        return false
    end
end

"""
    isdraw(board::Board)

Detects draws by stalemate, insufficient material, 50 move rule, and 3-fold repetition. Returns `true` if drawn.
"""
function isdraw(board::Board)
    isstalemate(board) || isdrawbymaterial(board) || is50moverule(board) || isrepetition(board)
end

"""
    isstalemate(board::Board)

Detects stalemate draws. Returns `true` if drawn.
"""
function isstalemate(board::Board)
    if ischeck(board)
        return false
    else
        ml = MoveStack(100)
        gen_quiet_moves!(ml, board)
        if !(length(ml) == 0)
            return false
        end
        clear!(ml)
        gen_noisy_moves!(ml, board)
        if !(length(ml) == 0)
            return false
        else
            return true
        end
    end
end

"""
    isdrawbymaterial(board::Board)

Detects draws by insufficient material. Returns `true` if the position is drawn.
"""
function isdrawbymaterial(board::Board)
    piece_count = count(white(board)) + count(black(board))
    if piece_count > 4
        return false
    elseif piece_count == 2
        return true
    elseif piece_count == 3
        if isone(bishops(board))
            return true
        elseif isone(knights(board))
            return true
        end
    elseif piece_count == 4
        if ismany(bishops(board) & LIGHT)
            return true
        elseif ismany(bishops(board) & DARK)
            return true
        end
    end
    return false
end

"""
    is50moverule(board::Board)

Detects draws by the 50 move rule. Returns `true` if the position is drawn.
"""
function is50moverule(board::Board)
    if board.halfmovecount > 99
        return true
    else
        return false
    end
end

"""
    isrepetition(board::Board)

Detects 3-fold repetition by analysing the hash history. Returns `true` if the position is drawn.
"""
function isrepetition(board::Board)
    if board.halfmovecount < 8
        return false
    end
    reps = 0
    for i in (board.movecount):-2:(board.movecount - board.halfmovecount)
        if (board.hash == board.history[i])
            reps += 1
        end
        if reps == 3
            return true
        end
    end
    return false
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
using MyStats
using Base.Test


text:
"""
$(SIGNATURES)

Scale each "row" (along dimension `d`) of array `x` to match the totals given in `totalV`. The result satisfies

`sum(selectdim(x, d, j)) == totalV[j]`
"""
function scale_array!(x :: AbstractArray{F1}, d :: Integer, totalV :: AbstractVector{F1}) where F1 <: AbstractFloat

    J = length(totalV);
    @assert J == size(x, d)  "Size mismatch"
    for j = 1 : J
        xd = selectdim(x, d, j);
        xd .*= (totalV[j] / sum(xd));
    end
    return nothing
end

"""
$(SIGNATURES)

Bracket an array between lower and upper bounds. In place.
"""
function bracket_array!(x :: F1, lb :: F1, ub :: F1) where F1 <: AbstractArray
    @assert size(x) == size(lb) == size(ub)
    for (k, v) in enumerate(x)
        if v > ub[k]
            x[k] = ub[k];
        elseif v < lb[k]
            x[k] = lb[k];
        end
    end
    return nothing
end

function bracket_array!(xM :: AbstractArray{F1}, lb :: F1, ub :: F1) where F1
    for (k, x) in enumerate(xM)
        if x > ub
            xM[k] = ub;
        end
        if x < lb
            xM[k] = lb;
        end
    end
end

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


text:
module MetaModuleTests
using  Base.Test

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

module A
using  ExpressionPatterns.Dispatch

@metamodule export @f,@g

@macromethod f(x) 1
@macromethod g(x) 1

end

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

module B
using  ExpressionPatterns.Dispatch

@metamodule import ..A.@f

@macromethod f(x+y) 2
@macromethod g(x+y) 2

end

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

module C
using  ExpressionPatterns.Dispatch

@metamodule importall ..A

@macromethod f(x-y) 3
@macromethod g(x-y) 3

end

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

@test A.@f(x+y) == B.@f(x+y)
@test A.@g(x+y) != B.@g(x+y)

@test A.@f(x-y) == C.@f(x-y)
@test A.@f(x-y) == C.@f(x-y)

end


text:
module TSPImages

using LightGraphs
using Images
using Random
using TravelingSalesmanHeuristics: solve_tsp
using Gadfly
using LinearAlgebra: norm

#=
Note: throughout this package I index pixels in images according to the indexing of normal
Julia matrices. That is, index (1,1) is the top left corner, (1,2) is directly to the right
of the top left corner, and (2,1) is directly below the top left corner. So pixel indices
are of the form (y,x) where increasing y corresponds to moving down. Note that this
indexing convention doesn't match that of many image editing programs.
=#

"Load an image to grayscale. If `invert`, invert the grayscale image. Return a tuple
(image_grayscale, image)."
function load_gray_image(path; invert = false)
    img = load(path)
    img_gray = Gray.(img)
    if invert
        img_gray .= oneunit(eltype(img_gray)) .- img_gray
    end
    return img_gray, img
end

"Sample `n` points according to a grayscale image. If the image has height and width
`(h,w)`, points returned will be in (0,h] x (0,w]. The probability of sampling a given
point is proportional to the darkness of the corresponding pixel in the input image.

Return a 2 by n matrix with one column per sampled point"
function sample_from_image(img::Matrix, n::Int)
    h, w = size(img)
    points = zeros(2,n)
    for pt_ind in 1:n
        # we sample with (1 - rand()) to get values in (0, 1]
        while true
            pt_y = h * (1 - rand())
            pt_x = w * (1 - rand())
            pt_pixel_val = img[ceil(Int, pt_y), ceil(Int, pt_x)]
            if rand() >= pt_pixel_val # accept this point, darker => more likely to accept
                points[:, pt_ind] .= [pt_y, pt_x]
                break
            end
        end
    end
    return points
end

"Partition a matrix of points into components. Two points in separate components are at
least `dist_cutoff` pixels from each other.

Internally, this function constructs a simple graph over the passed points. This graph
has an edge between any two points with distance less than `dist_cutoff`, and we then
return the connected components of this graph.

Return a Vector{Vector{Int}} where each `i` in `1:size(pts, 2)` appears exactly once."
function get_components(pts, dist_cutoff)
    graph = LightGraphs.euclidean_graph(pts; cutoff = dist_cutoff)[1]
    return LightGraphs.connected_components(graph)
end

"Order a matrix of points according to an approximate TSP solution. This will increase
the number of columns in the `pts` matrix by one as the path start/end will be duplicated.

Optional kwargs are `quality_factor` which is passed to `solve_tsp` and `dist_cutoff`
which specifies the maximum allowable edge distance. A value of nothing corresponds to
no maximum."
function order_points(pts; quality_factor = 50, dist_cutoff = nothing)
    n = size(pts, 2)
    distmat = [norm(pts[:,i] - pts[:,j]) for i in 1:n, j in 1:n]
    if !(isnothing(dist_cutoff))
        distmat[distmat .> dist_cutoff] .*= 100
    end
    path, _ = solve_tsp(distmat; quality_factor = quality_factor)
    return pts[:, path]
end

"Create a Gadfly plot layer from ordered points. `pts` should be a 2 by n matrix where
columns correspond to points and the order of columns defines a TSP path."
function points_to_layer(pts, colors; styling = (), frac = 1.0)
    n = size(pts, 2)
    nplot = 1 + round(Int, frac * (n-1))
    return layer(
        x = pts[2, 1:(nplot-1)],
        xend = pts[2, 2:nplot],
        y = -pts[1, 1:(nplot-1)],
        yend = -pts[1, 2:nplot],
        color = colors,
        Geom.segment,
        style(;styling...)
    )
end

"Take a collection of layers and build a plot which combines the layers."
function build_plot(layers, height, width)
    p = plot(layers[1], Coord.cartesian(xmin = 0, xmax = width, ymin = -height, ymax = 0))
    for i in 2:length(layers)
        push!(p, layers[i])
    end
    return p
end

"Take a plot, set the background to white, and remove labels/grids."
function _whiten_plot!(p)
    push!(p, style(background_color = colorant"white"))
    push!(p, Guide.xlabel(nothing))
    push!(p, Guide.ylabel(nothing))
    push!(p, Guide.xticks(;ticks = nothing))
    push!(p, Guide.yticks(;ticks = nothing))
    return nothing
end

"Load an image, sample points, solve TSP subproblems, and plot!

Positional arguments:
- `path`: path to image file
- `n::Int`: number of points to sample from the image
- `component_dist::Real`: distance in pixels between image components. If a group of sapmled
  points is separated from all other sampled points by at least `component_dist` pixels,
  then that group of points will form its own TSP path.

Keyword arguments:
- `layer_styling`
- `quality_factor::Real`: quality factor keyword passed to `solve_tsp`. See
  `TravelingSalesmanHeuristics.solve_tsp` for more information. Defaults to 50.
- `invert_image::Bool`: Whether to invert the grayscale input image so that light areas are
  sampled more than dark areas.
- `path_dist_cutoff::Union{Real, Nothing}`: maximum distance in pixels of an edge in any
  TSP path, or `nothing` for no limit.
- `verbose::Bool`: whether to print timing information during execution.

Return a NamedTuple with keys :plot, :height, and :width corresponding to
the plot object and dimensions of the input image."
function image_to_tsp_plot(
    path, n, component_dist;
    layer_styling = (:default_color => colorant"black",),
    quality_factor = 50,
    invert_image = false,
    path_dist_cutoff = nothing,
    single_color = nothing,
    verbose = true
)
    _time_load = @elapsed img_gray, img = load_gray_image(path; invert = invert_image)
    verbose && @show _time_load
    height, width = size(img)
    _time_sample = @elapsed pts = sample_from_image(img_gray, n)
    verbose && @show _time_sample
    _time_components = @elapsed components = get_components(pts, component_dist)
    verbose && @show _time_components
    _time_tsp = @elapsed layers = map(components) do comp
        pts_ordered = order_points(pts[:, comp];
                                   quality_factor = quality_factor,
                                   dist_cutoff = path_dist_cutoff)
        segment_midpoints = (pts_ordered[:,1:(end-1)] .+ pts_ordered[:,2:end]) / 2
        if isnothing(single_color)
            colors = [img[ceil(Int, col[1]), ceil(Int, col[2])] for col in eachcol(segment_midpoints)]
        else
            colors = [single_color]
        end
        return points_to_layer(pts_ordered, colors; styling = layer_styling)
    end
    verbose && @show _time_tsp
    p = build_plot(layers, height, width)
    _whiten_plot!(p)
    return (plot = p, width = width, height = height)
end

function save_plot(outfile, plot, width, height)
    draw(SVG(outfile, 1px * width, 1px * height), plot)
    return nothing
end

end # module


text:
const FrontMatter = Dict{String,Any}

"""
	frontmatter(nb::Notebook; raise::Bool=false)::Dict{String,Any}
	frontmatter(nb_path::String; raise::Bool=false)::Dict{String,Any}

Extract frontmatter from a notebook, which is extra meta-information that the author attaches to the notebook, often including *title*, *description*, *tags*, *author*, and more. Search for *frontmatter* online to learn more.

If `raise` is true, then parsing errors will be rethrown. If `false`, this function will always return a `Dict`.
"""
function frontmatter(nb::Notebook; raise::Bool=false)
	convert(FrontMatter,
		get(() -> FrontMatter(), get_metadata(nb), "frontmatter")
	)
end

function frontmatter(abs_path::String; raise::Bool=false)
    try
		# this will load the notebook to analyse, it won't run it
		frontmatter(load_notebook_nobackup(abs_path); raise)
	catch e
		if raise
			rethrow(e)
		else
			@error "Error reading notebook file." abs_path exception=(e,catch_backtrace())
			FrontMatter()
		end
	end
end


text:
function SciMLBase.solve(cache::LinearCache, alg::AbstractFactorization; kwargs...)
    if cache.isfresh
        fact = do_factorization(alg, cache.A, cache.b, cache.u)
        cache = set_cacheval(cache, fact)
    end

    y = ldiv!(cache.u, cache.cacheval, cache.b)
    SciMLBase.build_linear_solution(alg,y,nothing,cache)
end

# Bad fallback: will fail if `A` is just a stand-in
# This should instead just create the factorization type.
init_cacheval(alg::AbstractFactorization, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose) = do_factorization(alg, A, b, u)

## LU Factorizations

struct LUFactorization{P} <: AbstractFactorization
    pivot::P
end

function LUFactorization()
    pivot = @static if VERSION < v"1.7beta"
        Val(true)
    else
        RowMaximum()
    end
    LUFactorization(pivot)
end

function do_factorization(alg::LUFactorization, A, b, u)
    A isa Union{AbstractMatrix,AbstractDiffEqOperator} ||
        error("LU is not defined for $(typeof(A))")

    if A isa AbstractDiffEqOperator
        A = A.A
    end
    fact = lu!(A, alg.pivot)
    return fact
end

init_cacheval(alg::LUFactorization, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose) = ArrayInterface.lu_instance(A)

# This could be a GenericFactorization perhaps?
Base.@kwdef struct UMFPACKFactorization <: AbstractFactorization
    reuse_symbolic::Bool = true
end

function do_factorization(::UMFPACKFactorization, A, b, u)
    if A isa AbstractDiffEqOperator
        A = A.A
    end
    if A isa SparseMatrixCSC
        return lu(A)
    else
        error("Sparse LU is not defined for $(typeof(A))")
    end
end

function SciMLBase.solve(cache::LinearCache, alg::UMFPACKFactorization)
    A = cache.A
    if A isa AbstractDiffEqOperator
        A = A.A
    end
    if cache.isfresh
        if cache.cacheval !== nothing && alg.reuse_symbolic
            # If we have a cacheval already, run umfpack_symbolic to ensure the symbolic factorization exists
            # This won't recompute if it does.
            SuiteSparse.UMFPACK.umfpack_symbolic!(cache.cacheval)
            fact = lu!(cache.cacheval, A)
        else
            fact = do_factorization(alg, A, cache.b, cache.u)
        end
        cache = set_cacheval(cache, fact)
    end

    y = ldiv!(cache.u, cache.cacheval, cache.b)
    SciMLBase.build_linear_solution(alg,y,nothing,cache)
end

Base.@kwdef struct KLUFactorization <: AbstractFactorization
    reuse_symbolic::Bool = true
end

function do_factorization(::KLUFactorization, A, b, u)
    if A isa AbstractDiffEqOperator
        A = A.A
    end
    if A isa SparseMatrixCSC
        return klu(A)
    else
        error("KLU is not defined for $(typeof(A))")
    end
end

function SciMLBase.solve(cache::LinearCache, alg::KLUFactorization)
    A = cache.A
    if A isa AbstractDiffEqOperator
        A = A.A
    end
    if cache.isfresh
        if cache.cacheval !== nothing && alg.reuse_symbolic
            # If we have a cacheval already, run umfpack_symbolic to ensure the symbolic factorization exists
            # This won't recompute if it does.
            KLU.klu_analyze!(cache.cacheval)
            fact = klu!(cache.cacheval, A)
        else
            fact = do_factorization(alg, A, cache.b, cache.u)
        end
        cache = set_cacheval(cache, fact)
    end

    y = ldiv!(cache.u, cache.cacheval, cache.b)
    SciMLBase.build_linear_solution(alg,y,nothing,cache)
end

## QRFactorization

struct QRFactorization{P} <: AbstractFactorization
    pivot::P
    blocksize::Int
end

function QRFactorization()
    pivot = @static if VERSION < v"1.7beta"
        Val(false)
    else
        NoPivot()
    end
    QRFactorization(pivot, 16)
end

function do_factorization(alg::QRFactorization, A, b, u)
    A isa Union{AbstractMatrix,AbstractDiffEqOperator} ||
        error("QR is not defined for $(typeof(A))")

    if A isa AbstractDiffEqOperator
        A = A.A
    end
    fact = qr!(A, alg.pivot; blocksize = alg.blocksize)
    return fact
end

## SVDFactorization

struct SVDFactorization{A} <: AbstractFactorization
    full::Bool
    alg::A
end

SVDFactorization() = SVDFactorization(false, LinearAlgebra.DivideAndConquer())

function do_factorization(alg::SVDFactorization, A, b, u)
    A isa Union{AbstractMatrix,AbstractDiffEqOperator} ||
        error("SVD is not defined for $(typeof(A))")

    if A isa AbstractDiffEqOperator
        A = A.A
    end

    fact = svd!(A; full = alg.full, alg = alg.alg)
    return fact
end

## GenericFactorization

struct GenericFactorization{F} <: AbstractFactorization
    fact_alg::F
end

GenericFactorization(;fact_alg = LinearAlgebra.factorize) =
    GenericFactorization(fact_alg)

function do_factorization(alg::GenericFactorization, A, b, u)
    A isa Union{AbstractMatrix,AbstractDiffEqOperator} ||
        error("GenericFactorization is not defined for $(typeof(A))")

    if A isa AbstractDiffEqOperator
        A = A.A
    end
    fact = alg.fact_alg(A)
    return fact
end

## RFLUFactorization

RFLUFactorization() = GenericFactorization(;fact_alg=RecursiveFactorization.lu!)
init_cacheval(alg::GenericFactorization{typeof(RecursiveFactorization.lu!)}, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose) = ArrayInterface.lu_instance(A)


text:
module MOBC
using JuMP, vOptGeneric, PyPlot, NSGAII, Suppressor, NamedTuples, Nullables, Compat, ProgressMeter, OffsetArrays

export solve_stidsen, solve_parragh, parseGAP, parseGAP_orlib

include("Node.jl")
include("NonDomPoints.jl")
include("CoverCuts.jl")
include("ParserGAP.jl")
include("Stidsen.jl")
include("Parragh.jl")
include("Plots.jl")
include("RCNP.jl")

isonlybinary(m, obj) = all(m.colCat[map(x->getfield(x, :col), obj.aff.vars)] .== :Bin)
isbinary(x) = all(y-> y==1. || y==0., x)
dominates(::Type{Min}, x, y) = x[1] < y[1] && x[2] <= y[2] || x[1] <= y[1] && x[2] < y[2]
dominates(::Type{Max}, x, y) = x[1] > y[1] && x[2] >= y[2] || x[1] >= y[1] && x[2] > y[2]
weakly_dominates(::Type{Min}, x, y) = x[1] <= y[1] && x[2] <= y[2]
weakly_dominates(::Type{Max}, x, y) = x[1] >= y[1] && x[2] >= y[2]

evaluate(x, obj) = evaluate(x, obj.aff.constant, obj.aff.coeffs, map(x->getfield(x, :col), obj.aff.vars))
function evaluate(x, cst, coeffs, vars)::Float64
	res = cst
	for i = 1:length(coeffs)
	    @inbounds res += coeffs[i] * x[vars[i]]
	end
	res
end

# function integerbranch(n)
# 	newnode = copy(n, false)
# 	inds0 = find(x->x==0., n.x)
# 	inds1 = find(x->x==1., n.x)
# 	#NO-GOOD CSTR
# 	@constraint(newnode.m, sum(JuMP.Variable(newnode.m, j) for j in inds0) + sum(1. - JuMP.Variable(newnode.m, j) for j in inds1) >= 1.)
# 	return newnode
# end

end # module


text:
#
# Running this file generates the data needed to generate Figure 5 in
#
#   D. Neumann, T. Wiese, and W. Utschick, Learning the MMSE Channel Estimator,
#   IEEE Transactions on Signal Processing, 2018.
#

push!(LOAD_PATH,".")
using DataFrames
using CSV
import SCM3GPP; const scm = SCM3GPP
import CondNormalMMSE; const mmse = CondNormalMMSE
import CondNormalTF; const cntf = CondNormalTF

include("sim_helpers.jl")
include("more_estimators.jl") # OMP and Genie MMSE

verbose = true
#-------------------------------------
# Simulation parameters
#
write_file = true
filename   = "results/figure5.csv"
nBatches   = 100
nBatchSize = 100

#-------------------------------------
# Channel Model
#
snr        = 0 # [dB]
antennas   = [8,16,32,64,96]
AS         = 2.0 # standard deviation of Laplacian (angular spread)
nCoherence = 1
Channel    = scm.SCMMulti(pathAS=AS, nPaths=1)
# method that generates "nBatches" channel realizations
get_channel(nAntennas, nCoherence, nBatches) = scm.generate_channel(Channel, nAntennas, nCoherence=nCoherence, nBatches = nBatches)
# method that samples C_delta from delta prior
get_cov(nAntennas) = scm.toeplitzHe( scm.generate_channel(Channel, nAntennas, nCoherence=1)[2][:] )
# get circulant vector that generates all covariance matrices for arbitrary delta (here: use delta=0)
get_circ_cov_generator(nAntennas) = real(scm.best_circulant_approximation(scm.scm_channel([0.0],[1.0],nAntennas,AS=AS)[2]))

results = DataFrame()
for iAntenna in 1:length(antennas)
    nAntennas     = antennas[iAntenna]

    verbose && println("Simulating with ", nAntennas, " antennas")

    # Conditionally normal estimators
    cn_est = Dict{Symbol,Any}()
    cn_est[:FastMMSE]     = mmse.FastMMSE(snr, get_circ_cov_generator(nAntennas))
    cn_est[:CircMMSE]     = mmse.StructuredMMSE(snr, () -> get_cov(nAntennas), nSamples=16*nAntennas, transform = circ_trans)
    cn_est[:ToepMMSE]     = mmse.StructuredMMSE(snr, () -> get_cov(nAntennas), nSamples=16*nAntennas, transform = toep_trans)
    cn_est[:DiscreteMMSE] = mmse.DiscreteMMSE(snr,   () -> get_cov(nAntennas), nSamples=16*nAntennas)
    cn_est[:CircML]       = mmse.MLEst(snr, transform = circ_trans)

    algs = Dict{Symbol,Any}()
    algs[:GenieMMSE] = (y,h,h_cov) -> mmse_genie(y, h_cov, snr)
    algs[:GenieOMP]  = (y,h,h_cov) -> omp_genie(y, h)
    for (alg,cn) in cn_est
        algs[alg] = (y,h,h_cov) -> mmse.estimate(cn, y)
    end

    (errs,rates) = evaluate(algs, snr = snr, nBatches = nBatches, get_channel = () -> get_channel(nAntennas, nCoherence, nBatchSize), verbose = verbose)

    for alg in keys(algs)
        new_row = DataFrame(MSE        = errs[alg],
                            rate       = rates[alg],
                            Algorithm  = String(alg),
                            SNR        = snr,
                            nAntennas  = nAntennas,
                            nCoherence = nCoherence)

        if isempty(results)
            results = new_row
        else
            results = vcat(results,new_row)
        end
    end
end
if write_file
    CSV.write(filename, results)
end

max_stars_repo_name: tum-msv/learning-mmse-est


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>20.0, :alg=>GLGM06(δ=0.001)),
    "X0"   => BallInf(zeros(120), 1.0),
    "info" => Dict(
        "name" => "cdplayer",
        "dim" => 120,
        "linear" => true,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
##############
# Owner Type #
##############

type Owner <: GitHubType
    typ::Nullable{GitHubString}
    email::Nullable{GitHubString}
    name::Nullable{GitHubString}
    login::Nullable{GitHubString}
    bio::Nullable{GitHubString}
    company::Nullable{GitHubString}
    location::Nullable{GitHubString}
    gravatar_id::Nullable{GitHubString}
    id::Nullable{Int}
    public_repos::Nullable{Int}
    owned_private_repos::Nullable{Int}
    total_private_repos::Nullable{Int}
    public_gists::Nullable{Int}
    private_gists::Nullable{Int}
    followers::Nullable{Int}
    following::Nullable{Int}
    collaborators::Nullable{Int}
    blog::Nullable{HttpCommon.URI}
    url::Nullable{HttpCommon.URI}
    html_url::Nullable{HttpCommon.URI}
    updated_at::Nullable{Dates.DateTime}
    created_at::Nullable{Dates.DateTime}
    date::Nullable{Dates.DateTime}
    hireable::Nullable{Bool}
    site_admin::Nullable{Bool}
end

Owner(data::Dict) = json2github(Owner, data)
Owner(login::AbstractString, isorg = false) = Owner(Dict("login" => login, "typ" => isorg ? "User" : "Organization"))

namefield(owner::Owner) = owner.login

typprefix(isorg) = isorg ? "orgs" : "users"

#############
# Owner API #
#############

isorg(owner::Owner) = get(owner.typ, "") == "Organization"

owner(owner_obj::Owner; options...) = owner(name(owner_obj), isorg(owner_obj); options...)

function owner(owner_obj, isorg = false; options...)
    result = gh_get_json("/$(typprefix(isorg))/$(name(owner_obj))"; options...)
    return Owner(result)
end

function users(; options...)
    results, page_data = gh_get_paged_json("/users"; options...)
    return map(Owner, results), page_data
end

function check_membership(org, user; public_only = false, options...)
    scope = public_only ? "public_members" : "members"
    resp = gh_get("/orgs/$(name(org))/$scope/$(name(user))"; handle_error = false, allow_redirects = false,  options...)
    if resp.status == 204
        return true
    elseif resp.status == 404
        return false
    elseif resp.status == 302
        # For convenience, still check public membership. Otherwise, we don't know, so error
        @assert !public_only
        is_public_member = check_membership(org, user; public_only = true, options...)
        is_public_member && return true
        error("Enquiring about an Organization to which you do not have access.\n"*
              "Set `public_only=true` or provide authentication.")
    else
        handle_response_error(resp)
    end
end

function orgs(owner; options...)
    results, page_data = gh_get_paged_json("/users/$(name(owner))/orgs"; options...)
    return map(Owner, results), page_data
end

function followers(owner; options...)
    results, page_data = gh_get_paged_json("/users/$(name(owner))/followers"; options...)
    return map(Owner, results), page_data
end

function following(owner; options...)
    results, page_data = gh_get_paged_json("/users/$(name(owner))/following"; options...)
    return map(Owner, results), page_data
end

function pubkeys(owner; options...)
    results, page_data = gh_get_paged_json("/users/$(name(owner))/keys"; options...)
    output = Dict{Int,GitHubString}([(key["id"], key["key"]) for key in results])
    return output, page_data
end

repos(owner::Owner; options...) = repos(name(owner), isorg(owner); options...)

function repos(owner, isorg = false; options...)
    results, page_data = gh_get_paged_json("/$(typprefix(isorg))/$(name(owner))/repos"; options...)
    return map(Repo, results), page_data
end


text:
using GLPlot, GLAbstraction, GeometryTypes, FileIO, GLVisualize, Colors
##########################################################
# Image

a = [RGBA{N0f8}(i/512,j/512,0,1)for i=1:512, j=1:512]
# Without ImmutableArays, the color dimension is not known and you need to supply it
b = [Gray{Float32}((i*j)/512^2) for i=1:512, j=1:512]
c = loadasset("racoon.png")
d = loadasset("kitty.png")
layout!(SimpleRectangle(0f0, 0f0, 100f0, 100f0), glplot(a))
layout!(SimpleRectangle(0f0, 100f0, 100f0, 100f0), glplot(b))
layout!(SimpleRectangle(100f0, 0f0, 100f0, 100f0), glplot(c))
layout!(SimpleRectangle(100f0, 100f0, 100f0, 100f0), glplot(d))


text:
# this algorithm is modelled after the python networkx implementation:
# https://github.com/networkx/networkx/blob/master/networkx/generators/random_graphs.py#L655
function create_network(
    n::Int64, m0::Int64
)
    # check if n is smaller than m0
    if n >= m0
        # setup
        g = SimpleDiGraph(n)
        # set of nodes to connect to
        targets = collect(1:m0)
        # growing set of nodes for preferential attachment
        repeated_nodes = Array{Int64}(undef, 0)
        # initial source node
        source = m0 + 1
        # preferential attachment algorithm
        while source <= n
            for e in zip(targets, fill(source, m0))
                add_edge!(g, e[1], e[2])
            end
            append!(repeated_nodes, targets)
            targets = shuffle(repeated_nodes)[1:m0]
            source += 1
        end
        return g
    else
        error("n cannot be smaller than m0")
    end
end

function update_network!(
    state::Tuple{AbstractGraph, AbstractArray},
    config::Config
)
    graph, agent_list = state
    pref_attach_list = [src(e) for e in edges(graph) if agent_list[src(e)].active]
    for _ in 1:config.network.growth_rate
        push!(agent_list, Agent(generate_opinion(), generate_inclin_interact(), generate_check_regularity()))
        add_vertex!(graph)
        shuffle!(pref_attach_list)
        for i in 1:config.network.initial_follows
            add_edge!(graph, nv(graph), pref_attach_list[i])
        end
    end
end

# suppress output of include()
;

max_stars_repo_name: digitalemuendigkeit/twitter-abm

max_issues_repo_name: digitalemuendigkeit/twitter-abm

max_forks_repo_name: digitalemuendigkeit/twitter-abm


text:
"""
$(SIGNATURES)

Create customized distinguishable colormap for interior regions
"""
region_cmap(n)=distinguishable_colors(max(5,n),
                                      [RGB(0.85,0.6,0.6), RGB(0.6,0.85,0.6),RGB(0.6,0.6,0.85)],
                                      lchoices = range(70, stop=80, length=5),
                                      cchoices = range(25, stop=65, length=15),
                                      hchoices = range(20, stop=360, length=15)
                                      )

"""
$(SIGNATURES)

Create customized distinguishable colormap for boundary regions
"""
bregion_cmap(n)=distinguishable_colors(max(5,n),
                                      [RGB(1,0,0), RGB(0,1,0), RGB(0,0,1)],
                                      lchoices = range(50, stop=75, length=10),
                                      cchoices = range(75, stop=100, length=10),
                                      hchoices = range(20, stop=360, length=30)
                                      )

function Colors.RGB(c::String)
    c64=Colors.color_names[c]
    RGB(c64[1]/255,c64[2]/255, c64[3]/255)
end

Colors.RGB(c::Symbol)=Colors.RGB(String(c))
Colors.RGB(c::Tuple)=Colors.RGB(c...)
rgbtuple(c)=rgbtuple(Colors.RGB(c))
rgbtuple(c::RGB)=(red(c),green(c),blue(c))

"""
$(SIGNATURES)

Extract visible tetrahedra - those intersecting with the planes
`x=xyzcut[1]` or `y=xyzcut[2]`  or `z=xyzcut[3]`.

Return corresponding points and facets for each region for drawing as mesh (Makie,MeshCat)
or trisurf (pyplot)
"""
function extract_visible_cells3D(grid::ExtendableGrid,xyzcut; primepoints=zeros(0,0),Tp=SVector{3,Float32},Tf=SVector{3,Int32})
    coord=grid[Coordinates]
    cellnodes=grid[CellNodes]
    cellregions=grid[CellRegions]
    nregions=grid[NumCellRegions]

    function take(coord,simplex,xyzcut)
        all_lt=@MVector ones(Bool,3)
        all_gt=@MVector ones(Bool,3)
        for idim=1:3
            for inode=1:4
                c=coord[idim,simplex[inode]]-xyzcut[idim]
                all_lt[idim]=all_lt[idim] && (c<0.0)
                all_gt[idim]=all_gt[idim] && (c>0.0)
            end
        end
        tke=false
        tke=tke  ||   (!all_lt[1])  &&  (!all_gt[1]) && (!all_gt[2]) && (!all_gt[3])
        tke=tke  ||   (!all_lt[2])  &&  (!all_gt[2]) && (!all_gt[1]) && (!all_gt[3])
        tke=tke  ||   (!all_lt[3])  &&  (!all_gt[3]) && (!all_gt[1]) && (!all_gt[2])
    end

    faces=[Vector{Tf}(undef,0) for iregion=1:nregions]
    points=[Vector{Tp}(undef,0) for iregion=1:nregions]

    for iregion=1:nregions
        for iprime=1:size(primepoints,2)
            @views push!(points[iregion],Tp(primepoints[:,iprime]))
        end
    end
    tet=zeros(Int32,4)

    for itet=1:size(cellnodes,2)
        iregion=cellregions[itet]
        for i=1:4
            tet[i]=cellnodes[i,itet]
        end
        if take(coord,tet,xyzcut)
            npts=size(points[iregion],1)
            @views begin
                push!(points[iregion],coord[:,cellnodes[1,itet]])
                push!(points[iregion],coord[:,cellnodes[2,itet]])
                push!(points[iregion],coord[:,cellnodes[3,itet]])
                push!(points[iregion],coord[:,cellnodes[4,itet]])
                push!(faces[iregion],(npts+1,npts+2,npts+3))
                push!(faces[iregion],(npts+1,npts+2,npts+4))
                push!(faces[iregion],(npts+2,npts+3,npts+4))
                push!(faces[iregion],(npts+3,npts+1,npts+4))
            end
        end
    end
    points,faces
end

function extract_visible_bfaces3D(grid::ExtendableGrid,xyzcut; primepoints=zeros(0,0), Tp=SVector{3,Float32},Tf=SVector{3,Int32})
    cutcoord=zeros(3)

    function take(coord,simplex,xyzcut)
        for idim=1:3
            all_gt=true
            for inode=1:3
                c=coord[idim,simplex[inode]]-xyzcut[idim]
                all_gt= all_gt && c>0
            end
            if all_gt
                return false
            end
        end
        return true
    end

    coord=grid[Coordinates]
    nbfaces=num_bfaces(grid)
    bfacenodes=grid[BFaceNodes]
    bfaceregions=grid[BFaceRegions]
    nbregions=grid[NumBFaceRegions]

    faces=[Vector{Tf}(undef,0) for iregion=1:nbregions]
    points=[Vector{Tp}(undef,0) for iregion=1:nbregions]

    for iregion=1:nbregions
        for iprime=1:size(primepoints,2)
            @views push!(points[iregion],Tp(primepoints[:,iprime]))
        end
    end

    for i=1:nbfaces
        iregion=bfaceregions[i]
        tri=view(bfacenodes,:, i)
        if take(coord,tri,xyzcut)
            npts=size(points[iregion],1)
            @views push!(points[iregion],coord[:,bfacenodes[1,i]])
            @views push!(points[iregion],coord[:,bfacenodes[2,i]])
            @views push!(points[iregion],coord[:,bfacenodes[3,i]])
            push!(faces[iregion],(npts+1,npts+2,npts+3))
        end
    end
    points,faces
end

# old version with function values
function extract_visible_bfaces3D(grid::ExtendableGrid,func,xyzcut)
    cutcoord=zeros(3)

    function take(coord,simplex,xyzcut)
        for idim=1:3
            for inode=1:3
                cutcoord[inode]=coord[idim,simplex[inode]]-xyzcut[idim]
            end
            if !mapreduce(a->a<=0,*,cutcoord)
                return false
            end
        end
        return true
    end

    coord=grid[Coordinates]
    nbfaces=num_bfaces(grid)
    bfacenodes=grid[BFaceNodes]

    pmark=zeros(UInt32,size(coord,2))
    faces=ElasticArray{UInt32}(undef,3,0)
    npoints=0

    for i=1:nbfaces
        tri=view(bfacenodes,:, i)
        if take(coord,tri,xyzcut)
            for inode=1:3
                if pmark[tri[inode]]==0
                    npoints+=1
                    pmark[tri[inode]]=npoints
                end
            end
            tri=map(i->pmark[i],tri)
            append!(faces,tri)
        end
    end

    points=Array{Float32,2}(undef,3,npoints)
    values=Vector{Float32}(undef,npoints)
    for i=1:size(coord,2)
        if pmark[i]>0
            @views points[:,pmark[i]].=(coord[1,i],coord[2,i],coord[3,i])
            values[pmark[i]]=func[i]
        end
    end
    points,faces,values
end

"""
  $(SIGNATURES)
  Calculate intersections between tetrahedron with given linear
  fucntion data and plane

  Adapted from https://github.com/j-fu/gltools/blob/master/glm-3d.c#L341

  A non-empty intersection is either a triangle or a planar quadrilateral,
  define by either 3 or 4 intersection points between tetrahedron edges
  and the plane.

  Input:
  -       pointlist: 3xN array of grid point coordinates
  -    node_indices: 4 element array of node indices (pointing into pointlist and function_values)
  -   planeq_values: 4 element array of plane equation evaluated at the node coordinates
  - function_values: N element array of function values

  Mutates:
  -  ixcoord: 3x4 array of plane - tetedge intersection coordinates
  - ixvalues: 4 element array of fuction values at plane - tetdedge intersections

  Returns:
  - nxs,ixcoord,ixvalues

  This method can be used both for the evaluation of plane sections and for
  the evaluation of function isosurfaces.
"""
function tet_x_plane!(ixcoord,ixvalues,pointlist,node_indices,planeq_values,function_values; tol=0.0)

    # If all nodes lie on one side of the plane, no intersection
    if (mapreduce(a->a< -tol,*,planeq_values) || mapreduce(a->a>tol,*,planeq_values))
        return 0
    end
    # Interpolate coordinates and function_values according to
    # evaluation of the plane equation
    nxs=0
    for n1=1:3
        for n2=n1+1:4
            if planeq_values[n1]*planeq_values[n2]<tol
                nxs+=1
                t= planeq_values[n1]/(planeq_values[n1]-planeq_values[n2])
                for i=1:3
                    ixcoord[i,nxs]=pointlist[i,node_indices[n1]]+t*(pointlist[i,node_indices[n2]]-pointlist[i,node_indices[n1]])
                end
                ixvalues[nxs]=function_values[node_indices[n1]]+t*(function_values[node_indices[n2]]-function_values[node_indices[n1]])
            end
        end
    end
    return nxs
end

"""
 We should be able to parametrize this
 with a pushdata function which will remove one copy
 step for GeometryBasics.mesh creation - perhaps a meshcollector struct we
 can dispatch on.
 flevel could be flevels
 xyzcut could be a vector of plane data
 perhaps we can also collect isolines.
 Just an optional collector parameter, defaulting to somethig makie independent.

    Better yet:

 struct TetrahedronMarcher
  ...
 end
 tm=TetrahedronMarcher(planes,levels)

 foreach tet
       collect!(tm, tet_node_coord, node_function_values)
 end
 tm.colors=AbstractPlotting.interpolated_getindex.((cmap,), mcoll.vals, (fminmax,))
 mesh!(collect(mcoll),backlight=1f0)
"""

"""
   $(SIGNATURES)

"""
function marching_tetrahedra(grid::ExtendableGrid,func,planes,flevels; kwargs...)
    coord=grid[Coordinates]
    cellnodes=grid[CellNodes]
    marching_tetrahedra(coord,cellnodes,func,planes,flevels;kwargs...)
end

"""
   $(SIGNATURES)

"""
function marching_tetrahedra(coord,cellnodes,func,planes,flevels;
                             tol=0.0,
                             primepoints=zeros(0,0),
                             primevalues=zeros(0),
                             Tv=Float32,
                             Tp=SVector{3,Float32},
                             Tf=SVector{3,Int32})
    # We could rewrite this for Meshing.jl
    # CellNodes::Vector{Ttet}, Coord::Vector{Tpt}

    nplanes=length(planes)
    nlevels=length(flevels)

    all_ixfaces=Vector{Tf}(undef,0)
    all_ixcoord=Vector{Tp}(undef,0)
    all_ixvalues=Vector{Tv}(undef,0)

    @assert(length(primevalues)==size(primepoints,2))
    for iprime=1:size(primepoints,2)
        @views push!(all_ixcoord,primepoints[:,iprime])
        @views push!(all_ixvalues,primevalues[iprime])
    end

    planeq=zeros(4)
    ixcoord=zeros(3,6)
    ixvalues=zeros(6)
    cn=zeros(4)
    node_indices=zeros(Int32,4)

    plane_equation(plane,coord)=coord[1]*plane[1]+coord[2]*plane[2]+coord[3]*plane[3]+plane[4]

    function pushtris(ns,ixcoord,ixvalues)
        # number of intersection points can be 3 or 4
        if ns>=3
            last_i=length(all_ixvalues)
            for is=1:ns
                @views push!(all_ixcoord,ixcoord[:,is])
                push!(all_ixvalues,ixvalues[is])
            end
            push!(all_ixfaces,(last_i+1,last_i+2,last_i+3))
            if ns==4
                push!(all_ixfaces,(last_i+3,last_i+2,last_i+4))
            end
        end
    end

    # allocation free (besides push!)
    for itet=1:size(cellnodes,2)
        for i=1:4
            node_indices[i]=cellnodes[i,itet]
        end

        for iplane=1:nplanes
            @views map!(inode->plane_equation(planes[iplane],coord[:,inode]),planeq,node_indices)
            nxs=tet_x_plane!(ixcoord,ixvalues,coord,node_indices,planeq,func,tol=tol)
            pushtris(nxs,ixcoord,ixvalues)
        end

        for ilevel=1:nlevels
            map!(inode->(func[inode]-flevels[ilevel]),planeq,node_indices)
            nxs=tet_x_plane!(ixcoord,ixvalues,coord,node_indices,planeq,func)
            pushtris(nxs,ixcoord,ixvalues)
        end

    end
    all_ixcoord, all_ixfaces, all_ixvalues
end

########################################################################################
"""
    $(SIGNATURES)

    Collect isoline snippets on triangles ready for linesegments!

"""

function marching_triangles(grid::ExtendableGrid,func,levels)
    coord::Matrix{Float64}=grid[Coordinates]
    cellnodes::Matrix{Int32}=grid[CellNodes]
    marching_triangles(coord,cellnodes,func,levels)
end

function marching_triangles(coord,cellnodes,func,levels)
    points=Vector{Point2f0}(undef,0)
    function isect(nodes)
        (i1,i2,i3)=(1,2,3)

        f=(func[nodes[1]],func[nodes[2]],func[nodes[3]])

        f[1]  <= f[2]  ?  (i1,i2) = (1,2)   : (i1,i2) = (2,1)
        f[i2] <= f[3]  ?  i3=3              : (i2,i3) = (3,i2)
        f[i1] >  f[i2] ?  (i1,i2) = (i2,i1) : nothing

        (n1,n2,n3)=(nodes[i1],nodes[i2],nodes[i3])

        dx31=coord[1,n3]-coord[1,n1]
        dx21=coord[1,n2]-coord[1,n1]
        dx32=coord[1,n3]-coord[1,n2]

        dy31=coord[2,n3]-coord[2,n1]
        dy21=coord[2,n2]-coord[2,n1]
        dy32=coord[2,n3]-coord[2,n2]

        df31 = f[i3]!=f[i1] ? 1/(f[i3]-f[i1]) : 0.0
        df21 = f[i2]!=f[i1] ? 1/(f[i2]-f[i1]) : 0.0
        df32 = f[i3]!=f[i2] ? 1/(f[i3]-f[i2]) : 0.0

        for level ∈ levels
            if  (f[i1]<=level) && (level<f[i3])
	        α=(level-f[i1])*df31
	        x1=coord[1,n1]+α*dx31
	        y1=coord[2,n1]+α*dy31

	        if (level<f[i2])
	            α=(level-f[i1])*df21
	            x2=coord[1,n1]+α*dx21
		    y2=coord[2,n1]+α*dy21
                else
	            α=(level-f[i2])*df32
	            x2=coord[1,n2]+α*dx32
	            y2=coord[2,n2]+α*dy32
                end
                push!(points,Point2f0(x1,y1))
                push!(points,Point2f0(x2,y2))
            end
        end
    end

    for itri=1:size(cellnodes,2)
        @views isect(cellnodes[:,itri])
    end

    points
end

##############################################
# Create meshes from grid data
function regionmesh(grid,iregion)
    coord=grid[Coordinates]
    cn=grid[CellNodes]
    cr=grid[CellRegions]
    @views points=[Point2f0(coord[:,i]) for i=1:size(coord,2)]
    faces=Vector{GLTriangleFace}(undef,0)
    for i=1:length(cr)
        if cr[i]==iregion
            @views push!(faces,cn[:,i])
        end
    end
    Mesh(points,faces)
end

function bfacesegments(grid,ibreg)
    coord=grid[Coordinates]
    nbfaces=num_bfaces(grid)
    bfacenodes=grid[BFaceNodes]
    bfaceregions=grid[BFaceRegions]
    points=Vector{Point2f0}(undef,0)
    for ibface=1:nbfaces
        if bfaceregions[ibface]==ibreg
            push!(points,Point2f0(coord[1,bfacenodes[1,ibface]],coord[2,bfacenodes[1,ibface]]))
            push!(points,Point2f0(coord[1,bfacenodes[2,ibface]],coord[2,bfacenodes[2,ibface]]))
        end
    end
    points
end

function bfacesegments3(grid,ibreg)
    coord=grid[Coordinates]
    nbfaces=num_bfaces(grid)
    bfacenodes=grid[BFaceNodes]
    bfaceregions=grid[BFaceRegions]
    points=Vector{Point3f0}(undef,0)
    for ibface=1:nbfaces
        if bfaceregions[ibface]==ibreg
            push!(points,Point3f0(coord[1,bfacenodes[1,ibface]],coord[2,bfacenodes[1,ibface]],0.0))
            push!(points,Point3f0(coord[1,bfacenodes[2,ibface]],coord[2,bfacenodes[2,ibface]],0.0))
        end
    end
    points
end

############################################
"""
$(SIGNATURES)

Assume that `points` are nodes of a polyline.
Place `nmarkers` equidistant markers  at the polyline, under
the assumption that the points are transformed via the transformation
matrix M vor visualization.
"""
function markerpoints(points,nmarkers,transform)
    dist(p1,p2)=norm(transform*(p1-p2))

    llen=0.0
    for i=2:length(points)
       llen+=dist(points[i],points[i-1])
    end

    mdist=llen/(nmarkers-1)

    mpoints=[points[1]]

    i=2
    l=0.0
    lnext=l+mdist
    while i<length(points)
        d=dist(points[i],points[i-1])
        while l+d <= lnext && i<length(points)
            i=i+1
            l=l+d
            d=dist(points[i],points[i-1])
        end

        while lnext <=l+d &&  length(mpoints)<nmarkers-1
            α=(lnext-l)/d
            push!(mpoints,Point2f0(α*points[i]+ (1-α)*points[i-1]))
            lnext=lnext+mdist
        end
    end
    push!(mpoints,points[end])
end


text:
#=
UTF32Str type (UTF-32 encoding)

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones,
and other contributors to the Julia language
Licensed under MIT License, see LICENSE.md
Based in part on code for UTF32String that used to be in Julia
=#

# UTF-32 basic functions

is_unicode(str::MS_UTF32)  = true

is_ascii(str::MS_SubUTF32)   = false
is_latin(str::MS_SubUTF32)   = false
is_bmp(str::MS_SubUTF32)     = false
is_unicode(str::MS_SubUTF32) = true

# Speed this up by accessing 64 bits or more at a time
function _cnt_non_bmp(len, pnt::Ptr{UInt32})
    cnt = 0
    @inbounds for i = 1:len
        cnt += get_codeunit(pnt, i) > 0x0ffff
    end
    cnt
end

# Type _UTF32Str must have at least 1 character > 0xffff, so use other types as well
convert(::Type{<:Str{_UTF32CSE}}, ch::Unsigned) =
    (ch > 0xff
     ? (is_unicode(ch)
        ? (ch <= 0xffff ? _convert(_UCS2CSE, ch%UInt16) : _convert(_UTF32CSE, ch%UInt32))
        : strerror(StrErrors.INVALID, 0, ch))
     : _convert(_LatinCSE, ch%UInt8))

function convert(::Type{<:Str{UTF32CSE}}, str::AbstractString)
    is_empty(str) && return empty_utf32
    len, flags = unsafe_check_string(str)
    buf, pnt = _allocate(UInt32, len)
    @inbounds for ch in str
        set_codeunit!(pnt, UInt32(ch))
        pnt += sizeof(UInt32)
    end
    Str(UTF32CSE, buf)
end
convert(::Type{_UTF32Str}, str::AbstractString) = Str(str)

# This needs to handle the fact that the String type can contain invalid data
function convert(::Type{<:Str{UTF32CSE}}, str::MS_RawUTF8)
    # handle zero length string quickly
    (siz = sizeof(str)) == 0 && return empty_str(UTF32CSE)
    @preserve str begin
        pnt = pointer(str)
        # Validate UTF-8 encoding, and get number of words to create
        len, flags = fast_check_string(pnt, siz)
        # Optimize case where no characters > 0x7f, no invalid
        Str(UTF32CSE, flags == 0 ? _cvtsize(UInt32, pnt, len) : _encode_utf32(pnt, len))
    end
end

function convert(::Type{<:Str{UTF32CSE}}, str::MaybeSub{<:Str{Text2CSE}})
    # handle zero length string quickly
    (len = ncodeunits(str)) == 0 && return empty_utf32
    @preserve str begin
        pnt = pointer(str)
        # Validate and get number of characters to create
        cnt, flags = fast_check_string(pnt, len)
        Str(UTF32CSE,
            cnt == len ? _cvtsize(UInt32, pnt, cnt) : _transcode_utf16_to_utf32(pnt, cnt))
    end
end

@inline function get_cp(pnt)
    ch = get_codeunit(pnt)%UInt32
    # Handle ASCII characters
    if ch <= 0x7f
        ch, pnt + 1
    elseif ch < 0xe0 # Handle range 0x80-0x7ff
        ((ch & 0x1f) << 6) | (get_codeunit(pnt + 1) & 0x3f), pnt + 2
    elseif ch < 0xf0 # Handle range 0x800-0xffff
        ch = get_utf8_3byte(pnt += 2, ch)%UInt32
        # Handle surrogate pairs (should have been encoded in 4 bytes)
        (is_surrogate_lead(ch)
            # Build up 32-bit character from ch and trailing surrogate in next 3 bytes
         ? (get_supplementary(ch, get_utf8_3byte(pnt + 3, get_codeunit(pnt + 1)%UInt32)), pnt + 4)
         : (ch, pnt + 1))
    else # Handle range 0x10000-0x10ffff
        get_utf8_4byte(pnt + 3, ch), pnt + 4
    end
end

function _encode_utf32(pnt, len)
    buf, out = _allocate(UInt32, len)
    # has multi-byte UTF-8 sequences
    fin = bytoff(out, len)
    @inbounds while out < fin
        ch, pnt = get_cp(pnt)
        set_codeunit!(out, ch)
        out += sizeof(UInt32)
    end
    buf
end

# transcode to vector of UInt32 from validated UTF8
function _transcode_utf8_to_utf32(pnt, len)
    buf, out = _allocate(UInt32, len)
    fin = bytoff(out, len)
    while out < fin
        ch = get_codeunit(pnt)%UInt32
        # Handle ASCII characters
        if ch <= 0x7f
            set_codeunit!(out, ch)
        # Handle range 0x80-0x7ff
        elseif ch < 0xe0
            set_codeunit!(out, ((ch & 0x1f) << 6) | (get_codeunit(pnt += 1) & 0x3f))
        elseif ch < 0xf0
            # Handle range 0x800-0xffff
            set_codeunit!(out, get_utf8_3byte(pnt += 2, ch)%UInt32)
        else
            # Handle range 0x10000-0x10ffff
            set_codeunit!(out, get_utf8_4byte(pnt += 3, ch))
        end
        pnt += 1
        out += sizeof(UInt32)
    end
    buf
end

# transcode to vector of UInt32 from validated UTF16
function _transcode_utf16_to_utf32(pnt, len)
    buf, out = _allocate(UInt32, len)
    fin = bytoff(out, len)
    @inbounds while out < fin
        ch = get_codeunit(pnt)%UInt32
        pnt += sizeof(UInt16)
        is_surrogate_lead(ch) &&
            (ch = get_supplementary(ch, get_codeunit(pnt)) ; pnt += sizeof(UInt16))
        set_codeunit!(out, ch)
        out += sizeof(UInt32)
    end
    buf
end

# This can rely on the fact that a UTF8Str is always valid
function convert(::Type{<:Str{UTF32CSE}}, str::MS_UTF8)
    # handle zero length string quickly
    (len = ncodeunits(str)) == 0 && return empty_str(UTF32CSE)
    cnt = _length(MCU, str)
    # Optimize case where no characters > 0x7f, otherwise has multi-byte UTF-8 sequences
    Str(UTF32CSE, @preserve str (cnt == len
                                 ? _cvtsize(UInt32, pointer(str), cnt)
                                 : _transcode_utf8_to_utf32(pointer(str), cnt)))
end

# This can rely on the fact that an ASCIIStr, LatinStr, UCS2Str is always valid
convert(::Type{<:Str{UTF32CSE}},
        str::MaybeSub{Str{C}}) where {C<:Union{Binary_CSEs,ASCIICSE,Latin_CSEs,UCS2_CSEs}} =
    ((len = ncodeunits(str)) == 0
     ? empty_str(UTF32CSE)
     : @preserve str Str(UTF32CSE, _cvtsize(UInt32, pointer(str), len)))

# This can rely on the fact that a UTF16Str is always valid
function convert(::Type{<:Str{UTF32CSE}}, str::MS_UTF16)
    # handle zero length string quickly
    (len = ncodeunits(str)) == 0 && return empty_utf32
    # Get number of characters to create
    cnt = _length(MCU, str)
    # No surrogate pairs, do optimized copy
    Str(UTF32CSE, @preserve str (cnt == len
                                 ? _cvtsize(UInt32, pointer(str), cnt)
                                 : _transcode_utf16_to_utf32(pointer(str), cnt)))
end

function convert(::Type{<:Str{UTF16CSE}}, str::MaybeSub{Str{<:UTF32_CSEs}})
    # handle zero length string quickly
    (len = ncodeunits(str)) == 0 && return empty_utf16
    # get number of words to allocate
    # This can be faster just be checking how many > 0xffff
    @preserve str begin
        pnt = pointer(str)
        nonbmp = _cnt_non_bmp(len, pnt)
        # optimized path, no surrogates
        Str(UTF16CSE, nonbmp == 0 ? _cvtsize(UInt16, pnt, len) : _encode_utf16(pnt, len + nonbmp))
    end
end

function convert(::Type{<:Str{C}}, str::MaybeSub{Str{<:UTF32_CSEs}}) where {C<:UCS2_CSEs}
    # Might want to have an invalids_as argument
    # handle zero length string quickly
    (len = ncodeunits(str)) == 0 && return empty_str(C === UCS2CSE ? C : ASCIICSE)
    # Check if conversion is valid
    is_bmp(str) || strerror(StrErrors.INVALID_UCS2)
    @preserve str Str(UCS2CSE, _cvtsize(UInt16, pointer(str), len))
end

const UniRawChar = Union{UInt32, Int32, Text4Chr, Char}

function convert(::Type{<:Str{C}}, dat::AbstractArray{<:UniRawChar}) where {C<:UTF32_CSEs}
    (len = length(dat)) == 0 && return empty_str(C === UTF32CSE ? C : ASCIICSE)
    buf, pnt = _allocate(UInt32, len)
    C === _UTF32CSE && (msk = 0%UInt32)
    @preserve buf begin
        fin = bytoff(pnt, len)
        pos = 0
        while pnt < fin
            @inbounds ch = dat[pos += 1]%UInt32
            C === _UTF32CSE && (msk |= ch)
            set_codeunit!(pnt, check_valid(ch, pos))
            pnt += sizeof(UInt32)
        end
        (C === UTF32CSE || msk > 0xffff) && return Str(C, buf)
        msk > 0xff && return Str(_UCS2CSE, _cvtsize(UInt16, pnt, len))
        Str(msk > 0x7f ? _LatinCSE : ASCIICSE, _cvtsize(UInt8, pnt, len))
    end
end

function convert(::Type{<:Str{UTF32CSE}}, dat::Vector{<:Union{UInt32,Int32,Text4Chr}})
    is_valid(UTF32Str, dat) || strerror(StrErrors.INVALID)
    Str(UTF32CSE, _copysub(dat))
end

function is_valid(::Type{<:Str{UTF32CSE}}, dat::Vector{<:Union{UInt32,Int32,Text4Chr}})
    @preserve dat begin
        pnt = pointer(dat)
        fin = pnt + sizeof(dat)
        while pnt < fin
            ch = get_codeunit(pnt)
            (!is_surrogate_codeunit(ch) && ch <= 0x10ffff) || return false
            pnt += 4
        end
    end
    true
end

# Not sure this is valid anymore, want to avoid type piracy
convert(::Type{T}, v::AbstractArray{<:UniRawChar}) where {T<:AbstractString} =
    convert(T, convert(UTF32Str, v))

# To do, make this more generic, add a function that can create a vector & fill it from an Str
function convert(::Type{Vector{UInt32}}, str::MaybeSub{<:Str{<:Quad_CSEs}})
    len = ncodeunits(str)
    vec = create_vector(UInt32, len)
    @preserve str _memcpy(pointer(vec), pointer(str), len)
    vec
end

convert(::Type{<:Str{Text4CSE}}, str::MaybeSub{<:Str{C}}
        ) where {C<:Union{ASCIICSE,Text1CSE,BinaryCSE,Latin_CSEs, UCS2_CSEs}} =
            @preserve str Str(Text4CSE, _cvtsize(UInt32, pointer(str), ncodeunits(str)))

convert(::Type{<:Str{Text4CSE}}, str::MaybeSub{<:Str{<:UTF32_CSEs}}) = Str(Text4CSE, _copysub(str))

# Should check for 0xxxxxfeff and 0xfffexxxx as well, might be 16-bit encoded
_convert(pnt::Ptr{T}, len, T1) where {T<:Union{UInt32,UInt32_U,UInt32_S,UInt32_US}} =
    ((ch = unsafe_load(pnt)) == 0xfffe0000
     ? _convert(reinterpret(Ptr{T1}, pnt + 4), len - 1)
     : (ch == 0x0feff ? _convert(pnt + 4, len - 1) : _convert(pnt, len)))

function convert(::Type{T}, bytes::AbstractArray{UInt8}) where {C<:UTF32_CSEs,T<:Str{C}}
    is_empty(bytes) && return empty_utf32
    # Note: there are much better ways of detecting what the likely encoding is here,
    # this only deals with big or little-ending UTF-32
    # It really should detect at a minimum UTF-8, UTF-16 big and little
    len = length(bytes)
    len & 3 == 0 || strerror(StrErrors.ODD_BYTES_32, len, 0)
    len >>>= 2
    @preserve bytes begin
        pnt = pointer(bytes)
        if reinterpret(UInt, pnt) & 3 == 0
            buf, out = _convert(reinterpret(Ptr{UInt32}, pnt), len, swappedtype(UInt32))
        else
            buf, out = _convert(reinterpret(Ptr{UInt32_U}, pnt), len, swappedtype(UInt32_U))
        end
        # Todo, this needs better handling
        is_valid(T, out, len) || strerror(StrErrors.INVALID)
        Str(C, buf)
    end
end

function is_valid(::Type{<:Str{UTF32CSE}}, str::AbstractArray{<:UniRawChar})
    @inbounds for c in str
        ch = c%UInt32
        (!is_surrogate_codeunit(ch) && ch <= 0x10ffff) || return false
    end
    true
end
function is_valid(::Type{<:Str{_UTF32CSE}}, str::AbstractArray{<:UniRawChar})
    msk = 0%UInt32
    @inbounds for c in str
        ch = c%UInt32
        msk |= ch
        (!is_surrogate_codeunit(ch) && ch <= 0x10ffff) || return false
    end
    msk > 0xffff
end

is_valid(vec::AbstractArray{Char}) = is_unicode(vec)

function map(fun, str::MS_UTF32)
    @preserve str buf begin
        len = ncodeunits(str)
        buf, out = _allocate(UInt32, len)
        pnt = pointer(str)
        fin = bytoff(pnt, len)
        while pnt < fin
            set_codeunit!(out, check_valid(fun(get_codeunit(pnt))))
            pnt += sizeof(UInt32)
            out += sizeof(UInt32)
        end
        Str(UTF32CSE, buf)
    end
end

# Make sure the result actually still has at least one character > 0xffff
function map(fun, str::MS_SubUTF32)
    @preserve str buf begin
        msk = 0%UInt32
        len = ncodeunits(str)
        buf, out = _allocate(UInt32, len)
        pnt = pointer(str)
        fin = bytoff(pnt, len)
        while pnt < fin
            ch = check_valid(fun(get_codeunit(pnt)))
            msk |= ch
            set_codeunit!(out, ch)
            pnt += sizeof(UInt32)
            out += sizeof(UInt32)
        end
        Str(ifelse(msk < 0xffff, UTF32CSE, _UTF32CSE), buf)
    end
end


text:
abstract type Table <: AbstractMetaData end

"""
    Table()

general interface to accessing the gamedata.
"""
function Table()
    tablelist = string.(xlsxfilenames(CACHE["config"]))

    menu = RadioMenu(tablelist, pagesize=8)
    choice = request("Choose a Table to load (press 'q' to cancel) ", menu)

    if choice == -1 || choice > length(tablelist)
        println("Menu canceled.")
        return
    end
    fname = tablelist[choice]
    return Table(fname)
end
function Table(fname)
    loadtable(fname)
end

"""
    XLSXTable

wrapper around JSONWorkbook

# Note
consider providing option to convert `JSONWorksheet` to `IndexedTables.jl` because the Unreal Engine **DataTable** requires `key` column
"""
mutable struct XLSXTable{FileName} <: Table
    data::Union{JSONWorkbook, String}
    localizedata::Dict{String, Any}
    out::Dict{String, String} # output filename
    schemas::Dict{String, Any} # JSONSchema per sheet
    localize_key::Dict{String, Any}
    kwargs::Dict{String, Any}
    mtime::Float64
end
function XLSXTable(file, config)
    FileName = splitext(basename(file))[1] |> Symbol

    out = Dict{String, String}()
    schema = Dict{String, Any}()
    localize_key = Dict{String, Any}()
    localizedata = Dict{String, Any}()
    kwargs = Dict{String, Any}()

    for row in config["workSheets"]
        sheetname = row["name"]
        out[sheetname] = row["out"]
        kwargs[sheetname] = begin
            haskey(row, "kwargs") ? namedtuple(row["kwargs"]) : namedtuple(Dict{String,Any}())
        end

        # localizer
        localize_key[sheetname] = begin
            loc = get(row, "localize", missing)
            # default: uses row number for localize key
            if !ismissing(loc)
                loc = get(loc, "keycolumn", "")
            end
            loc
        end
        localizedata[sheetname] = missing
        # JSONSchema
        schema[sheetname] = lookfor_jsonschema(row["out"])
    end
    XLSXTable{FileName}(file, localizedata, out, schema, localize_key, kwargs, 0.)
end

isloaded(x::XLSXTable) = isa(x.data, JSONWorkbook)

function loaddata!(tb::XLSXTable)
    if ismodified(tb)
        kwarg_per_sheet = tb.kwargs
        tb.data = JSONWorkbook(xlsxpath(tb), keys(kwarg_per_sheet), kwarg_per_sheet)
        tb.mtime = mtime(xlsxpath(tb))
        localize!(tb)
    end
    return tb
end

function Base.getindex(tb::XLSXTable, i)
    if !isloaded(tb)
        loaddata!(tb)
    end
    getindex(tb.data, i)
end

function lookfor_jsonschema(filename)
    schema = missing
    if endswith(filename, ".json")
        if haskey(GAMEENV, "JSONSCHEMA")
            schemafile = joinpath(GAMEENV["JSONSCHEMA"], filename)
            if isfile(schemafile)
                schema = SchemaData(schemafile)
            end
        end
    end
    return schema
end

loadtable(fname::AbstractString) = loadtable(Symbol(fname))
function loadtable(fname::Symbol)
    config = CACHE["config"]
    if !in(fname, xlsxfilenames(config))
        hay = string.(xlsxfilenames(config))
        needle = string(fname)
        for (i, h) in enumerate(hay)
            if lowercase(h) == lowercase(needle)
                fname = Symbol(h)
                break
            end
            if i == length(hay)
                throw_fuzzylookupname(hay, needle)
            end
        end
    end
    loaddata!(config.tables[fname])
end

# fallback function
Base.basename(xgd::XLSXTable) = basename(xlsxpath(xgd))
Base.dirname(xgd::XLSXTable) = dirname(xlsxpath(xgd))

_filename(xgd::XLSXTable{NAME}) where {NAME} = NAME

function index(x::XLSXTable)
    if isloaded(tb)
        x.data.sheetindex
    else
        throw(AssertionError("xlsx data is not loaded"))
    end
end
function XLSXasJSON.sheetnames(xgd::XLSXTable)
    collect(keys(xgd.out))
end
filepath(tb::XLSXTable) = xlsxpath(tb)
function XLSXasJSON.xlsxpath(tb::XLSXTable)::String
    if isloaded(tb)
        xlsxpath(tb.data)
    else
        tb.data
    end
end

"""
    xlookup(value, jws::JSONWorksheet, lookup_col::JSONPointer, return_col::JSONPointer;
                find_mode = findfirst, lt=<comparison>)

https://support.office.com/en-us/article/xlookup-function-b7fd680e-6d10-43e6-84f9-88eae8bf5929

## Arguements
- lt: equality operator. you can use `==`, `<=`, `>=`
- find_mode: `findfirst`, `findlast`, `findall`

## Examples
- xlookup(5, Table("Items")["Equipment"], j"/Key", j"/Text/\$Name")
- xlookup("HP", Table("Items")["Consumable"], j"/Type", :; find_mode = findlast)
- xlookup(10, Table("Ability")["Data"], j"/Value", j"/\$Name"; lt = >=, find_mode = findall)
"""
@inline function xlookup(value, jws::JSONWorksheet, lookup_col, return_col; kwargs...)
    xlookup(
        value,
        jws,
        JSONPointer.Pointer(lookup_col),
        JSONPointer.Pointer(return_col);
        kwargs...,
    )
end
@inline function xlookup(
    value,
    jws::JSONWorksheet,
    lookup_col::JSONPointer.Pointer,
    return_col;
    find_mode::Function=findfirst,
    lt::Function=isequal,
)

    if !haskey(jws, lookup_col)
        throw(ArgumentError("$(lookup_col) does not exist"))
    end
    if isa(return_col, JSONPointer.Pointer)
        if !haskey(jws, return_col)
            throw(ArgumentError("$(return_col) does not exist"))
        end
    end

    idx = _xlookup_findindex(value, jws, lookup_col, find_mode, lt)

    if isnothing(idx)
        r = nothing
    elseif isempty(idx)
        r = Any[]
    else
        if isa(return_col, Array)
            col_indicies = map(return_col) do this
                i = findfirst(el -> el.tokens == this.tokens, keys(jws))
                if isa(i, Nothing)
                    throw(KeyError(this))
                end
                i
            end
            r = jws[idx, col_indicies]
        else
            r = jws[idx, return_col]
        end
    end
    return r
end

@inline function _xlookup_findindex(value, jws, lookup_col, find_mode, lt)
    find_mode(el -> lt(el[lookup_col], value), jws.data)
end


text:
module Utils
using jInv.Mesh
using SparseArrays
using Distributed
using LinearAlgebra
include("synthModels.jl")
include("jointStuff.jl")

export MATBased,MATFree,RBFBasedSimple1,RBFBasedSimple2,RBF10BasedSimple1, RBF10BasedSimple2,RBFBased,RBF10Based,RBF5Based;
MATBased        = "matrix";
MATFree         = "matfree";
RBFBasedSimple1 = "RBFsSimple1";
RBFBasedSimple2 = "RBFsSimple2";
RBF10BasedSimple1 = "RBFs10Simple1";
RBF10BasedSimple2 = "RBFs10Simple2";
RBFBased        = "RBFs";
RBF10Based        = "RBFs10";
RBF5Based        = "RBFs5";
export arrangeRemoteCallPCIntoLocalData
export getDivergenceMatrix2
export arrangeRemoteCallDataIntoLocalData,getIndicesOfKthWorker,divideDataToWorkersData, divideDirectDataToWorkersData
export dividePointCloudDataToWorkersData
export PCFun
function arrangeRemoteCallDataIntoLocalData(DremoteCall::Array{Future,1})
N = length(DremoteCall);
DdivideLocal = Array{Array{Float64,2}}(undef,N);
for k=1:N
	DdivideLocal[k] = fetch(DremoteCall[k]);
end
l = 0;
for k=1:N
	l+=size(DdivideLocal[k],2);
end
Dlocal = zeros(eltype(DdivideLocal[1]),size(DdivideLocal[1],1),l)
for k=1:N
	Dlocal[:,k:N:end] = DdivideLocal[k];
end
return Dlocal;
end

function arrangeRemoteCallPCIntoLocalData(DremoteCall::Array{Future,1})
N = length(DremoteCall);
println("length of dremote:",N);
println("size of dremote[1]:",size(DremoteCall)[1])
DdivideLocal = Array{Array{Float64}}(undef,N);
for k=1:N
println("size of fetch",size(fetch(DremoteCall[k])));
	DdivideLocal[k] = fetch(DremoteCall[k]);
end
l = 0;
for k=1:N
	l+=size(DdivideLocal[k],2);
end
Dlocal = zeros(eltype(DdivideLocal[1]),size(DdivideLocal[1],1),l)
for k=1:N
	Dlocal[:,k:N:end] = DdivideLocal[k];
end
return Dlocal;
end

function divideDataToWorkersData(nworkers::Int64,Dlocal)
ndips = size(Dlocal,2);
dobs 	= Array{Array{Float32}}(undef,nworkers);
for i=1:nworkers
	I_i = getIndicesOfKthWorker(nworkers,i,ndips);
	dobs[i] = Dlocal[:,I_i];
end
return dobs;
end

function divideDirectDataToWorkersData(nworkers::Int64,Dlocal)
nslices = size(Dlocal,3);
dobs 	= Array{Array{Float32,3}}(undef,nworkers);
for i=1:nworkers
	I_i = getDirectIndicesOfKthWorker(nworkers,i,nslices);
	dobs[i] = Dlocal[:,:,I_i];
end
return dobs;
end

function dividePointCloudDataToWorkersData(nworkers::Int64,Dlocal)
nslices = size(Dlocal,1);
dobs 	= Array{Array{Float64,1}}(undef,nworkers);
for i=1:nworkers
	I_i = getDirectIndicesOfKthWorker(nworkers,i,nslices);
	dobs[i] = Dlocal;
end
return dobs;
end

function getIndicesOfKthWorker(nworkers::Int64,k::Int64,nsrc::Int64)
return k:nworkers:nsrc;
end

function getDirectIndicesOfKthWorker(nworkers::Int64,k::Int64,nsrc::Int64)
return k:nworkers:nsrc;
end

"""
	mis,dmis,d2mis = SSDFun(dc,dobs,Wd)

	Input:

		dc::Array   -  simulated data
		dobs::Array -  measured data
		Wd::Array   -  diagonal weighting

	Output:

		mis::Real   -  misfit, 0.5*|dc-dobs|_Wd^2
		dmis        -  gradient
		d2mis       -  diagonal of Hessian

"""
function PCFun(dc::Union{Array{Float64},Array{Float32}},dobs::Union{Array{Float64},Array{Float32}},Wd::Union{Array{Float64},Array{Float32}})
	#println("dobs size:",size(dobs));
	#println("dc size:",size(dc))
	res   = vec(dc)-vec(dobs) # predicted - observed data
	#res   = vec(dc) # predicted - observed data
	Wd    = vec(Wd)
	mis   = .5*real(dot(Wd.*res,Wd.*res))  # data misfit
	dmis  = 1*Wd.*(Wd.*res)
	d2mis = Wd.*Wd

	# res   = vec(dc) # predicted - observed data
	# #res   = vec(dc) # predicted - observed data
	# Wd    = vec(Wd)
	# mis   = real(dot(-1*(res),ones(size(res))))/length(res);  # data misfit
	# dmis  = ones(size(res))
	# d2mis = 0*(res)

	# eps = 10^-3
	# res   = log.(vec(dc) .+ eps) # predicted - observed data
	# #res   = vec(dc) # predicted - observed data
	# Wd    = vec(Wd)
	# mis   = -1*real(dot(res,ones(size(res))))
	# dmis  = -1 ./ (dc .+eps);

	# d2mis = 1 ./ (dc.^2 .+eps);

	return mis, dmis, d2mis
end # function SSDFun

function getDivergenceMatrix2(Mesh::AbstractTensorMesh;saveMat = true)
# Mesh.Div = getDivergenceMatrix(Mesh::AbstractTensorMesh) builds face-to-cc divergence operator
	if isempty(Mesh.Div)
		if Mesh.dim==2
			D1 = kron(sparse(1.0I, Mesh.n[2], Mesh.n[2]),ddx(Mesh.n[1]))
			D2 = kron(ddx(Mesh.n[2]),sparse(1.0I, Mesh.n[1], Mesh.n[1]))
			Div = [D1 D2]
		elseif Mesh.dim==3
			D1 = kron(sparse(1.0I, Mesh.n[3], Mesh.n[3]),kron(sparse(1.0I, Mesh.n[2], Mesh.n[2]),ddx(Mesh.n[1])))
			D2 = kron(sparse(1.0I, Mesh.n[3], Mesh.n[3]),kron(ddx(Mesh.n[2]),sparse(1.0I, Mesh.n[1], Mesh.n[1])))
			D3 = kron(ddx(Mesh.n[3]),kron(sparse(1.0I, Mesh.n[2], Mesh.n[2]),sparse(1.0I, Mesh.n[1], Mesh.n[1])))
			Div = [D1 D2 D3]
		end
		println("size o Div:",size(Div));
		Vi = getVolumeInv(Mesh)
		F  = getFaceArea(Mesh)
		println("size of Vi:",size(Vi)); println("size of F:",size(F));
		Div = Vi*(Div*F);
		if saveMat
			Mesh.Div = Div;
		end
		return Div;
	else
		return Mesh.Div;
	end
end

end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
# small scrip to compute the solution of Lipman Schinwer equation
# In this case we use the donwsampled function using linear interpolation

using PyPlot
using IterativeSolvers
using SpecialFunctions

include("../src/FastConvolutionDual.jl")
include("../src/FastConvolution.jl")
include("../src/Preconditioner.jl")
include("../src/FastConvolutionSlowDual.jl")
include("../src/FastConvolutionDualDownSampled.jl")

FFTW.set_num_threads(4);
BLAS.set_num_threads(4)

#Defining Omega
h = 0.00025
k = 1/h

# size of box
a = 1;
b = 0.125;
x = collect(-a/2:h:a/2)
y = collect(-b/2:h:b/2)
(n,m) = length(x), length(y)
N = n*m
X = repmat(x, 1, m)[:];
Y = repmat(y', n,1)[:];
# we solve \triangle u + k^2(1 + nu(x))u = 0

# We use the modified quadrature in Duan and Rohklin
(ppw,D) = referenceValsTrapRule();
D0 = D[1];

window(y,alpha) = 1*(abs.(y).<=alpha/2) + (abs.(y).>alpha/2).*(abs.(y).<alpha).*
                     exp.(2*exp.(-0.5*alpha./(abs.(y)-alpha/2))./
                         ((abs.(y)-alpha/2)./(0.5*alpha)-1) )

window(y,alpha, beta) = 1*(abs.(y).<=beta) + (abs.(y).>beta).*(abs.(y).<alpha).*
                           exp.(2*exp.(-(alpha- beta)./(abs.(y)-beta))./
                               ((abs.(y)-beta)./(alpha- beta)-1) )

# Defining the smooth perturbation of the slowness
nu(x,y) = -0.05*(sin(4*pi*x/(0.96))).*
          window(y,0.96*b/2, 0.48*b/2).*
          window(x,0.96*0.5, 0.3);

figure(1); clf();
imshow(real(reshape( (1 + nu(X,Y)),n,m)),
       extent=[y[1], y[end], x[1], x[end]]); cb =  colorbar();

fastconvSlowDual = buildFastConvolutionSlowDual(x,y,h,k,nu, [1.0, 0.0],
                                                quadRule ="Greengard_Vico");

rhsSlowDual = -k^2*nu(X,Y) + zeros(Complex128,N);

# allocating the solution
sigmaSlow = zeros(Complex128,N);

# solving the system using GMRES
@time info =  gmres!(sigmaSlow, fastconvSlowDual, rhsSlowDual, maxiter = 10 )
# println(info[2].residuals[:])

using Interpolations

downsampleX = 16
downsampleY = 2

fastconvSlowDualDown = FastDualDownSampled(fastconvSlowDual,
                                           downsampleX, downsampleY)

index1 = 1:downsampleX:n
index2 = 1:downsampleY:m

Sindex = spzeros(n,m)
for ii = 0:round(Integer,(n-1)/downsampleX)
    for jj = 0:round(Integer,(m-1)/downsampleY)
        Sindex[1+ii*downsampleX,1+jj*downsampleY] = 1
    end
end

index = find(Sindex[:])
S = speye(n*m, n*m)
Sampling = S[index,:];

sigmaSlowDown = zeros(Complex128, size(Sampling)[1])

@time info =  gmres!(sigmaSlowDown, fastconvSlowDualDown, Sampling*rhsSlowDual, maxiter = 10 )
# println(info[2].residuals[:])

sigmaDownsampled = reshape(sigmaSlowDown, round(Integer,(n-1)/downsampleX) +1,
                          round(Integer,(m-1)/downsampleY) +1)

figure(17);
clf();
imshow(real(sigmaDownsampled)); colorbar();
title("sigma Slow in coarse mesh")

knots = (collect(1:downsampleX:n), collect(1:downsampleY:m))

itp_real = interpolate(knots, real(sigmaDownsampled), Gridded(Linear()))
itp_imag = interpolate(knots, imag(sigmaDownsampled), Gridded(Linear()))

interpsigma = itp_real[collect(1:n),collect(1:m) ]
              + 1im*itp_imag[collect(1:n),collect(1:m) ]

figure(15);
clf();
imshow(real(interpsigma)); colorbar();
title("interpolated sigma slow")

normSigmaSlow = maximum(abs.(sigmaSlow))

figure(15);
clf();
imshow(real(interpsigma - reshape(sigmaSlow,n,m))/normSigmaSlow)
colorbar()
title("relative error");

max_stars_repo_name: Forgotten/Lippmann-Schingwer

max_issues_repo_name: Forgotten/Lippmann-Schingwer

max_forks_repo_name: Forgotten/Lippmann-Schingwer


text:
# Test file for Julia

# Comment

#= Multi-line
   comment =#

#BEGIN

#END

"""
 Multi-line string
"""
```
 Multi-line command
```
raw"string\a\\a"

let z = zip(1:2, 3:4, 5:6)
    @test size(z) == (2,)
    @test collect(z) == [(1,3,5), (2,4,6)]
    @test eltype(z) == Tuple{Int,Int,Int}
end

@testset "generic conversion from Integer" begin
    x = rand(Int128)
    @test BigInt(x) % Int128 === x
    y = rand(UInt128)
    @test BigInt(y) % UInt128 === y
end

@testset "show" begin
    @test sprint(show, BitSet()) == "BitSet([])"
    @test sprint(show, BitSet([1,2,3])) == "BitSet([1, 2, 3])"
    show(IOBuffer(), BitSet())
end

cmd = `$(Base.julia_cmd()) --depwarn=error --startup-file=no boundscheck_exec.jl`
if !success(pipeline(cmd; stdout=stdout, stderr=stderr))
    error("boundscheck test failed, cmd : $cmd")
end

@test iterate(I, CartesianIndex(3, typemax(Int)))[1] == CartesianIndex(4,typemax(Int))
@test iterate(I, CartesianIndex(4, typemax(Int)))    === nothing
@test_throws MethodError write(IOBuffer(), ASCIIChar('x'))
@test_throws MethodError read(IOBuffer('x'), ASCIIChar)

let header = "julia [switches] -- [programfile] [args...]"
    @test startswith(read(`$exename -h`, String), header)
    @test startswith(read(`$exename --help`, String), header)
end

@test isequal(exp(complex( Inf, NaN)), complex(-Inf, NaN))
@test isequal(exp(complex( Inf, Inf)), complex(-Inf, NaN))

# Numbers
0b10
0o01_70
0x00000000000000001111__22_2233334444
-0x0002im
1.0e10
0.00025f0
-1.5__5f0
0xdeadbeefim
0x1.8p3
0x.4p-1

chars = ['0', '1', '2', '3', 'a', 'b', 'c', 'd', 'e', 'X', 'Y', 'Z',
         '𠀋', '𠀌', '𠀍', '🂠', '🂡', '🂢', '🂣', '🂮']

@test docstrings_equal(@doc(ModuleMacroDoc), doc"I am a module")
match(r"^\s*(?:#|$)", "# a comment")

max_stars_repo_name: dawidsowa/syntax-highlighting

max_issues_repo_name: dawidsowa/syntax-highlighting

max_forks_repo_name: dawidsowa/syntax-highlighting


text:
function evalPsiFromRectangleCoil(r0,z0,dr,dz,Ic,quadRule_r,quadRule_z, xr,xz; derivative=0)

	surf=dr*dz
	jc=Ic/surf

	#define quadrature points and weigths in the rectangle
	xgr=r0 .+ quadRule_r.x *0.5 *dr
	xgz=z0 .+ quadRule_z.x *0.5 *dz

	xgr_rep=repeat(xgr,size(xgz,1))
	xgz_rep=reshape(repeat(xgz,1,size(xgr,1))',size(xgr,1)*size(xgz,1))

	wgr=quadRule_r.w *0.5 *dr
	wgz=quadRule_z.w *0.5 *dz

	wgr_rep=repeat(wgr,size(wgz,1))
	wgz_rep=reshape(repeat(wgz,1,size(wgr,1))',size(wgr,1)*size(wgz,1))

	#
	G = @.Green(xr,xz,xgr_rep,xgz_rep, derivative=derivative)
	G_PXPG = getindex.(G,:G_xy)

	val=wgr_rep.*wgz_rep.*G_PXPG*jc
	val=sum(val)

	if(derivative==0)
		return (psi=val,)
	end
	if(derivative>0)
		dxrG = getindex.(G,:dxrG_xy)
		dxzG = getindex.(G,:dxzG_xy)

		dxrval=wgr_rep.*wgz_rep.*dxrG*jc
		dxzval=wgr_rep.*wgz_rep.*dxzG*jc

		dxrval=sum(dxrval)
		dxzval=sum(dxzval)

		return (psi=val, dxrpsi=dxrval, dxzpsi=dxzval)

	end

end

function assemblyRhsRectangleCoil(r0,z0,dr,dz,Ic,quadRule_r, quadRule_z, BEMesh, quadRule)

	# dimension
	n=size(BEMesh.edges,1)

	#
	rhs=zeros(n)

	#quadrature rule on BEMesh
	xGauss=quadRule.x
	wGauss=quadRule.w
	nQuad=size(quadRule.w,1)

	#in order to remap the quadrature we need the midpoints and length of edges
	e=collect(1:n)
	edgeToNode=BEMesh.edges
	h=BEMesh.length

    P=0.5*(BEMesh.coordinates[edgeToNode[:,1],:]+BEMesh.coordinates[edgeToNode[:,2],:])
	#@show P

	#quadrature points on all edges
	PG=zeros(size(P,1),2,nQuad)
	for i in 1:nQuad
		PG[:,:,i]=P + xGauss[i]*0.5*(BEMesh.coordinates[edgeToNode[:,2],:]-BEMesh.coordinates[edgeToNode[:,1],:])

		evalC=@.evalPsiFromRectangleCoil((r0,),(z0,),(dr,),(dz,),(Ic,),(quadRule_r,),(quadRule_z,), PG[:,1,i], PG[:,2,i])

		psiC = getindex.(evalC,:psi)

		rhs=rhs+0.5*wGauss[i]*h.*psiC
	end

	return rhs

end