text:
# Convient function for rearranging contexts

# Create a new context containing the given contexts stacked horizontally.
#
# Args:
#  x0: X-position of the new root context
#  y0: Y-position of the new root context
#  height: Height of the root context.
#  aligned_contexts: One or more canveses accompanied with a vertical alignment
#                    specifier, giving the vertical positioning of the context.
#
function hstack(x0, y0, height, aligned_contexts::(Context, VAlignment)...)
    if isempty(aligned_contexts)
        return context(x0, y0, 0, height)
    end

    # To get the expected results, we scale width units, so that everything
    # fits.
    total_width_units = sum(Float64[context.box.width.cw
                                    for (context, _) in aligned_contexts])
    width = sum(Measure[context.box.width
                        for (context, _) in aligned_contexts])
    width = Measure(width,
                    cw=total_width_units > 0.0 ?
                        width.cw / total_width_units : 0.0)

    height = y_measure(height)

    root = context(x0, y0, width, height)
    x = Measure()
    for (context, aln) in aligned_contexts
        context = copy(context)
        context.box = copy(context.box)

        if context.box.width.cw != 0.0
            context.box =
                BoundingBox(context.box,
                    width=Measure(context.box.width,
                        cw=context.box.width.cw / total_width_units))
        end

        # Should we interpret vbottom to mean 0?
        context.box = BoundingBox(context.box, x0=x)
        if aln == vtop
            context.box = BoundingBox(context.box, y0=Measure{T}())
        elseif aln == vcenter
            context.box = BoundingBox(context.box,
                            y0=(height / 2) - (context.box.height / 2))
        elseif aln == vbottom
            context.box = BoundingBox(context.box,
                            y0=height - context.box.height)
        end

        root = compose!(root, context)
        x += context.box.width
    end

    return root
end

hstack() = context()

# Create a new context containing the given contexts stacked horizontally.
#
# This is the simple version of hstack. The root context will be placed on 0cx,
# 0cy, and its height will be the maximum of the contexts it contains. All
# contexts will be centered vertically.
#
function hstack(contexts::Context...; x0::MeasureOrNumber=0,
                y0::MeasureOrNumber=0, height=0)
    if height == 0
        height = maximum([context.box.height for context in contexts])
    end
    return hstack(x0, y0, height, [(context, vcenter) for context in contexts]...)
end

# Create a new context containing the given contexts stacked vertically.
#
# Args:
#  x0: X-position of the new root context
#  y0: Y-position of the new root context
#  width: Height of the root context.
#  aligned_contexts: One or more canveses accompanied with a horizontal alignment
#                    specifier, giving the horizontal positioning of the context.
#
function vstack(x0, y0, width, aligned_contexts::(Context, HAlignment)...)

    if isempty(aligned_contexts)
        return context(x0, y0, width, 0)
    end

    # Scale height units
    total_height_units = sum(Float64[context.box.height.ch
                                     for (context, _) in aligned_contexts])
    height = sum(Measure[context.box.height
                         for (context, _) in aligned_contexts])
    height = Measure(height,
                     ch=total_height_units > 0.0 ?
                          height.ch / total_height_units : 0.0)

    width = x_measure(width)

    root = context(x0, y0, width, height)
    y = Measure()
    for (context, aln) in aligned_contexts
        context = copy(context)
        context.box = copy(context.box)

        if context.box.height.ch != 0.0
            context.box =
                BoundingBox(context.box,
                    height=Measure(context.box.height,
                        ch=context.box.height.ch / total_height_units))
        end

        context.box = BoundingBox(context.box, y0=y)
        if aln == hleft
            context.box = BoundingBox(context.box, x0=Measure{T}())
        elseif aln == hcenter
            context.box = BoundingBox(context.box,
                            x0=(width / 2) - (context.box.width / 2))
        elseif aln == hright
            context.box = BoundingBox(context.box,
                            x0 = width - context.box.width)
        end

        root = compose(root, context)
        y += context.box.height
    end

    return root
end

vstack() = context()

# Create a new context containing the given contexts stacked horizontally.
#
# The simple version of vstack. The root context will be placed on 0cx, 0cy, and
# its width will be the maximum of the contexts it contains. All contexts will
# be centered horizontally..
#
function vstack(contexts::Context...; x0::MeasureOrNumber=0,
                y0::MeasureOrNumber=0, width::MeasureOrNumber=0)
    if width == 0
        width = max([context.box.width for context in contexts]...)
    end
    return vstack(x0, y0, width, [(context, hcenter) for context in contexts]...)
end

ext: jl


text:
 using IJulia

println("The browser will open, please navigate to .ipynb in the
 project directory")
 notebook()

max_stars_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_issues_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_forks_repo_name: jmulloka/Breast-Cancer-Text-Mining


text:
#=##############################################################################
# # DESCRIPTION
#     Examples of VTK formatting and usage of VTKtools
# # AUTHORSHIP
#   * Author    : Eduardo J Alvarez
#   * Email     : Edo.AlvarezR@gmail.com
#   * Created   : Nov 2017
#   * License   : MIT License
=###############################################################################

"""
  Example of using lines2vtk with generateVTK. Generates a simple 6-sided box.
"""
function simple_box(; prompt=true)
  file_name1 = "temp_vtk_example"
  file_name2 = "temp_vtk_example_w_data"
  # Defining points
  p0 = [0,0,0]
  p1 = [1,0,0]
  p2 = [1,1,0]
  p3 = [0,1,0]
  p4 = [0,0,1]
  p5 = [1,0,1]
  p6 = [1,1,1]
  p7 = [0,1,1]

  # Defining sides (cells)
  c0 = [p0, p3, p2, p1] # Bottom
  c1 = [p4, p5, p6, p7] # Top
  c2 = [p0, p4, p7, p3] # Left
  c3 = [p1, p2, p6, p5] # Right
  c4 = [p0, p1, p5, p4] # Front
  c5 = [p2, p3, p7, p6] # Back

  # Formats the cells into vtk format
  points, vtk_cells = gt.lines2vtk([c0,c1,c2,c3,c4,c5])

  # Generates the vtk file
  gt.generateVTK(file_name1, points; cells=vtk_cells)

  # Calls paraview
  run(`paraview --data="$(file_name1).vtk;"`)

  # ------ Example of outputting a data field ----------------------------------
  # Data for each point in each cell
  d0 = [0, 3, 2, 1] # Data of cell 0
  d1 = [4, 5, 6, 7]
  d2 = [0, 4, 7, 3]
  d3 = [1, 2, 6, 5]
  d4 = [0, 1, 5, 4]
  d5 = [2, 3, 7, 6]

  # Formats the cells into vtk format
  points, vtk_cells, vtk_values = gt.lines2vtk([c0,c1,c2,c3,c4,c5];
                                             values=[d0,d1,d2,d3,d4,d5])

  # Formats the data for generateVTK
  data = []
  push!(data, Dict(
                  "field_name" => "Point_index",
                  "field_type" => "scalar",
                  "field_data" => vtk_values
                  )
       )
  # Generates the vtk file
  gt.generateVTK(file_name2, points;
               cells=vtk_cells, point_data=data)
  # Calls paraview
  run(`paraview --data="$(file_name2).vtk;"`)

  # Deletes files
  if prompt
    print("Delete vtk files? ([y]/n) ")
    y = readline()
  else
    y = "y"
  end
  if y=="y"; run(`rm -f $(file_name1).vtk $(file_name2).vtk`); end;

end

"""
  Example of using a parametrically meshed surface using `discretize()`,
  `lines2vtkcells()`, and `generateVTK()`.
"""
function parametric_mesh(; prompt=true)
  file_name = "temp_mesh_example"

  # Helicoid function
  revs = 2                        # Number of revolutions
  rho1 = 8                        # Outer radius
  rho2 = rho1/4                   # Inner radius
  k = 1                           # Period
  f1(z) = [rho1*cos(k*z), rho1*sin(k*z), z] # Outer edge
  f2(z) = [rho2*cos(k*z), rho2*sin(k*z), z] # Inner edge

  # Discretize edges
  xlow, xhigh = 0, revs*2*pi      # Bounds of parametric edges
  n = 100                         # Number of cells
  r = 1.0                         # Expansion ratio
  line1 = gt.discretize(f1, xlow, xhigh, n, r)
  line2 = gt.discretize(f2, xlow, xhigh, n, r)

  # Dummy point data for good looking visuals
  pd1 = [i for i in 1:size(line1)[1]]
  pd2 = size(line1)[1] .+ [i for i in 1:size(line2)[1]]

  # Generates cells in VTK Legacy format
  out = gt.lines2vtkcells(line1, line2; point_data1=pd1, point_data2=pd2)
  points, vtk_cells, point_data = out

  # Formats the point data for generateVTK
  data = []
  push!(data, Dict(
                  "field_name" => "Point_index",
                  "field_type" => "scalar",
                  "field_data" => point_data
                  )
       )

   # Generates the vtk file
   gt.generateVTK(file_name, points; cells=vtk_cells, point_data=data)

   # Calls paraview
   run(`paraview --data="$(file_name).vtk;"`)

   # Deletes files
   if prompt
     print("Delete vtk files? ([y]/n) ")
     y = readline()
   else
     y = "y"
   end
   if y=="y"; run(`rm -f $(file_name).vtk`); end;
end

"""
  Example of using a parametrically meshed surface using `discretize()`,
  `lines2vtkmulticells()`, and `generateVTK()`.
"""
function parametric_mesh2(; prompt=true)
  file_name = "temp_mesh_example2"

  # Helicoid function
  revs = 2                        # Number of revolutions
  rho1 = 8                        # Outer radius
  rho2 = rho1/4                   # Inner radius
  k = 1                           # Period
  f1(z) = [rho1*cos(k*z), rho1*sin(k*z), z] # Outer edge
  f2(z) = [rho2*cos(k*z), rho2*sin(k*z), z] # Inner edge

  # Discretize edges
  xlow, xhigh = 0, revs*2*pi      # Bounds of parametric edges
  n = 100                         # Number of cells
  r = 1.0                         # Expansion ratio
  line1 = gt.discretize(f1, xlow, xhigh, n, r)
  line2 = gt.discretize(f2, xlow, xhigh, n, r)

  # Discretize width between lines
  nwidth = 25                     # Number of rows between edges
  rwidth = 10.0                   # Expansion ratio
  central = true                  # Expands about center between edges
  sections = [(1.0, nwidth, rwidth, central)]   # Discretization sections

  # Dummy point data for good looking visuals
  pd1 = [i for i in 1:size(line1)[1]]
  pd2 = size(line1)[1] .+ [i for i in 1:size(line2)[1]]

  # Generates cells in VTK Legacy format
  out = gt.lines2vtkmulticells(line1, line2, sections;
                                            point_data1=pd1, point_data2=pd2)
  points, vtk_cells, point_data = out

  # Formats the point data for generateVTK
  data = []
  push!(data, Dict(
                  "field_name" => "Point_index",
                  "field_type" => "scalar",
                  "field_data" => point_data
                  )
       )

   # Generates the vtk file
   gt.generateVTK(file_name, points; cells=vtk_cells, point_data=data)

   # Calls paraview
   run(`paraview --data="$(file_name).vtk;"`)

   # Deletes files
   if prompt
     print("Delete vtk files? ([y]/n) ")
     y = readline()
   else
     y = "y"
   end
   if y=="y"; run(`rm -f $(file_name).vtk`); end;
end

"""
  Example of using lines2vtk with generateVTK. Generates a simple 6-sided box.
"""
function simple_box2(; prompt=true)
  file_name1 = "temp_vtk_example"
  # Defining points
  p0 = [0,0,0]
  p1 = [1,0,0]
  p2 = [1,1,0]
  p3 = [0,1,0]
  p4 = [0,0,1]
  p5 = [1,0,1]
  p6 = [1,1,1]
  p7 = [0,1,1]

  points = [p0, p1, p2, p3, p4, p5, p6 ,p7]

  # Data for each point in each cell
  c0 = [0, 3, 2, 1] # Data of cell 0
  c1 = [4, 5, 6, 7]
  c2 = [0, 4, 7, 3]
  c3 = [1, 2, 6, 5]
  c4 = [0, 1, 5, 4]
  c5 = [2, 3, 7, 6]

  cells = [c0, c1, c2, c3, c4, c5]

  cell_index = [1,2,3,4,5,6]

  cell_index_field = Dict(
                            "field_name" => "CellInd",
                            "field_type" => "scalar",
                            "field_data" => cell_index
                          )

  cell_data = [cell_index_field]

  # Generates the vtk file
  gt.generateVTK(file_name1, points; cells=cells, cell_data=cell_data)

  # Calls paraview
  run(`paraview --data="$(file_name1).vtk;"`)

end


text:
@doc raw"""
    Rule{T<:AbstractFloat}

Quadrature rule.

Approximation of an integral over `[a, b]` by a sum over discrete points `x` with weights `w`:
```math
    ∫ f(x) ω(x) dx ≈ ∑_i f(x_i) w_i
```
where we generally have superexponential convergence for smooth ``f(x)`` in
the number of quadrature points.
"""
struct Rule{T}
    x::Vector{T}
    w::Vector{T}
    a::T
    b::T

    x_forward::Vector{T}
    x_backward::Vector{T}

    function Rule(x, w, a=-1, b=1, x_forward=x .- a, x_backward=b .- x)
        a ≤ b || error("a must be ≤ b")
        all(≤(b), x) || error("x must be ≤ b")
        all(≥(a), x) || error("x must be ≥ a")
        issorted(x) || error("x must be strictly increasing")
        length(x) == length(w) ||
            throw(DimensionMismatch("x and w must have the same length"))
        return new{eltype(x)}(x, w, a, b, x_forward, x_backward)
    end
end

"""
    quadrature(rule, f)

Approximate `f`'s integral.
"""
function quadrature(rule, f)
    return sum(rule.w .* f.(rule.x))
end

"""
    reseat(rule, a, b)

Reseat quadrature rule to new domain.
"""
function reseat(rule::Rule, a, b)
    scaling = (b - a) / (rule.b - rule.a)
    x = (rule.x .- (rule.a + rule.b) / 2) * scaling .+ (a + b) / 2
    w = rule.w * scaling
    x_forward = rule.x_forward * scaling
    x_backward = rule.x_backward * scaling
    return Rule(x, w, a, b, x_forward, x_backward)
end

"""
    scale(rule, factor)

Scale weights by `factor`.
"""
function scale(rule, factor)
    return Rule(rule.x, rule.w * factor, rule.a, rule.b, rule.x_forward, rule.x_backward)
end

"""
    piecewise(rule, edges)

Piecewise quadrature with the same quadrature rule, but scaled.
"""
function piecewise(rule, edges::Vector)
    start = edges[begin:(end - 1)]
    stop = edges[(begin + 1):end]
    issorted(edges) || error("edges must be monotonically increasing")
    return joinrules([reseat(rule, a, b) for (a, b) in zip(start, stop)])
end

"""
    joinrules(rules)

Join multiple Gauss quadratures together.
"""
function joinrules(rules::AbstractVector{Rule{T}}) where {T<:AbstractFloat}
    for i in Iterators.drop(eachindex(rules), 1)
        rules[i - 1].b == rules[i].a || error("rules must be contiguous")
    end

    x = reduce(vcat, rule.x for rule in rules)
    w = reduce(vcat, rule.w for rule in rules)
    a = first(rules).a
    b = last(rules).b

    x_forward = reduce(vcat, rule.x_forward .+ (rule.a - a) for rule in rules)
    x_backward = reduce(vcat, rule.x_backward .+ (b - rule.b) for rule in rules)

    return Rule(x, w, a, b, x_forward, x_backward)
end

"""
    legendre(n[, T])

Gauss-Legendre quadrature with `n` points on [-1, 1].
"""
legendre(n) = Rule(gauss(n)...)
legendre(n, T) = Rule(gauss(T, n)...)

"""
    legendre_collocation(rule, n=length(rule.x))

Generate collocation matrix from Gauss-Legendre rule.
"""
function legendre_collocation(rule, n=length(rule.x))
    res = (legvander(rule.x, n - 1) .* rule.w)'
    invnorm = range(0.5; length=n)
    res .*= invnorm
    return res
end

function Base.convert(::Type{Rule{T}}, rule::Rule) where {T}
    return Rule(
        T.(rule.x), T.(rule.w),
        T(rule.a), T(rule.b),
        T.(rule.x_forward), T.(rule.x_backward),
    )
end

"""
    NestedRule{T}

Nested quadrature rule.
"""
struct NestedRule{T}
    rule::Rule{T}
    v::Vector{T}
end

function reseat(nrule::NestedRule, a, b)
    newrule = reseat(nrule.rule, a, b)
    newv = (b - a) / (rule.b - rule.a) * nrule.v
    return NestedRule(newrule, newv)
end

function kronrod_21()
    x, w, v = kronrod(10)
    append!(x, -reverse(x)[2:end])
    append!(w, reverse(w)[2:end])
    append!(v, reverse(v))
    return NestedRule(Rule(x, w), v)
end


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

using Documenter, ErrorTypes

DocMeta.setdocmeta!(ErrorTypes, :DocTestSetup, :(using ErrorTypes); recursive=true)

makedocs(;
    modules=[ErrorTypes],
    format=Documenter.HTML(),

    pages = [
        "Home" => "index.md",
        "Usage" => "usage.md",
        "Why use ErrorTypes?" => "motivation.md",
    ],
    repo="https://github.com/jakobnissen/ErrorTypes.jl/blob/{commit}{path}#L{line}",
    sitename="ErrorTypes.jl",
    authors="Jakob Nybo Nissen",
    assets=String[],
)

deploydocs(;
    repo="github.com/jakobnissen/ErrorTypes.jl.git",
)


text:
@testset "Arithmetic" begin
    x = rand(Complex{Double64})
    y = rand(Double64)

    @test x + y isa Complex{Double64}
    @test y + x isa Complex{Double64}
    @test x - y isa Complex{Double64}
    @test y - x isa Complex{Double64}
    @test x * y isa Complex{Double64}
    @test y * x isa Complex{Double64}
    @test x / y isa Complex{Double64}
    @test y / x isa Complex{Double64}
    @test x ^ y isa Complex{Double64}
    @test y ^ x isa Complex{Double64}
    @test x ^ x isa Complex{Double64}
    @test y ^ y isa Double64

    x, y = rand(2)
    @test add2(x, y) == (Double64(x) + Double64(y))
    @test sub2(x, y) == (Double64(x) - Double64(y))
    @test mul2(x, y) == (Double64(x) * Double64(y))
    @test div2(x, y) == (Double64(x) / Double64(y))

    @test x ⊕ y == (Double64(x) + Double64(y))
    @test x ⊖ y == (Double64(x) - Double64(y))
    @test x ⊗ y == (Double64(x) * Double64(y))
    @test x ⊘ y == (Double64(x) / Double64(y))
end

@testset "Trig functions" begin

    x = rand(Double64)
    sinx, cosx = sincos(x)
    @test sinx ≈ sin(x)
    @test cosx ≈ cos(x)
end


text:
using Documenter
using Literate
using BioMedQuery

# compile all examples in BioMedQuery/examples/literate_src into markdown and jupyter notebooks for documentation
for (root, dirs, files) in walkdir("examples/literate_src")
    for file in files
        Literate.notebook(joinpath(root,file), joinpath(@__DIR__, "src", "notebooks"), execute=false, documenter=true)
    end
end

for (root, dirs, files) in walkdir("examples/literate_src")
    for file in files
        Literate.markdown(joinpath(root,file), joinpath(@__DIR__, "src", "examples"), documenter=true)
    end
end

# make assets dir if doesn't exist
assets_dir = joinpath(@__DIR__,"src","assets")
if !isdir(assets_dir)
    mkdir(assets_dir)
end

base_url = "https://raw.githubusercontent.com/bcbi/code_style_guide/master/assets/"

# get/replace favicon
favicon_url = base_url*"favicon.ico"
favicon_path = joinpath(assets_dir,"favicon.ico")
run(`curl -g -L -f -o $favicon_path $favicon_url`)

# get/replace css
css_url = base_url*"bcbi.css"
css_path = joinpath(assets_dir,"bcbi.css")
run(`curl -g -L -f -o $css_path $css_url`)

# get/replace logo
logo_url = base_url*"bcbi-white-v.png"
logo_path = joinpath(assets_dir,"logo.png")
run(`curl -g -L -f -o $logo_path $logo_url`)

makedocs(
    modules = [ BioMedQuery],
    assets = [
        "assets/favicon.ico",
        "assets/bcbi.css",
        "assets/logo.png"
        ],
    sitename = "BioMedQuery.jl",
    debug = true,
    pages = [
        "Home" => "index.md",
        "Examples" => Any[
            "Overview" => "examples.md",
            "Pubmed Search and Save" => "examples/1_pubmed_search_and_save.md",
            "MeSH/UMLS Map and Filtering" => "examples/2_umls_map_and_filter.md",
            "Export to Citations" => "examples/4_pubmed_export_citations.md",
            "Load MEDLINE" => "examples/5_load_medline.md",
        ],
        "Manual" => Any[
            "Processes/Workflows" => "processes.md",
            "PubMed" => "pubmed.md",
            "Clinical Trials" => "ct.md",
            "Database Utilities" => "dbutils.md",
            "Library" => "library.md"
            ]
        ]
    )

deploydocs(
    repo = "github.com/bcbi/BioMedQuery.jl.git",
)

max_stars_repo_name: bcbi/NLM


text:
module TerminalLoggers

@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), "```julia" => "```jldoctest")
end TerminalLoggers

using Logging:
    AbstractLogger,
    LogLevel, BelowMinLevel, Debug, Info, Warn, Error, AboveMaxLevel

import Logging:
    handle_message, shouldlog, min_enabled_level, catch_exceptions

using LeftChildRightSiblingTrees: Node, addchild, isroot, prunebranch!
using ProgressLogging: asprogress
using UUIDs: UUID

export TerminalLogger

const ProgressLevel = LogLevel(-1)

include("ProgressMeter/ProgressMeter.jl")
using .ProgressMeter:
    BarGlyphs

include("StickyMessages.jl")
include("ProgressBar.jl")
include("TerminalLogger.jl")

end # module


text:
a = @queue 1 2 3 4
b = @stack 1 2 3 4

values(a)
values(b)

max_stars_repo_name: frankill/ftest

max_issues_repo_name: frankill/ftest

max_forks_repo_name: frankill/ftest


text:
"""
    ℤd{T <: Integer, d} <: Integer

The  following immutable struct stores a value of type `T <: Integer` modulo `d`.

The type `T` is customizable, but should be choosen as the smallest `Integer` type that fits `d`.
For normal use cases, `Int8` is enough.

Example:
```julia
ℤd{Int8, 2}(1) # stores a 1
ℤd{Int8, 2}(2) # stores mod(2,2) = 0
```

For convenience it is recommended to use type aliases if you (mostly) work with fixed `d`:

```julia
# Define a type alias for the integers modulo two
const ℤ₂ = ℤd{Bool, 2}
# Define a type alias for the integers modulo three
const ℤ₃ = ℤd{Int8, 3}
```
Another method is to use the provided type aliases via `import`:

    ℤ₂ = Zd_Arithmetics.Z2
    ℤ₂ = Zd_Arithmetics.ℤ₂
or

    import Zd_Arithmetics.Z2
    import Zd_Arithmetics.ℤ₂

"""
struct ℤd{T <: Integer, d} <: Integer
    value::T
    ℤd{T, d}(x::Integer) where {T <: Integer, d} = new(T(mod(x, d)))
end

"""
    ℤd{T, d}(x::ℤd{T, d′}; force=false)::ℤd{T,d} where {T <: Integer, d, d′}

Upwards converting from `ℤd′` to `ℤd` for `d′ < d` (possible, because ``\\mathbb{Z}_{d'} \\subset \\mathbb{Z}_d`` for ``d' < d``).

For `d == d′` the original `x` will be returned.

For `d > d′` an error is thrown if the optional parameter `force` is not `true` (default: `false`).
If the `force` flag is enabled, then the value of `x` will be mapped to ``x \\mod d`` regardless of `d` or `d′` (possible loss of information).
"""
function ℤd{T, d}(x::ℤd{T, d′}; force=false) where {T <: Integer, d, d′}
    if d < d′ && !force
        # do not allow conversion from Zd′ → Zd
        throw("Modulus mismatch: can only convert from one subset to a superset! ℤ$d ⊂ ℤ$d′. Use force=true if you want to convert anyways (possible loss of information)")
    elseif d == d′
        # no conversion needed
        return x
    else
        # convert
        return ℤd{T, d}(x.value)
    end
end

# some type aliases, that are not exported by default
"""
    ℤ₂ = ℤd{Bool, 2}
    Z2 = ℤ₂

Type aliases for integers modulo 2.

The type `T` is set to `Bool` instead of `Int8` for performance-reasons.
The simplified operations for the ℤ₂ case can be found at the bottom of operations.jl.
"""
const ℤ₂ = ℤd{Bool, 2}
const Z2 = ℤ₂

"""
    ℤ₃ = ℤd{Int8, 3}
    Z3 = ℤ₃

Type aliases for integers modulo 3.
"""
const ℤ₃ = ℤd{Int8, 3}
const Z3 = ℤ₃

"""
    ℤ₄ = ℤd{Int8, 4}
    Z4 = ℤ₄

Type aliases for integers modulo 4.
"""
const ℤ₄ = ℤd{Int8, 4}
const Z4 = ℤ₄

"""
    ℤ₅ = ℤd{Int8, 5}
    Z5 = ℤ₅

Type aliases for integers modulo 5.
"""
const ℤ₅ = ℤd{Int8, 5}
const Z5 = ℤ₅

"""
    ℤ₆ = ℤd{Int8, 6}
    Z6 = ℤ₆

Type aliases for integers modulo 6.
"""
const ℤ₆ = ℤd{Int8, 6}
const Z6 = ℤ₆


text:
# numbers

function to_N(N::Type{NUM}, v::NUM) where {NUM<:Real}
    return v
end

function to_N(N::Type{<:Real}, v::Real)
    return N(v)
end

function to_N(N::Type{<:Rational{NUM}}, v::Real) where {NUM}
    return rationalize(NUM, v)
end

# arrays

function to_N(N::Type{NUM}, a::AbstractVector{NUM}) where {NUM<:Real}
    return a
end

function to_N(N::Type{<:Real}, a::AbstractVector{<:Real})
    r = Vector{N}(undef, length(a))
    for i in eachindex(a)
        r[i] = to_N(N, a[i])
    end
    return r
end

# matrices

function to_N(N::Type{NUM}, m::AbstractMatrix{NUM}) where {NUM<:Real}
    return m
end

function to_N(N::Type{<:Real}, m::AbstractMatrix{<:Real})
    r = Matrix{N}(undef, size(m))
    for j = 1:size(m,2)
        for i = 1:size(m,1)
            r[i, j] = to_N(N, m[i, j])
        end
    end
    return r
end

# arrays of arrays

function to_N(N::Type{NUM}, a::AbstractVector{<:AbstractVector{NUM}}) where {NUM<:Real}
    return a
end

function to_N(N::Type{<:Real}, a::AbstractVector{<:AbstractVector{<:Real}})
    r = Vector{Vector{N}}(undef, length(a))
    for i in eachindex(a)
        v = a[i]
        b = Vector{N}(undef, length(v))
        for j in eachindex(v)
            b[j] = to_N(N, v[j])
        end
        r[i] = b
    end
    return r
end


text:
function parse_config(path::String)
    cf = init_config()
    f = open(path, "r")
    for i in eachline(f, keep = true)
        if first(i) == '['
            continue
        end
        idx = something(findfirst(isequal('='), i), 0)
        var = rstrip(i[1:idx-1])
        val = strip(i[idx+1:end])
        cf[var] = val
    end
    close(f)
    cf
end

# Defaults
function init_config()
    a = Dict{String, String}()

    a["version"] = "unknown"
    a["connect_four_neighbors_only"] =  "False"
    a["connect_using_avg_resistances"] = "False"
    a["use_polygons"] =  "False"
    a["polygon_file"] = "(Browse for a short-circuit region file)"
    a["source_file"] = "(Browse for a current source file)"
    a["ground_file"]  = "(Browse for a ground point file)"
    a["ground_file_is_resistances"] = "True"
    a["use_unit_currents"] = "False"
    a["use_direct_grounds"] = "False"
    a["remove_src_or_gnd"]  = "keepall"
    a["mask_file"] =  "None"
    a["use_mask"] = "False"
    a["preemptive_memory_release"] = "False"
    a["low_memory_mode"] = "False"
    a["parallelize"] =  "False"
    a["max_parallel"] =  "0"
    a["print_timings"] = "False"
    a["print_rusages"] = "False"
    a["solver"] =  "cg+amg"
    a["use_variable_source_strengths"] = "False"
    a["variable_source_file"] =  "None"
    a["set_null_currents_to_nodata"] = "False"
    a["output_file"] = "(Choose a base name for output files)"
    a["write_cum_cur_map_only"] =  "False"
    a["log_transform_maps"] = "False"
    a["write_max_cur_maps"] = "False"
    a["compress_grids"] = "False"
    a["set_null_voltages_to_nodata"] = "False"
    a["set_focal_node_currents_to_zero"] = "False"
    a["write_volt_maps"] = "False"
    a["write_cur_maps"] = "False"
    a["habitat_map_is_resistances"] = "True"
    a["habitat_file"] = "(Browse for a resistance file)"
    a["scenario"] = "not entered"
    a["data_type"] = "raster"
    a["use_included_pairs"] = "False"
    a["included_pairs_file"] = "(Browse for a file with pairs to include or exclude)"
    a["point_file"] = "(Browse for file with locations of focal points or regions)"
    a["use_reclass_table"] = "False"
    a["reclass_file"] = "(Browse for file with reclassification data)"
    a["profiler_log_file"]  = "None"
    a["log_file"] = "None"
    a["log_level"] = "INFO"
    a["screenprint_log"] = "False"
    a["precision"] = "Double"
    a["log_file"] = "None"
    a["log_level"] = "INFO"
    a["cholmod_batch_size"] = "1000"
    a["use_64bit_indexing"] = "true"

    a
end

function update!(cfg, new)
    for (key,val) in new
        cfg[key] = val
    end
end

function write_config(cfg)
    open(cfg["output_file"], "w") do f
        write(f, """
        [Circuitscape Mode]
        data_type = $(cfg["data_type"])
        scenario = $(cfg["scenario"])

        [Version]
        version = 5.0.0

        [Habitat raster or graph]
        habitat_file = $(cfg["habitat_file"])
        habitat_map_is_resistances = $(cfg["habitat_map_is_resistances"])

        [Connection Scheme for raster habitat data]
        connect_four_neighbors_only = $(cfg["connect_four_neighbors_only"] in TRUELIST)
        connect_using_avg_resistances = $(cfg["connect_using_avg_resistances"] in TRUELIST)

        [Short circuit regions (aka polygons)]
        use_polygons = $(cfg["use_polygons"] in TRUELIST)
        polygon_file = $(cfg["use_polygons"])

        [Options for advanced mode]
        ground_file_is_resistances = $(cfg["ground_file_is_resistances"] in TRUELIST)
        source_file = $(cfg["source_file"])
        remove_src_or_gnd = $(cfg["remove_src_or_gnd"])
        ground_file = $(cfg["ground_file"])
        use_unit_currents = $(cfg["use_unit_currents"] in TRUELIST)
        use_direct_grounds = $(cfg["use_direct_grounds"] in TRUELIST)

        [Mask file]
        use_mask = $(cfg["use_mask"] in TRUELIST)
        mask_file = $(cfg["mask_file"])

        [Options for one-to-all and all-to-one modes]
        use_variable_source_strengths = $(cfg["use_variable_source_strengths"] in TRUELIST)
        variable_source_file = $(cfg["variable_source_file"])

        [Options for pairwise and one-to-all and all-to-one modes]
        included_pairs_file = $(cfg["included_pairs_file"])
        use_included_pairs = $(cfg["use_included_pairs"] in TRUELIST)
        point_file = $(cfg["point_file"])

        [Calculation options]
        solver = cg+amg

        [Output options]
        write_cum_cur_map_only = $(cfg["write_cum_cur_map_only"] in TRUELIST)
        log_transform_maps = $(cfg["log_transform_maps"] in TRUELIST)
        output_file = $(cfg["output_file"])
        write_max_cur_maps = $(cfg["write_max_cur_maps"] in TRUELIST)
        write_volt_maps = $(cfg["write_volt_maps"] in TRUELIST)
        set_null_currents_to_nodata = $(cfg["set_null_currents_to_nodata"] in TRUELIST)
        set_null_voltages_to_nodata = $(cfg["set_null_voltages_to_nodata"] in TRUELIST)
        compress_grids = $(cfg["compress_grids"] in TRUELIST)
        write_cur_maps = $(cfg["write_cur_maps"] in TRUELIST)
        """)
    end
end


text:
module BinaryIO

using Formatting

export readBinary!, writeBinary!

function readBinary!(
    filename::AbstractString,
    arr::AbstractArray{Float64},
    buffer::AbstractArray{UInt8};
    endianess::Symbol=:little_endian,
    delete::Bool=false
)
    local func, need_swap, actual_nbread

    if length(arr) * 8 != length(buffer)
        throw(ErrorException("Length of buffer should be exactly 8 times of data array's."))
    end

    nbread = length(buffer)

    open(filename, "r") do io
        actual_nbread = readbytes!(io, buffer)
    end
    if nbread != actual_nbread
        throw(ErrorException(format("Number of bytes read: {:d} does not equal to {:d}", actual_nbread, nbread)))
    end

    arr[:] = reinterpret(Float64, buffer)

    if endianess == :little_endian
        func = ltoh
    elseif endianess == :big_endian
        func = ntoh
    else
        throw(ErrorException("Unknown symbol: " * String(endianess)))
    end

    arr[:] = func.(arr)

    if delete
        rm(filename, force=true)
    end
end

function writeBinary!(
    filename::AbstractString,
    arr::AbstractArray{Float64},
    buffer::AbstractArray{UInt8};
    endianess::Symbol=:little_endian
)
    local func

    if length(arr) * 8 != length(buffer)
        throw(ErrorException("Length of buffer should be exactly 8 times of data array's."))
    end

    nbwrite = length(buffer)

    if endianess == :little_endian
        func = htol
    elseif endianess == :big_endian
        func = hton
    else
        throw(ErrorException("Unknown symbol: " * String(endianess)))
    end

    buffer[:] = reinterpret(UInt8, func.(arr))
    actual_nbwrite = write(filename, buffer)

    if nbwrite != actual_nbwrite
        throw(ErrorException(format("Number of bytes write: {:d} does not equal to {:d}", actual_nbwrite, nbwrite)))
    end

end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Test
using ModelingToolkit, OrdinaryDiffEq

include("../examples/rc_model.jl")

sys = structural_simplify(rc_model)
@test !isempty(ModelingToolkit.defaults(sys))
u0 = [
      capacitor.v => 0.0
      capacitor.p.i => 0.0
      resistor.v => 0.0
     ]
prob = ODEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Rodas4())

@test sol[resistor.p.i] == sol[capacitor.p.i]
@test sol[resistor.n.i] == -sol[capacitor.p.i]
@test sol[capacitor.n.i] == -sol[capacitor.p.i]
@test iszero(sol[ground.g.i])
@test iszero(sol[ground.g.v])
@test sol[resistor.v] == sol[source.p.v] - sol[capacitor.p.v]

u0 = [
      capacitor.v => 0.0
     ]
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())

@test sol[resistor.p.i] == sol[capacitor.p.i]
@test sol[resistor.n.i] == -sol[capacitor.p.i]
@test sol[capacitor.n.i] == -sol[capacitor.p.i]
@test iszero(sol[ground.g.i])
@test iszero(sol[ground.g.v])
@test sol[resistor.v] == sol[source.p.v] - sol[capacitor.p.v]
#using Plots
#plot(sol)

include("../examples/serial_inductor.jl")
sys = structural_simplify(ll_model)
u0 = [
      inductor1.i => 0.0
      inductor2.i => 0.0
      inductor2.v => 0.0
     ]
prob = ODEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Rodas4())

prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())

@variables t x1(t) x2(t) x3(t) x4(t)
D = Differential(t)
@named sys1_inner = ODESystem([D(x1) ~ x1], t)
@named sys1_partial = compose(ODESystem([D(x2) ~ x2], t; name=:foo), sys1_inner)
@named sys1 = extend(ODESystem([D(x3) ~ x3], t; name=:foo), sys1_partial)
@named sys2 = compose(ODESystem([D(x4) ~ x4], t; name=:foo), sys1)
@test_nowarn sys2.sys1.sys1_inner.x1 # test the correct nesting

# compose tests
@parameters t

function record_fun(;name)
    pars = @parameters a=10 b=100
    ODESystem(Equation[], t, [], pars; name)
end

function first_model(;name)
    @named foo=record_fun()

    defs = Dict()
    defs[foo.a] = 3
    defs[foo.b] = 300
    pars = @parameters x=2 y=20
    compose(ODESystem(Equation[], t, [], pars; name, defaults=defs), foo)
end
@named goo = first_model()
@unpack foo = goo
@test ModelingToolkit.defaults(goo)[foo.a] == 3
@test ModelingToolkit.defaults(goo)[foo.b] == 300


text:
mutable struct officeRent <: Rent

end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
lattice = load_lattice(joinpath(@__DIR__, "./Mn2RuGa.in"))
atoms   = load_atoms(  joinpath(@__DIR__, "./Mn2RuGa.in"))
atoms   = attach_pseudos(atoms, Mn="hgh/pbe/mn-q15.hgh",
                                Ru="hgh/pbe/ru-q16.hgh",
                                Ga="hgh/pbe/ga-q3.hgh")
magnetic_moments = build_magnetic_moments(atoms; Ga=0.0, Mn=5.0, Ru=5.0)

config[@__DIR__] = (
    directory        =  @__DIR__,
    lattice          = lattice,
    atoms            = atoms,
    magnetic_moments = magnetic_moments,
    #
    kwargs_model     = (smearing=Smearing.Gaussian(), temperature=0.01),
    Ecut             = 35,
    kgrid            = [13, 13, 13],
    supersampling    = 2.0,
    potential        = nothing,
    kwargs_scf       = (; ),
    mixing           = KerkerMixing(),
    cases            = standard_cases(),
)

max_stars_repo_name: mfherbst/supporting-adaptive-damping

max_issues_repo_name: mfherbst/supporting-adaptive-damping

max_forks_repo_name: mfherbst/supporting-adaptive-damping


text:
using Statistics
using StatsBase:sample, sample!
using EvoTrees
using BenchmarkTools
using CUDA

# prepare a dataset
features = rand(Int(1.25e6), 100)
# features = rand(100, 10)
X = features
Y = rand(size(X, 1))
𝑖 = collect(1:size(X, 1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace=false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1)) + 1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

###################################################
# GPU
###################################################
params_g = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:none,
    nrounds=100,
    λ=1.0, γ=0.1, η=0.1,
    max_depth=2, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=64);

model_g, cache_g = EvoTrees.init_evotree_gpu(params_g, X_train, Y_train);

params_g = model_g.params;
X_size = size(cache_g.X_bin);

# select random rows and cols
𝑖c = cache_g.𝑖_[sample(params_g.rng, cache_g.𝑖_, ceil(Int, params_g.rowsample * X_size[1]), replace=false, ordered=true)]
𝑖 = CuVector(𝑖c)
𝑗c = cache_g.𝑗_[sample(params_g.rng, cache_g.𝑗_, ceil(Int, params_g.colsample * X_size[2]), replace=false, ordered=true)]
𝑗 = CuVector(𝑗c)

cache_g.nodes[1].𝑖 = 𝑖
cache_g.𝑗 .= 𝑗
# build a new tree
# 144.600 μs (23 allocations: 896 bytes) - 5-6 X time faster on GPU
@time CUDA.@sync EvoTrees.update_grads_gpu!(params_g.loss, cache_g.δ𝑤, cache_g.pred_gpu, cache_g.Y_gpu)
# @btime CUDA.@sync EvoTrees.update_grads_gpu!($params_g.loss, $cache_g.δ𝑤, $cache_g.pred_gpu, $cache_g.Y_gpu)
# sum Gradients of each of the K parameters and bring to CPU

# 33.447 ms (6813 allocations: 307.27 KiB)
tree = EvoTrees.TreeGPU(params_g.max_depth, model_g.K, params_g.λ)
CUDA.@time EvoTrees.grow_tree_gpu!(tree, cache_g.nodes, params_g, cache_g.δ𝑤, cache_g.edges, cache_g.𝑗, cache_g.out, cache_g.left, cache_g.right, cache_g.X_bin, cache_g.K)
CUDA.@time EvoTrees.grow_tree_gpu!(tree, cache_g.nodes, params_g, cache_g.δ𝑤, cache_g.edges, cache_g.𝑗, cache_g.out, cache_g.left, cache_g.right, cache_g.X_bin, cache_g.K)

# CUDA.@time EvoTrees.grow_tree_gpu!(tree, params_g, cache_g.δ, cache_g.hist, cache_g.histL, cache_g.histR, cache_g.gains, cache_g.edges, 𝑖, 𝑗, 𝑛, cache_g.X_bin);
# @btime EvoTrees.grow_tree_gpu!(EvoTrees.TreeGPU(UInt32($params_g.max_depth), $model_g.K, $params_g.λ), $params_g, $cache_g.δ, $cache_g.hist, $cache_g.histL, $cache_g.histR, $cache_g.gains, $cache_g.edges, $𝑖, $𝑗, $𝑛, $cache_g.X_bin);
# @code_warntype EvoTrees.grow_tree_gpu!(EvoTrees.TreeGPU(params_g.max_depth, model_g.K, params_g.λ), params_g, cache_g.δ, cache_g.hist, cache_g.histL, cache_g.histR, cache_g.gains, cache_g.edges, 𝑖, 𝑗, 𝑛, cache_g.X_bin);

# push!(model_g.trees, tree);
# # 2.736 ms (93 allocations: 13.98 KiB)
# @time CUDA.@sync EvoTrees.predict_gpu!(cache_g.pred_gpu, tree, cache_g.X_bin)
# @btime CUDA.@sync EvoTrees.predict_gpu!($cache_g.pred_gpu, $tree, $cache_g.X_bin)

###########################
# Tree GPU
###########################
δ𝑤, K, edges, X_bin, nodes, out, left, right = cache_g.δ𝑤, cache_g.K, cache_g.edges, cache_g.X_bin, cache_g.nodes, cache_g.out, cache_g.left, cache_g.right;

# 9.613 ms (81 allocations: 13.55 KiB)
# 𝑗2 = CuArray(sample(UInt32.(1:100), 50, replace=false, ordered=true))
# @time EvoTrees.update_hist_gpu!(params_g.loss, nodes[1].h, δ𝑤, X_bin, 𝑖, 𝑗, K)
# println(nodes[1].h)
CUDA.@time EvoTrees.update_hist_gpu!(params_g.loss, nodes[1].h, δ𝑤, X_bin, 𝑖, 𝑗, K)
CUDA.@time EvoTrees.update_hist_gpu!(params_g.loss, nodes[1].h, δ𝑤, X_bin, 𝑖, 𝑗, K)

# @btime EvoTrees.update_hist_gpu!($params_g.loss, $nodes[1].h, $δ𝑤, $X_bin, $𝑖, $𝑗, $K)
# @btime EvoTrees.update_hist_gpu!($nodes[1].h, $δ𝑤, $X_bin, $nodes[1].𝑖, $𝑗)
# @code_warntype EvoTrees.update_hist_gpu!(hist, δ, X_bin, 𝑖, 𝑗, 𝑛)

# depth=1
# nid = 2^(depth - 1):2^(depth) - 1
# # 97.000 μs (159 allocations: 13.09 KiB)
# @time CUDA.@sync EvoTrees.update_gains_gpu!(gains::AbstractArray{T,3}, hist::AbstractArray{T,4}, histL::AbstractArray{T,4}, histR::AbstractArray{T,4}, 𝑗::AbstractVector{S}, params_g, nid, depth);
# @btime CUDA.@sync EvoTrees.update_gains_gpu!(gains::AbstractArray{T,3}, hist::AbstractArray{T,4}, histL::AbstractArray{T,4}, histR::AbstractArray{T,4}, 𝑗::AbstractVector{S}, params_g, nid, depth);
# gains[:,:,1]

# tree = EvoTrees.TreeGPU(UInt32(params_g.max_depth), model_g.K, params_g.λ)
# n = 1
# best = findmax(view(gains, :,:,n))
# if best[2][1] != params_g.nbins && best[1] > -Inf
#     tree.gain[n] = best[1]
#     tree.feat[n] = best[2][2]
#     tree.cond_bin[n] = best[2][1]
#     tree.cond_float[n] = edges[tree.feat[n]][tree.cond_bin[n]]
# end
# tree.split[n] = tree.cond_bin[n] != 0

# # 673.900 μs (600 allocations: 29.39 KiB)
# @time CUDA.@sync EvoTrees.update_set_gpu!(𝑛, 𝑖, X_bin, tree.feat, tree.cond_bin, params_g.nbins)
# @btime CUDA.@sync EvoTrees.update_set_gpu!($𝑛, $𝑖, $X_bin, $tree.feat, $tree.cond_bin, $params_g.nbins)


text:
module TestDeprecated

using Test, DataFrames, CategoricalArrays

const ≅ = isequal

@testset "DataFrame!" begin
    x = [1,2,3]
    y = [4,5,6]
    @test DataFrame!(x=x, y=y, copycols=true) == DataFrame(x=x,y=y)
    df1 = DataFrame(x=x, y=y)
    df2 = DataFrame!(df1)
    @test df1 == df2
    @test df1.x === df2.x
    @test df1.y === df2.y

    a=[1,2,3]
    df = DataFrame!(:a=>a, :b=>1, :c=>1:3)
    @test propertynames(df) == [:a, :b, :c]
    @test df.a === a

    df = DataFrame!("a"=>a, "b"=>1, "c"=>1:3)
    @test propertynames(df) == [:a, :b, :c]
    @test df."a" === a

    df = DataFrame!(Dict(:a=>a, :b=>1, :c=>1:3))
    @test propertynames(df) == [:a, :b, :c]
    @test df.a === a

    df = DataFrame!(Dict("a"=>a, "b"=>1, "c"=>1:3))
    @test propertynames(df) == [:a, :b, :c]
    @test df."a" === a

    df = DataFrame!((x, y))
    @test propertynames(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame!((x, y), (:x1, :x2))
    @test propertynames(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame!((x, y), ("x1", "x2"))
    @test names(df) == ["x1", "x2"]
    @test df."x1" === x
    @test df."x2" === y

    @test_throws MethodError DataFrame!([1 2; 3 4], copycols=false)
    @test_throws MethodError DataFrame!([1 2; 3 4])
    @test_throws MethodError DataFrame!([Union{Int, Missing}, Union{Float64, Missing}],
                                        [:x1, :x2], 2)
end

@testset "test categorical" begin
    df = DataFrame(x=["a", "b", "c"],
                   y=["a", "b", missing],
                   z=[1,2,3])
    for x in [df, view(df, :, :)]
        y = categorical(x)
        @test y isa DataFrame
        @test x ≅ y
        @test x.x !== y.x
        @test x.y !== y.y
        @test x.z !== y.z
        @test y.x isa CategoricalVector{String}
        @test y.y isa CategoricalVector{Union{Missing, String}}
        @test y.z isa Vector{Int}

        y = categorical(x, Int)
        @test y isa DataFrame
        @test x ≅ y
        @test x.x !== y.x
        @test x.y !== y.y
        @test x.z !== y.z
        @test y.x isa Vector{String}
        @test y.y isa Vector{Union{Missing, String}}
        @test y.z isa CategoricalVector{Int}

        for colsel in [:, names(x), [1,2,3], [true,true,true], r"", Not(r"a")]
            y = categorical(x, colsel)
            @test y isa DataFrame
            @test x ≅ y
            @test x.x !== y.x
            @test x.y !== y.y
            @test x.z !== y.z
            @test y.x isa CategoricalVector{String}
            @test y.y isa CategoricalVector{Union{Missing, String}}
            @test y.z isa CategoricalVector{Int}
        end

        for colsel in [:x, "x", 1, [:x], ["x"], [1], [true, false, false], r"x", Not(2:3)]
            y = categorical(x, colsel)
            @test y isa DataFrame
            @test x ≅ y
            @test x.x !== y.x
            @test x.y !== y.y
            @test x.z !== y.z
            @test y.x isa CategoricalVector{String}
            @test y.y isa Vector{Union{Missing, String}}
            @test y.z isa Vector{Int}
        end

        for colsel in [:z, "z", 3, [:z], ["z"], [3], [false, false, true], r"z", Not(1:2)]
            y = categorical(x, colsel)
            @test y isa DataFrame
            @test x ≅ y
            @test x.x !== y.x
            @test x.y !== y.y
            @test x.z !== y.z
            @test y.x isa Vector{String}
            @test y.y isa Vector{Union{Missing, String}}
            @test y.z isa CategoricalVector{Int}
        end

        for colsel in [Int[], Symbol[], [false, false, false], r"a", Not(:)]
            y = categorical(x, colsel)
            @test y isa DataFrame
            @test x ≅ y
            @test x.x !== y.x
            @test x.y !== y.y
            @test x.z !== y.z
            @test y.x isa Vector{String}
            @test y.y isa Vector{Union{Missing, String}}
            @test y.z isa Vector{Int}
        end
    end
end

@testset "test categorical!" begin
    using DataFrames: _columns
    df = DataFrame(A = Vector{Union{Int, Missing}}(1:3), B = Vector{Union{Int, Missing}}(4:6))
    DRT = CategoricalArrays.DefaultRefType
    @test all(c -> isa(c, Vector{Union{Int, Missing}}), eachcol(categorical!(deepcopy(df))))
    @test all(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
              eachcol(categorical!(deepcopy(df), [1,2])))
    @test all(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
              eachcol(categorical!(deepcopy(df), [:A,:B])))
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), [:A]))) == 1
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), :A))) == 1
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), [1]))) == 1
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), 1))) == 1

    @test all(c -> isa(c, Vector{Union{Int, Missing}}), eachcol(categorical!(deepcopy(df))))
    @test all(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
              eachcol(categorical!(deepcopy(df), Not(Not([1,2])))))
    @test all(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
              eachcol(categorical!(deepcopy(df), Not(Not([:A,:B])))))
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), Not(Not([:A]))))) == 1
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), Not(Not(:A))))) == 1
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), Not(Not([1]))))) == 1
    @test findfirst(c -> typeof(c) <: CategoricalVector{Union{Int, Missing}},
                    _columns(categorical!(deepcopy(df), Not(Not(1))))) == 1
end

@testset "categorical!" begin
    df = DataFrame([["a", "b"], ['a', 'b'], [true, false], 1:2, ["x", "y"]])
    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df)))),
                  [CategoricalArrays.CategoricalValue{String,UInt32},
                   Char, Bool, Int,
                   CategoricalArrays.CategoricalValue{String,UInt32}]))
    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df), :))),
                  [CategoricalArrays.CategoricalValue{String,UInt32},
                   CategoricalArrays.CategoricalValue{Char,UInt32},
                   CategoricalArrays.CategoricalValue{Bool,UInt32},
                   CategoricalArrays.CategoricalValue{Int,UInt32},
                   CategoricalArrays.CategoricalValue{String,UInt32}]))
    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df), compress=true))),
                  [CategoricalArrays.CategoricalValue{String,UInt8},
                   Char, Bool, Int,
                   CategoricalArrays.CategoricalValue{String,UInt8}]))
    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df), names(df)))),
                  [CategoricalArrays.CategoricalValue{String,UInt32},
                   CategoricalArrays.CategoricalValue{Char,UInt32},
                   CategoricalArrays.CategoricalValue{Bool,UInt32},
                   CategoricalArrays.CategoricalValue{Int,UInt32},
                   CategoricalArrays.CategoricalValue{String,UInt32}]))
    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df), names(df), compress=true))),
                  [CategoricalArrays.CategoricalValue{String,UInt8},
                   CategoricalArrays.CategoricalValue{Char,UInt8},
                   CategoricalArrays.CategoricalValue{Bool,UInt8},
                   CategoricalArrays.CategoricalValue{Int,UInt8},
                   CategoricalArrays.CategoricalValue{String,UInt8}]))
    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df), Not(1:0)))),
                  [CategoricalArrays.CategoricalValue{String,UInt32},
                   CategoricalArrays.CategoricalValue{Char,UInt32},
                   CategoricalArrays.CategoricalValue{Bool,UInt32},
                   CategoricalArrays.CategoricalValue{Int,UInt32},
                   CategoricalArrays.CategoricalValue{String,UInt32}]))
    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df), Not(1:0), compress=true))),
                  [CategoricalArrays.CategoricalValue{String,UInt8},
                   CategoricalArrays.CategoricalValue{Char,UInt8},
                   CategoricalArrays.CategoricalValue{Bool,UInt8},
                   CategoricalArrays.CategoricalValue{Int,UInt8},
                   CategoricalArrays.CategoricalValue{String,UInt8}]))

    @test all(map(<:, eltype.(eachcol(categorical!(deepcopy(df), Integer))),
                  [String, Char,
                   CategoricalArrays.CategoricalValue{Bool,UInt32},
                   CategoricalArrays.CategoricalValue{Int,UInt32},
                   String]))

    df = DataFrame([["a", missing]])
    categorical!(df)
    @test df.x1 isa CategoricalVector{Union{Missing, String}}

    df = DataFrame(x1=[1, 2])
    categorical!(df)
    @test df.x1 isa Vector{Int}
    categorical!(df, :)
    @test df.x1 isa CategoricalVector{Int}
end

@testset "categorical with All and Between" begin
    df = DataFrame(x1=["a", "b"], y=[2, 3])
    categorical(df, All())
    categorical(df, Between(1, 2))
    categorical!(df, All())
    categorical!(df, Between(1,2))
end

@testset "deprecated describe syntax" begin
    @test describe(DataFrame(a=[1,2]), cols = :a, :min, :min2 => minimum, "max2" => maximum, :max) ==
          DataFrame(variable=:a, min=1, min2=1, max2=2, max=2)
end

end # module


text:
using FastArrayOps
using BenchmarkLite
using ArrayViews

blas_set_num_threads(1)
type BenchCase{Op, FT} <: Proc end

abstract MethodType
abstract FuncType
type BM_TEST end
type BM_JBase1 end
type BM_JBase2 end
type BM_BLAS end
type BM_Broadcast end
type BM_Forloop end
type BM_FAO end

# procedure name
Base.string{Op, FT}(::BenchCase{Op, FT}) = string("$Op"[4:end])
Base.string{FT<:FuncType}(::FT) = lowercase(string("$FT"))
# procedure codes
Base.length(p::BenchCase, n::Int) = n
Base.isvalid(p::BenchCase, n::Int) = (n > 0)
Base.done(p::BenchCase, n, s) = nothing

function runFAO(procs::Vector{Proc}, cfgs::Vector;
             nruns::Int=0, verbose::Int=2, allowgc::Bool=true, logger::IO=STDOUT)

    bt = BenchmarkTable(cfgs, procs)
    m = length(cfgs)
    n = length(procs)
    for (j, p) in enumerate(procs)
        procname = string(p)
        verbose >= 1 && println(logger, "Benchmarking $procname ...")

        for (i, cfg) in enumerate(cfgs)
            cfgname = string(cfg)
            (nr, et) = run(p, cfg; nruns=2, allowgc=allowgc) #warm up
            runs = max(div(nruns, cfg), 4)
            (nr, et) = run(p, cfg; nruns=runs, allowgc=allowgc)

            bt.nruns[i, j] = nr
            bt.etime[i, j] = et

            verbose >= 2 &&
                println(logger, "  $procname with cfg = $cfgname: nruns = $nr, elapsed = $et secs")
        end
    end
    return bt
end


text:
# http://web.mit.edu/puzzle/www/2007/puzzles/1_1_1/
@testset "1-1=1" begin
    f = best_feature(wordlist(IOBuffer("""
        STRIFE
        SEAMAN
        NIX
        ETCH
        POST
        QUEERART
        FOO
        TALKS
        REPAYS
        STU
        HUMF
        UNDERHID
        SIXTEENS
        BOWMEN
        """)), 1)
    @test f.description == "has a 1-letter transdeletion"
end


text:
using DbgMacro

using Suppressor: @capture_out

code = quote
    function foo()
        m = [1 2; 3 4]
        @dumpct asd(1 + x, :y, $(Expr(:$, :z)), k=v)
        42
    end
end

out = @capture_out eval(code)

expected_out = """
Expr
  head: Symbol call
  args: Array{Any}((5,))
    1: Symbol asd
    2: Expr
      head: Symbol call
      args: Array{Any}((3,))
        1: Symbol +
        2: Int64 1
        3: Symbol x
    3: QuoteNode
      value: Symbol y
    4: Expr
      head: Symbol \$
      args: Array{Any}((1,))
        1: Symbol z
    5: Expr
      head: Symbol kw
      args: Array{Any}((2,))
        1: Symbol k
        2: Symbol v
"""

@test out == expected_out


text:
# Rows grouping.
# Maps row contents to the indices of all the equal rows.
# Used by groupby(), join(), nonunique()
immutable RowGroupDict{T<:AbstractDataTable}
    "source data table"
    dt::T
    "number of groups"
    ngroups::Int
    "row hashes"
    rhashes::Vector{UInt}
    "hashindex -> index of group-representative row"
    gslots::Vector{Int}
    "group index for each row"
    groups::Vector{Int}
    "permutation of row indices that sorts them by groups"
    rperm::Vector{Int}
    "starts of ranges in rperm for each group"
    starts::Vector{Int}
    "stops of ranges in rperm for each group"
    stops::Vector{Int}
end

# "kernel" functions for hashrows()
# adjust row hashes by the hashes of column elements
function hashrows_col!(h::Vector{UInt}, v::AbstractVector)
    @inbounds for i in eachindex(h)
        h[i] = hash(v[i], h[i])
    end
    h
end

if !isdefined(Base, :unsafe_get)
    unsafe_get(x::Nullable) = x.value
    unsafe_get(x::Any) = x
end

function hashrows_col!{T<:Nullable}(h::Vector{UInt}, v::AbstractVector{T})
    @inbounds for i in eachindex(h)
        h[i] = isnull(v[i]) ?
               h[i] + Base.nullablehash_seed :
               hash(unsafe_get(v[i]), h[i])
    end
    h
end

# should give the same hash as AbstractVector{T}
function hashrows_col!{T}(h::Vector{UInt}, v::AbstractCategoricalVector{T})
    # TODO is it possible to optimize by hashing the pool values once?
    @inbounds for (i, ref) in enumerate(v.refs)
        h[i] = hash(CategoricalArrays.index(v.pool)[ref], h[i])
    end
    h
end

# should give the same hash as AbstractNullableVector{T}
# enables efficient sequential memory access pattern
function hashrows_col!{T}(h::Vector{UInt}, v::AbstractNullableCategoricalVector{T})
    # TODO is it possible to optimize by hashing the pool values once?
    @inbounds for (i, ref) in enumerate(v.refs)
        h[i] = ref == 0 ?
               h[i] + Base.nullablehash_seed :
               hash(CategoricalArrays.index(v.pool)[ref], h[i])
    end
    h
end

# Calculate the vector of `dt` rows hash values.
function hashrows(dt::AbstractDataTable)
    res = zeros(UInt, nrow(dt))
    for col in columns(dt)
        hashrows_col!(res, col)
    end
    return res
end

# Helper function for RowGroupDict.
# Returns a tuple:
# 1) the number of row groups in a data table
# 2) vector of row hashes
# 3) slot array for a hash map, non-zero values are
#    the indices of the first row in a group
# Optional group vector is set to the group indices of each row
row_group_slots(dt::AbstractDataTable, groups::Union{Vector{Int}, Void} = nothing) =
    row_group_slots(ntuple(i -> dt[i], ncol(dt)), hashrows(dt), groups)

function row_group_slots(cols::Tuple{Vararg{AbstractVector}},
                         rhashes::AbstractVector{UInt},
                         groups::Union{Vector{Int}, Void} = nothing)
    @assert groups === nothing || length(groups) == length(cols[1])
    # inspired by Dict code from base cf. https://github.com/JuliaData/DataTables.jl/pull/17#discussion_r102481481
    sz = Base._tablesz(length(rhashes))
    @assert sz >= length(rhashes)
    szm1 = sz-1
    gslots = zeros(Int, sz)
    ngroups = 0
    @inbounds for i in eachindex(rhashes)
        # find the slot and group index for a row
        slotix = rhashes[i] & szm1 + 1
        gix = 0
        probe = 0
        while true
            g_row = gslots[slotix]
            if g_row == 0 # unoccupied slot, current row starts a new group
                gslots[slotix] = i
                gix = ngroups += 1
                break
            elseif rhashes[i] == rhashes[g_row] # occupied slot, check if miss or hit
                if isequal_row(cols, i, g_row) # hit
                    gix = groups !== nothing ? groups[g_row] : 0
                end
                break
            end
            slotix = slotix & szm1 + 1 # check the next slot
            probe += 1
            @assert probe < sz
        end
        if groups !== nothing
            groups[i] = gix
        end
    end
    return ngroups, rhashes, gslots
end

# Builds RowGroupDict for a given datatable.
# Partly uses the code of Wes McKinney's groupsort_indexer in pandas (file: src/groupby.pyx).
function group_rows(dt::AbstractDataTable)
    groups = Vector{Int}(nrow(dt))
    ngroups, rhashes, gslots = row_group_slots(dt, groups)

    # count elements in each group
    stops = zeros(Int, ngroups)
    @inbounds for g_ix in groups
        stops[g_ix] += 1
    end

    # group start positions in a sorted table
    starts = Vector{Int}(ngroups)
    if !isempty(starts)
        starts[1] = 1
        @inbounds for i in 1:(ngroups-1)
            starts[i+1] = starts[i] + stops[i]
        end
    end

    # define row permutation that sorts them into groups
    rperm = Vector{Int}(length(groups))
    copy!(stops, starts)
    @inbounds for (i, gix) in enumerate(groups)
        rperm[stops[gix]] = i
        stops[gix] += 1
    end
    stops .-= 1
    return RowGroupDict(dt, ngroups, rhashes, gslots, groups, rperm, starts, stops)
end

# Find index of a row in gd that matches given row by content, 0 if not found
function findrow(gd::RowGroupDict,
                 dt::DataTable,
                 gd_cols::Tuple{Vararg{AbstractVector}},
                 dt_cols::Tuple{Vararg{AbstractVector}},
                 row::Int)
    (gd.dt === dt) && return row # same table, return itself
    # different tables, content matching required
    rhash = rowhash(dt_cols, row)
    szm1 = length(gd.gslots)-1
    slotix = ini_slotix = rhash & szm1 + 1
    while true
        g_row = gd.gslots[slotix]
        if g_row == 0 || # not found
            (rhash == gd.rhashes[g_row] &&
            isequal_row(gd_cols, g_row, dt_cols, row)) # found
            return g_row
        end
        slotix = (slotix & szm1) + 1 # miss, try the next slot
        (slotix == ini_slotix) && break
    end
    return 0 # not found
end

# Find indices of rows in 'gd' that match given row by content.
# return empty set if no row matches
function findrows(gd::RowGroupDict,
                  dt::DataTable,
                  gd_cols::Tuple{Vararg{AbstractVector}},
                  dt_cols::Tuple{Vararg{AbstractVector}},
                  row::Int)
    g_row = findrow(gd, dt, gd_cols, dt_cols, row)
    (g_row == 0) && return view(gd.rperm, 0:-1)
    gix = gd.groups[g_row]
    return view(gd.rperm, gd.starts[gix]:gd.stops[gix])
end

function Base.getindex(gd::RowGroupDict, dtr::DataTableRow)
    g_row = findrow(gd, dtr.dt, ntuple(i -> gd.dt[i], ncol(gd.dt)),
                    ntuple(i -> dtr.dt[i], ncol(dtr.dt)), dtr.row)
    (g_row == 0) && throw(KeyError(dtr))
    gix = gd.groups[g_row]
    return view(gd.rperm, gd.starts[gix]:gd.stops[gix])
end


text:
export dbscan

using mlpack._Internal.io

import mlpack_jll
const dbscanLibrary = mlpack_jll.libmlpack_julia_dbscan

# Call the C binding of the mlpack dbscan binding.
function dbscan_mlpackMain()
  success = ccall((:dbscan, dbscanLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module dbscan_internal
  import ..dbscanLibrary

end # module

"""
    dbscan(input; [epsilon, min_size, naive, selection_type, single_mode, tree_type, verbose])

This program implements the DBSCAN algorithm for clustering using accelerated
tree-based range search.  The type of tree that is used may be parameterized, or
brute-force range search may also be used.

The input dataset to be clustered may be specified with the `input` parameter;
the radius of each range search may be specified with the `epsilon` parameters,
and the minimum number of points in a cluster may be specified with the
`min_size` parameter.

The `assignments` and `centroids` output parameters may be used to save the
output of the clustering. `assignments` contains the cluster assignments of each
point, and `centroids` contains the centroids of each cluster.

The range search may be controlled with the `tree_type`, `single_mode`, and
`naive` parameters.  `tree_type` can control the type of tree used for range
search; this can take a variety of values: 'kd', 'r', 'r-star', 'x',
'hilbert-r', 'r-plus', 'r-plus-plus', 'cover', 'ball'. The `single_mode`
parameter will force single-tree search (as opposed to the default dual-tree
search), and '`naive` will force brute-force range search.

An example usage to run DBSCAN on the dataset in `input` with a radius of 0.5
and a minimum cluster size of 5 is given below:

```julia
julia> using CSV
julia> input = CSV.read("input.csv")
julia> _, _ = dbscan(input; epsilon=0.5, min_size=5)
```

# Arguments

 - `input::Array{Float64, 2}`: Input dataset to cluster.
 - `epsilon::Float64`: Radius of each range search.  Default value `1`.

 - `min_size::Int`: Minimum number of points for a cluster.  Default value
      `5`.

 - `naive::Bool`: If set, brute-force range search (not tree-based) will
      be used.  Default value `false`.

 - `selection_type::String`: If using point selection policy, the type of
      selection to use ('ordered', 'random').  Default value `ordered`.

 - `single_mode::Bool`: If set, single-tree range search (not dual-tree)
      will be used.  Default value `false`.

 - `tree_type::String`: If using single-tree or dual-tree search, the type
      of tree to use ('kd', 'r', 'r-star', 'x', 'hilbert-r', 'r-plus',
      'r-plus-plus', 'cover', 'ball').  Default value `kd`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `assignments::Array{Int, 1}`: Output matrix for assignments of each
      point.
 - `centroids::Array{Float64, 2}`: Matrix to save output centroids to.

"""
function dbscan(input;
                epsilon::Union{Float64, Missing} = missing,
                min_size::Union{Int, Missing} = missing,
                naive::Union{Bool, Missing} = missing,
                selection_type::Union{String, Missing} = missing,
                single_mode::Union{Bool, Missing} = missing,
                tree_type::Union{String, Missing} = missing,
                verbose::Union{Bool, Missing} = missing,
                points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, dbscanLibrary), Nothing, ());

  IORestoreSettings("DBSCAN clustering")

  # Process each input argument before calling mlpackMain().
  IOSetParamMat("input", input, points_are_rows)
  if !ismissing(epsilon)
    IOSetParam("epsilon", convert(Float64, epsilon))
  end
  if !ismissing(min_size)
    IOSetParam("min_size", convert(Int, min_size))
  end
  if !ismissing(naive)
    IOSetParam("naive", convert(Bool, naive))
  end
  if !ismissing(selection_type)
    IOSetParam("selection_type", convert(String, selection_type))
  end
  if !ismissing(single_mode)
    IOSetParam("single_mode", convert(Bool, single_mode))
  end
  if !ismissing(tree_type)
    IOSetParam("tree_type", convert(String, tree_type))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("assignments")
  IOSetPassed("centroids")
  # Call the program.
  dbscan_mlpackMain()

  return IOGetParamURow("assignments"),
         IOGetParamMat("centroids", points_are_rows)
end


text:
#!/usr/bin/env julia

clones_dir = "clones"
static_dir = "static"
get_old_package_artifacts = false

import Pkg
import Pkg.TOML
import Pkg.Artifacts: download_artifact, artifact_path
import LibGit2

mkpath(clones_dir)
mkpath(static_dir)

const tar_opts = ```
    --format=posix
    --numeric-owner
    --owner=0
    --group=0
    --mode=go-w,+X
    --mtime=1970-01-01
    --pax-option=exthdr.name=%d/PaxHeaders/%f,delete=atime,delete=ctime,delete=mtime
    --no-recursion
```
# reproducible tarball options based on
# http://h2.jaguarpaw.co.uk/posts/reproducible-tar/

const compress = `gzip -9`
const decompress = `gzcat`

function make_tarball(
    tarball::AbstractString,
    tree_path::AbstractString,
)
    paths = String[]
    for (root, dirs, files) in walkdir(tree_path)
        path = root != tree_path ? relpath(root, tree_path) : ""
        for file in [dirs; files]
            push!(paths, joinpath(path, file))
        end
    end
    sort!(paths)
    mktemp() do paths_file, io
        for path in paths
            print(io, "$path\0")
        end
        close(io)
        open(tarball, write=true) do io
            tar_cmd = `gtar $tar_opts -cf - -C $tree_path --null -T $paths_file`
            run(pipeline(tar_cmd, compress, io))
        end
    end
    return
end

function create_git_tarball(
    tarball::AbstractString,
    repo_path::AbstractString,
    tree_hash::AbstractString,
)
    repo = LibGit2.GitRepo(repo_path)
    tree = LibGit2.GitObject(repo, tree_hash)
    tree_path = mktempdir()
    opts = LibGit2.CheckoutOptions(
        checkout_strategy = LibGit2.Consts.CHECKOUT_FORCE,
        target_directory = Base.unsafe_convert(Cstring, tree_path)
    )
    LibGit2.checkout_tree(repo, tree, options=opts)
    make_tarball(tarball, tree_path)
    verify_tarball_hash(tarball, tree_hash)
    return
end

function verify_tarball_hash(
    tarball::AbstractString,
    tree_hash::AbstractString,
)
    local hash
    mktempdir() do tmp_dir
        run(pipeline(`$decompress $tarball`, `tar -C $tmp_dir -x -`))
        hash = bytes2hex(Pkg.GitTools.tree_hash(tmp_dir))
        chmod(tmp_dir, 0o777, recursive=true)
    end
    hash == tree_hash || error("""
        tree hash mismatch:
        - expected: $tree_hash
        - computed: $hash
        """)
    return
end

function process_artifact(info::Dict)
    local tree_hash, tree_path, tarball
    try
        tree_hash = info["git-tree-sha1"]
        tree_sha1 = Pkg.Types.SHA1(tree_hash)
        tarball = joinpath(static_dir, "artifact", tree_hash)
        isfile(tarball) && return
        downloads = info["download"]
        downloads isa Array || (downloads = [downloads])
        for download in downloads
            url = download["url"]
            hash = download["sha256"]
            download_artifact(tree_sha1, url, hash, verbose=true) && break
        end
        tree_path = artifact_path(tree_sha1, honor_overrides=false)
        isdir(tree_path) || error("artifact install failed")
    catch err
        @warn err
    end
    mkpath(dirname(tarball))
    make_tarball(tarball, tree_path)
    verify_tarball_hash(tarball, tree_hash)
    return
end

registries = Dict{String,String}()

for depot in DEPOT_PATH
    depot_regs = joinpath(depot, "registries")
    isdir(depot_regs) || continue
    for reg_dir in readdir(depot_regs, join=true)
        isdir(reg_dir) || continue
        reg_file = joinpath(reg_dir, "Registry.toml")
        isfile(reg_file) || continue
        reg_data = TOML.parsefile(reg_file)
        # generate registry tarball
        let tree_hash = readchomp(`git -C $reg_dir rev-parse 'HEAD^{tree}'`)
            uuid = reg_data["uuid"]
            tarball = joinpath(static_dir, "registry", uuid, tree_hash)
            mkpath(dirname(tarball))
            create_git_tarball(tarball, reg_dir, tree_hash)
            registries[reg_data["uuid"]] = tree_hash
        end
        for (uuid, info) in reg_data["packages"]
            name = info["name"]
            path = info["path"]
            pkg_dir = joinpath(reg_dir, path)
            pkg_info = TOML.parsefile(joinpath(pkg_dir, "Package.toml"))
            versions = TOML.parsefile(joinpath(pkg_dir, "Versions.toml"))
            # generate archive of each version
            static_pkg_dir = joinpath(static_dir, "package", uuid)
            mkpath(static_pkg_dir)
            clone_dir = joinpath(clones_dir, uuid)
            updated = false
            if !isdir(clone_dir)
                try run(`git clone --mirror $pkg_repo $clone_dir`)
                catch err
                    println(stderr, "Cannot clone $name [$uuid]")
                    break
                end
                updated = true
            end
            for (ver, info) in versions
                tree_hash = info["git-tree-sha1"]
                tarball = joinpath(static_pkg_dir, tree_hash)
                if (is_new_tarball = !isfile(tarball))
                    pkg_repo = pkg_info["repo"]
                    @label again
                    try create_git_tarball(tarball, clone_dir, tree_hash)
                    catch err
                        if updated
                            println(stderr, "Cannot checkout $name [$uuid] $tree_hash")
                            rm(tarball, force=true)
                            continue
                        end
                        updated = true
                        try run(`git -C $clone_dir remote update`)
                        catch err
                            println(stderr, "Cannot update $name [$uuid]")
                            continue
                        end
                        @goto again
                    end
                    isfile(tarball) || continue
                end
                is_new_tarball || get_old_package_artifacts || continue
                # look for artifact files
                for path in eachline(pipeline(`$decompress $tarball`, `gtar -t`))
                    # NOTE: the above can't handle paths with newlines
                    # doesn't seem to be a way to get tar to use \0 instead
                    basename(path) in Pkg.Artifacts.artifact_names || continue
                    extract = pipeline(`$decompress $tarball`, `gtar -x -O $path`)
                    artifacts = TOML.parse(read(extract, String))
                    for (key, val) in artifacts
                        if val isa Dict
                            process_artifact(val)
                        elseif val isa Vector
                            foreach(process_artifact, val)
                        else
                            @warn "invalid artifact file entry: $val"
                        end
                    end
                end
            end
            isempty(readdir(static_pkg_dir)) && rm(static_pkg_dir)
        end
    end
end

# generate current registries file
open(joinpath(static_dir, "registries"), write=true) do io
    for (uuid, tree_hash) in sort!(collect(registries))
        println(io, "/registry/$uuid/$tree_hash")
    end
end


text:
# Note: Terminology overall taken from:
# https://en.wikipedia.org/wiki/Interval_(mathematics)#Terminology

"""
    Bound <: Any

Abstract type representing all possible endpoint classifications (e.g. open, closed,
unbounded).
"""
:Bound

"""
    Bounded <: Bound

Abstract type indicating that the endpoint of an interval is not unbounded (e.g. open or
closed).
"""
:Bounded

"""
    Closed <: Bounded <: Bound

Type indicating that the endpoint of an interval is closed (the endpoint value is *included*
in the interval).
"""
:Closed

"""
    Open <: Bounded <: Bound

Type indicating that the endpoint of an interval is open (the endpoint value is *not
included* in the interval).
"""
:Open

"""
    Unbounded <: Bound

Type indicating that the endpoint of an interval is unbounded (the endpoint value is
effectively infinite).
"""
:Unbounded

"""
    LeftEndpoint <: Endpoint

Represents the lower endpoint of an `AbstractInterval`. Useful for comparing two endpoints
to each other.

### Examples

```jldoctest; setup = :(using Intervals; using Intervals: LeftEndpoint)
julia> LeftEndpoint(Interval(0.0, 1.0))
Intervals.Endpoint{Float64, Intervals.Direction{:Left}(), Closed}(0.0)

julia> LeftEndpoint{Closed}(1.0)
Intervals.Endpoint{Float64, Intervals.Direction{:Left}(), Closed}(1.0)

julia> LeftEndpoint{Closed}(1) < LeftEndpoint{Closed}(2)
true

julia> LeftEndpoint{Closed}(0) < LeftEndpoint{Open}(0)
true

julia> LeftEndpoint{Open}(0) <= LeftEndpoint{Closed}(0)
false
```

See also: [`RightEndpoint`](@ref)
"""
:LeftEndpoint

"""
    RightEndpoint <: Endpoint

Represents the upper endpoint of an `AbstractInterval`. Useful for comparing two endpoints
to each other.

### Examples

```jldoctest; setup = :(using Intervals; using Intervals: RightEndpoint)
julia> RightEndpoint(Interval(0.0, 1.0))
Intervals.Endpoint{Float64, Intervals.Direction{:Right}(), Closed}(1.0)

julia> RightEndpoint{Closed}(1.0)
Intervals.Endpoint{Float64, Intervals.Direction{:Right}(), Closed}(1.0)

julia> RightEndpoint{Closed}(1) < RightEndpoint{Closed}(2)
true

julia> RightEndpoint{Open}(0) < RightEndpoint{Closed}(0)
true

julia> RightEndpoint{Closed}(0) <= RightEndpoint{Open}(0)
false
```

See also: [`LeftEndpoint`](@ref)
"""
:RightEndpoint

"""
    first(interval::AbstractInterval{T}) -> Union{T,Nothing}

The value of the lower endpoint. When the lower endpoint is unbounded `nothing` will be
returned.
"""
Base.first(::AbstractInterval)

"""
    last(interval::AbstractInterval{T}) -> Union{T,Nothing}

The value of the upper endpoint. When the upper endpoint is unbounded `nothing` will be
returned.
"""
Base.last(::AbstractInterval)

"""
    span(interval::AbstractInterval) -> Any

The delta between the upper and lower endpoints. For bounded intervals returns a
non-negative value while intervals with any unbounded endpoints will throw an
`ArgumentError`.

To avoid having to capture the exception use the pattern:
```julia
Intervals.isbounded(interval) ? span(interval) : infinity
```
Where `infinity` is a variable representing the value you wish to use to represent an
unbounded, or infinite, span.
"""
span(::AbstractInterval)

"""
    isclosed(interval) -> Bool

Is a closed-interval: includes both of its endpoints.
"""
isclosed(::AbstractInterval)

"""
    isopen(interval) -> Bool

Is an open-interval: excludes both of its endpoints.
"""
Base.isopen(::AbstractInterval)

"""
    isunbounded(interval) -> Bool

Is an unbounded-interval: unbounded at both ends.
"""
isunbounded(::AbstractInterval)

"""
    isbounded(interval) -> Bool

Is a bounded-interval: either open, closed, left-closed/right-open, or
left-open/right-closed.

Note using `!isbounded` is commonly used to determine if any end of the interval is
unbounded.
"""
isbounded(::AbstractInterval)

"""
    minimum(interval::AbstractInterval{T}; [increment]) -> T

The minimum value contained within the `interval`.

If left-closed, returns `first(interval)`.
If left-open, returns `first(interval) + eps(first(interval))`
If left-unbounded, returns minimum value possible for type `T`.

A `BoundsError` is thrown for empty intervals or when the increment results in a minimum value
not-contained by the interval.
"""
minimum(::AbstractInterval; increment)

"""
    maximum(interval::AbstractInterval{T}; [increment]) -> T

The maximum value contained within the `interval`.

If right-closed, returns `last(interval)`.
If right-open, returns `first(interval) + eps(first(interval))`
If right-unbounded, returns maximum value possible for type `T`.

A `BoundsError` is thrown for empty intervals or when the increment results in a maximum value
not-contained by the interval.
"""
maximum(::AbstractInterval; increment)


text:
# test/runtests.jl

using AVSfldIO
using Test: @test, @testset, detect_ambiguities

list = [
"fld-read.jl"
"fld-write.jl"
"test-io.jl"
]

for file in list
    @testset "$file" begin
        include(file)
    end
end

@test length(detect_ambiguities(AVSfldIO)) == 0

max_stars_repo_name: JeffFessler/AVSfldIO

max_issues_repo_name: JeffFessler/AVSfldIO


text:
using consonance

push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/")
push!(LOAD_PATH, "/n/home07/iczekala/DiskJockey/")
push!(LOAD_PATH, "/pool/scout0/DiskJockey/")

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    # "--opt1"
    # help = "an option with an argument"
    "--run_index", "-r"
    help = "Output run index"
    arg_type = Int
end

parsed_args = parse_args(ARGS, s)
run_index = parsed_args["run_index"]

imgdir = scratchdir * @sprintf("d%02d/", run_index)

# make the directory where RADMC will make all its files exist
if !ispath(imgdir)
    println("Creating ", imgdir)
    mkdir(imgdir)
end

cd(imgdir)
println("Now in ", pwd())

files = ["dustopac.inp", "dustkappa_silicate.inp", "stars.inp", "lines.inp", "molecule_co.inp", "wavelength_micron.inp"]
for file in files
    src = homedir * file
    dst = imgdir * file
    cp(src, dst)
end

cp(homedir * "radmc3d.inp.dust", imgdir * "radmc3d.inp")

using constants
using image
using model

function make_image(pars, id::Int)

    write_dust(pars, "", grid)
    write_model(pars, "", grid)

    # RADMC conventions for inclination and PA
    incl = pars.incl # [deg]
    PA = pars.PA # [deg] Position angle runs counter clockwise

    lam0 = cc/230.538e9 * 1e4 # [microns]

    run(`radmc3d mctherm` |> DevNull)
    run(`radmc3d image incl $incl posang $PA npix $npix lambda $lam0` |> DevNull)

    src = "image.out"
    dst = outdir * @sprintf("dimage%04d.out", id)

    cp(src, dst)

end

const global npix = 256 # number of pixels
const global grid = Grid(100, 32, 0.5, 800, true)
write_grid("", grid)

start = run_index * nframes_per_proc + 1
ids = Int[i for i=start:(start + nframes_per_proc)]

for i in ids
    if i <= nframes
        make_image(pars[i], i)
    else
        break
    end
end

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
const _RXNS_VAR_KEY = :rxnv
const _STOI_CON_KEY = :stoi_conv
const _LB_CON_KEY = :lb_conv
const _UB_CON_KEY = :ub_conv

_length(lp_model) = length(_get_vars(lp_model))

function _up_con_rhs!(lp_model::JuMP.Model, conv_key::Symbol, vals::Vector, cidxs = 1:_length(lp_model))
    @assert length(vals) == length(cidxs)
    conv = lp_model[conv_key]
    for (ci, val) in zip(cidxs, vals)
        JuMP.set_normalized_rhs(conv[ci], val)
    end
    return lp_model
end
_up_con_rhs!(lp_model::JuMP.Model, conv_key::Symbol, val::Number, cidx::Integer) = _up_con_rhs!(lp_model, conv_key, [val], [cidx])

_get_con_rhs(lp_model::JuMP.Model, conv_key::Symbol) = JuMP.normalized_rhs.(lp_model[conv_key])

function _delete!(lp_model, con)
    try
        JuMP.delete(lp_model, con)
        catch err; !(err isa JuMP.MOI.InvalidIndex) && rethrow(err)
    end
    return lp_model
end
_delete!(lp_model, sym::Symbol) = (_delete!(lp_model, lp_model[sym]); JuMP.unregister(lp_model, sym))

_get_vars(lp_model::JuMP.Model) = lp_model[_RXNS_VAR_KEY]
_set_vars!(lp_model::JuMP.Model, N) = (lp_model[_RXNS_VAR_KEY] = @JuMP.variable(lp_model, [1:N]); lp_model)

function build_lp_model(S, b, lb, ub, solver)
    lp_model = JuMP.Model(solver)
    JuMP.set_silent(lp_model)
    _set_vars!(lp_model, size(S, 2))
    set_lb_con!(lp_model, lb)
    set_ub_con!(lp_model, ub)
    set_stoi_con!(lp_model, S, b)
    return lp_model
end
build_lp_model(net::MetNets.MetNet, solver) = build_lp_model(net.S, net.b, net.lb, net.ub, solver)
build_lp_model(net::MetNets.MetNet; solver = Clp.Optimizer) = build_lp_model(net, solver)

function set_start_value(lp_model, vals)
    xs = _get_vars(lp_model)
    for (x, val) in zip(xs, vals)
        JuMP.set_start_value(x, val)
    end
    return lp_model
end

function set_stoi_con!(lp_model, S, b)
    haskey(lp_model, _STOI_CON_KEY) && _delete!(lp_model, _STOI_CON_KEY)
    x = _get_vars(lp_model)
    lp_model[_STOI_CON_KEY] = @JuMP.constraint(lp_model, S * x .- b .== 0.0)
    return lp_model
end

function set_lb_con!(lp_model, cidxs, lb)
    if haskey(lp_model, _LB_CON_KEY)
        _up_con_rhs!(lp_model, _LB_CON_KEY, lb, cidxs)
    else
        x = _get_vars(lp_model)
        lp_model[_LB_CON_KEY] = @JuMP.constraint(lp_model, x .>= lb, base_name = string(_LB_CON_KEY))
    end
    return lp_model
end
set_lb_con!(lp_model, lb) = set_lb_con!(lp_model, 1:_length(lp_model), lb)

function set_ub_con!(lp_model, cidxs, ub)
    if haskey(lp_model, _UB_CON_KEY)
        _up_con_rhs!(lp_model, _UB_CON_KEY, ub, cidxs)
    else
        x = _get_vars(lp_model)
        lp_model[_UB_CON_KEY] = @JuMP.constraint(lp_model, x .<= ub, base_name = string(_UB_CON_KEY))
    end
    return lp_model
end
set_ub_con!(lp_model, lb) = set_ub_con!(lp_model, 1:_length(lp_model), lb)

function optimize!(lp_model, obj_idx; sense = JuMP.MOI.MAX_SENSE)
    x = _get_vars(lp_model)
    @JuMP.objective(lp_model, sense, x[obj_idx])
    JuMP.optimize!(lp_model)
    return lp_model
end

function FBAOut(lp_model::JuMP.Model, obj_idx; drop_LPsol = true)
    N = length(_get_vars(lp_model))
    status = JuMP.termination_status(lp_model)
    LPsol = drop_LPsol ? nothing : lp_model
    if status != JuMP.MOI.OPTIMAL
        return FBAOut(fill(NaN, N), NaN, obj_idx, LPsol)
    end
    x = _get_vars(lp_model)
    v = JuMP.value.(x)
    return FBAOut(v, v[obj_idx], obj_idx, LPsol)
end


text:
## The scope of this test set is to check that the helper functions in core/utils.jl and core/start_values_methods.jl don't break
@testset "utils and start value methods" begin

    data = _PMD.parse_file(joinpath(BASE_DIR, "test/data/extra/networks/case3_unbalanced.dss"); data_model=_PMD.MATHEMATICAL)

    @testset "dimension reduction" begin
        #test 1: three active connections, do nothing
        case3_data = deepcopy(data)
        reduce_single_phase_loadbuses!(case3_data, exclude = [])
        @test length(case3_data["bus"]["3"]["terminals"]) == 3
        #test 2: two active connections, do nothing
        delete!(case3_data["load"], "2")
        reduce_single_phase_loadbuses!(case3_data, exclude = [])
        @test length(case3_data["bus"]["3"]["terminals"]) == 3
        #test 3: one active connection, but exclude, do nothing
        delete!(case3_data["load"], "3")
        reduce_single_phase_loadbuses!(case3_data, exclude = [3])
        @test length(case3_data["bus"]["3"]["terminals"]) == 3
        #test 4: one active connection, reduce!
        delete!(case3_data["load"], "3")
        reduce_single_phase_loadbuses!(case3_data, exclude = [])
        @test case3_data["bus"]["3"]["terminals"] == [2]

        pf_result= _PMD.solve_mc_pf(case3_data, _PMD.ACPUPowerModel, ipopt_solver)
        @test pf_result["termination_status"] == LOCALLY_SOLVED
    end

    msr_path = joinpath(BASE_DIR, "test/data/extra/measurements/case3_meas.csv")
    _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
    pf_result= _PMD.solve_mc_pf(data, _PMD.ACPUPowerModel, ipopt_solver)

    @testset "other utils" begin

        vn = 1.0
        pg = 1.0
        qg = 0.3*pg
        pd = 0.05
        qd = 0.3*pd

        update_all_bounds!(data; v_min=0.9*vn, v_max=1.1*vn, pg_min=-pg, pg_max=pg, qg_min=-qg, qg_max=qg, pd_min=-pd, pd_max=pd, qd_min=-qd, qd_max=qd)

        @test data["bus"]["1"]["vmin"][1] == data["bus"]["1"]["vmin"][2] == 0.9*vn
        @test data["bus"]["1"]["vmax"][1] == data["bus"]["1"]["vmax"][2] == 1.1*vn
        @test data["gen"]["1"]["pmin"][1] == data["gen"]["1"]["pmin"][2] == -pg
        @test data["gen"]["1"]["pmax"][1] == data["gen"]["1"]["pmax"][2] == pg
        @test data["load"]["1"]["pmin"][1] == data["load"]["1"]["pmin"][2] == -pd
        @test data["load"]["1"]["pmax"][1] == data["load"]["1"]["pmax"][2] == pd
        @test data["gen"]["1"]["qmin"][1] == data["gen"]["1"]["qmin"][2] == -qg
        @test data["gen"]["1"]["qmax"][1] == data["gen"]["1"]["qmax"][2] == qg
        @test data["load"]["1"]["qmin"][1] == data["load"]["1"]["qmin"][2] == -qd
        @test data["load"]["1"]["qmax"][1] == data["load"]["1"]["qmax"][2] == qd

        _PMDSE.assign_basic_individual_criteria!(data; chosen_criterion="rwls")

        @test data["meas"]["2"]["crit"] == "rwls"
        @test data["meas"]["10"]["crit"] == "rwls"

        for (m, meas) in data["meas"]
            if meas["var"] ∈ [:pd, :qd]
               _PMDSE.assign_basic_individual_criteria!(data["meas"][m]; chosen_criterion="mle")
           else
              _PMDSE.assign_basic_individual_criteria!(data["meas"][m]; chosen_criterion="rwlav")
           end
        end

        @test data["meas"]["7"]["crit"] == "mle"
        @test data["meas"]["10"]["crit"] == "mle"
        @test data["meas"]["1"]["crit"] == "rwlav"
        @test data["meas"]["2"]["crit"] == "rwlav"
        @test data["meas"]["6"]["crit"] == "rwlav"

        n_meas = length(data["meas"])
        add_measurement!(data, :p, :branch, 1, [0.0012], [0.0003])
        @test length(data["meas"]) == n_meas+1
        @test data["meas"]["$(n_meas+1)"]["var"] == :p
    end

    @testset "start values" begin
        assign_start_to_variables!(data)
        @test isapprox(data["bus"]["4"]["vm_start"][1], 0.9959; atol = 1e-7)
        @test isapprox(data["load"]["1"]["qd_start"][1], 0.006; atol = 1e-7)

        assign_start_to_variables!(data, pf_result)
        @test isapprox(data["bus"]["3"]["vm_start"][2], 0.981757; atol = 1e-7)
        @test isapprox(data["load"]["3"]["pd_start"][1], 0.018; atol = 1e-7)

        _PMDSE.assign_residual_ub!(data)
        @test data["meas"]["7"]["res_max"] == 100

        data["se_settings"] = Dict{String, Any}("rescaler" => 50)
        _PMDSE.assign_residual_ub!(data, chosen_upper_bound=10.0, rescale = true)
        @test data["meas"]["7"]["res_max"] == 500
    end
end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module AbstractEncounterModelImpl

export AbstractEncounterModel

abstract AbstractEncounterModel

end

module AbstractEncounterModelInterfaces

export
       # common
       generateEncounter,
       getInitialState,
       getNextCommand,
       getTrajectory,

       # CorrAEM
       generateEncountersToFile,
       setInitialDistributions,
       validate,
       getInitialSample

# common

function generateEncounter()
end

function getInitialState()
end

function getNextCommand()
end

function getTrajectory()
end

# CorrAEM

function generateEncountersToFile()
end

function setInitialDistributions()
end

function validate()
end

function getInitialSample()
end

end