text:
include("node.jl")

using PyPlot
using DataStructures

ϵ= 1e-12

type Tree
    dict::Dict{Int, Node}
    idMax
    root
    nSource
    nTarget

    rank
    maxLevel

    sourceTree
    normalTree
    targetTree

    center
    radius
    Tree() = new(Dict(), 0, -1, 0, 0, 0, 0, [], [], [], [Inf, Inf], [Inf, Inf])
end

function populate!(t::Tree, source, normal, target, nSource::Int, nTarget::Int, rank::Int, maxLevel::Int)
    t.rank = rank
    t.nSource = nSource
    t.nTarget = nTarget
    t.maxLevel = 0
    t.sourceTree = source
    t.normalTree = normal
    t.targetTree = target

    getCenterRadius!(t, t.sourceTree)
    # set root's level and nodeIndex
    t.root = 1; t.dict[1] = Node(0, 0)
    t.idMax = t.root
    t.dict[t.root].nSource = t.nSource
    t.dict[t.root].nTarget = t.nTarget
    t.dict[t.root].center = t.center
    t.dict[t.root].radius = t.radius
    t.dict[t.root].sourceIndex = Array(1:nSource)
    t.dict[t.root].targetIndex = Array(1:nTarget)

    assignChildren!(t, t.root, maxLevel)
    buildTree!(t)
    checkTarget(t)

end

function getCenterRadius!(t::Tree, source)
    maxX = maximum(source[:, 1])
    minX = minimum(source[:, 1])
    maxY = maximum(source[:, 2])
    minY = minimum(source[:, 2])
    t.center = ([minX, minY] + [maxX, maxY])/2.0
    t.radius = ([maxX, maxY] - [minX, minY])/2.0
end

function assignChildren!(t::Tree, id::Int, maxLevel::Int)
    @assert(length(t.dict) >= id)
    @assert(t.root != -1)
    node = t.dict[id]
    # get target points
    for k = 1:node.nTarget
        push!(node.target, t.targetTree[node.targetIndex[k], :])
    end

    # no source points
    if (node.nSource == 0)
        node.isLeaf = true
        node.isEmpty = true
    end

    if (node.nSource != 0)
        for k = 1:node.nSource
            push!(node.source, t.sourceTree[node.sourceIndex[k], :])
            push!(node.normal, t.normalTree[node.sourceIndex[k], :])
        end

        # divide into subtrees
        if (node.nSource <= t.rank || node.nLevel == maxLevel)
            # at leaf
            node.isLeaf = true
            if t.maxLevel < node.nLevel
                t.maxLevel = node.nLevel
            end
        else
            # not a leaf, cut into children
            for k = 1:4
                t.idMax += 1
                push!(node.child, t.idMax)
                t.dict[t.idMax] = Node(node.nLevel + 1, k)
                t.dict[t.idMax].parent = id
                t.dict[t.idMax].center[1] = node.center[1] + (((k-1) & 1) - 0.5) * node.radius[1]
                t.dict[t.idMax].center[2] = node.center[2] + (((k-1) >> 1) - 0.5) * node.radius[2]
                t.dict[t.idMax].radius[1] = node.radius[1] * 0.5
                t.dict[t.idMax].radius[2] = node.radius[2] * 0.5
                t.dict[t.idMax].nSource = 0 # pre-allocate
                t.dict[t.idMax].nTarget = 0 # pre-allocate
            end

            # distribute all the particles, complete
            # sourceIndex, targetIndex, nSource, nTarget
            for k = 1:node.nSource
                # id is parent.
                parentSourceIndex = node.sourceIndex[k]
                firstBit = t.sourceTree[parentSourceIndex, 2] < node.center[2] ? 0:1
                secondBit = t.sourceTree[parentSourceIndex, 1] < node.center[1] ? 0:1
                childIndex = 2 * firstBit + secondBit
                cid = node.child[childIndex + 1]
                push!(t.dict[cid].sourceIndex, parentSourceIndex)
                t.dict[cid].nSource += 1
            end

            for k = 1:node.nTarget
                parentTargetIndex = node.targetIndex[k]
                firstBit = t.targetTree[parentTargetIndex, 2] < node.center[2] ? 0:1
                secondBit = t.targetTree[parentTargetIndex, 1] < node.center[1] ? 0:1
                childIndex = 2 * firstBit + secondBit
                cid = node.child[childIndex + 1]
                push!(t.dict[cid].targetIndex, parentTargetIndex)
                t.dict[cid].nTarget += 1
            end

            for k = 1: 4
                assignChildren!(t, node.child[k], maxLevel)
            end
        end
    end
end

function buildTree!(t::Tree)
    rootId = t.root
    d = t.dict
    rootCenter = d[rootId].center
    rootRadius = d[rootId].radius
    xMin, xMax = rootCenter[1] - rootRadius[1], rootCenter[1] + rootRadius[1]
    yMin, yMax = rootCenter[2] - rootRadius[2], rootCenter[2] + rootRadius[2]

    nodeQueue = Queue(Int)
    enqueue!(nodeQueue, rootId)
    while(!isempty(nodeQueue))
        frontId = dequeue!(nodeQueue)
        buildNode!(t, frontId, xMin, xMax, yMin, yMax)
        if length(d[frontId].child) != 0
            for k = 1:4
                enqueue!(nodeQueue, d[frontId].child[k])
            end
        end
    end
end

function buildNode!(t::Tree, id::Int, xMin, xMax, yMin, yMax)
    uSet = Set{Int}()
    vSet = Set{Int}()
    wSet = Set{Int}()
    xSet = Set{Int}()
    d = t.dict
    node = d[id]

    if (node.parent != -1)
        parentNode = t.dict[node.parent]
        dx = node.radius[1]
        dy = node.radius[2]
        xStart = parentNode.center[1] - dx
        yStart = parentNode.center[2] - dy

        for col = -2:3
            for row = -2:3
                # in 6 x 6 grid
                currentX = xStart + 2 * col * dx
                currentY = yStart + 2 * row * dy
                if  (currentX >= xMin &&  currentX <= xMax &&
                     currentY >= yMin && currentY <= yMax &&
                    !(abs(node.center[1]-currentX) <ϵ  && abs(node.center[2] - currentY) < ϵ))
                    # within bounds and not equal to itself.
                    currentId = findNode(t, 1, [currentX, currentY])
                    adjacent = isAdjacent(t, id, currentId)

                    currentNode = d[currentId]
                    if (currentNode.nLevel < node.nLevel)
                        # from coarse mesh
                        if (adjacent)
                            if currentNode.isLeaf
                                push!(uSet, currentId)
                            end
                        else
                            push!(xSet, currentId)
                        end
                    end

                    if (currentNode.nLevel == node.nLevel)
                        if (!adjacent)
                            push!(vSet, currentId)
                        else
                            if (node.isLeaf)
                                rest = Queue(Int)
                                enqueue!(rest, currentId)
                                while(!isempty(rest))
                                    frontId = dequeue!(rest)
                                    frontNode = d[frontId]
                                    if (!isAdjacent(t, frontId, id))
                                        push!(wSet, frontId)
                                    else
                                        if (frontNode.isLeaf)
                                            push!(uSet, frontId)
                                        else

                                            for k = 1:4
                                                enqueue!(rest, frontNode.child[k])
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    if (node.isLeaf)
        push!(uSet, id)
    end

    node.uList = Int[]
    node.vList = Int[]
    node.wList = Int[]
    node.xList = Int[]

    for u in uSet
        push!(node.uList, u)
    end

    for v in vSet
        push!(node.vList, v)
    end

    for w in wSet
        push!(node.wList, w)
    end

    for x in xSet
        push!(node.xList, x)
    end

    node.nUList = length(node.uList)
    node.nVList = length(node.vList)
    node.nWList = length(node.wList)
    node.nXList = length(node.xList)
end

function findNode(t::Tree, id::Int, center)
    # tree traversal
    node = t.dict[id]
    if norm(node.center - center, 1) < ϵ
        return id
    else
        if node.isLeaf
            return id
        else
            col = node.center[1] > center[1] ? 0 : 1
            row = node.center[2] > center[2] ? 0 : 1
            id = 2 * row + col
            return findNode(t, node.child[id + 1], center)
        end
    end
end

function isAdjacent(t::Tree, aId, bId)
    nodeA = t.dict[aId]
    nodeB = t.dict[bId]
    diff = abs(nodeA.center - nodeB.center)
    r = nodeA.radius + nodeB.radius

    xAdjacent = (abs(diff[1]- r[1]) < ϵ) && (r[2] >= diff[2] - ϵ) # not far on the other direction
    yAdjacent = (abs(diff[2]- r[2]) < ϵ) && (r[1] >= diff[1] - ϵ) # not far on the other direction

    return xAdjacent || yAdjacent
end

function display(t::Tree, checknode=[])

    #scatter(t.sourceTree[:, 1], t.sourceTree[:, 2])
    lines = Any[]
    d = t.dict
    for i = 1 : length(d)
        cx = d[i].center[1]
        cy = d[i].center[2]
        rx = d[i].radius[1]
        ry = d[i].radius[2]
        xs = [cx - rx, cx + rx, cx + rx, cx - rx]
        ys = [cy - ry, cy - ry, cy + ry, cy + ry]
        push!(lines, collect(zip(xs, ys)))
    end
    line_seg = matplotlib[:collections][:LineCollection](lines)
    fig = figure("quadtree")
    ax = axes()
    ax[:add_collection](line_seg)
    #scatter(t.sourceTree[:, 1], t.sourceTree[:, 2], s = 3.8, alpha=0.5)
    #scatter(t.targetTree[:, 1], t.targetTree[:, 2], s = 5, alpha = 0.9)
    axis("image")
    if length(checknode) != 0
        for id in checknode
            checkList(t, id)
        end
    end

end

function checkList(t::Tree, nodeId)
    node = t.dict[nodeId]
    ucenterX = []
    ucenterY = []
    vcenterX = []
    vcenterY = []
    wcenterX = []
    wcenterY = []
    xcenterX = []
    xcenterY = []
    for id in node.uList
        nd = t.dict[id]
        push!(ucenterX, nd.center[1])
        push!(ucenterY, nd.center[2])
    end

    for id in node.vList
        nd = t.dict[id]
        push!(vcenterX, nd.center[1])
        push!(vcenterY, nd.center[2])
    end
    for id in node.wList
        nd = t.dict[id]
        push!(wcenterX, nd.center[1])
        push!(wcenterY, nd.center[2])
    end
    for id in node.xList
        nd = t.dict[id]
        push!(xcenterX, nd.center[1])
        push!(xcenterY, nd.center[2])
    end

    scatter(ucenterX, ucenterY, color = [1,0,0], s = 16.0, alpha = 0.7)
    scatter(vcenterX, vcenterY, color = [0,1.0,0.], s = 16.0, alpha = 0.7)
    scatter(wcenterX, wcenterY, color = [0,0,1], s = 16.0, alpha = 0.7)
    scatter(xcenterX, xcenterY, color = [0.5,0.5,0.5], s = 16.0, alpha = 0.7)
    scatter(node.center[1], node.center[2], color=[0,0,0], s = 16.0, alpha = 0.7)
end

function checkTarget(t::Tree)
    d, l = t.dict, length(t.dict)
    for i = 1:l
        node = d[i]
        ru = node.center + node.radius
        ld = node.center - node.radius
        nt = node.nTarget
        ti = node.targetIndex
        for j = 1: nt
            coordinate = t.targetTree[ti[j], :]
            if !((coordinate[1] <= ru[1] + ϵ && ϵ + coordinate[1] >= ld[1])
                && (coordinate[2] <= ru[2] + ϵ && ϵ + coordinate[2] >= ld[2]))
                @printf("(%6.4f, %6.4f) is outside of source box.\n", coordinate[1], coordinate[2])
            end
        end
    end
end


text:
cd(@__DIR__); include("setups/grid7x3.jl");

## TSD (T = 0.1)
D = natural_eigdist(𝚽, 𝛌, Q; T = 0.1, distance = :TSD)
E = transform(fit(MDS, D, maxoutdim=2, distances=true))
# E[1, :] .*= sign(E[1, 2]); E[2, :] .*= -sign(E[1, 2])
plt = grid7x3_mds_heatmaps(E, 𝚽; backend = :pyplot)
savefig(plt, "../figs/Grid7x3_MDS_TSD_T01.png")
display(plt)

ext: jl


text:
function sum(dv::AbstractDistributionVector{<:LogNormal}, ms::MissingStrategy=PassMissing();
    isgapfilled::AbstractVector{Bool} = Falses(length(dv)))
    length(dv) == length(isgapfilled) || error(
        "argument gapfilled must have the same length as dv ($(length(dv))" *
        "but was $(length(isgapfilled)).")
    if isa(ms, HandleMissingStrategy)
        nonmissing = findall(.!ismissing.(dv))
        #return( dv[nonmissing], isgapfilled[nonmissing])
        if !isempty(nonmissing)
            return(sum(@inbounds(dv[nonmissing]),
                isgapfilled = @inbounds(isgapfilled[nonmissing])))
        end
    end
    # uncorrelated, only sum diagonal
    Ssum = s = Ssumnonfilled = zero(eltype(nonmissingtype(eltype(dv))))#zero(T)
    nterm = 0
    for (i,d) in enumerate(dv)
        μ,σ = params(d)
        Si = exp(μ + abs2(σ)/2)
        Ssum += Si
        if !isgapfilled[i]
            Ssumnonfilled += Si
            s += abs2(σ) * abs2(Si)
            nterm += 1
        end
    end
    nterm > 0 || error("Expected at least one nonmissing term, but mu = $μ")
    σ2eff = s/abs2(Ssumnonfilled)
    μ_sum = log(Ssum) - σ2eff/2
    LogNormal(μ_sum, √σ2eff)
end

function sum(dv::AbstractDistributionVector{D},
    acf::AutoCorrelationFunction, ms::MissingStrategy=PassMissing();
    isgapfilled::AbstractVector{Bool} = Falses(length(dv)),
    storage::AbstractVector{Union{Missing,ST}} =
        Vector{Union{Missing,eltype(D)}}(undef, length(dv)),
    method::Val{M} = Val(:vector)) where
    {D<:LogNormal, ST<:eltype(D), M}
    #storage = Vector{Union{Missing,eltype(D)}}(undef, length(dv))
    if M == :vector
        return(sum_lognormals(
            dv, acf, ms, isgapfilled=isgapfilled, storage = storage))
    end
    if M == :bandedmatrix
        corrM = Symmetric(cormatrix_for_acf(length(dv), coef(acf)))
        return(sum_lognormals(
            dv, corrM, ms; isgapfilled=isgapfilled, storage = storage))
    end
    error("Unknown method $method. Supported are Val(:vector) and Val(:bandedmatrix)")
end

function sum_lognormals(dv::AbstractDistributionVector{D},
    acf::AutoCorrelationFunction, ms::MissingStrategy=PassMissing();
    isgapfilled::AbstractVector{Bool} = Falses(length(dv)),
    storage::AbstractVector{Union{Missing,DS}} =
       Vector{Union{Missing,eltype(D)}}(undef, length(dv))) where
    #{D<:LogNormal, ST<:, B}
    {D<:LogNormal, DS<:eltype(D)}
    #details<< Implements estimation according to
    # Messica A(2016) A simple low-computation-intensity model for approximating
    # the distribution function of a sum of non-identical lognormals for
    # financial applications. 10.1063/1.4964963
    μ = params(dv, Val(1))
    σ = params(dv, Val(2))
    coef_acf1 = @view coef(acf)[2:end]
    corrlength = length(coef_acf1)
    #acfm = vcat(reverse(coef_acf), 1, coef_acf)
    #use OffsetArrays so that index corresponds to lag: acfm[0] == 1
    acfm = OffsetArray(vcat(reverse(coef_acf1), 1, coef_acf1), -length(coef_acf1)-1)
    n = length(μ)
    @. storage = exp(μ + abs2(σ)/2)
    nmissing = count(ismissing.(storage))
    !(isa(ms, HandleMissingStrategy)) && nmissing != 0 && error(
        "Found missing values. Use argument 'SkipMissing()' " *
        "to sum over nonmissing.")
    # 0 in storage has the effect of not contributing to Ssum nor s
    # For excluding terms of gapfilled records, must make
    # to either set storage[gapfilled] to zero or check in the product
    Spure = disallowmissing(replace(storage, missing => 0.0))
    σpure = disallowmissing(replace(σ, missing => 0.0))
    Ssum = sum(Spure)
    # after computing Ssum, can also set gapfilled to zero
    Spure[isgapfilled] .= 0.0
    s = Ssumunfilled = zero(eltype(σpure))
    for i in 1:n
        iszero(Spure[i]) && continue # nothing added
        Ssumunfilled += Spure[i]
        jstart = max(1, i - corrlength)
        jend = min(n, i + corrlength)
        for j in jstart:jend
            #acf_ind = (j-i + corrlength +1)
            # sij will be zero if sigma or storage is missing (replaced by zero)
            # Sj moved to start to help multiplication by early zero
            sij = Spure[j] * acfm[j-i] * σpure[i] * σpure[j] * Spure[i]
            s += sij
        end
    end
    σ2eff = s/abs2(Ssumunfilled)
    μ_sum = log(Ssum) - σ2eff/2
    tmp = μ_sum, √σ2eff
    #@show Ssum, s, n - nmissing
    LogNormal(μ_sum, √σ2eff)
end

function sum(dv::AbstractDistributionVector{D},
    corr::Symmetric, ms::MissingStrategy=PassMissing();
    isgapfilled::AbstractArray{Bool,1}=Falses(length(dv)),
    storage::AbstractVector{Union{Missing,DS}} =
       Vector{Union{Missing,eltype(D)}}(undef, length(dv))) where
    {D<:LogNormal, DS<:eltype(D)}
    sum_lognormals(
        dv, corr, ms; isgapfilled=isgapfilled, storage = storage)
end

function sum_lognormals(dv::AbstractDistributionVector{D},
    corr::Symmetric, ms::MissingStrategy=PassMissing();
    isgapfilled::AbstractVector{Bool} = Falses(length(dv)),
    storage::AbstractVector{Union{Missing,DS}} =
        Vector{Union{Missing,eltype(D)}}(undef, length(dv))) where
    {D<:LogNormal, DS<:eltype(D)}
    μ = params(dv, Val(1))
    σ = params(dv, Val(2))
    # storage = allowmissing(similar(μ))
    @. storage = exp(μ + abs2(σ)/2)
    nmissing = count(ismissing, storage)
    anymissing = nmissing != 0
    !(isa(ms, HandleMissingStrategy)) && anymissing && error(
         "Found missing values. Use argument 'SkipMissing()' " *
         "to sum over nonmissing.")
    Ssum::nonmissingtype(eltype(storage)) = sum(skipmissing(storage))
    # gapfilled records only used for Ssum, can set the to 0 now
    # so they do not contribute to s and Ssumfilled for computation of σ2eff
    storage[isgapfilled] .= 0
    Ssumunfilled::nonmissingtype(eltype(storage)) = sum(skipmissing(storage))
    @. storage = σ * storage  # do only after Ssum
    # setting storage to zero results in summing zero for missing records
    # which is the same as filtering both storage and corr
    anymissing && replace!(storage, missing => 0.0)
    #s = transpose(disallowmissing(storage)) * corr * disallowmissing(storage)
    #Spure = view_nonmissing(storage) # non-allocating
    Spure = disallowmissing(storage) # allocating - tested: is faster than the view
    s = transpose(Spure) * corr * Spure
    σ2eff = s/abs2(Ssumunfilled)
    μ_sum = log(Ssum) - σ2eff/2
    #@show Ssum, s, length(storage) - nmissing
    LogNormal(μ_sum, √σ2eff)
end

mean(dv::AbstractDistributionVector{<:LogNormal},
    ms::MissingStrategy=PassMissing(); kwargs...) =
    mean_lognormals(dv, ms; kwargs...)
mean(dv::AbstractDistributionVector{<:LogNormal}, corr::Symmetric,
    ms::MissingStrategy=PassMissing(); kwargs...) =
    mean_lognormals(dv, corr, ms; kwargs...)
mean(dv::AbstractDistributionVector{<:LogNormal}, acf::AutoCorrelationFunction,
    ms::MissingStrategy=PassMissing();
    kwargs...) =
    mean_lognormals(dv, acf, ms; kwargs...)

function mean_lognormals(dv::AbstractDistributionVector{<:LogNormal}, x...;
    kwargs...)
    ds = sum(dv, x...; kwargs...)
    n = count(x -> !ismissing(x), dv)
    # multiplicative uncertainty stays the same
    LogNormal(ds.μ/n, ds.σ)
end


text:
export Circulant, Toeplitz, embed

struct Toeplitz{T} <: AbstractArray{T, 2}
	c :: Vector{T}
end

getindex(T::Toeplitz, i::Int, j::Int) = T.c[i-j+div(length(T.c)+1,2)]
isassigned(T::Toeplitz, i::Int, j::Int) = isassigned(T.c,i-j+div(length(T.c)+1,2))

size(T::Toeplitz, r::Int) = (r==1 || r==2) ? 1 + div(length(T.c),2) :
    throw(ArgumentError("Invalid dimension $r"))
size(T::Toeplitz) = size(T,1), size(T,2)

# Fast matrix x vector multiplication via embedding Toeplitz() into Circulant()
function *(A::Toeplitz{T},x::Vector{T}) where T
    n=length(A.c)
    k=div(n+1,2)
    C=Circulant([A.c[k:n];A.c[1:k-1]])
    (C*[x;zeros(T,k-1)])[1:k]
end

function mul!(y::StridedVector{T},A::Toeplitz{T},x::StridedVector{T}) where T
    n=length(A.c)
    k=div(n+1,2)
    C=Circulant([A.c[k:n];A.c[1:k-1]])
    xx=[x;zeros(T,k-1)]
    yy=mul!(similar(xx),C,xx)
    copy!(y, 1, yy, 1, length(y))
    return y
end

function Matrix(To::Toeplitz{T}) where T
	n=size(To, 1)
	M=Array{T}(undef,n,n)
	for i=1:n
		M[i:n,i] = To.c[n:2n-i]
		M[1:i-1,i] = To.c[n-i+1:n-1]
	end
	M
end

struct Circulant{T} <: AbstractArray{T, 2}
	c :: Vector{T}
end

getindex(C::Circulant, i::Int, j::Int) = C.c[mod(i-j,length(C.c))+1]
isassigned(C::Circulant, i::Int, j::Int) = isassigned(C.c,mod(i-j,length(C.c))+1)
size(C::Circulant, r::Int) = (r==1 || r==2) ? length(C.c) :
    throw(ArgumentError("Invalid dimension $r"))
size(C::Circulant) = size(C,1), size(C,2)

# Fast matrix x vector via fft()
# see Golub, van Loan, Matrix Computations, John Hopkins, Baltimore, 1996, p. 202
function *(C::Circulant{T},x::Vector{T}) where T
    xt=fft(x)
    vt=fft(C.c)
    yt=vt.*xt
    typeof(x[1])==Int ? map(Int,round.(real(ifft(yt)))) : ( (T <: Real) ? map(T,real(ifft(yt))) : ifft(yt))
end

function mul!(y::StridedVector{T},C::Circulant{T},x::StridedVector{T}) where T
    xt=fft(x)
    vt=fft(C.c)
    yt=ifft(vt.*xt)
    if T<: Int
        map!(round,y,yt)
    elseif T<: Real
        map!(real,y,yt)
    else
        copy!(y,yt)
    end
    return y
end

function Matrix(C::Circulant{T}) where T
	n=size(C, 1)
	M=Array{T}(undef,n,n)
	for i=1:n
		M[i:n,i] = C.c[1:n-i+1]
		M[1:i-1,i] = C.c[n-i+2:n]
	end
	M
end

function embed(To::Toeplitz{T}) where T
    return Circulant([To.c[div(end+1,2):end];To.c[1:div(end-1,2)]])
end


text:
# Julia wrapper for header: /usr/local/include/libavutil/mem.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_malloc,
    av_malloc_array,
    av_realloc,
    av_reallocp,
    av_realloc_array,
    av_reallocp_array,
    av_free,
    av_mallocz,
    av_mallocz_array,
    av_strdup,
    av_freep,
    av_memcpy_backptr,
    av_fast_realloc,
    av_fast_malloc

function av_malloc(size::Csize_t)
    ccall((:av_malloc,libavutil),Ptr{Void},(Csize_t,),size)
end

function av_malloc_array(nmemb::Csize_t,size::Csize_t)
    ccall((:av_malloc_array,libavutil),Ptr{Void},(Csize_t,Csize_t),nmemb,size)
end

function av_realloc(ptr,size::Csize_t)
    ccall((:av_realloc,libavutil),Ptr{Void},(Ptr{Void},Csize_t),ptr,size)
end

function av_reallocp(ptr,size::Csize_t)
    ccall((:av_reallocp,libavutil),Cint,(Ptr{Void},Csize_t),ptr,size)
end

function av_realloc_array(ptr,nmemb::Csize_t,size::Csize_t)
    ccall((:av_realloc_array,libavutil),Ptr{Void},(Ptr{Void},Csize_t,Csize_t),ptr,nmemb,size)
end

function av_reallocp_array(ptr,nmemb::Csize_t,size::Csize_t)
    ccall((:av_reallocp_array,libavutil),Cint,(Ptr{Void},Csize_t,Csize_t),ptr,nmemb,size)
end

function av_free(ptr)
    ccall((:av_free,libavutil),Void,(Ptr{Void},),ptr)
end

function av_mallocz(size::Csize_t)
    ccall((:av_mallocz,libavutil),Ptr{Void},(Csize_t,),size)
end

function av_mallocz_array(nmemb::Csize_t,size::Csize_t)
    ccall((:av_mallocz_array,libavutil),Ptr{Void},(Csize_t,Csize_t),nmemb,size)
end

function av_strdup(s)
    ccall((:av_strdup,libavutil),Ptr{UInt8},(Ptr{UInt8},),s)
end

function av_freep(ptr)
    ccall((:av_freep,libavutil),Void,(Ptr{Void},),ptr)
end

function av_memcpy_backptr(dst,back::Integer,cnt::Integer)
    ccall((:av_memcpy_backptr,libavutil),Void,(Ptr{UInt8},Cint,Cint),dst,back,cnt)
end

function av_fast_realloc(ptr,size,min_size::Csize_t)
    ccall((:av_fast_realloc,libavutil),Ptr{Void},(Ptr{Void},Ptr{UInt32},Csize_t),ptr,size,min_size)
end

function av_fast_malloc(ptr,size,min_size::Csize_t)
    ccall((:av_fast_malloc,libavutil),Void,(Ptr{Void},Ptr{UInt32},Csize_t),ptr,size,min_size)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
import Base: Expr

function Expr(x::EXPR)
    ret = Expr(:call)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{TopLevel})
    ret = Expr(:toplevel)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end
# Expr(x::HEAD{Tokens.LBRACE}) = :cell1d
Expr(x::HEAD{T}) where {T} = Symbol(lowercase(string(T)))
Expr(x::KEYWORD{T}) where {T} = Symbol(lowercase(string(T)))

Expr(x::EXPR{IDENTIFIER}) = Symbol(x.val)

function Expr(x::EXPR{OPERATOR{O,K,dot}}) where {O, K, dot}
    if dot
        Symbol(:., UNICODE_OPS_REVERSE[K])
    else
        UNICODE_OPS_REVERSE[K]
    end
end

Expr(x::EXPR{LITERAL{Tokens.TRUE}}) = true
Expr(x::EXPR{LITERAL{Tokens.FALSE}}) = false
function Expr(x::EXPR{HEAD{:nothing}}) end

Expr(x::EXPR{LITERAL{T}}) where {T} = Base.parse(x.val)
Expr(x::EXPR{LITERAL{Tokens.FLOAT}}) = Base.parse(x.val)
Expr(x::EXPR{LITERAL{Tokens.MACRO}}) = Symbol(x.val)
Expr(x::EXPR{LITERAL{Tokens.STRING}}) = x.val
Expr(x::EXPR{LITERAL{Tokens.TRIPLE_STRING}}) = x.val

# cross compatability for line number insertion in macrocalls
@static if VERSION < v"0.7.0-DEV.357"
    Expr(x::EXPR{LITERAL{Tokens.CMD}}) = Expr(:macrocall, Symbol("@cmd"), x.val[2:end-1])

    function Expr(x::EXPR{x_Str})
        if x.args[1] isa EXPR{BinarySyntaxOpCall}
            mname = Expr(x.args[1])
            mname.args[2] = QuoteNode(Symbol("@", mname.args[2].value, "_str"))
            ret = Expr(:macrocall, mname)
        else
            ret = Expr(:macrocall, Symbol("@", x.args[1].val, "_str"))
        end
        for i = 2:length(x.args)
            push!(ret.args, x.args[i].val)
        end
        return ret
    end

    function Expr(x::EXPR{x_Cmd})
        ret = Expr(:macrocall, Symbol("@", x.args[1].val, "_cmd"))
        for i = 2:length(x.args)
            push!(ret.args, x.args[i].val)
        end
        return ret
    end

    function Expr(x::EXPR{MacroCall})
        ret = Expr(:macrocall)
        for a in x.args
            if !(a isa EXPR{P} where P <: PUNCTUATION)
                push!(ret.args, Expr(a))
            end
        end
        ret
    end

    """
        remlineinfo!(x)
    Removes line info expressions. (i.e. Expr(:line, 1))
    """
    function remlineinfo!(x)
        if isa(x, Expr)
            id = find(map(x -> (isa(x, Expr) && x.head == :line) || (isdefined(:LineNumberNode) && x isa LineNumberNode), x.args))
            deleteat!(x.args, id)
            for j in x.args
                remlineinfo!(j)
            end
        end
        x
    end
else
    Expr(x::EXPR{LITERAL{Tokens.CMD}}) = Expr(:macrocall, Symbol("@cmd"), nothing, x.val[2:end-1])

    function Expr(x::EXPR{x_Str})
        if x.args[1] isa EXPR{BinarySyntaxOpCall}
            mname = Expr(x.args[1])
            mname.args[2] = QuoteNode(Symbol("@", mname.args[2].value, "_str"))
            ret = Expr(:macrocall, mname, nothing)
        else
            ret = Expr(:macrocall, Symbol("@", x.args[1].val, "_str"), nothing)
        end
        for i = 2:length(x.args)
            push!(ret.args, x.args[i].val)
        end
        return ret
    end

    function Expr(x::EXPR{x_Cmd})
        ret = Expr(:macrocall, Symbol("@", x.args[1].val, "_cmd"), nothing)
        for i = 2:length(x.args)
            push!(ret.args, x.args[i].val)
        end
        return ret
    end

    function Expr(x::EXPR{MacroCall})
        ret = Expr(:macrocall)
        for a in x.args
            if !(a isa EXPR{P} where P <: PUNCTUATION)
                push!(ret.args, Expr(a))
            end
        end
        insert!(ret.args, 2, nothing)
        ret
    end
    """
        remlineinfo!(x)
    Removes line info expressions. (i.e. Expr(:line, 1))
    """
    function remlineinfo!(x)
        if isa(x, Expr)
            if x.head == :macrocall && x.args[2] != nothing
                id = find(map(x -> (isa(x, Expr) && x.head == :line) || (isdefined(:LineNumberNode) && x isa LineNumberNode), x.args))
                deleteat!(x.args, id)
                for j in x.args
                    remlineinfo!(j)
                end
                insert!(x.args, 2, nothing)
            else
                id = find(map(x -> (isa(x, Expr) && x.head == :line) || (isdefined(:LineNumberNode) && x isa LineNumberNode), x.args))
                deleteat!(x.args, id)
                for j in x.args
                    remlineinfo!(j)
                end
            end
        end
        x
    end
end

Expr(x::EXPR{PUNCTUATION{K}}) where {K} = string(K)

Expr(x::EXPR{Quotenode}) = QuoteNode(Expr(x.args[end]))
Expr(x::EXPR{KEYWORD{t}}) where t = Symbol(lowercase(string(t)))

function Expr(x::EXPR{Call})
    ret = Expr(:call)
    for a in x.args
        if a isa EXPR{Parameters}
            insert!(ret.args, 2, Expr(a))
        elseif !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{Comparison})
    ret = Expr(:comparison)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{ChainOpCall})
    ret = Expr(:call, Expr(x.args[2]))
    for i = 1:length(x.args)
        if isodd(i)
            push!(ret.args, Expr(x.args[i]))
        end
    end
    ret
end

Expr(x::EXPR{BinaryOpCall}) = Expr(:call, Expr(x.args[2]), Expr(x.args[1]), Expr(x.args[3]))

function Expr(x::EXPR{BinarySyntaxOpCall})
    if x.args[2] isa EXPR{OPERATOR{WhereOp,Tokens.WHERE,false}}
        ret = Expr(:where, Expr(x.args[1]))
        for i = 3:length(x.args)
            a = x.args[i]
            if !(a isa EXPR{P} where P <: PUNCTUATION || a isa EXPR{K} where K <: KEYWORD)
                push!(ret.args, Expr(a))
            end
        end
    else
        ret = Expr(Expr(x.args[2]), Expr(x.args[1]), Expr(x.args[3]))
    end
    return ret
end

Expr(x::EXPR{ConditionalOpCall}) = Expr(:if, Expr(x.args[1]), Expr(x.args[3]), Expr(x.args[5]))

Expr(x::EXPR{ColonOpCall}) = Expr(:(:), Expr(x.args[1]), Expr(x.args[3]), Expr(x.args[5]))

function Expr(x::EXPR{UnarySyntaxOpCall})
    if x.args[1] isa EXPR{OP} where OP <: OPERATOR
        return Expr(Expr(x.args[1]), Expr(x.args[2]))
    else
        return Expr(Expr(x.args[2]), Expr(x.args[1]))
    end
end

function Expr(x::EXPR{UnaryOpCall})
    return Expr(:call, Expr(x.args[1]), Expr(x.args[2]))
end

Expr(x::EXPR{Struct}) = Expr(:type, false, Expr(x.args[2]), Expr(x.args[3]))

Expr(x::EXPR{Mutable}) = length(x.args) == 4 ? Expr(:type, true, Expr(x.args[2]), Expr(x.args[3])) : Expr(:type, true, Expr(x.args[3]), Expr(x.args[4]))

Expr(x::EXPR{Abstract}) = length(x.args) == 2 ? Expr(:abstract, Expr(x.args[2])) : Expr(:abstract, Expr(x.args[3]))
Expr(x::EXPR{Bitstype}) = Expr(:bitstype, Expr(x.args[2]), Expr(x.args[3]))
Expr(x::EXPR{Primitive}) = Expr(:bitstype, Expr(x.args[4]), Expr(x.args[3]))
Expr(x::EXPR{TypeAlias}) = Expr(:typealias, Expr(x.args[2]), Expr(x.args[3]))

function Expr(x::EXPR{Block})
    ret = Expr(:block)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    return ret
end

function Expr(x::EXPR{TupleH})
    ret = Expr(:tuple)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    return ret
end

Expr(x::EXPR{Kw}) = Expr(:kw, Expr(x.args[1]), Expr(x.args[3]))

function Expr(x::EXPR{Parameters})
    ret = Expr(:parameters)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    return ret
end

Expr(x::EXPR{InvisBrackets}) = Expr(x.args[2])
Expr(x::EXPR{Begin}) = Expr(x.args[2])

function Expr(x::EXPR{Quote})
    if x.args[2] isa EXPR{InvisBrackets} && (x.args[2].args[2] isa EXPR{OP} where OP <: OPERATOR|| x.args[2].args[2] isa EXPR{L} where L <: LITERAL || x.args[2].args[2] isa EXPR{IDENTIFIER})
        return QuoteNode(Expr(x.args[2]))
    else
        return Expr(:quote, Expr(x.args[2]))
    end
end

function Expr(x::EXPR{If})
    ret = Expr(:if)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION || a isa EXPR{K} where K <: KEYWORD)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{FunctionDef})
    ret = Expr(:function)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION || a isa EXPR{K} where K <: KEYWORD)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{Try})
    ret = Expr(:try)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION || a isa EXPR{K} where K <: KEYWORD)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{Let})
    ret = Expr(:let, Expr(x.args[end - 1]))
    for i = 1:length(x.args) - 2
        a = x.args[i]
        if !(a isa EXPR{P} where P <: PUNCTUATION || a isa EXPR{K} where K <: KEYWORD)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{Do})
    ret = Expr(x.args[1])
    insert!(ret.args, 2, Expr(:->, Expr(x.args[3]), Expr(x.args[4])))
    ret
end

fix_range(a) = Expr(a)
function fix_range(a::EXPR{BinaryOpCall})
    if (a.args[2] isa EXPR{OPERATOR{ComparisonOp,Tokens.IN,false}} || a.args[2] isa EXPR{OPERATOR{ComparisonOp,Tokens.ELEMENT_OF,false}})
        Expr(:(=), Expr(a.args[1]), Expr(a.args[3]))
    else
        Expr(a)
    end
end

function Expr(x::EXPR{For})
    ret = Expr(:for)
    if x.args[2] isa EXPR{Block}
        arg = Expr(:block)
        for a in x.args[2].args
            if !(a isa EXPR{P} where P <: PUNCTUATION)
                push!(arg.args, fix_range(a))
            end
        end
        push!(ret.args, arg)
    else
        push!(ret.args, fix_range(x.args[2]))
    end
    push!(ret.args, Expr(x.args[3]))
    ret
end

function Expr(x::EXPR{While})
    ret = Expr(:while)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION || a isa EXPR{K} where K <: KEYWORD)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{Return})
    ret = Expr(:return)
    for i = 2:length(x.args)
        a = x.args[i]
        push!(ret.args, Expr(a))
    end
    ret
end

function Expr(x::EXPR{Global})
    ret = Expr(:global)
    if x.args[2] isa EXPR{Const}
        ret = Expr(:const, Expr(:global, Expr(x.args[2].args[2])))
    elseif length(x.args) == 2 && x.args[2] isa EXPR{TupleH}
        for a in x.args[2].args
            if !(a isa EXPR{P} where P <: PUNCTUATION)
                push!(ret.args, Expr(a))
            end
        end
    else
        for i = 2:length(x.args)
            a = x.args[i]
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{Local})
    ret = Expr(:local)
    if x.args[2] isa EXPR{Const}
        ret = Expr(:const, Expr(:global, Expr(x.args[2].args[2])))
    elseif length(x.args) == 2 && x.args[2] isa EXPR{TupleH}
        for a in x.args[2].args
            if !(a isa EXPR{P} where P <: PUNCTUATION)
                push!(ret.args, Expr(a))
            end
        end
    else
        for i = 2:length(x.args)
            a = x.args[i]
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{Const})
    ret = Expr(:const)
    for i = 2:length(x.args)
        a = x.args[i]
        push!(ret.args, Expr(a))
    end
    ret
end
Expr(x::EXPR{Break}) = Expr(:break)
Expr(x::EXPR{Continue}) = Expr(:continue)

function Expr(x::EXPR{Curly})
    ret = Expr(:curly)
    for a in x.args
        if a isa EXPR{Parameters}
            insert!(ret.args, 2, Expr(a))
        elseif !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{Vect})
    ret = Expr(:vect)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{Macro})
    Expr(:macro, Expr(x.args[2]), Expr(x.args[3]))
end

Expr(x::EXPR{HEAD{:globalrefdoc}}) = GlobalRef(Core, Symbol("@doc"))

function Expr(x::EXPR{Row})
    ret = Expr(:row)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{Hcat})
    ret = Expr(:hcat)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{Vcat})
    ret = Expr(:vcat)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{Ref})
    ret = Expr(:ref)
    for a in x.args
        if a isa EXPR{Parameters}
            insert!(ret.args, 2, Expr(a))
        elseif !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{TypedHcat})
    ret = Expr(:typed_hcat)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{TypedVcat})
    ret = Expr(:typed_vcat)

    for a in x.args
        if a isa EXPR{Parameters}
            insert!(ret.args, 2, Expr(a))
        elseif !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{Comprehension})
    ret = Expr(:comprehension)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end
# NEEDS FIX (only works for 2 levels )
function Expr(x::EXPR{Flatten})
    ret = Expr(:flatten, Expr(x.args[1]))

    gen1 = ret.args[1].args[1].args[2]
    gen2 = ret.args[1].args[2]
    ret.args[1].args[1].args[2] = gen2
    ret.args[1].args[2] = gen1
    return ret
end

function Expr(x::EXPR{Generator})
    ret = Expr(:generator, Expr(x.args[1]))
    for i = 3:length(x.args)
        a = x.args[i]
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            if a isa EXPR{BinaryOpCall} && (a.args[2] isa EXPR{OPERATOR{ComparisonOp,Tokens.IN,false}} || a.args[2] isa EXPR{OPERATOR{ComparisonOp,Tokens.ELEMENT_OF,false}})
                push!(ret.args, Expr(:(=), Expr(a.args[1]), Expr(a.args[3])))
            else
                push!(ret.args, Expr(a))
            end
        end
    end
    ret
end

function Expr(x::EXPR{Filter})
    ret = Expr(:filter)
    for a in x.args
        if !(a isa EXPR{KEYWORD{Tokens.IF}})
            if a isa EXPR{BinaryOpCall} && (a.args[2] isa EXPR{OPERATOR{ComparisonOp,Tokens.IN,false}} || a.args[2] isa EXPR{OPERATOR{ComparisonOp,Tokens.ELEMENT_OF,false}})
                push!(ret.args, Expr(:(=), Expr(a.args[1]), Expr(a.args[3])))
            else
                push!(ret.args, Expr(a))
            end
        end
    end
    ret
end

function Expr(x::EXPR{TypedComprehension})
    ret = Expr(:typed_comprehension)
    for a in x.args
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

function Expr(x::EXPR{Export})
    ret = Expr(:export)
    for i = 2:length(x.args)
        a = x.args[i]
        if !(a isa EXPR{P} where P <: PUNCTUATION)
            push!(ret.args, Expr(a))
        end
    end
    ret
end

Expr(x::EXPR{ModuleH}) = Expr(:module, true, Expr(x.args[2]), Expr(x.args[3]))
Expr(x::EXPR{BareModule}) = Expr(:module, false, Expr(x.args[2]), Expr(x.args[3]))

function _get_import_block(x, i, ret)
    while x.args[i + 1] isa EXPR{OPERATOR{DotOp,Tokens.DOT,false}}
        i += 1
        push!(ret.args, :.)
    end
    while i < length(x.args) && !(x.args[i + 1] isa EXPR{PUNCTUATION{Tokens.COMMA}})
        i += 1
        a = x.args[i]
        if !(a isa EXPR{P} where P <: PUNCTUATION) && !(a isa EXPR{OPERATOR{DotOp,Tokens.DOT,false}} || a isa EXPR{OPERATOR{ColonOp,Tokens.COLON,false}})
            push!(ret.args, Expr(a))
        end
    end

    return i
end

Expr(x::EXPR{Import}) = expr_import(x, :import)
Expr(x::EXPR{ImportAll}) = expr_import(x, :importall)
Expr(x::EXPR{Using}) = expr_import(x, :using)

function expr_import(x, kw)
    col = find(a isa EXPR{o} where o <: OPERATOR{ColonOp} for a in x.args)
    comma = find(a isa EXPR{PUNCTUATION{Tokens.COMMA}} for a in x.args)
    if isempty(comma)
        ret = Expr(kw)
        i = 1
        _get_import_block(x, i, ret)
    elseif isempty(col)
        ret = Expr(:toplevel)
        i = 1
        while i < length(x.args)
            nextarg = Expr(kw)
            i = _get_import_block(x, i, nextarg)
            if i < length(x.args) && (x.args[i + 1] isa EXPR{PUNCTUATION{Tokens.COMMA}})
                i += 1
            end
            push!(ret.args, nextarg)
        end
    else
        ret = Expr(:toplevel)
        top = Expr(kw)
        i = 1
        while x.args[i + 1] isa EXPR{OPERATOR{DotOp,Tokens.DOT,false}}
            i += 1
            push!(top.args, :.)
        end
        while i < length(x.args) && !(x.args[i + 1] isa EXPR{o} where o <: OPERATOR{ColonOp})
            i += 1
            a = x.args[i]
            if !(a isa EXPR{P} where P <: PUNCTUATION) && !(a isa EXPR{OPERATOR{DotOp,Tokens.DOT,false}} || a isa EXPR{OPERATOR{ColonOp,Tokens.COLON,false}})
                push!(top.args, Expr(a))
            end
        end
        while i < length(x.args)
            nextarg = Expr(kw, top.args...)
            i = _get_import_block(x, i, nextarg)
            if i < length(x.args) && (x.args[i + 1] isa EXPR{PUNCTUATION{Tokens.COMMA}})
                i += 1
            end
            push!(ret.args, nextarg)
        end
    end
    return ret
end

function Expr(x::EXPR{FileH})
    ret = Expr(:file)
    for a in x.args
        push!(ret.args, Expr(a))
    end
    ret
end

function Expr(x::EXPR{StringH})
    ret = Expr(:string)
    for a in x.args
        push!(ret.args, Expr(a))
    end
    ret
end

UNICODE_OPS_REVERSE = Dict{Tokenize.Tokens.Kind,Symbol}()
for (k, v) in Tokenize.Tokens.UNICODE_OPS
    UNICODE_OPS_REVERSE[v] = Symbol(k)
end

UNICODE_OPS_REVERSE[Tokens.EQ] = :(=)
UNICODE_OPS_REVERSE[Tokens.PLUS_EQ] = :(+=)
UNICODE_OPS_REVERSE[Tokens.MINUS_EQ] = :(-=)
UNICODE_OPS_REVERSE[Tokens.STAR_EQ] = :(*=)
UNICODE_OPS_REVERSE[Tokens.FWD_SLASH_EQ] = :(/=)
UNICODE_OPS_REVERSE[Tokens.FWDFWD_SLASH_EQ] = :(//=)
UNICODE_OPS_REVERSE[Tokens.OR_EQ] = :(|=)
UNICODE_OPS_REVERSE[Tokens.CIRCUMFLEX_EQ] = :(^=)
UNICODE_OPS_REVERSE[Tokens.DIVISION_EQ] = :(÷=)
UNICODE_OPS_REVERSE[Tokens.REM_EQ] = :(%=)
UNICODE_OPS_REVERSE[Tokens.LBITSHIFT_EQ] = :(<<=)
UNICODE_OPS_REVERSE[Tokens.RBITSHIFT_EQ] = :(>>=)
UNICODE_OPS_REVERSE[Tokens.LBITSHIFT] = :(<<)
UNICODE_OPS_REVERSE[Tokens.RBITSHIFT] = :(>>)
UNICODE_OPS_REVERSE[Tokens.UNSIGNED_BITSHIFT] = :(>>>)
UNICODE_OPS_REVERSE[Tokens.UNSIGNED_BITSHIFT_EQ] = :(>>>=)
UNICODE_OPS_REVERSE[Tokens.BACKSLASH_EQ] = :(\=)
UNICODE_OPS_REVERSE[Tokens.AND_EQ] = :(&=)
UNICODE_OPS_REVERSE[Tokens.COLON_EQ] = :(:=)
UNICODE_OPS_REVERSE[Tokens.PAIR_ARROW] = :(=>)
UNICODE_OPS_REVERSE[Tokens.APPROX] = :(~)
UNICODE_OPS_REVERSE[Tokens.EX_OR_EQ] = :($=)
UNICODE_OPS_REVERSE[Tokens.XOR_EQ] = :(⊻=)
UNICODE_OPS_REVERSE[Tokens.RIGHT_ARROW] = :(-->)
UNICODE_OPS_REVERSE[Tokens.LAZY_OR] = :(||)
UNICODE_OPS_REVERSE[Tokens.LAZY_AND] = :(&&)
UNICODE_OPS_REVERSE[Tokens.ISSUBTYPE] = :(<:)
UNICODE_OPS_REVERSE[Tokens.ISSUPERTYPE] = :(>:)
UNICODE_OPS_REVERSE[Tokens.GREATER] = :(>)
UNICODE_OPS_REVERSE[Tokens.LESS] = :(<)
UNICODE_OPS_REVERSE[Tokens.GREATER_EQ] = :(>=)
UNICODE_OPS_REVERSE[Tokens.GREATER_THAN_OR_EQUAL_TO] = :(≥)
UNICODE_OPS_REVERSE[Tokens.LESS_EQ] = :(<=)
UNICODE_OPS_REVERSE[Tokens.LESS_THAN_OR_EQUAL_TO] = :(≤)
UNICODE_OPS_REVERSE[Tokens.EQEQ] = :(==)
UNICODE_OPS_REVERSE[Tokens.EQEQEQ] = :(===)
UNICODE_OPS_REVERSE[Tokens.IDENTICAL_TO] = :(≡)
UNICODE_OPS_REVERSE[Tokens.NOT_EQ] = :(!=)
UNICODE_OPS_REVERSE[Tokens.NOT_EQUAL_TO] = :(≠)
UNICODE_OPS_REVERSE[Tokens.NOT_IS] = :(!==)
UNICODE_OPS_REVERSE[Tokens.NOT_IDENTICAL_TO] = :(≢)
UNICODE_OPS_REVERSE[Tokens.IN] = :(in)
UNICODE_OPS_REVERSE[Tokens.ISA] = :(isa)
UNICODE_OPS_REVERSE[Tokens.LPIPE] = :(<|)
UNICODE_OPS_REVERSE[Tokens.RPIPE] = :(|>)
UNICODE_OPS_REVERSE[Tokens.COLON] = :(:)
UNICODE_OPS_REVERSE[Tokens.DDOT] = :(..)
UNICODE_OPS_REVERSE[Tokens.EX_OR] = :($)
UNICODE_OPS_REVERSE[Tokens.PLUS] = :(+)
UNICODE_OPS_REVERSE[Tokens.MINUS] = :(-)
UNICODE_OPS_REVERSE[Tokens.PLUSPLUS] = :(++)
UNICODE_OPS_REVERSE[Tokens.OR] = :(|)
UNICODE_OPS_REVERSE[Tokens.STAR] = :(*)
UNICODE_OPS_REVERSE[Tokens.FWD_SLASH] = :(/)
UNICODE_OPS_REVERSE[Tokens.REM] = :(%)
UNICODE_OPS_REVERSE[Tokens.BACKSLASH] = :(\)
UNICODE_OPS_REVERSE[Tokens.AND] = :(&)
UNICODE_OPS_REVERSE[Tokens.FWDFWD_SLASH] = :(//)
UNICODE_OPS_REVERSE[Tokens.CIRCUMFLEX_ACCENT] = :(^)
UNICODE_OPS_REVERSE[Tokens.DECLARATION] = :(::)
UNICODE_OPS_REVERSE[Tokens.CONDITIONAL] = :(?)
UNICODE_OPS_REVERSE[Tokens.DOT] = :(.)
UNICODE_OPS_REVERSE[Tokens.NOT] = :(!)
UNICODE_OPS_REVERSE[Tokens.PRIME] = Symbol(''')
UNICODE_OPS_REVERSE[Tokens.DDDOT] = :(...)
UNICODE_OPS_REVERSE[Tokens.TRANSPOSE] = Symbol(".'")
UNICODE_OPS_REVERSE[Tokens.ANON_FUNC] = :(->)
UNICODE_OPS_REVERSE[Tokens.WHERE] = :where


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 9e3f6282-3575-11eb-35e8-69e4cf57b6f4
using LinearAlgebra, Plots, PlutoUI

# ╔═╡ 16df5992-3574-11eb-3682-0981f8882895
md"""
# Generalized optimal transport

Traditional optimal transport find the optimal transportation map given fixed marginals. For many applications, this requirement is too stringent, for example, because the marginals are not balanced. So, one can define **generalized optimal transport** as:

$$\max_Q\, \langle M, Q \rangle + \epsilon H(Q) + \alpha D_\text{KL}(\mathbf{a}\mid Q \mathbb{1}) + \beta D_\text{KL}(\mathbf{b}\mid Q^\intercal \mathbb{1})\,.$$

One notes that the **Kullback-Leibler** is now used as a soft-constraints to drive the marginals of $Q$ to match $\mathbf{a}$ and $\mathbf{b}$, these themselfs don't have to be normalized.

One can solve this problem using the **generalized Sinkhorn algorithm**.

"""

# ╔═╡ 244cb164-3575-11eb-3b9d-837acaa0e9e9
function gen_sinkhorn(M, a, b; ϵ=1.0, α=1.0, β=1.0)
	K = exp.(M ./ ϵ)
	u, v = ones(length(a)), ones(length(b))
	for i in 1:10
		u .= (a ./ K * v) .^ (α / (α + ϵ))
		v .= (b ./ K' * u) .^ (β / (β + ϵ))
	end
	K .*= u
	K .*= v'
	K ./= sum(K)
	return K
end

# ╔═╡ 62aa2bf2-3576-11eb-0cc6-e31f7ad771d0
n, m = 20, 10

# ╔═╡ 3486c424-3576-11eb-35a3-0bcc30b65445
M = randn(n, m)

# ╔═╡ ac95ae3a-357b-11eb-1e16-b56d76851314
normalize(x) = x ./ sum(x)

# ╔═╡ 5ecc17ca-3576-11eb-00ce-2fe7c5d8f55b
a = rand(n) #|> normalize

# ╔═╡ 72fc2c12-3576-11eb-23d1-c1d14da6de89
b = rand(m) #|> normalize

# ╔═╡ 3d0844ac-3577-11eb-3086-112eb6425f79
plot(bar(a), bar(b), layout=(2,1))

# ╔═╡ d9ffdbac-3576-11eb-24f1-55d257632108
@bind α Slider(0:0.01:10, default=1)

# ╔═╡ f22bc95c-3576-11eb-2a14-9d92a4497f68
@bind β Slider(0:0.01:10, default=1)

# ╔═╡ cd231944-357b-11eb-24b6-dbbb08cb40d3
@bind ϵ Slider(0:0.01:10, default=1)

# ╔═╡ 7702f944-3576-11eb-08e4-0f1f61b1d3bd
Q = gen_sinkhorn(M, a, b; ϵ, α, β)

# ╔═╡ c80cf7c2-3576-11eb-0e81-05121d909271
sum(Q)

# ╔═╡ d0a700bc-3576-11eb-2d6c-2354903c6a7d
heatmap(Q)

# ╔═╡ 84a3ccea-357b-11eb-0924-b9f6fe8c481f
ã = sum(Q, dims=2)[:]

# ╔═╡ 8fd9f8dc-357b-11eb-0fa8-812f0ea106e5
b̃ = sum(Q, dims=1)[:]

# ╔═╡ 7c05eabe-357b-11eb-0f67-4fa981d4d725
plot(bar(ã), bar(b̃), layout=(2,1))

# ╔═╡ Cell order:
# ╠═16df5992-3574-11eb-3682-0981f8882895
# ╠═9e3f6282-3575-11eb-35e8-69e4cf57b6f4
# ╠═244cb164-3575-11eb-3b9d-837acaa0e9e9
# ╠═62aa2bf2-3576-11eb-0cc6-e31f7ad771d0
# ╠═3486c424-3576-11eb-35a3-0bcc30b65445
# ╠═ac95ae3a-357b-11eb-1e16-b56d76851314
# ╠═5ecc17ca-3576-11eb-00ce-2fe7c5d8f55b
# ╠═72fc2c12-3576-11eb-23d1-c1d14da6de89
# ╠═3d0844ac-3577-11eb-3086-112eb6425f79
# ╠═d9ffdbac-3576-11eb-24f1-55d257632108
# ╠═f22bc95c-3576-11eb-2a14-9d92a4497f68
# ╠═cd231944-357b-11eb-24b6-dbbb08cb40d3
# ╠═7702f944-3576-11eb-08e4-0f1f61b1d3bd
# ╠═c80cf7c2-3576-11eb-0e81-05121d909271
# ╠═d0a700bc-3576-11eb-2d6c-2354903c6a7d
# ╠═84a3ccea-357b-11eb-0924-b9f6fe8c481f
# ╠═8fd9f8dc-357b-11eb-0fa8-812f0ea106e5
# ╠═7c05eabe-357b-11eb-0f67-4fa981d4d725

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
using ReservoirComputing

#model parameters
const approx_res_size = 30
const radius = 1.2
const activation = tanh
const degree = 6
const sigma = 0.1
const beta = 0.0
const alpha = 1.0
const nla_type = NLADefault()
const in_size = 3
const out_size = 3
W_in = ReservoirComputing.init_input_layer(approx_res_size, in_size, sigma)
W = ReservoirComputing.init_reservoir_givendeg(approx_res_size, radius, degree)
const extended_states = false
const h_steps = 2

const train_len = 50
const predict_len = 12
data = ones(Float64, in_size, 100)
train = data[:, 1:1+train_len-1]
test = data[:, train_len:train_len+predict_len-1]

#constructor 1
esn = ESN(approx_res_size,
    train,
    degree,
    radius,
    activation = activation,
    sigma = sigma,
    alpha = alpha,
    nla_type = nla_type,
    extended_states = extended_states)

#test constructor
@test isequal(Integer(floor(approx_res_size/in_size)*in_size), esn.res_size)
@test isequal(train, esn.train_data)
@test isequal(alpha, esn.alpha)
@test isequal(activation, esn.activation)
@test isequal(nla_type, esn.nla_type)
@test size(esn.W) == (esn.res_size, esn.res_size)
@test size(esn.W_in) == (esn.res_size, esn.in_size)
@test size(esn.states) == (esn.res_size, train_len)

#test train
W_out = ESNtrain(esn, beta)
@test size(W_out) == (out_size, esn.res_size)
#test predict
output = ESNpredict(esn, predict_len, W_out)
@test size(output) == (out_size, predict_len)
#test h steps predict
output = ESNpredict_h_steps(esn, predict_len, h_steps, test, W_out)
@test size(output) == (out_size, predict_len)

#constructor 2
esn = ESN(W,
    train,
    activation = activation,
    sigma = sigma,
    alpha = alpha,
    nla_type = nla_type,
    extended_states = extended_states)

#test constructor
@test isequal(approx_res_size, esn.res_size)
@test isequal(train, esn.train_data)
@test isequal(alpha, esn.alpha)
@test isequal(activation, esn.activation)
@test isequal(nla_type, esn.nla_type)
@test size(esn.W) == (esn.res_size, esn.res_size)
@test size(esn.W_in) == (esn.res_size, esn.in_size)
@test size(esn.states) == (esn.res_size, train_len)

#test train
W_out = ESNtrain(esn, beta)
@test size(W_out) == (out_size, esn.res_size)
#test predict
output = ESNpredict(esn, predict_len, W_out)
@test size(output) == (out_size, predict_len)

#constructor 3
esn = ESN(approx_res_size,
    train,
    degree,
    radius,
    W_in,
    activation = activation,
    alpha = alpha,
    nla_type = nla_type,
    extended_states = extended_states)

#test constructor
@test isequal(Integer(floor(approx_res_size/in_size)*in_size), esn.res_size)
@test isequal(train, esn.train_data)
@test isequal(alpha, esn.alpha)
@test isequal(activation, esn.activation)
@test isequal(nla_type, esn.nla_type)
@test size(esn.W) == (esn.res_size, esn.res_size)
@test size(esn.W_in) == (esn.res_size, esn.in_size)
@test size(esn.states) == (esn.res_size, train_len)

#test train
W_out = ESNtrain(esn, beta)
@test size(W_out) == (out_size, esn.res_size)
#test predict
output = ESNpredict(esn, predict_len, W_out)
@test size(output) == (out_size, predict_len)

#constructor 4
esn = ESN(W,
    train,
    W_in,
    activation = activation,
    alpha = alpha,
    nla_type = nla_type,
    extended_states = extended_states)

#test constructor
@test isequal(approx_res_size, esn.res_size)
@test isequal(train, esn.train_data)
@test isequal(alpha, esn.alpha)
@test isequal(activation, esn.activation)
@test isequal(nla_type, esn.nla_type)
@test size(esn.W) == (esn.res_size, esn.res_size)
@test size(esn.W_in) == (esn.res_size, esn.in_size)
@test size(esn.states) == (esn.res_size, train_len)

#test train
W_out = ESNtrain(esn, beta)
@test size(W_out) == (out_size, esn.res_size)
#test predict
output = ESNpredict(esn, predict_len, W_out)
@test size(output) == (out_size, predict_len)


text:
module GetSupporting

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export mid_get_supporting

function mid_get_supporting(m)

    @expression(m, sum_of_boxes_below_box_i[i=1:length(boxes)],
        sum([m[:if_on_top][i, j] for j in 1:length(boxes) if i!=j])
    )

    @variable(m, two_or_more_below[i=1:length(boxes)], Bin)
    @constraint(m, two_or_more_below_c1[i=1:length(boxes)],
        m[:sum_of_boxes_below_box_i][i] + length(boxes) * (1 - m[:two_or_more_below][i]) >= 2
    )

    #min(flbx_front_left_bot_x[j]) * if_j_below_i
    #min(flby_front_left_bot_y[j]) * if_j_below_i

    @variable(m, min_flbx_of_supporting_ij[i=1:length(boxes), j=1:length(boxes); i!=j] >= 0)
    @constraint(m, get_min_flbx_of_supporting_ij[i=1:length(boxes), j=1:length(boxes); i!=j], (
        m[:min_flbx_of_supporting_ij][i, j] .>=
        [m[:if_on_top][i, j] * m[:flbx_front_left_bot_x][j], (1 - m[:if_on_top][i, j]) * m[:flbx_front_left_bot_x][i]]
        )
    )

    @variable(m, min_flby_of_supporting_ij[i=1:length(boxes), j=1:length(boxes); i!=j] >= 0)
    @constraint(m, get_min_flby_of_supporting_ij[i=1:length(boxes), j=1:length(boxes); i!=j], (
        m[:min_flby_of_supporting_ij][i, j] .>=
        [m[:if_on_top][i, j] * m[:flby_front_left_bot_y][j], (1 - m[:if_on_top][i, j]) * m[:flby_front_left_bot_y][i]]
        )
    )

    #max(flbx_front_left_bot_x[j] + box_x_axis_cover[j]) * if_j_below_i
    #max(flby_front_left_bot_y[j] + box_y_axis_cover[j]) * if_j_below_i

    @variable(m, max_flbx_of_supporting_ij[i=1:length(boxes), j=1:length(boxes); i!=j] >= 0)
    @constraint(m, get_max_flbx_of_supporting_ij[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:max_flbx_of_supporting_ij][i, j] .>=
        m[:if_on_top][i, j] * (m[:flbx_front_left_bot_x][j] + m[:box_x_axis_cover][j])
    )

    @variable(m, max_flby_of_supporting_ij[i=1:length(boxes), j=1:length(boxes); i!=j] >= 0)
    @constraint(m, get_max_flby_of_supporting_ij[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:max_flby_of_supporting_ij][i, j] .>=
        m[:if_on_top][i, j] * (m[:flby_front_left_bot_y][j] + m[:box_y_axis_cover][j])
    )

    #Change expression so that we have one value for each box
    @variable(m, min_flbx_of_supporting_i[i=1:length(boxes)] >= 0)
    @variable(m, min_flby_of_supporting_i[i=1:length(boxes)] >= 0)
    @variable(m, max_flbx_of_supporting_i[i=1:length(boxes)] >= 0)
    @variable(m, max_flby_of_supporting_i[i=1:length(boxes)] >= 0)

    @constraint(m, get_min_flbx_of_supporting_i[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:min_flbx_of_supporting_ij][i, j] >= m[:min_flbx_of_supporting_i][i]
    )
    @constraint(m, get_min_flby_of_supporting_i[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:min_flby_of_supporting_ij][i, j] >= m[:min_flby_of_supporting_i][i]
    )
    @constraint(m, get_max_flbx_of_supporting_i[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:max_flbx_of_supporting_i][i] >= m[:max_flbx_of_supporting_ij][i, j]
    )
    @constraint(m, get_max_flby_of_supporting_i[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:max_flby_of_supporting_i][i] >= m[:max_flby_of_supporting_ij][i, j]
    )

    #want to minimise this in the objective: minimise max, and maximise min
    @expression(m, ds_middle_minmax_objective[i=1:length(boxes), j=1:length(boxes); i!=j],
        10 * (1 - m[:two_or_more_below][i]) +
        m[:max_flbx_of_supporting_ij][i, j] +
        m[:max_flby_of_supporting_ij][i, j] +
        m[:min_flbx_of_supporting_ij][i, j] +
        m[:min_flby_of_supporting_ij][i, j] +
        m[:max_flbx_of_supporting_i][i] +
        m[:max_flby_of_supporting_i][i] +
        (L - m[:min_flbx_of_supporting_i][i]) +
        (W - m[:min_flby_of_supporting_i][i])
    )

    return m
end

end


text:
module NIfTI

using CodecZlib, Mmap, MappedArrays, TranscodingStreams

import Base.getindex, Base.size, Base.ndims, Base.length, Base.write, Base64
export NIVolume, niread, niwrite, voxel_size, time_step, vox, getaffine, setaffine

include("parsers.jl")
include("extensions.jl")
include("volume.jl")
include("headers.jl")

"""
    NIVolume{T<:Number,N,R} <: AbstractArray{T,N}
An `N`-dimensional NIfTI volume, with raw data of type
`R`. Note that if `R <: Number`, it will be converted to `Float32`. Additionally, the header is automatically
updated to be consistent with the raw volume.

# Members
- `header`: a `NIfTI1Header`
- `extensions`: a Vector of `NIfTIExtension`s
- `raw`: Raw data of type `R`
"""
struct NIVolume{T<:Number,N,R} <: AbstractArray{T,N}
    header::NIfTI1Header
    extensions::Vector{NIfTIExtension}
    raw::R
end
NIVolume(header::NIfTI1Header, extensions::Vector{NIfTIExtension}, raw::R) where {R}=
    niupdate(new(header, extensions, raw))

NIVolume(header::NIfTI1Header, extensions::Vector{NIfTIExtension}, raw::AbstractArray{T,N}) where {T<:Number,N} =
    NIVolume{typeof(one(T)*1f0+1f0),N,typeof(raw)}(header, extensions, raw)
NIVolume(header::NIfTI1Header, raw::AbstractArray{T,N}) where {T<:Number,N} =
    NIVolume{typeof(one(T)*1f0+1f0),N,typeof(raw)}(header, NIfTIExtension[], raw)
NIVolume(header::NIfTI1Header, extensions::Vector{NIfTIExtension}, raw::AbstractArray{Bool,N}) where {N} =
    NIVolume{Bool,N,typeof(raw)}(header, extensions, raw)
NIVolume(header::NIfTI1Header, raw::AbstractArray{Bool,N}) where {N} =
    NIVolume{Bool,N,typeof(raw)}(header, NIfTIExtension[], raw)

header(x::NIVolume) = getfield(x, :header)

include("coordinates.jl")

"""
    voxel_size(header::NIfTI1Header)

Get the voxel size **in mm** from a `NIfTI1Header`.
"""
voxel_size(header::NIfTI1Header) =
    [header.pixdim[i] * SPATIAL_UNIT_MULTIPLIERS[header.xyzt_units & Int8(3)] for i = 2:min(header.dim[1], 3)+1]

# Always in ms
"""
    time_step(header::NIfTI1Header)

Get the TR **in ms** from a `NIfTI1Header`.
"""
time_step(header::NIfTI1Header) =
    header.pixdim[5] * TIME_UNIT_MULTIPLIERS[header.xyzt_units >> 3]

function to_dim_info(dim_info::Tuple{Integer,Integer,Integer})
    if dim_info[1] > 3 || dim_info[1] < 0
        error("Invalid frequency dimension $(dim_info[1])")
    elseif dim_info[2] > 3 || dim_info[2] < 0
        error("Invalid phase dimension $(dim_info[2])")
    elseif dim_info[3] > 3 || dim_info[3] < 0
        error("Invalid slice dimension $(dim_info[3])")
    end

    return Int8(dim_info[1] | (dim_info[2] << 2) | (dim_info[3] << 4))
end

# Returns or sets dim_info as a tuple whose values are the frequency, phase, and slice dimensions
function dim_info(header::NIfTI1Header)
    return (
        header.dim_info & int8(3),
        (header.dim_info >> 2) & int8(3),
        (header.dim_info >> 4) & int8(3)
    )
end
function dim_info(header::NIfTI1Header, dim_info::Tuple{T, T, T}) where {T<:Integer}
    header.dim_info = to_dim_info(dim_info)
end

# Gets the size of a type in bits
nibitpix(t::Type) = Int16(sizeof(t)*8)
nibitpix(::Type{Bool}) = Int16(1)

# Constructor
function NIVolume(
    # Optional MRI volume; if not given, an empty volume is used
    raw::AbstractArray{T}=Int16[],
    extensions::Union{Vector{NIfTIExtension},Nothing}=nothing;

    # Fields specified as UNUSED in NIfTI1 spec
    data_type::AbstractString="", db_name::AbstractString="", extents::Integer=Int32(0),
    session_error::Integer=Int16(0), regular::Integer=Int8(0), glmax::Integer=Int32(0),
    glmin::Integer=Int16(0),

    # The frequency encoding, phase encoding and slice dimensions.
    dim_info::NTuple{3, Integer}=(0, 0, 0),
    # Describes data contained in the file; for valid values, see
    # http://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/group__NIfTI1__INTENT__CODES.html
    intent_p1::Real=0f0, intent_p2::Real=0f0, intent_p3::Real=0f0,
    intent_code::Integer=Int16(0), intent_name::AbstractString="",
    # Information about which slices were acquired, in case the volume has been padded
    slice_start::Integer=Int16(0), slice_end::Integer=Int16(0), slice_code::Int8=Int8(0),
    # The size of each voxel and the time step. These are formulated in mm unless xyzt_units is
    # explicitly specified
    voxel_size::NTuple{3, Real}=(1f0, 1f0, 1f0), time_step::Real=0f0, xyzt_units::Int8=Int8(18),
    # Slope and intercept by which volume shoudl be scaled
    scl_slope::Real=1f0, scl_inter::Real=0f0,
    # These describe how data should be scaled when displayed on the screen. They are probably
    # rarely used
    cal_max::Real=0f0, cal_min::Real=0f0,
    # The amount of time taken to acquire a slice
    slice_duration::Real=0f0,
    # Indicates a non-zero start point for time axis
    toffset::Real=0f0,

    # "Any text you like"
    descrip::AbstractString="",
    # Name of auxiliary file
    aux_file::AbstractString="",

    # Parameters for Method 2. See the NIfTI spec
    qfac::Float32=0f0, quatern_b::Real=0f0, quatern_c::Real=0f0, quatern_d::Real=0f0,
    qoffset_x::Real=0f0, qoffset_y::Real=0f0, qoffset_z::Real=0f0,
    # Orientation matrix for Method 3
    orientation::Union{Matrix{Float32},Nothing}=nothing) where {T <: Number}
    local t
    if isempty(raw)
        raw = Int16[]
        t = Int16
    else
        t = T
    end

    if extensions == nothing
        extensions = NIfTIExtension[]
    end

    method2 = qfac != 0 || quatern_b != 0 || quatern_c != 0 || quatern_d != 0 || qoffset_x != 0 ||
        qoffset_y != 0 || qoffset_z != 0
    method3 = orientation != nothing

    if method2 && method3
        error("Orientation parameters for Method 2 and Method 3 are mutually exclusive")
    end

    if method3
        if size(orientation) != (3, 4)
            error("Orientation matrix must be of dimensions (3, 4)")
        end
    else
        orientation = zeros(Float32, 3, 4)
    end

    if slice_start == 0 && slice_end == 0 && dim_info[3] != 0
        slice_start = 0
        slice_end = size(raw, dim_info[3]) - 1
    end

    NIVolume(NIfTI1Header(SIZEOF_HDR1, string_tuple(data_type, 10), string_tuple(db_name, 18), extents, session_error,
        regular, to_dim_info(dim_info), to_dim_i16(size(raw)), intent_p1, intent_p2,
        intent_p3, intent_code, eltype_to_int16(t), nibitpix(t),
        slice_start, (qfac, voxel_size..., time_step, 0, 0, 0), 352,
        scl_slope, scl_inter, slice_end, slice_code,
        xyzt_units, cal_max, cal_min, slice_duration,
        toffset, glmax, glmin, string_tuple(descrip, 80), string_tuple(aux_file, 24), (method2 || method3),
        method3, quatern_b, quatern_c, quatern_d,
        qoffset_x, qoffset_y, qoffset_z, (orientation[1, :]...,),
        (orientation[2, :]...,), (orientation[3, :]...,), string_tuple(intent_name, 16), NP1_MAGIC), extensions, raw)
end

# Validates the header of a volume and updates it to match the volume's contents
function niupdate(vol::NIVolume{T}) where {T}
    vol.header.sizeof_hdr = SIZEOF_HDR1
    vol.header.dim = to_dim_i16(size(vol.raw))
    vol.header.datatype = eltype_to_int16(T)
    vol.header.bitpix = nibitpix(T)
    vol.header.vox_offset = isempty(vol.extensions) ? Int32(352) :
        Int32(mapreduce(esize, +, vol.extensions) + SIZEOF_HDR1)
    vol
end

# Avoid method ambiguity
write(io::Base64.Base64EncodePipe, vol::NIVolume{UInt8,1}) = invoke(write, (IO, NIVolume{UInt8,1}), io, vol)

# Write a NIfTI file
function write(io::IO, vol::NIVolume)
    write(io, niupdate(vol).header)
    if isempty(vol.extensions)
        write(io, Int32(0))
    else
        for ex in vol.extensions
            sz = esize(ex)
            write(io, Int32(sz))
            write(io, Int32(ex.ecode))
            write(io, ex.edata)
            write(io, zeros(UInt8, sz - length(ex.edata)))
        end
    end
    if eltype(vol.raw) == Bool
        write(io, BitArray(vol.raw))
    else
        write(io, vol.raw)
    end
end

"""
    niwrite(path::AbstractString, vol::NIVolume)

Write a NIVolume to a file specified by `path`.
"""
function niwrite(path::AbstractString, vol::NIVolume)
    if split(path,".")[end] == "gz"
        io = open(path, "w")
        stream = GzipCompressorStream(io)
        write(stream, vol)
        close(stream)
        close(io)
    else
        io = open(path, "w")
        write(io, vol)
        close(io)
    end
end

# Read header from a NIfTI file
function read_header(io::IO)
    header, swapped = read(io, NIfTI1Header)
    if header.sizeof_hdr != SIZEOF_HDR1
        error("This is not a NIfTI-1 file")
    end
    header, swapped
end

# Look for a gzip header in an IOStream
function isgz(io::IO)
    try
        ret = read(io, UInt8) == 0x1F && read(io, UInt8) == 0x8B
        seekstart(io)
        ret
    catch err
        if isa(err, EOFError)
            @debug "reading the file resulted in an EOF error and \nthe end of the file was read.\nNo more data was available to read from the filestream.\nIt is likely that the file was corrupted or is empty (0 bytes)."
            rethrow(err)
        end
    end
end

"""
    niread(file; mmap=false, mode="r")

Read a NIfTI file to a NIVolume. Set `mmap=true` to memory map the volume.
"""
function niread(file::AbstractString; mmap::Bool=false, mode::AbstractString="r")
    io = niopen(file, mode)
    hdr, swapped = read_header(io)
    ex = read_extensions(io, hdr.vox_offset - 352)

    if hdr.magic === NP1_MAGIC
        vol = read_volume(io, to_eltype(hdr.datatype), to_dimensions(hdr.dim), mmap)
    else
        vol = read_volume(niopen(hdr_to_img(file), mode), to_eltype(hdr.datatype), to_dimensions(hdr.dim), mmap)
    end

    if swapped && sizeof(eltype(vol)) > 1
        return NIVolume(hdr, ex, mappedarray(ntoh, hton, vol))
    else
        return NIVolume(hdr, ex, vol)
    end
end

# Allow file to be indexed like an array, but with indices yielding scaled data
@inline getindex(f::NIVolume{T}, idx::Vararg{Int}) where {T} =
    getindex(f.raw, idx...,) * f.header.scl_slope + f.header.scl_inter

add1(x::Union{AbstractArray{T},T}) where {T<:Integer} = x + 1
add1(::Colon) = Colon()

"""
    vox(f::NIVolume, args...,)

Get the value of a voxel from volume `f`, scaled by slope and intercept given in header, with 0-based indexing.
`args` are the voxel indices and the length should be the number of dimensions in `f`.
"""
@inline vox(f::NIVolume, args...,) = getindex(f, map(add1, args)...,)

size(f::NIVolume) = size(f.raw)
size(f::NIVolume, d) = size(f.raw, d)
ndims(f::NIVolume) = ndims(f.raw)
length(f::NIVolume) = length(f.raw)
lastindex(f::NIVolume) = lastindex(f.raw)

end


text:
using StrideArrays, LinearAlgebra, Aqua, ArrayInterface
using Test

import InteractiveUtils
InteractiveUtils.versioninfo(stdout; verbose = true)

@show StrideArrays.VectorizationBase.register_count()
const START_TIME = time()

@time @testset "StrideArrays.jl" begin
  @time Aqua.test_all(StrideArrays, ambiguities = false, project_toml_formatting = false)
  # Currently, there is one method ambiguity:
  # - map(f::F, A::AbstractStrideArray, args::Vararg{Any, K}) where {F, K} in StrideArrays at StrideArrays/src/miscellaneous.jl:22
  # - map(f, a1::AbstractArray, a2::StaticArrays.StaticArray, as::AbstractArray...) in StaticArrays at StaticArrays/0bweZ/src/mapreduce.jl:33
  @time @test length(Test.detect_ambiguities(StrideArrays)) <= 1
  # @test isempty(detect_unbound_args(StrideArrays))
  @time include("matmul_tests.jl")
  @time include("misc.jl")
  @time include("broadcast_tests.jl")
end

const ELAPSED_MINUTES = (time() - START_TIME) / 60
@test ELAPSED_MINUTES < 180


text:
struct Equation{DT <: Number, FT <: Function, FPRE <: Union{Function,Nothing}, FPOST <: Union{Function,Nothing}}
    f::FT
    x₀::Array{DT,2}
    f_pre::FPRE
    f_post::FPOST

    function Equation(f::FT, x₀::AbstractArray{DT,2}; f_preproc::FPRE=nothing, f_postproc::FPOST=nothing) where {DT, FT, FPRE, FPOST}
        new{DT, FT, FPRE, FPOST}(f, Array(x₀), f_preproc, f_postproc)
    end
end

ndims(equ::Equation) = length(axes(equ.x₀,1))
nsamples(equ::Equation) = length(axes(equ.x₀,2))

haspreprocessing(::Equation{<:Number, <:Function, Nothing}) = false
haspreprocessing(::Equation{<:Number, <:Function, <:Function}) = true

haspostprocessing(::Equation{<:Number, <:Function, <:Union{Function,Nothing}, Nothing}) = false
haspostprocessing(::Equation{<:Number, <:Function, <:Union{Function,Nothing}, <:Function}) = true

function preprocessing(equ::Equation, x₀::AbstractArray)
    if haspreprocessing(equ)
        equ.f_pre(x₀)
    end
end

function postprocessing(equ::Equation, x₁::AbstractArray)
    if haspostprocessing(equ)
        equ.f_post(x₁)
    end
end

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
# Handers for execute_request and related messages, which are
# the core of the Jupyter protocol: execution of Julia code and
# returning results.

import Base.Libc: flush_cstdio
import Pkg

Base.showable(a::AbstractVector{<:MIME}, x) = any(m -> showable(m, x), a)

"""
A vector of MIME types (or vectors of MIME types) that IJulia will try to
render. IJulia will try to render every MIME type specified in the first level
of the vector. If a vector of MIME types is specified, IJulia will include only
the first MIME type that is renderable (this allows for the expression of
priority and exclusion of redundant data).

For example, since "text/plain" is specified as a first-child of the array,
IJulia will always try to include a "text/plain" representation of anything that
is displayed. Since markdown and latex are specified within a sub-vector, IJulia
will always try to render "text/markdown", and will only try to render
"text/latex" if markdown isn't possible.
"""
const ijulia_mime_types = Vector{Union{MIME, AbstractVector{MIME}}}([
    MIME("text/plain"),
    MIME("image/svg+xml"),
    [MIME("image/png"),MIME("image/jpeg")],
    [
        MIME("text/markdown"),
        MIME("text/html"),
        MIME("text/latex"), # Jupyter expects this
        MIME("application/x-latex"), # but this is more standard?
    ],
])

"""
MIME types that when rendered (via stringmime) return JSON data. See
`ijulia_mime_types` for a description of how MIME types are selected.

This is necessary to embed the JSON as is in the displaydata bundle (rather than
as stringify'd JSON).
"""
const ijulia_jsonmime_types = Vector{Union{MIME, Vector{MIME}}}([
    [MIME("application/vnd.vegalite.v2+json"), MIME("application/vnd.vega.v3+json")],
    MIME("application/vnd.dataresource+json"),
])

register_mime(x::Union{MIME, Vector{MIME}})= push!(ijulia_mime_types, x)
register_mime(x::AbstractVector{<:MIME}) = push!(ijulia_mime_types, Vector{Mime}(x))
register_jsonmime(x::Union{MIME, Vector{MIME}}) = push!(ijulia_jsonmime_types, x)
register_jsonmime(x::AbstractVector{<:MIME}) = push!(ijulia_jsonmime_types, Vector{Mime}(x))

include("magics.jl")

# return a String=>Any dictionary to attach as metadata
# in Jupyter display_data and pyout messages
metadata(x) = Dict()

"""
Generate the preferred MIME representation of x.

Returns a tuple with the selected MIME type and the representation of the data
using that MIME type.
"""
function display_mimestring(mime_array::Vector{MIME}, x)
    for m in mime_array
        if showable(m, x)
            return display_mimestring(m, x)
        end
    end
    error("No displayable MIME types in mime array.")
end

display_mimestring(m::MIME, x) = (m, limitstringmime(m, x))

"""
Generate the preferred json-MIME representation of x.

Returns a tuple with the selected MIME type and the representation of the data
using that MIME type (as a `JSONText`).
"""
function display_mimejson(mime_array::Vector{MIME}, x)
    for m in mime_array
        if showable(m, x)
            return display_mimejson(m, x)
        end
    end
    error("No displayable MIME types in mime array.")
end

display_mimejson(m::MIME, x) = (m, JSON.JSONText(limitstringmime(m, x)))

"""
Generate a dictionary of `mime_type => data` pairs for all registered MIME
types. This is the format that Jupyter expects in display_data and
execute_result messages.
"""
function display_dict(x)
    data = Dict{String, Union{String, JSONText}}()
    for m in ijulia_mime_types
        if showable(m, x)
            mime, mime_repr = display_mimestring(m, x)
            data[string(mime)] = mime_repr
        end
    end

    for m in ijulia_jsonmime_types
        if showable(m, x)
            mime, mime_repr = display_mimejson(m, x)
            data[string(mime)] = mime_repr
        end
    end

    return data

end

# queue of objects to display at end of cell execution
const displayqueue = Any[]

# remove x from the display queue
function undisplay(x)
    i = findfirst(isequal(x), displayqueue)
    i !== nothing && splice!(displayqueue, i)
    return x
end

import Base: ip_matches_func

function show_bt(io::IO, top_func::Symbol, t, set)
    # follow PR #17570 code in removing top_func from backtrace
    eval_ind = findlast(addr->ip_matches_func(addr, top_func), t)
    eval_ind !== nothing && (t = t[1:eval_ind-1])
    Base.show_backtrace(io, t)
end

# wrapper for showerror(..., backtrace=false) since invokelatest
# doesn't support keyword arguments.
showerror_nobt(io, e, bt) = showerror(io, e, bt, backtrace=false)

# return the content of a pyerr message for exception e
function error_content(e, bt=catch_backtrace();
                       backtrace_top::Symbol=SOFTSCOPE[] ? :softscope_include_string : :include_string,
                       msg::AbstractString="")
    tb = map(x->String(x), split(sprint(show_bt,
                                        backtrace_top,
                                        bt, 1:typemax(Int)),
                                  "\n", keepempty=true))

    ename = string(typeof(e))
    evalue = try
        # Peel away one LoadError layer that comes from running include_string on the cell
        isa(e, LoadError) && (e = e.error)
        sprint((io, e, bt) -> invokelatest(showerror_nobt, io, e, bt), e, bt)
    catch
        "SYSTEM: show(lasterr) caused an error"
    end
    pushfirst!(tb, evalue) # fperez says this needs to be in traceback too
    if !isempty(msg)
        pushfirst!(tb, msg)
    end
    Dict("ename" => ename, "evalue" => evalue,
                 "traceback" => tb)
end

#######################################################################

# global variable so that display can be done in the correct Msg context
execute_msg = Msg(["julia"], Dict("username"=>"jlkernel", "session"=>uuid4()), Dict())
# global variable tracking the number of bytes written in the current execution
# request
const stdio_bytes = Ref(0)

import REPL: helpmode

# use a global array to accumulate "payloads" for the execute_reply message
const execute_payloads = Dict[]

const stdout_name = isdefined(Base, :stdout) ? "stdout" : "STDOUT"

function execute_request(socket, msg)
    code = msg.content["code"]
    @vprintln("EXECUTING ", code)
    global execute_msg = msg
    global n, In, Out, ans
    stdio_bytes[] = 0
    silent = msg.content["silent"]
    store_history = get(msg.content, "store_history", !silent)
    empty!(execute_payloads)

    if !silent
        n += 1
        send_ipython(publish[],
                     msg_pub(msg, "execute_input",
                             Dict("execution_count" => n,
                                          "code" => code)))
    end

    silent = silent || REPL.ends_with_semicolon(code)
    if store_history
        In[n] = code
    end

    # "; ..." cells are interpreted as shell commands for run
    code = replace(code, r"^\s*;.*$" =>
                   m -> string(replace(m, r"^\s*;" => "Base.repl_cmd(`"),
                               "`, ", stdout_name, ")"))

    # "] ..." cells are interpreted as pkg shell commands
    if occursin(r"^\].*$", code)
        code = "IJulia.Pkg.REPLMode.do_cmd(IJulia.minirepl[], \"" *
            escape_string(code[2:end]) * "\"; do_rethrow=true)"
    end

    # a cell beginning with "? ..." is interpreted as a help request
    hcode = replace(code, r"^\s*\?" => "")

    try
        for hook in preexecute_hooks
            invokelatest(hook)
        end

        ans = result = if hcode != code # help request
            Core.eval(Main, helpmode(hcode))
        else
            #run the code!
            occursin(magics_regex, code) ? magics_help(code) :
                SOFTSCOPE[] ? softscope_include_string(current_module[], code, "In[$n]") :
                include_string(current_module[], code, "In[$n]")
        end

        if silent
            result = nothing
        elseif result !== nothing
            if store_history
                Out[n] = result
            end
        end

        user_expressions = Dict()
        for (v,ex) in msg.content["user_expressions"]
            user_expressions[v] = invokelatest(parse, ex)
        end

        for hook in postexecute_hooks
            invokelatest(hook)
        end

        # flush pending stdio
        flush_all()

        undisplay(result) # dequeue if needed, since we display result in pyout
        invokelatest(display) # flush pending display requests

        if result !== nothing
            result_metadata = invokelatest(metadata, result)
            result_data = invokelatest(display_dict, result)
            send_ipython(publish[],
                         msg_pub(msg, "execute_result",
                                 Dict("execution_count" => n,
                                              "metadata" => result_metadata,
                                              "data" => result_data)))

        end
        send_ipython(requests[],
                     msg_reply(msg, "execute_reply",
                               Dict("status" => "ok",
                                            "payload" => execute_payloads,
                                            "execution_count" => n,
                                            "user_expressions" => user_expressions)))
        empty!(execute_payloads)
    catch e
        bt = catch_backtrace()
        try
            # flush pending stdio
            flush_all()
            for hook in posterror_hooks
                invokelatest(hook)
            end
        catch
        end
        empty!(displayqueue) # discard pending display requests on an error
        content = error_content(e,bt)
        send_ipython(publish[], msg_pub(msg, "error", content))
        content["status"] = "error"
        content["execution_count"] = n
        send_ipython(requests[], msg_reply(msg, "execute_reply", content))
    end
end


text:
"""
RANSAC algorithm used when there are more points than algorithm requires.

# Arguments:

- `sample_selector`: Function that given set of ids of the size `n_samples`,
    returns input for the `kernel` method.
- `kernel`: Function that calculates candidates from the sampled data.
- `rank`: Function that evaluates candidates, computed by `kernel`.
    It should return `n_inliers, model`, where `model` can be of any type.
    And is not used by the ransac itself. It only uses `n_inliers` to select
    which model is better.
- `n_points`: Total number of points in the data.
- `n_samples`: Number of points, that `kernel` function accepts.
    E.g. `five_point_ransac` sets this value to `5`.
- `iterations`: Maximum number of iterations to attempt. Default is `100`.
- `confidence`: Confidence in `[0, 1]` range. The higher the value, the more
    certain algorithm is that the picked model does not contain outliers.
    Number of iterations performed by RANSAC depends on this value as well.

# Returns:

Number of inliers and the selected model in the format returned by `kernel`.
"""
function ransac(
    sample_selector, kernel, rank, n_points, n_samples;
    iterations = 100, confidence = 0.99,
)
    best_M = nothing
    best_n_inliers = 0
    n = iterations

    ϵ = eps()
    current_iteration = 0
    while n > current_iteration
        ids = randperm(n_points)[1:n_samples]
        M = kernel(sample_selector(ids)...)
        n_inliers, M = rank(M; sample_ids=ids)

        if n_inliers > best_n_inliers
            best_M = M
            best_n_inliers = n_inliers
            # Update `n` estimate, which is the number of trials to ensure
            # a dataset with no outliers is picked with given `confidence`.
            p_no_outliers = 1.0 - (n_inliers / n_points) ^ n_samples
            p_no_outliers = clamp(p_no_outliers, ϵ, 1.0 - ϵ)
            n = log(1.0 - confidence) / log(p_no_outliers)
        end

        current_iteration += 1
        current_iteration > iterations && break
    end

    best_M ≡ nothing && return nothing
    best_n_inliers, best_M
end


text:
module ASTIR

include("utils.jl")
include("frontend.jl")
include("unbox.jl")
include("controlflow.jl")
include("dummy.jl")

end


text:
### A Pluto.jl notebook ###
# v0.14.4

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 3d4d751a-b972-4e88-b568-1f19f0f20564
begin
        const DIR = "$(ENV["HOME"])/Projects/ML_DL/Notebooks/julia-notebooks/Julia_Investigation/"

        using Pkg
        Pkg.activate(DIR)
end

# ╔═╡ e95db99e-a9ec-437c-9557-4dec4f42926f
begin
        using PlutoUI
        using CSV, DataFrames, TextAnalysis
        using Pipe: @pipe
        using Plots
end

# ╔═╡ 35b8ad6d-7c5e-43bf-b014-ce24010a959a
using TextAnalysis: NaiveBayesClassifier, fit!, predict

# ╔═╡ cbad2e04-9d92-11eb-3d84-9da6f9a3c818
md"""

## Sentiment analysis with Julia

Based on [HoJ - Tom Kwong](https://github.com/Humans-of-Julia/Challenges/tree/main/Week2-TextAnalysis.jl/tk3369) and the presentation [Sentiment Analysis in Julia with TextAnalysis.jl by Tom Kwong](https://www.youtube.com/watch?v=vZr1SOxjDms)

This notebook explores the [Twitter US Airline Sentiment](https://www.kaggle.com/crowdflower/twitter-airline-sentiment) data set from Kaggle.

We will use [TextAnalysis.jl](https://github.com/JuliaText/TextAnalysis.jl) as the primary tool for analyzing textual data.

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ 26957ae3-1a74-4be7-b96e-d3292991a82a
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ 225ebbac-2232-47ff-b7dc-e2c59814f196
md"""
### Loading the Data
"""

# ╔═╡ 3ab74396-3e42-4545-af0b-f5746345324d
begin
        cd(DIR)
        df = CSV.File("./data/Tweets.csv") |> DataFrame
end

# ╔═╡ 13a43945-1a9c-4526-9f63-03a306b8d459
md"""
### Data Wrangling

Let us take at look at the data to get some understanding about what is going on.
"""

# ╔═╡ 7e812f27-4129-48ba-b4f5-a49ddc4169ca
describe(df, :eltype, :nmissing, first => :first, last => :last)

# ╔═╡ 3a0b11b9-a51e-47f1-833c-8c9a2e4d73ce
md"""
It looks like the column `:airline_sentiment` is our label.
"""

# ╔═╡ 5f3e9cc0-fb07-4f2d-be7c-590f13bbd03d
let gp = combine(groupby(df, :airline_sentiment), nrow)
        bar(gp.airline_sentiment, gp.nrow,
                title="Airline Sentiment",
                label="",
                # legend=:topleft
        )
end

# ╔═╡ 1c0d74d4-aac6-4fa4-8957-84457c73107e
md"""
OK, so this dataset is unbalanced, there are more negative sentiment than positivie one (maybe not surprising - we write to complain about...)
"""

# ╔═╡ 668b9aac-454a-43a6-8c40-44a0c023129a
histogram(df.airline_sentiment_confidence;
        legend=nothing, title="Airline Sentiment Confidence")

# ╔═╡ 23c79d73-3431-44c3-a4ef-4dd8337f7931
let gp = groupby(dropmissing(df, :negativereason), :negativereason) |>
                df_ -> combine(df_, nrow)
        bar(gp.negativereason, gp.nrow;
                title="Negative Reasons", label=nothing, xrotation=45)
end

# ╔═╡ ca585d7b-d7f1-48f6-ac88-10aa582b6dbe
md"""
The CSV file contains over 14,000 tweets. Let's quickly examine some individual data.

Before we go further, it would be nice to display a single record in table format. We can define a `table` function that converts an indexable object into Markdown format, which can be displayed in this Pluto notebook.
"""

# ╔═╡ 1289f7cf-6c2a-4b67-a817-3de3926d46ec
function table(nt)
        io = IOBuffer()   ## the link with the bind macro
        println(io, "|name|value|")
        println(io, "|---:|:----|")
        for k ∈ keys(nt)
                println(io, "|`", k, "`|", nt[k], "|")
        end
        return Markdown.parse(String(take!(io)))
end;

# ╔═╡ dbc871c5-9ca9-432c-98f8-af386ebdb22e
md"""
Let us define a variable called `row` and bind it to a slider for quick experimentation.
"""

# ╔═╡ e4c7dc20-f703-4d56-9835-af4fdcb32668
@bind row html"""
<input type="range" min="1" max ="100" value="36" />
"""

# ╔═╡ 6d7d638a-93a5-4c40-8d60-899ccafb0cdc
"Current Record: $row"

# ╔═╡ 4a649215-b35b-4e5c-96ab-ee9808d2de8e
table(df[row, :])

# ╔═╡ 9e89132b-ae35-4d85-ae98-15b1ea326991
md"""
As an example, record #36 has the tweet text as:
```
Nice RT @VirginAmerica: Vibe with the moodlight from takeoff to touchdown. #MoodlitMonday #ScienceBehindTheExperience http://t.co/Y7O0uNxTQP
```
This is a tricky one because it contains all of the followings:
- mention (`@VirginAmeria`)
- hash tag (`#MoodlitMonday` and `#ScienceBehindTheExperience`)
- URL (`http://t.co/Y7O0uNxTQP`)
Technically `RT` is a shorthand for "retweet" so perhaps it should be expanded but let's not worry about that for now.
"""

# ╔═╡ de978888-c378-46e9-8158-2cdbc1a6fe2f
md"""
### Handling mentions and hastags

If we just ignore these problems then it can be a disaster.
"""

# ╔═╡ e31a7663-6705-450e-9100-346b8fe51e15
let s = df[36, :text]
        sd = StringDocument(lowercase(s))
        op = 0x00
        op |= strip_punctuation
        op |= strip_stopwords
        op |= strip_html_tags
        prepare!(sd, op)
        stem!(sd)
        table(ngrams(sd))
end

# ╔═╡ 8db4042d-fba1-4bb7-ba43-0acc718e1389
md"""
We can see some problems here. It seems that when we stripped punctuations, it also took the `@` and `#` signs away. The URL also became weird. OK, that is what stripping punctuation means...
"""

# ╔═╡ 5fc69267-810a-479c-acdf-6026ca901f64
md"""
#### Extracting mentions, hash tags, and URL's.

This neat idea came from José Bayoán Santiago Calderón when I asked the question on Slack. Let's define some functions using regular expressions.
"""

# ╔═╡ 50e09727-5fb9-4f02-9a71-f2fa65eb19b9
const Regexp = Dict{Symbol, Regex}(
        :mention => r"@\w+",
        :hashtag => r"#\w+",
        :url => r"http[s]?://(?:[a-zA-Z]|[0-9]|[$/\-_@\.&\+]|[!\*\(\),]|(?:%[0-9a-fA-F]{2}))+",
)

# ╔═╡ 9ea51080-e3f8-4d70-a139-aa29ddf6d13b
extract_tokens(s::AbstractString, token_type::Symbol) =
        collect(x.match for x in eachmatch(Regexp[token_type], s))

# ╔═╡ d55d2b3c-82cf-4328-9aee-279e839d0a92
function remove_tokens(s::AbstractString)
        for re ∈ values(Regexp)
                s = replace(s, re => "")
        end
        s
end

# ╔═╡ 69c18f2f-7e22-4188-9464-9705f4aed7e3
remove_tokens("Nice RT @VirginAmerica: Vibe with the moodlight from takeoff to touchdown. #MoodlitMonday #ScienceBehindTheExperience http://t.co/Y7O0uNxTQP")

# ╔═╡ e4186b6a-2ee3-4fd7-b1a3-06a499aeaf07
values(Regexp) .=> ""

# ╔═╡ 7f5d2d98-e170-46ff-a470-e7708cbe3fb3
md"""
#### Create a new DataFrame with extrated and clean textfields
"""

# ╔═╡ 10a4ddbc-66c4-472d-999e-99838da6652d
begin
        d = Dict(
                :airline_sentiment => df.airline_sentiment,
                :text => df.text,
                :mentions => extract_tokens.(df.text, :mention),
                :hashtags => extract_tokens.(df.text, :hashtag),
                :urls => extract_tokens.(df.text, :url),
                :clean_text => remove_tokens.(df.text) .|> strip .|> lowercase
        )
        df₂ = DataFrame(d);
end

# ╔═╡ a65e9da9-0e8b-4356-820f-266eaba857c8
table(df₂[36, :])

# ╔═╡ 09441471-d246-4511-8d8e-e941b1166159
md"""
As we can see, the mentions/hashtags/urls are extracted into separate columns in the DataFrame. The `clean_text` field contains the cleaned version of `text`.
"""

# ╔═╡ 5855ee27-4d34-41e4-9881-70ad4b5f064e
md"""
### Using Naive Bayes Classifier

"""

# ╔═╡ 819af8d7-c242-4f7b-a61f-121263bdb54c
md"""
In our DataFrame, we already have a column `x_string_doc` with `StringDocuments` values. So we can just fit them to the classifier.
"""

# ╔═╡ 41a18e5c-7e01-4828-9df0-53d446da9eb1
function create_string_doc(s::AbstractString)::StringDocument
        sd = StringDocument(s)
        op = 0x00
        op |= strip_punctuation
        op |= strip_stopwords
        op |= strip_html_tags
        prepare!(sd, op)
        stem!(sd)
        sd
end

# ╔═╡ 3ac25d84-b396-4341-956b-c550054038c0
model = let
        labels = df₂.airline_sentiment .|> strip |>
                unique ## ≡ labels = ["neutral", "positive", "negative"]
        nbc = NaiveBayesClassifier(labels)

        for (clean_txt, y) ∈ zip(df₂.clean_text, df₂.airline_sentiment)
                x = create_string_doc(clean_txt)
                fit!(nbc, x, y |> strip)
        end

        nbc
end;

# ╔═╡ cd495081-805b-4e71-814a-10add9a65d03
typeof(model)

# ╔═╡ bcae0e95-22fc-4526-9c71-a186e20d5d8d
typeof(NaiveBayesClassifier)

# ╔═╡ 032ddbf5-5995-4ff9-be43-3cfa16eaac91
function test_model(model, tweets::Vector{<: AbstractString})
        df = DataFrame(text=tweets)
        df.doc = remove_tokens.(tweets) .|>
                create_string_doc .|>
                TextAnalysis.text
        # TextAnalysis.text.(create_string_doc.(remove_tokens.(tweets)))
        df.analysis = predict.(Ref(model), df.doc)

        df.positive = getindex.(df.analysis, "positive")
        df.negative = getindex.(df.analysis, "negative")
        df.neutral = getindex.(df.analysis, "neutral")

        select!(df, Not(:analysis))
        df
end

# ╔═╡ 6887b10b-0c44-4a1a-bbe0-e08032ac1dff
let
        tweets = [
                "whatever airline sucks!",
                "i love @american service :-)",
                "just ok",
                "hello world"]
        test_model(model, tweets)
end

# ╔═╡ 0c7ad5a2-e8b1-4189-82f5-d39cbc122d2d
md"""
#### Determining accuracy

OK, so how well does the Naive Bayes Classifier perform?

As the `predict` fucntion returns a `Dict` object with the probabilities assigned to each class, we need to choose the best one. So let us write a function for that:
"""

# ╔═╡ 2c853e55-843c-43c8-87b4-9ef6810dc03d
function best_pred(model::NaiveBayesClassifier, sd::StringDocument)
        predict(model, sd) |> argmax
end

# ╔═╡ 2484eca4-198c-4178-8b78-26ee246f786b
md"""
Now le us make predictions over all 14K tweets.
"""

# ╔═╡ 902b8634-7cde-4808-9e77-c13caac19c8f
ŷ = let x = remove_tokens.(df₂.text) .|> lowercase .|> create_string_doc
        best_pred.(Ref(model), x)
end;

# ╔═╡ 53c34fd9-f9c4-4321-94f3-2e4a6c3cd70a
hits = count(df₂.airline_sentiment .== ŷ)

# ╔═╡ e3ff7e20-c376-4bf7-b6da-0769551d5561
misses = length(ŷ) - hits

# ╔═╡ 92378e5b-c5e9-4a75-ae1e-216d67f4fd60
way_off = count(
        (df₂.airline_sentiment .≠ ŷ) .&
        (df₂.airline_sentiment .≠ "neutral") .&
        (ŷ .≠ "neutral")
)

## predicted positive and actually negative and conversely.

# ╔═╡ 240e86a4-4787-41ca-91ad-27893db438cc
begin
        accuracy_perc = hits / (hits + misses) * 100.
        round(accuracy_perc; digits=2)
end

# ╔═╡ f9e6c472-3dfb-4617-a0ef-99d1cbd34427
begin
        slightly_off_perc = (misses - way_off) / (hits + misses) * 100.
        round(slightly_off_perc; digits=2)
end

# ╔═╡ 05e7f1ba-ac19-4160-b158-0ef83af00e89
begin
        way_off_perc = way_off / (hits + misses) * 100.
        round(way_off_perc; digits=2)
end

# ╔═╡ fab2cd7c-cecc-41a8-ab28-b71ba00fec95
bar(["hits", "misses slightly", "misses way-off"],
        [accuracy_perc, slightly_off_perc , way_off_perc];
        legend = :none)

# ╔═╡ 52eef2d7-27b5-418f-b432-64f8561472b3
md"""
### Analyze some random tweets

Using the tweet API (ok, requires an account)...
"""

# ╔═╡ 8102896d-c25a-4210-81e5-94195b44f9e3
## using HTTP, JSON3
# token = readline("...<internal_dir>/twitter-bearer")

# ╔═╡ acbe6150-4732-48d5-a3cd-b695d3e0887c
# resp = HTTP.get("https://api.twitter.com/1.1/search/tweets.json?q=lang%3Aen%20flight", ["authorization" => "Bearer $token"]);

# ╔═╡ fdfabb92-3be3-4644-906a-eae308bdda0e
# begin
#       data = JSON3.read(response.body);
#       data[:statuses][1]
# end

# ╔═╡ 7e4bcd5e-0378-4c18-941d-feddba9ae872
# tweets = [x.text for x in data.statuses]

# ╔═╡ c63aa4f9-79b3-4651-85b9-50cb97b6c7c5
# result = test_model(model, tweets)

# ╔═╡ c5c2dd33-ad6e-451a-bbcc-12e84d0786b9
# table(result[9,:])

# ╔═╡ b10fe2c1-b451-48c4-8088-c03c85d31ed1
html"""

"""

# ╔═╡ Cell order:
# ╟─cbad2e04-9d92-11eb-3d84-9da6f9a3c818
# ╠═3d4d751a-b972-4e88-b568-1f19f0f20564
# ╠═e95db99e-a9ec-437c-9557-4dec4f42926f
# ╟─26957ae3-1a74-4be7-b96e-d3292991a82a
# ╟─225ebbac-2232-47ff-b7dc-e2c59814f196
# ╠═3ab74396-3e42-4545-af0b-f5746345324d
# ╟─13a43945-1a9c-4526-9f63-03a306b8d459
# ╠═7e812f27-4129-48ba-b4f5-a49ddc4169ca
# ╟─3a0b11b9-a51e-47f1-833c-8c9a2e4d73ce
# ╠═5f3e9cc0-fb07-4f2d-be7c-590f13bbd03d
# ╟─1c0d74d4-aac6-4fa4-8957-84457c73107e
# ╠═668b9aac-454a-43a6-8c40-44a0c023129a
# ╠═23c79d73-3431-44c3-a4ef-4dd8337f7931
# ╟─ca585d7b-d7f1-48f6-ac88-10aa582b6dbe
# ╠═1289f7cf-6c2a-4b67-a817-3de3926d46ec
# ╟─dbc871c5-9ca9-432c-98f8-af386ebdb22e
# ╠═e4c7dc20-f703-4d56-9835-af4fdcb32668
# ╟─6d7d638a-93a5-4c40-8d60-899ccafb0cdc
# ╠═4a649215-b35b-4e5c-96ab-ee9808d2de8e
# ╟─9e89132b-ae35-4d85-ae98-15b1ea326991
# ╟─de978888-c378-46e9-8158-2cdbc1a6fe2f
# ╠═e31a7663-6705-450e-9100-346b8fe51e15
# ╟─8db4042d-fba1-4bb7-ba43-0acc718e1389
# ╟─5fc69267-810a-479c-acdf-6026ca901f64
# ╠═50e09727-5fb9-4f02-9a71-f2fa65eb19b9
# ╠═9ea51080-e3f8-4d70-a139-aa29ddf6d13b
# ╠═d55d2b3c-82cf-4328-9aee-279e839d0a92
# ╠═69c18f2f-7e22-4188-9464-9705f4aed7e3
# ╠═e4186b6a-2ee3-4fd7-b1a3-06a499aeaf07
# ╟─7f5d2d98-e170-46ff-a470-e7708cbe3fb3
# ╠═10a4ddbc-66c4-472d-999e-99838da6652d
# ╠═a65e9da9-0e8b-4356-820f-266eaba857c8
# ╟─09441471-d246-4511-8d8e-e941b1166159
# ╟─5855ee27-4d34-41e4-9881-70ad4b5f064e
# ╟─819af8d7-c242-4f7b-a61f-121263bdb54c
# ╠═35b8ad6d-7c5e-43bf-b014-ce24010a959a
# ╠═41a18e5c-7e01-4828-9df0-53d446da9eb1
# ╠═3ac25d84-b396-4341-956b-c550054038c0
# ╠═cd495081-805b-4e71-814a-10add9a65d03
# ╠═bcae0e95-22fc-4526-9c71-a186e20d5d8d
# ╠═032ddbf5-5995-4ff9-be43-3cfa16eaac91
# ╠═6887b10b-0c44-4a1a-bbe0-e08032ac1dff
# ╟─0c7ad5a2-e8b1-4189-82f5-d39cbc122d2d
# ╠═2c853e55-843c-43c8-87b4-9ef6810dc03d
# ╟─2484eca4-198c-4178-8b78-26ee246f786b
# ╠═902b8634-7cde-4808-9e77-c13caac19c8f
# ╠═53c34fd9-f9c4-4321-94f3-2e4a6c3cd70a
# ╠═e3ff7e20-c376-4bf7-b6da-0769551d5561
# ╠═92378e5b-c5e9-4a75-ae1e-216d67f4fd60
# ╠═240e86a4-4787-41ca-91ad-27893db438cc
# ╠═f9e6c472-3dfb-4617-a0ef-99d1cbd34427
# ╠═05e7f1ba-ac19-4160-b158-0ef83af00e89
# ╠═fab2cd7c-cecc-41a8-ab28-b71ba00fec95
# ╟─52eef2d7-27b5-418f-b432-64f8561472b3
# ╠═8102896d-c25a-4210-81e5-94195b44f9e3
# ╠═acbe6150-4732-48d5-a3cd-b695d3e0887c
# ╠═fdfabb92-3be3-4644-906a-eae308bdda0e
# ╠═7e4bcd5e-0378-4c18-941d-feddba9ae872
# ╠═c63aa4f9-79b3-4651-85b9-50cb97b6c7c5
# ╠═c5c2dd33-ad6e-451a-bbcc-12e84d0786b9
# ╟─b10fe2c1-b451-48c4-8088-c03c85d31ed1

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    UniformSampler(size, replace=false)

A method for uniform sampling from spatial objects that produces
samples of given `size` with or without replacement depending on
the option `replace`.
"""
struct UniformSampler <: AbstractSampler
  size::Int
  replace::Bool
end

UniformSampler(size::Int) = UniformSampler(size, false)

function sample(object, sampler::UniformSampler)
  n = nelms(object)
  s = sampler.size
  r = sampler.replace
  if s > n && r == false
    @error "invalid sample size for sampling without replacement"
  end
  view(object, sample(1:n, s, replace=r))
end


text:
@testset "CheckNumElementsUnderRank" begin
    src = Shape([2, 3, 5, 7, 11, 13])
    @test 1 == lower_volume(src, 0)
    @test 2 == lower_volume(src, 1)
    @test 2 * 3 == lower_volume(src, 2)
    @test 2 * 3 * 5 == lower_volume(src, 3)
    @test 2 * 3 * 5 * 7 == lower_volume(src, 4)
    @test 2 * 3 * 5 * 7 * 11 == lower_volume(src, 5)
    @test 2 * 3 * 5 * 7 * 11 * 13 == lower_volume(src, 6)
end


text:
function f = fTr(x,Param)
switch lower(Param.ProfTr)
  case 'bickley'
    f=Param.RhoTheta*sin(Param.k*x(2));
  case 'sphericalharmonics'
    [lon,lat,r]=cart2sphere(x(1),x(2),x(3));
    %f1=real(harmonicY(Param.lHar,Param.mHar,pi/2-lat,lon,r));
    f=harmonicS(Param.lHar,Param.mHar,lat,lon,r);
end
end

function intG=integrandG(tau,RadEarth)
global Omega uM lat0G lat1G eN
f=2.0*Omega*sin(tau);
if (tau<=lat0G) || (tau>=lat1G)
  uStart=0.0;
else
  uStart=uM/eN*exp(1.0/((tau-lat0G)*(tau-lat1G)));
end
if abs(tau)<0.5*pi
  intG=(RadEarth*f+uStart*tan(tau))*uStart;
else
  intG=0.0;
end
end


text:
# Linear Regression with L1 penalty
PlotResults=true

#useproc("GPU") # GPU is significantly faster for larger systems
useproc("CPU")

StartCode()
X=ADnode()
w=ADvariable()
Y=ADnode()
loss=meanSquareLoss(X*w,Y)+0.1*meanAbs(w)+0.05*meanSquare(w) # mixed regularisation term
net=EndCode() # defines the graph

# instantiate parameter nodes and inputs:

D=100 # dimension of weight vector
N=500 # number of datapoints
truew=randn(D,1).*(rand(D,1).>0.8)
net.value[w]=truew
net.value[X]=randn(N,D)
net.value[Y]=net.value[X]*net.value[w] # make a realisable problem
net.value[w]=randn(D,1) # start with a random w

net=compile(net) # compile and preallocate memory

@gpu CUDArt.init([0]) # let the user do device management
@gpu net=convert(net,"GPU")
#gradcheck(net)

# Training:
println("Training: using $(net.gpu==true? "GPU" : "CPU") ")
ParsToUpdate=Parameters(net)
nupdates=500
LearningRate=0.1
error=Array(Float64,0)
velo=NesterovInit(net) # Nesterov velocity
for i=1:nupdates
    ADforward!(net)
    ADbackward!(net)
    push!(error,extract(net.value[net.FunctionNode])[1])
    printover("iteration $i: training loss = $(error[i])")
    for par in ParsToUpdate
        #GradientDescentUpdate!(net.value[par],net.gradient[par],LearningRate)
        NesterovGradientDescentUpdate!(net.value[par],net.gradient[par],velo[par],LearningRate,i)
    end
end

@gpu CUDArt.device_synchronize()

println("\ntrue weights   ",round(100*truew')/100)
println("learned weights",round(100*extract(net.value[w])')/100)

if PlotResults
    figure(1);plot(truew);title("true weights")
    figure(2);plot(extract(net.value[w])); title("learned weights")
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff