text:
abstract type Accumulator{A,B} end

intype(::Type{Accumulator{A}}) where A = A
intype(acc::Accumulator) = intype(typeof(acc))

outtype(::Type{Accumulator{A,B}}) where {A,B} = B
outtype(acc::Accumulator) = outtype(typeof(acc))

"""
"""
add!(::Accumulator{A}, ::A) where A = nothing

"""
"""
summarize(::Accumulator) = nothing

"""
"""
reset!(::Accumulator) = nothing

abstract type SimpleAccumulator{A} <: Accumulator{A,A} end


text:
module PrivateRyan

using Lazy

export @___!!checked!!___

const private = symbol("@__!!private!!1__")
const public = symbol("@free")

secretname(n) = symbol("#__$(n)__")

mod(T) = error("T is not checked.")

macro ___!!checked!!___(ex)
  if isexpr(ex, :type)
    name = namify(ex.args[2])
    for (i, field) in enumerate(ex.args[3].args)
      isexpr(field, :line, LineNumberNode) && continue
      (isexpr(field, :macrocall) && field.args[1] in [public, private]) ||
        error("Field must be labelled as public or private")
      ispub = field.args[1] == public
      fieldname = namify(field.args[2])
      !ispub && (field.args[2] = replace(field.args[2], fieldname, secretname(fieldname)))
      ex.args[3].args[i] = field.args[2]
    end
    quote
      $(esc(ex))
      PrivateRyan.mod(::$(esc(name))) = $(current_module())
      $(esc(name))
    end

  elseif isexpr(ex, :.)
    T = esc(ex.args[1])
    f = secretname(ex.args[2].args[1])
    mod = current_module()
    quote
      $mod == mod($T) ||
        error("$($(Expr(:quote, ex))) can only be accessed from $(mod($T))")
      $T.$f
    end

  elseif isexpr(ex, :(=)) && isexpr(ex.args[1], :.)
    T = esc(ex.args[1].args[1])
    f = secretname(ex.args[1].args[2].args[1])
    mod = current_module()
    quote
      $mod == mod($T) ||
        error("$($(Expr(:quote, ex))) can only be accessed from $(mod($T))")
      $T.$f = $(esc(ex.args[2]))
    end

  else
    error("Unsupported syntax")
  end
end

end # module

ext: jl


text:
"""
    RanefInfo

Information on random effects conditional modes/means, variances, etc.

Used for creating caterpillar plots.

!!! note
    This functionality may be moved upstream into MixedModels.jl in the near future.
"""
struct RanefInfo{T<:AbstractFloat}
    cnames::Vector{String}
    levels::Vector
    ranef::Matrix{T}
    stddev::Matrix{T}
end

"""
    ranefinfo(m::MixedModel)

Return a `NamedTuple{fnames(m), NTuple(k, RanefInfo)}` from model `m`
"""
function ranefinfo(m::MixedModel{T}) where {T}
    fn = fnames(m)
    val = sizehint!(RanefInfo[], length(fn))
    re = ranef(m)
    for grp in fn
        push!(val, ranefinfo(m, grp, re))
    end
    return NamedTuple{fn}((val...,))
end

"""
    ranefinfo(m::MixedModel, gf::Symbol)

Return a `RanefInfo` corresponding to the grouping variable `gf` model `m`.
"""
function ranefinfo(m::LinearMixedModel, gf::Symbol, re=ranef(m))
    idx = findfirst(==(gf), fnames(m))
    isnothing(idx) &&
        throw(ArgumentError("$gf is not the name of a grouping variable in the model"))

    # XXX replace ranef(m)[idx] with ranef(m, gf) when that becomes available upstream
    re, eff, cv = m.reterms[idx], re[idx], condVar(m, gf)
    return RanefInfo(
        re.cnames,
        re.levels,
        Matrix(adjoint(eff)),
        Matrix(adjoint(dropdims(sqrt.(mapslices(diag, cv; dims=1:2)); dims=2))),
    )
end

ranefinfo(m::GeneralizedLinearMixedModel, args...; kwargs...) = ranefinfo(m.LMM, args...; kwargs...)

"""
    ranefinfotable(ri::RanefInfo)

Return the information in `ri` as a column table (`NamedTuple` of `Vector`s)

The columns are

- `name`: name of the random effect
- `level`: level of the grouping factor
- `cmode`: conditional mode of the random effect
- `cstddev`: conditional standard deviation of the random effect

"""
function ranefinfotable(ri::RanefInfo)
    cnames, levels = ri.cnames, ri.levels
    k = length(cnames)
    l = length(levels)
    (;
        name=repeat(cnames; inner=l),
        level=repeat(levels; outer=k),
        cmode=vec(ri.ranef),
        cstddev=vec(ri.stddev),
    )
end

"""
    caterpillar!(f::Figure, r::RanefInfo; orderby=1)

Add Axes of a caterpillar plot from `r` to `f`.

The order of the levels on the vertical axes is increasing `orderby` column
of `r.ranef`, usually the `(Intercept)` random effects.
Setting `orderby=nothing` will disable sorting, i.e. return the levels in the
order they are stored in.

!!! note
    Even when not sorting the levels, they might have already been sorted during
    model matrix construction. If you want impose a particular ordering on the
    levels, then you must sort the relevant fields in the `RanefInfo` object before
    calling `caterpillar!`.
"""
function caterpillar!(f::Figure, r::RanefInfo; orderby=1)
    rr = r.ranef
    y = axes(rr, 1)
    ord = isnothing(orderby) ? y : sortperm(view(rr, :, orderby))
    cn = r.cnames
    axs = [Axis(f[1, j]) for j in axes(rr, 2)]
    linkyaxes!(axs...)
    for (j, ax) in enumerate(axs)
        xvals = view(rr, ord, j)
        scatter!(ax, xvals, y; color=(:red, 0.2))
        errorbars!(ax, xvals, y, 1.960 * view(r.stddev, ord, j); direction=:x)
        ax.xlabel = cn[j]
        ax.yticks = y
        j > 1 && hideydecorations!(ax; grid=false)
    end
    axs[1].yticks = (y, string.(r.levels[ord]))
    return f
end

"""
    caterpillar(m::LinearMixedModel, gf::Symbol)

Returns a `Figure` of a "caterpillar plot" of the random-effects means and prediction intervals

A "caterpillar plot" is a horizontal error-bar plot of conditional means and standard deviations
of the random effects.
"""
function caterpillar(m::MixedModel, gf::Symbol=first(fnames(m)))
    return caterpillar!(Figure(; resolution=(1000, 800)), ranefinfo(m)[gf])
end

"""
    qqcaterpillar!(f::Figure, r::RanefInfo; cols=axes(r.cnames, 1))

Update the figure with a caterpillar plot with the vertical axis on the Normal() quantile scale.

The display can be restricted to a subset of random effects associated with a grouping variable by specifying `cols`, either by indices or term names.
"""
function qqcaterpillar!(f::Figure, r::RanefInfo; cols=axes(r.cnames, 1))
    cols = _cols_to_idx(r, cols)
    cn, rr = r.cnames, r.ranef
    y = zquantile.(ppoints(size(rr, 1)))
    axs = [Axis(f[1, j]) for j in axes(cols, 1)]
    linkyaxes!(axs...)
    for (j, k) in enumerate(cols)
        ax = axs[j]
        xvals = rr[:, k]
        ord = sortperm(xvals)
        xvals = xvals[ord]
        scatter!(ax, xvals, y; color=(:red, 0.2))
        errorbars!(ax, xvals, y, 1.960 * view(r.stddev, ord, k); direction=:x)
        ax.xlabel = string(cn[k])
        j > 1 && hideydecorations!(ax; grid=false)
    end
    return f
end

_cols_to_idx(r, cols) = cols
_cols_to_idx(r, cols::AbstractVector{<:Symbol}) = _cols_to_idx(r, string.(cols))
function _cols_to_idx(r, cols::Vector{<:AbstractString})
    return [i for (i, c) in enumerate(r.cnames) if c in cols]
end

"""
    qqcaterpillar(m::LinearMixedModel, gf::Symbol=first(fnames(m)); cols=nothing)

Returns a `Figure` of a "qq-caterpillar plot" of the random-effects means and prediction intervals.

The display can be restricted to a subset of random effects associated with a grouping variable by specifying `cols`, either by indices or term names.
"""
function qqcaterpillar(m::MixedModel, gf::Symbol=first(fnames(m)); cols=nothing)
    reinfo = ranefinfo(m, gf)
    cols = something(cols, axes(reinfo.cnames, 1))
    return qqcaterpillar!(Figure(; resolution=(1000, 800)), reinfo; cols=cols)
end


text:
"""
    MatModule

Module for abstract material.
"""
module MatModule

__precompile__(true)

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict

"""
    AbstractMat

Abstract type of material.
"""
abstract type AbstractMat; end

"""
    massdensity(self::AbstractMat)

Return mass density.
"""
function massdensity(self::AbstractMat)
	return self.mass_density
end

end


text:
function smolyak_derivative(weights::Array{T,1},node::Array{T,1},multi_index::Array{S,2},pos::S) where {T<:AbstractFloat,S<:Integer}

  unique_multi_index = sort(unique(multi_index))
  unique_orders = m_i(unique_multi_index).-1

  # Here we construct the base polynomials

  base_polynomials = Array{Array{T,2},1}(undef,length(unique_orders))
  base_polynomial_derivatives = Array{Array{T,2},1}(undef,length(unique_orders))
  for i = 1:length(unique_orders)
    base_polynomials[i] = chebyshev_polynomial(unique_orders[i],node)
    base_polynomial_derivatives[i] = chebyshev_polynomial_derivative(unique_orders[i],node)
  end

  # Compute the unique polynomial terms from the base polynomials

  unique_base_polynomials = Array{Array{T,2},1}(undef,length(unique_orders))
  unique_base_polynomial_derivatives = Array{Array{T,2},1}(undef,length(unique_orders))
  for i = length(unique_orders):-1:2
    unique_base_polynomials[i] = base_polynomials[i][:,size(base_polynomials[i-1],2)+1:end]
    unique_base_polynomial_derivatives[i] = base_polynomial_derivatives[i][:,size(base_polynomial_derivatives[i-1],2)+1:end]
  end
  unique_base_polynomials[1] = base_polynomials[1]
  unique_base_polynomial_derivatives[1] = base_polynomial_derivatives[1]

  # Construct the first row of the interplation matrix

  polynomials = Array{T,1}(undef,length(weights))

  # Iterate over nodes, doing the above three steps at each iteration

  l = 1
  @inbounds for j = 1:size(multi_index,1)
    if pos == 1
      new_polynomials = unique_base_polynomial_derivatives[multi_index[j,1]][1,:]
    else
      new_polynomials = unique_base_polynomials[multi_index[j,1]][1,:]
    end
    for i = 2:size(multi_index,2)
      if pos == i
        new_polynomials = kron(new_polynomials,unique_base_polynomial_derivatives[multi_index[j,i]][i,:])
      else
        new_polynomials = kron(new_polynomials,unique_base_polynomials[multi_index[j,i]][i,:])
      end
    end
    m = length(new_polynomials)
    polynomials[l:l+m-1] = new_polynomials
    l += m
  end

  evaluated_derivative = zero(T)

  for i = 1:length(polynomials)
    evaluated_derivative += polynomials[i]*weights[i]
  end

  return evaluated_derivative

end

function smolyak_derivative(weights::Array{T,1},node::Array{T,1},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}},pos::S) where {T<:AbstractFloat,S<:Integer}

  nodes = copy(node)

  if size(domain,2) != length(node)
    error("domain is inconsistent with the number of dimensions")
  end

  d = length(node)
  for i = 1:d
    nodes[i] = normalize_node(node[i],domain[:,i])
  end

  evaluated_derivative = smolyak_derivative(weights,nodes,multi_index,pos)

  return evaluated_derivative

end

function smolyak_gradient(weights::Array{T,1},node::Array{T,1},multi_index::Array{S,2}) where {T<:AbstractFloat,S<:Integer}

  d = length(node)
  gradient = Array{T,2}(undef,1,d)

  for i = 1:d
    gradient[i] = smolyak_derivative(weights,node,multi_index,i)
  end

  return gradient

end

function smolyak_gradient(weights::Array{T,1},node::Array{T,1},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,S<:Integer}

  d = length(node)
  gradient = Array{T,2}(undef,1,d)

  for i = 1:d
    gradient[i] = smolyak_derivative(weights,node,multi_index,domain,i)
  end

  return gradient

end

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

function smolyak_gradient(weights::Array{T,1},multi_index::Array{S,2}) where {T<:AbstractFloat,S<:Integer}

  function goo(x::Array{T,1}) where {T <: AbstractFloat}

    return smolyak_gradient(weights,x,multi_index)

  end

  return goo

end

function smolyak_gradient(weights::Array{T,1},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,S<:Integer}

  function goo(x::Array{T,1}) where {T <: AbstractFloat}

    return smolyak_gradient(weights,x,multi_index,domain)

  end

  return goo

end

function smolyak_derivative(weights::Array{T,1},multi_index::Array{S,2},pos::S) where {T<:AbstractFloat,S<:Integer}

  function goo(x::Array{T,1}) where {T <: AbstractFloat}

    return smolyak_derivative(weights,x,multi_index,pos)

  end

  return goo

end

function smolyak_derivative(weights::Array{T,1},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}},pos::S) where {T<:AbstractFloat,S<:Integer}

  function goo(x::Array{T,1}) where {T <: AbstractFloat}

    return smolyak_derivative(weights,x,multi_index,domain,pos)

  end

  return goo

end


text:
using Pkg
Pkg.activate(".")

using Distributed
@everywhere using ParallelDataTransfer

@everywhere function rule30(ca::Array{Bool})
    lastv = ca[1]
    for i in 2:(length(ca)-1)
        current = ca[i]
        ca[i] = xor(lastv, ca[i] || ca[i+1])
        lastv = current
    end
end

@everywhere function getsetborder(ca::Array{Bool},
                                  neighbours::Tuple{Int64,Int64})
    ca[1] = (@fetchfrom neighbours[1] caa[end-1])
    ca[end] = (@fetchfrom neighbours[2] caa[2])
end

function printsimdist(workers::Array{Int})
    for w in workers
        dat = @fetchfrom w caa
        for b in dat[2:end-1]
            print(b ? "#" : " ")
        end
    end
    println()
end

function runca(steps::Int, visualize::Bool)
    @sync for w in workers()
        @async @fetchfrom w fill!(caa, false)
    end
    @fetchfrom wks[Int(nwks/2)+1] caa[2]=true
    visualize && printsimdist(workers())
    for i in 1:steps
        @sync for w in workers()
            @async @fetchfrom w getsetborder(caa, neighbours)
        end
        @sync for w in workers()
            @async @fetchfrom w rule30(caa)
        end
        visualize && printsimdist(workers())
    end
end

wks = workers()
nwks = length(wks)
for i in 1:nwks
    sendto(wks[i],neighbours = (i==1 ? wks[nwks] : wks[i-1],
                                i==nwks ? wks[1] : wks[i+1]))
    fetch(@defineat wks[i] const caa = zeros(Bool, 15+2));
end

runca(20,true)


text:
module DH32StatsHelpers

# package code goes here
export cMDS,
Cohen_d,
DataFrame_to_distmat,
zscore_log10_cols

#cMDS_of_DataFrame,
#kmeans_of_cMDS,

using DataFrames
using Distances
using Distributions
using MultivariateStats
using Clustering

using Distributions
using DataFrames

"""
Compute effect size measure Cohen's d

Input:
  - float array x1
  - float array x2
  - confidence level (default 0.95)

Output:
  - DataFrame with 1 row and 3 columns:
  d, lower_ci, upper_ci
"""
function Cohen_d(x1::Array{Float64,1}, x2::Array{Float64,1}, conf_level::Float64=0.95)

    n1 = length(x1)
    n2 = length(x2)
    dof = n1 + n2 - 2.
    d = (mean(x1) - mean(x2))/sqrt((((n1-1.)*var(x1)+(n2-1.)*var(x2))/(n1+n2-2.)))
    S_d = sqrt(((n1 + n2)/(n1 * n2) + 0.5 * d^2/dof) * ((n1 + n2)/dof))
    Z = -quantile(TDist(dof), (1.0-conf_level)/2.0)
    lower_ci = d - Z * S_d
    upper_ci = d + Z * S_d

    return DataFrame(d=d, lower_ci=lower_ci, upper_ci=upper_ci)
end

"""
Modify columns of a data frame of count data or other positive data by:

- Replacing 0 by 1 (to avoid NA after log)
- log10-transformation
- z-scores of log10-transformed

Output: transformed data frame
"""
function zscore_log10_cols(df::DataFrame)

    df_new = copy(df)

    for i in 1:length(df)

        #replace all zeros by 1 to avoid log(0)
        df_new[df_new[:,i] .== 0,i] = 1

        #take the log10 of all columns
        df_new[:,i] = map(log10,df[:,i])

        #take z-scores
        df_new[:,i] = zscore(convert(Array{Float64},df_new[:,i]))
    end

    return df_new
end

"""
Translates a data frame into a distance matrix.

Input:

   - data frame

   - distance (from Distances.jl), e.g. HellingerDist()

   - row_instances: true or false. Are rows of data frame the instances for which distances should be evaluated (true), or the columns (false)

Output: distance matrix (symmetric matrix with zero diagonal)

"""
function DataFrame_to_distmat(df::DataFrame, dist::Any, row_instances::Bool=true)

    x = convert(Array, df)

    if (row_instances == true)
        return pairwise(dist, x')
    end

    return pairwise(dist, x)

end

"""
Carry out a classical multi-dimensional scaling (cMDS).

Input:

    - distances (symmetric matrix)

    - dimension of cMDS output (default: 2)

    - should eigenvalues be computed? (true or false)

Output:

    - cMDS array with instances as columns

    - if eigenvalues == true: also 1-dim array of eigenvalues

"""
function cMDS(d::Array{Float64,2}, dim::Int64=2, eigenvalues::Bool=true)

    #do classical MDS:
    cMDS = classical_mds(d, dim)

    if eigenvalues == true
        #compute eigenvalues
        G = dmat2gram(d)
        E = eigfact!(Symmetric(G))
        #return cMDS array and eigenvalue vector descending from largest
        return cMDS, sort(E.values[:,1],rev=true)
    end

    return cMDS

end

end # module


text:
function u_test1(x, y)
    cos(x) * cos(y)
end

function u_test2(x, y)
    cos(x) - cos(y)
end

function mfunc_one!(u, m)
    m .= 1
end

function hfunc_u!(u, ϕ, h)
    h .= u
end

function hfunc_ϕ!(u, ϕ, h)
    h .= ϕ
end

# function testFokkerPlanckOperator()
#     M  = 64
#     N  = 64
#
#     grid = Grid2d(M,N)
#     u    = get_field(grid)
#     divJ = get_field(grid)
#     ft   = FourierTransform(grid)
#     op   = FokkerPlanckOperator(grid, ft)
#
#     evaluate_function_on_grid(grid, u_test, u)
#
#     collision_operator!(op, u, divJ)
#
# end

function testNullspace(utest, mfunc, hfunc)
    M  = 64
    N  = 64

    grid = Grid2d(M,N)
    u    = get_field(grid)
    divJ = get_field(grid)
    ft   = FourierTransform(grid)
    op   = FokkerPlanckOperator(grid, ft; mfunc=mfunc, hfunc=hfunc)

    evaluate_function_on_grid(grid, utest, u)

    collision_operator!(op, u, divJ)

    # println(maximum(abs.(divJ)))
    @test maximum(abs.(divJ)) ≈ zero(eltype(divJ)) atol=1E-9
    # @test maximum(abs.(divJ)) ≈ zero(eltype(divJ)) atol=1E-10

end

# testFokkerPlanckOperator()
testNullspace(u_test1, mfunc_one!, hfunc_ϕ!)
testNullspace(u_test2, mfunc_one!, hfunc_ϕ!)
testNullspace(u_test,  mfunc_one!, hfunc_u!)


text:
#=

  From BLOG example/poisson-ball.blog
  """
  Model file for balls in an urn, allowing observation errors.
  This version uses a Poisson prior for the number of balls.
  """

  Summary Statistics
     parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
         Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

   ...
       numBalls    5.8985    1.8454     0.0185    0.0695    677.0798    1.0002      125.6178
   ...
        numTrue    5.0069    1.5379     0.0154    0.0547    816.5438    1.0006      151.4924
       obsPost1    1.3249    0.4684     0.0047    0.0088   2952.0720    1.0000      547.6942
       obsPost2    1.6475    0.4778     0.0048    0.0096   2963.4386    1.0006      549.8031

  Distributions of variable numBalls (num:0)
  5.00000 =>    1996  (0.199600)
  8.00000 =>    1678  (0.167800)
  7.00000 =>    1588  (0.158800)
  6.00000 =>    1548  (0.154800)
  4.00000 =>    1371  (0.137100)
  3.00000 =>     827  (0.082700)
  9.00000 =>     752  (0.075200)
  2.00000 =>     204  (0.020400)
  1.00000 =>      36  (0.003600)

  Distributions of variable numTrue (num:0)
  5.00000 =>    2486  (0.248600)
  4.00000 =>    2180  (0.218000)
  6.00000 =>    2040  (0.204000)
  7.00000 =>    1216  (0.121600)
  3.00000 =>    1205  (0.120500)
  8.00000 =>     368  (0.036800)
  2.00000 =>     261  (0.026100)
  1.00000 =>     130  (0.013000)
  9.00000 =>      94  (0.009400)
  10.00000 =>      11  (0.001100)
  0.00000 =>       9  (0.000900)

  Distributions of variable obsPost1
  blue       =>    6751  (0.675100)
  green      =>    3249  (0.324900)

  Distributions of variable obsPost2
  green      =>    6475  (0.647500)
  blue       =>    3525  (0.352500)

  Cf ~/webppl/poisson_Ball.wppl

=#
using Turing, StatsPlots
include("jl_utils.jl")

@model function poisson_ball(obs=[1,2,1,2,1,2,1,2,1,2])
    n = length(obs)
    Blue = 1
    Green = 2;
    colors = [Blue,Green]

    # True color for each ball
    TrueColor ~ filldist(Categorical([0.5,0.5]),n) # colors

    # Number of balls
    numBalls ~ truncated(Poisson(6),0,n-1)

    # Which balls is drawn
    BallDrawn ~ filldist(DiscreteUniform(1,numBalls),n)

    # We observe Blue,Green,Blue,Green,Blue,Green,Blue,Green,Blue,Green
    for d in 1:n
        trueColor = TrueColor[BallDrawn[d]]
        if trueColor == Blue
            obs[d] ~ Categorical([0.8,0.2])
        else
            obs[d] ~ Categorical([0.2,0.8])
        end
    end

    # Posterior, first 2 balls
    obs_post = tzeros(2)
    for d in 1:2
        trueColor = TrueColor[BallDrawn[d]]
        if trueColor == Blue
            obs_post[d] ~ Categorical([0.8,0.2])
        else
            obs_post[d] ~ Categorical([0.2,0.8])
        end
    end

    numTrue ~ Dirac(sum([obs[d] == TrueColor[d] ? 1 : 0 for d in 1:n]))

    obsPost1 ~ Dirac(obs_post[1])
    obsPost2 ~ Dirac(obs_post[2])

end

model = poisson_ball()
# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(5), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:numBalls)
show_var_dist_pct(chns,:numTrue)
show_var_dist_pct(chns,:obsPost1,["blue","green"])
show_var_dist_pct(chns,:obsPost2,["blue","green"])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
function CountHorizontalBorders(r::Robot)
    counter = 0
    B = Array{Int, 1}(undef, 1000)
    B, counter = movetothecorner(r, West, B)

    num_horizontal_borders = snake(r)

    movements(r, Sud)
    movements(r, West)
    returnfromthecorner(r, inverse(West), B, counter)

    return num_horizontal_borders - 1
end

function snake(r::Robot)
    side = Ost
    num_horizontal_borders = movementsandcountborder(r, side)
    while isborder(r, Nord) == false
        move!(r, Nord)
        side = inverse(side)
        num_horizontal_borders += movementsandcountborder(r, side)
    end
    return num_horizontal_borders
end

function movementsandcountborder(r,side)
    num_borders = 0
    while !isborder(r,side)
        move!(r,side)
        if isborder(r,Nord)
            num_borders += 1
        end
        while isborder(r,Nord)
            if isborder(r,side)
                break
            end
            move!(r,side)
        end
    end
    return num_borders
end

function inverse(side)
    return HorizonSide((Int(side) + 2) % 4)
end

function movements(r, side)
    while isborder(r, side) == false
        move!(r, side)
    end
end

function movetothecorner(r, side, B)
    counter = 1
    while ((!isborder(r, side)) || (!isborder(r, HorizonSide((Int(side) + 1) % 4))))
        if !isborder(r, side)
            move!(r, side)
            B[counter] = 1
            counter += 1
        end
        if !isborder(r, HorizonSide((Int(side) + 1) % 4))
            move!(r, HorizonSide((Int(side) + 1) % 4))
            B[counter] = 2
            counter += 1
        end
    end
    return B, counter
end

function returnfromthecorner(r, side, B, counter)
    for i in 1:(counter - 1)
        if B[counter - i] == 1
            move!(r, side)
        end
        if B[counter - i] == 2
            move!(r, HorizonSide((Int(side) + 1) % 4))
        end
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule MariaDB_Connector_C_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("MariaDB_Connector_C")
JLLWrappers.@generate_main_file("MariaDB_Connector_C", UUID("aabc7e14-95f1-5e66-9f32-aea603782360"))
end  # module MariaDB_Connector_C_jll


text:
using Random
using Statistics, ZigZagBoomerang, LinearAlgebra, Test, SparseArrays
using ZigZagBoomerang: Partition
@testset "Partition" begin
    d = 6
    partition = ZigZagBoomerang.Partition(2, d)
    for i in 1:d
        a,b = partition(i)
        @test i == partition(a,b)
        @test a == (i > 3) + 1
        @test b == mod1(i, 3)
    end

    d = 6
    partition = ZigZagBoomerang.Partition(3, d)
    for i in 1:d
        a,b = partition(i)
        @test i == partition(a,b)
    end
end
@testset "Parallel ZigZag" begin
    Random.seed!(1)
    d = 20
    K = 2
    T = 1000.0
    Δ = 0.05

    d2 = d÷K

    Γ = sparse(SymTridiagonal(1.0ones(d), -0.4ones(d-1)))
    Γ2 = copy(Γ)
    for k1 in 0:K-1
        for k2 in 0:K-1
            k2 == k1 && continue
            Γ2[(k1*d2 + 1):(k1+1)*d2,(k2*d2 + 1):(k2+1)*d2] .= 0
        end
    end

    dropzeros!(Γ2)

    partition = ZigZagBoomerang.Partition(K, d)

    ∇ϕ(x, i, Γ) = ZigZagBoomerang.idot(Γ, i, x) # sparse computation

    t0 = 0.0
    x0 = 0.1randn(d)
    θ0 = rand([-1.0,1.0], d)

    G = [i => rowvals(Γ)[nzrange(Γ, i)] for i in eachindex(θ0)]

    c = 5*[norm(Γ[:, i], 2) for i in 1:d]

    dt = 0.5

    Z = ZigZag(Γ2, x0*0)
    println("Multithreaded: (", Threads.nthreads(), " cores)")
    tr, (t, x, θ), (acc, num) = @time ZigZagBoomerang.parallel_spdmp(partition, ∇ϕ, t0, x0, θ0, T, c, G, Z, Γ; Δ=Δ)
    @test 0.1/sqrt(T) < mean(abs.(mean(tr))) < 4/sqrt(T)
    @show mean(abs.(mean(tr)))

    if d < 100
         ts, xs = ZigZagBoomerang.sep(collect(discretize(tr, dt)))
         @test mean(abs.(cov(xs) - inv(Matrix(Γ)))) < 4/sqrt(T)
    end

    println("Single thread:")
    Z = ZigZag(Γ, x0*0)
    tr, _ = @time ZigZagBoomerang.spdmp(∇ϕ, t0, x0, θ0, T, c, Z, Γ)
    @test 0.1/sqrt(T) < mean(abs.(mean(tr))) < 4/sqrt(T)
    if d <= 100
        ts, xs = ZigZagBoomerang.sep(collect(discretize(tr, dt)))
        @test mean(abs.(cov(xs) - inv(Matrix(Γ)))) <4/sqrt(T)
    end

end


text:
process_value(::Val{:stencil}, key, val, global_p, current_p) =
    process_value(key, val, global_p, current_p, ())

process_value(::Val{:laser}, key, val, global_p, current_p) =
    process_value(key, val, global_p, current_p, ())

process_value(::Val{:species}, key, val, global_p, current_p) =
    process_value(key, val, global_p, current_p, ())

process_value(::Val{:window}, key, val, global_p, current_p) =
    process_value(key, val, global_p, current_p, ())


text:
# Benchmark the parsing of a PDB file given as an argument

using Bio.Structure

pdb_filepath = ARGS[1]

# Run to JIT compile
read(pdb_filepath, PDB)

elapsed = @elapsed struc = read(pdb_filepath, PDB)

println(elapsed)

max_stars_repo_name: chemlove/pdb-benchmarks

max_issues_repo_name: chemlove/pdb-benchmarks

max_forks_repo_name: chemlove/pdb-benchmarks


text:
using PyPlot
using JLD2
using DelimitedFiles
using ADSeismic

n = 200
@load "data/dat0-1.jld2" U
U0 = U

Uk = zeros(size(U0,1), n, n)
@load "data/dat0-4.jld2" U
Uk[:,1:n÷2,1:n÷2] = U

@load "data/dat1-4.jld2" U
Uk[:,1:n÷2,n÷2+1:end] = U

@load "data/dat2-4.jld2" U
Uk[:,n÷2+1:end, 1:n÷2] = U

@load "data/dat3-4.jld2" U
Uk[:,n÷2+1:end, n÷2+1:end] = U

@show maximum(abs.(U0 - Uk))


text:
using RadioTools
using Base.Test

λ = cm2m(fre2wl(ghz2hz(12)))
D = 2

@test angulares(λ,D) == 3122.8380843749997


text:
# Important packages
using JLD

# Loading data...
println("== Loading data...")
d = load("backup_matrices/eigenvalues.jld")["values"]
v = load("backup_matrices/eigenvectors.jld")["vectors"]

# Defining the size of the word vectors
println("== Defining word vectors size...")
m = 50

# Sorting eigenvalues by its absolute value (BubbleSort)
println("== BubbleSorting eigenvalues by its absolute values...")
for i = 1:length(d)
    swapped = false
    for j = 1:length(d) - i - 1
        if (abs(d[j]) < abs(d[j + 1]))
            aux = d[j]
            d[j] = d[j+1]
            d[j + 1] = aux
            swapped = true
        end
    end
    if swapped == false
        break
    end
end

# Filling a trunkated matrix (MxN) w/ most significant eigenvalues
println("== Filling the trunkated matrix (MxN)...")
d2 = zeros(m, length(d))
for i = 1:m
    d2[i, i] += sqrt(d[i])
end

# Reconstructing matrix
println("== Constructing W matrix with all the word vectors...")
word_vectors = d2 * v.'

# Saving reconstructed matrix
println("== Exporting word vectors matrix...")
save("backup_matrices/wordvectors.jld", "wordvectors", word_vectors)