text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""

    AbstractMenu

The supertype for all Menu types.

# Functions

The following functions can be called on all <:AbstractMenu types.
Details can be found in

## Exported

  - `request(m::AbstractMenu)`
  - `request(msg::AbstractString, m::AbstractMenu)`

## Hidden

  - `printmenu(m::AbstractMenu, cursor::Int; init::Bool=false, oldstate=nothing)`

# Subtypes

All subtypes must be mutable, and must contain the fields `pagesize::Int` and
`pageoffset::Int`. They must also implement the following functions.

## Necessary Functions

These functions must be implemented for all subtypes of AbstractMenu.

  - `pick(m::AbstractMenu, cursor::Int)`
  - `cancel(m::AbstractMenu)`
  - `options(m::AbstractMenu)`   # `numoptions` is an alternative
  - `writeline(buf::IO, m::AbstractMenu, idx::Int, iscursor::Bool)`

If `m` does not have a field called `selected`, then you must also implement `selected(m)`.

## Optional Functions

These functions do not need to be implemented for all AbstractMenu
subtypes.

  - `header(m::AbstractMenu)`
  - `keypress(m::AbstractMenu, i::UInt32)`
  - `numoptions(m::AbstractMenu)`
  - `selected(m::AbstractMenu)`

"""
abstract type AbstractMenu end

function getproperty(m::AbstractMenu, name::Symbol)
    if name === :pagesize
        return getfield(m, :pagesize)::Int
    elseif name === :pageoffset
        return getfield(m, :pageoffset)::Int
    end
    return getfield(m, name)
end

# TODO Julia2.0: get rid of parametric intermediate, making it just
#   abstract type ConfiguredMenu <: AbstractMenu end
# Or perhaps just make all menus ConfiguredMenus
# Also consider making `cursor` a mandatory field in the Menu structs
# instead of going via the RefValue in `request`.
abstract type _ConfiguredMenu{C} <: AbstractMenu end
const ConfiguredMenu = _ConfiguredMenu{<:AbstractConfig}

# NECESSARY FUNCTIONS
# These functions must be implemented for all subtypes of AbstractMenu
######################################################################

"""
    pick(m::AbstractMenu, cursor::Int)

Defines what happens when a user presses the Enter key while the menu is open.
If `true` is returned, `request()` will exit.
`cursor` indexes the position of the selection.
"""
pick(m::AbstractMenu, cursor::Int) = error("unimplemented")

"""
    cancel(m::AbstractMenu)

Define what happens when a user cancels ('q' or ctrl-c) a menu.
`request()` will always exit after calling this function.
"""
cancel(m::AbstractMenu) = error("unimplemented")

"""
    options(m::AbstractMenu)

Return a list of strings to be displayed as options in the current page.

Alternatively, implement `numoptions`, in which case `options` is not needed.
"""
options(m::AbstractMenu) = error("unimplemented")

"""
    writeline(buf::IO, m::AbstractMenu, idx::Int, iscursor::Bool)

Write the option at index `idx` to `buf`. `iscursor`, if `true`, indicates that this
item is at the current cursor position (the one that will be selected by hitting "Enter").

If `m` is a `ConfiguredMenu`, `TerminalMenus` will print the cursor indicator.
Otherwise the callee is expected to handle such printing.

!!! compat "Julia 1.6"
    `writeline` requires Julia 1.6 or higher.

    On older versions of Julia, this was
        `writeLine(buf::IO, m::AbstractMenu, idx, iscursor::Bool)`
    and `m` is assumed to be unconfigured. The selection and cursor indicators can be
    obtained from `TerminalMenus.CONFIG`.

    This older function is supported on all Julia 1.x versions but will be dropped in Julia 2.0.
"""
function writeline(buf::IO, m::AbstractMenu, idx::Int, iscursor::Bool)
    # error("unimplemented")    # TODO: use this in Julia 2.0
    writeLine(buf, m, idx, iscursor)
end

# OPTIONAL FUNCTIONS
# These functions do not need to be implemented for all menu types
##################################################################

"""
    header(m::AbstractMenu) -> String

Returns a header string to be printed above the menu.
Defaults to "".
"""
header(m::AbstractMenu) = ""

"""
    keypress(m::AbstractMenu, i::UInt32) -> Bool

Handle any non-standard keypress event.
If `true` is returned, [`TerminalMenus.request`](@ref) will exit.
Defaults to `false`.
"""
keypress(m::AbstractMenu, i::UInt32) = false

"""
    numoptions(m::AbstractMenu) -> Int

Return the number of options in menu `m`. Defaults to `length(options(m))`.
"""
numoptions(m::AbstractMenu) = length(options(m))

"""
    selected(m::AbstractMenu)

Return information about the user-selected option.
By default it returns `m.selected`.
"""
selected(m::AbstractMenu) = m.selected

"""
    request(m::AbstractMenu; cursor=1)

Display the menu and enter interactive mode. `cursor` indicates the item
number used for the initial cursor position. `cursor` can be either an
`Int` or a `RefValue{Int}`. The latter is useful for observation and
control of the cursor position from the outside.

Returns `selected(m)`.
"""
request(m::AbstractMenu; kwargs...) = request(terminal, m; kwargs...)

function request(term::REPL.Terminals.TTYTerminal, m::AbstractMenu; cursor::Union{Int, Base.RefValue{Int}}=1, suppress_output=false)
    if cursor isa Int
        cursor = Ref(cursor)
    end
    menu_header = header(m)
    !suppress_output && !isempty(menu_header) && println(term.out_stream, menu_header)

    state = nothing
    if !suppress_output
        state = printmenu(term.out_stream, m, cursor[], init=true)
    end

    raw_mode_enabled = try
        REPL.Terminals.raw!(term, true)
        true
    catch err
        suppress_output || @warn("TerminalMenus: Unable to enter raw mode: $err")
        false
    end
    # hide the cursor
    raw_mode_enabled && !suppress_output && print(term.out_stream, "\x1b[?25l")

    try
        while true
            lastoption = numoptions(m)
            c = readkey(term.in_stream)

            if c == Int(ARROW_UP)
                cursor[] = move_up!(m, cursor[], lastoption)
            elseif c == Int(ARROW_DOWN)
                cursor[] = move_down!(m, cursor[], lastoption)
            elseif c == Int(PAGE_UP)
                cursor[] = page_up!(m, cursor[], lastoption)
            elseif c == Int(PAGE_DOWN)
                cursor[] = page_down!(m, cursor[], lastoption)
            elseif c == Int(HOME_KEY)
                cursor[] = 1
                m.pageoffset = 0
            elseif c == Int(END_KEY)
                cursor[] = lastoption
                m.pageoffset = lastoption - m.pagesize
            elseif c == 13 # <enter>
                # will break if pick returns true
                pick(m, cursor[]) && break
            elseif c == UInt32('q')
                cancel(m)
                break
            elseif c == 3 # ctrl-c
                cancel(m)
                ctrl_c_interrupt(m) ? throw(InterruptException()) : break
            else
                # will break if keypress returns true
                keypress(m, c) && break
            end

            if !suppress_output
                state = printmenu(term.out_stream, m, cursor[], oldstate=state)
            end
        end
    finally # always disable raw mode
        if raw_mode_enabled
            !suppress_output && print(term.out_stream, "\x1b[?25h") # unhide cursor
            REPL.Terminals.raw!(term, false)
        end
    end
    !suppress_output && println(term.out_stream)

    return selected(m)
end

"""
    request([term,] msg::AbstractString, m::AbstractMenu)

Shorthand for `println(msg); request(m)`.
"""
request(msg::AbstractString, m::AbstractMenu; kwargs...) = request(terminal, msg, m; kwargs...)

function request(term::REPL.Terminals.TTYTerminal, msg::AbstractString, m::AbstractMenu; kwargs...)
    println(term.out_stream, msg)
    request(term, m; kwargs...)
end

function move_up!(m::AbstractMenu, cursor::Int, lastoption::Int=numoptions(m))
    if cursor > 1
        cursor -= 1 # move selection up
        if cursor < (2+m.pageoffset) && m.pageoffset > 0
            m.pageoffset -= 1 # scroll page up
        end
    elseif scroll_wrap(m)
        # wrap to bottom
        cursor = lastoption
        m.pageoffset = lastoption - m.pagesize
    end
    return cursor
end

function move_down!(m::AbstractMenu, cursor::Int, lastoption::Int=numoptions(m))
    if cursor < lastoption
        cursor += 1 # move selection down
        pagepos = m.pagesize + m.pageoffset
        if pagepos <= cursor && pagepos < lastoption
            m.pageoffset += 1 # scroll page down
        end
    elseif scroll_wrap(m)
        # wrap to top
        cursor = 1
        m.pageoffset = 0
    end
    return cursor
end

function page_up!(m::AbstractMenu, cursor::Int, lastoption::Int=numoptions(m))
    # If we're at the bottom, move the page 1 less to move the cursor up from
    # the bottom entry, since we try to avoid putting the cursor at bounds.
    m.pageoffset -= m.pagesize - (cursor == lastoption ? 1 : 0)
    m.pageoffset = max(m.pageoffset, 0)
    return max(cursor - m.pagesize, 1)
end

function page_down!(m::AbstractMenu, cursor::Int, lastoption::Int=numoptions(m))
    m.pageoffset += m.pagesize - (cursor == 1 ? 1 : 0)
    m.pageoffset = min(m.pageoffset, lastoption - m.pagesize)
    return min(cursor + m.pagesize, lastoption)
end

"""
    printmenu(out, m::AbstractMenu, cursoridx::Int; init::Bool=false, oldstate=nothing) -> newstate

Display the state of a menu. `init=true` causes `m.pageoffset` to be initialized to start printing at
or just above the current cursor location; when `init` is false, the terminal will
preserve the current setting of `m.pageoffset` and overwrite the previous display.
Returns `newstate`, which can be passed in as `oldstate` on the next call to allow accurate
overwriting of the previous display.

!!! compat "Julia 1.6"
    `printmenu` requires Julia 1.6 or higher.

    On older versions of Julia, this was called `printMenu` and it lacked the `state` argument/return value.
    This older function is supported on all Julia 1.x versions but will be dropped in Julia 2.0.
"""
function printmenu(out::IO, m::AbstractMenu, cursoridx::Int; oldstate=nothing, init::Bool=false)
    # TODO Julia 2.0?: get rid of `init` and just use `oldstate`
    buf = IOBuffer()
    lastoption = numoptions(m)::Int
    ncleared = oldstate === nothing ? m.pagesize-1 : oldstate

    if init
        # like clamp, except this takes the min if max < min
        m.pageoffset = max(0, min(cursoridx - m.pagesize ÷ 2, lastoption - m.pagesize))
    else
        print(buf, "\x1b[999D\x1b[$(ncleared)A")   # move left 999 spaces and up `ncleared` lines
    end

    firstline = m.pageoffset+1
    lastline = min(m.pagesize+m.pageoffset, lastoption)

    for i in firstline:lastline
        # clearline
        print(buf, "\x1b[2K")

        upscrollable = i == firstline && m.pageoffset > 0
        downscrollable = i == lastline && i != lastoption

        if upscrollable && downscrollable
            print(buf, updown_arrow(m)::Union{Char,String})
        elseif upscrollable
            print(buf, up_arrow(m)::Union{Char,String})
        elseif downscrollable
            print(buf, down_arrow(m)::Union{Char,String})
        else
            print(buf, ' ')
        end

        printcursor(buf, m, i == cursoridx)
        writeline(buf, m, i, i == cursoridx)

        (firstline == lastline || i != lastline) && print(buf, "\r\n")
    end

    newstate = lastline-firstline  # final line doesn't have `\n`
    if newstate < ncleared && oldstate !== nothing
        # we printed fewer lines than last time. Erase the leftovers.
        for i = newstate+1:ncleared
            print(buf, "\r\n\x1b[2K")
        end
        print(buf, "\x1b[$(ncleared-newstate)A")
    end

    print(out, String(take!(buf)))

    return newstate
end

scroll_wrap(m::ConfiguredMenu) = scroll_wrap(m.config)
scroll_wrap(c::AbstractConfig) = scroll_wrap(c.config)
scroll_wrap(c::Config) = c.scroll_wrap
scroll_wrap(::AbstractMenu) = CONFIG[:scroll_wrap]

ctrl_c_interrupt(m::ConfiguredMenu) = ctrl_c_interrupt(m.config)
ctrl_c_interrupt(c::AbstractConfig) = ctrl_c_interrupt(c.config)
ctrl_c_interrupt(c::Config) = c.ctrl_c_interrupt
ctrl_c_interrupt(::AbstractMenu) = CONFIG[:ctrl_c_interrupt]

up_arrow(m::ConfiguredMenu) = up_arrow(m.config)
up_arrow(c::AbstractConfig) = up_arrow(c.config)
up_arrow(c::Config) = c.up_arrow
up_arrow(::AbstractMenu) = CONFIG[:up_arrow]

down_arrow(m::ConfiguredMenu) = down_arrow(m.config)
down_arrow(c::AbstractConfig) = down_arrow(c.config)
down_arrow(c::Config) = c.down_arrow
down_arrow(::AbstractMenu) = CONFIG[:down_arrow]

updown_arrow(m::ConfiguredMenu) = updown_arrow(m.config)
updown_arrow(c::AbstractConfig) = updown_arrow(c.config)
updown_arrow(c::Config) = c.updown_arrow
updown_arrow(::AbstractMenu) = CONFIG[:updown_arrow]

printcursor(buf, m::ConfiguredMenu, iscursor::Bool) = print(buf, iscursor ? cursor(m.config) : ' ', ' ')
cursor(c::AbstractConfig) = cursor(c.config)
cursor(c::Config) = c.cursor
printcursor(buf, ::AbstractMenu, ::Bool) = nothing   # `writeLine` is expected to do the printing (get from CONFIG[:cursor])

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
"""
Initialize model parameters (in-place).
"""
init!(x::IRTmodel, q, u) = _init!(x, q, u)
init(x::IRTmodel, q, u) = _init(x, q, u)

# graded
function _init!(x::graded, q, u)
    # u is a vector of item response
    freq = map(j -> count(i -> i == j, u), sort(unique(skipmissing(u))))
    prop = freq ./ sum(freq)
    b = @. quantile(Normal(), cumsum(prop))
    d = -b[begin:1:end-1]
    a = fill(1.5, length(q))
    a[.!q] .= 0.0
    x.a = a
    x.d = d
    x.group = [1]
    x.fixed = false
end

function _init(x::graded, q, u)
    # u is a vector of item response
    freq = map(j -> count(i -> i == j, u), sort(unique(skipmissing(u))))
    prop = freq ./ sum(freq)
    b = @. quantile(Normal(), cumsum(prop))
    d = -b[begin:1:end-1]
    a = fill(1.5, length(q))
    a[.!q] .= 0.0
    return graded(a = a, d = d, group = [1], fixed = false)
end

"""
    setgroupindex!(x::IRTmodel, u, g)
Set the group index to the struct::IRTmodel.
"""
function setgroupindex!(x::IRTmodel, u, g)
    group = unique(g[u .!== missing])
    x.group = group
end

# guessing

function _init!(x::guessing, q, u)
    prop = mean(u)
    d = [-quantile(Normal(), prop)]
    a = fill(1.5, count(q))
    x.a = a
    x.d = d
    x.c = [c]
    x.group = [1]
    x.fixed = false
    x.estc = true
end

function init(x::guessing, q, u)
    prop = mean(u)
    d = [-quantile(Normal(), prop)]
    a = fill(1.5, count(q))
    return guessing(a = a, d = d, c = [0.2], group = [1], fixed = false, estc = true)
end

# LPE
function _init!(x::LPE, q, u)
    # u is a vector of item response
    freq = map(j -> count(i -> i == j, u), sort(unique(skipmissing(u))))
    prop = freq ./ sum(freq)
    b = @. quantile(Normal(), cumsum(prop))
    d = -b[begin:1:end-1]
    a = fill(1.5, length(q))
    a[.!q] .= 0.0
    x.a = a
    x.d = d
    x.ξ = [log(1.0)]
    x.group = [1]
    x.fixed = false
end

function _init(x::LPE, q, u)
    # u is a vector of item response
    freq = map(j -> count(i -> i == j, u), sort(unique(skipmissing(u))))
    prop = freq ./ sum(freq)
    b = @. quantile(Normal(), cumsum(prop))
    d = -b[begin:1:end-1]
    a = fill(1.5, length(q))
    a[.!q] .= 0.0
    return LPE(a = a, d = d, ξ = [log(1.0)], group = [1], fixed = false)
end

ext: jl


text:
## Display string vector on fixed with screen
function show_string_vector(sV :: Vector{T1},  width :: T2 = 80;
    io :: IO = stdout)  where {T1 <: AbstractString,  T2 <: Integer}

    n = length(sV);
    if n < 1
        return nothing
    end

    iCol = 0;
    for s in sV
        len1 = length(s);
        if len1 > 0
            if iCol + len1 > width
                println(io, " ")
                iCol = 0;
            end
            print(io, s);
            iCol = iCol + len1;
            print(io, "    ");
            iCol = iCol + 4;
        end
    end
    println(io, " ")
end

function formatted_value(v :: AbstractFloat)
    return string(round(v, digits = 3))
end

function formatted_value(v :: Vector{T1}) where T1 <: AbstractFloat
    vStr = "";
    for j = 1 : length(v)
        vStr = vStr * formatted_value(v[j]);
        if j < length(v)
            vStr = vStr * " | ";
        end
    end
    return vStr
end

function formatted_value(v :: Array{T1}) where T1 <: AbstractFloat
    vStr = formatted_value(v[1]) * " ... " * formatted_value(v[end])
end

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


text:
using FileIO, Test, NomnomlJS, VisualRegressionTests

@testset "NomnomlJS" begin
    dir = joinpath(@__DIR__, "data")
    file = joinpath(dir, "test.noml")
    string = read(file, String)
    d1 = read(file, Diagram)
    d2 = Diagram(string)
    @test d1.src == d2.src
    # Check that the non-SVG results are similar to our reference. Testing the
    # SVG isn't important since they all get generated from that to begin with.
    for ext in ("png",) # TODO "pdf", "eps")
        gen = joinpath(dir, "generated.$ext")
        write(gen, d1)
        func = fname -> save(fname, (load(gen)))
        @visualtest func joinpath(dir, "reference.png") false 5.0
    end
    for ext in ("pdf", "eps")
        @test write(joinpath(dir, "generated.$ext"), d1) > 0
    end
end


text:
mutable struct PlanktonModel
    arch::Architecture          # architecture on which models will run
    t::Int64                    # time in second
    iteration::Int64            # model interation
    individuals::individuals    # initial individuals generated by `setup_agents`
    nutrients::NamedTuple       # initial nutrient fields
    grid::AbstractGrid          # grid information
    bgc_params::Dict            # biogeochemical parameter set
    timestepper::timestepper    # operating Tuples and arrays for timestep
    mode::AbstractMode          # Carbon, Quota, or MacroMolecular
end

"""
    PlanktonModel(arch::Architecture, grid::AbstractGrid;
                  mode = QuotaMode(),
                  N_species = 1,
                  N_individual = 1024,
                  max_individuals = 8*1024,
                  bgc_params = bgc_params_default(),
                  phyt_params = phyt_params_default(N_species, mode),
                  nut_initial = default_nut_init(),
                  t = 0.0,
                  )

Generate a `PlanktonModel` data structure.

Keyword Arguments (Required)
============================
- `arch` : `CPU()` or `GPU()`. Computer architecture being used to run the model.
- `grid` : a `AbstractGrid` structure. Discrete grid for the model (resolution and geometry).

Keyword Arguments (Optional)
============================
- `mode` : Phytoplankton physiology mode, choose among CarbonMode(), QuotaMode(), or MacroMolecularMode().
- `N_species` : Number of species.
- `N_individual` : Number of individuals per species.
- `max_individuals` : Maximum number of individuals per species the model can hold.
- `bgc_params` : Parameter set for biogeochemical processes modeled in the model.
- `phyt_params` : Parameter set for physiological processes of individuals modeled in the model.
- `nut_initial` : The source of initial conditions of nutrient fields, should be either a `NamedTuple`
                           or a `Dict` containing the file paths pointing to the files of nutrient initial conditions.
- `t` : Model time, start from 0 by default, in second.
"""
function PlanktonModel(arch::Architecture, grid::AbstractGrid;
                       mode = QuotaMode(),
                       N_species::Int64 = 1,
                       N_individual::Int64 = 1024,
                       max_individuals::Int64 = 8*1024,
                       bgc_params = bgc_params_default(),
                       phyt_params = phyt_params_default(N_species, mode),
                       nut_initial = default_nut_init(),
                       t::Int64 = 0,
                       )

    @assert N_individual ≤ max_individuals

    if arch == GPU() && !has_cuda()
        throw(ArgumentError("Cannot create a GPU model. No CUDA-enabled GPU was detected!"))
    end

    grid_d = replace_grid_storage(arch, grid)

    inds = generate_individuals(phyt_params, arch, N_species, N_individual, max_individuals, grid_d, mode)

    nutrients = generate_nutrients(arch, grid_d, nut_initial)

    ts = timestepper(arch, grid_d, max_individuals)

    iteration  = 0

    model = PlanktonModel(arch, t, iteration, inds, nutrients, grid_d, bgc_params, ts, mode)

    return model
end

function show(io::IO, model::PlanktonModel)
    Nsp = length(model.individuals.phytos)
    N = Int(dot(model.individuals.phytos.sp1.data.ac,model.individuals.phytos.sp1.data.ac))
    cap = length(model.individuals.phytos.sp1.data.ac)
    print(io, "PlanktonModel:\n",
              "├── grid: $(short_show(model.grid))\n",
              "├── $(model.mode) is selected for phytoplankton physiology\n",
              "├── individuals: $(Nsp) phytoplankton species each with $(N) individuals\n",
              "└── maximum number of individuals: $(cap) per species")
end


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

name = "MPFR"
version = v"4.1.0"

# Collection of sources required to build MPFR
sources = [
    ArchiveSource("https://www.mpfr.org/mpfr-$(version)/mpfr-$(version).tar.xz",
                  "0c98a3f1732ff6ca4ea690552079da9c597872d30e96ec28414ee23c95558a7f"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/mpfr-*
./configure --prefix=$prefix --host=$target --enable-shared --disable-static --with-gmp=${prefix} --enable-thread-safe --enable-shared-cache --disable-float128 --disable-decimal-float
make -j${nproc}
make install

# On Windows, make sure non-versioned filename exists...
if [[ ${target} == *mingw* ]]; then
    cp -v ${prefix}/bin/libmpfr-*.dll ${prefix}/bin/libmpfr.dll
fi
"""

# We enable experimental platforms as this is a core Julia dependency
platforms = supported_platforms(;experimental=true)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    # We explicitly ask for GMP v6.1.2 so that it is compatible with both GMP v6.1.2 and v6.2.0+
    Dependency("GMP_jll", v"6.1.2"),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"5")


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

import .Zygote

Zygote.@nograd warn_kmm_julialib
Zygote.@nograd safe_diagm


text:
# 2D magnetic flux function.

"""
    compute_flux_function(b, Δ, nG=2)

Calculate the 2D magnetic flux function ψ from the magnetic field `b` and discrete steps
`Δ`. `nG` is the number of ghost cells along each dimension in the vector field.
ψ is defined as
``\\psi = \\int B_x dz = - \\int B_z dx`` from ``\\mathbf{B} = \\hat{y}\\times\\nabla\\psi``
in the X-Z plane and Y is the out-of-plane direction. This is strictly true if the guide
field By is zero and B is divergence-free. However, numerically there will be errors.
The current implementation calculates ψ by integrating along -z boundary first, and then
going along z.
Reference: [Flux function](https://doi.org/10.1063/1.3657424)
"""
function compute_flux_function(b::AbstractArray{T,N}, Δ, nG=2) where {T,N}
   dx, dz = Δ
   flux = zeros(T, size(b,2), size(b,3))

   # Find ψ along the -z boundary cells
   @inbounds for i in axes(b,2)[1+nG:end-nG]
      bz = b[3,i,1+nG]
      flux[i,1+nG] = flux[i-1,1+nG] - bz * dx
   end

   # For each row, integrate in z-direction
   @inbounds for j in axes(b,3)[2+nG:end-nG], i in axes(b,2)[1+nG:end-nG]
      bx = b[1,i,j]
      flux[i,j] = flux[i,j-1] + bx * dz
   end

   flux
end

"""
    find_reconnection_points(ψ; retol=1e-1, method=1) -> indices_x, indices_o

Find X-point and O-point indices in 2D magnetic field topology from flux function `ψ`.
The current implementation does not work for the 3 layers near the boundary.

# Optional arguments
- `retol=1e-1`: determines the relative tolerance of the ratio w.r.t. |∇ψ|² to accept a
gradient as 0.
- `method=1`: method 1 compute the cell-centered 1st and 2nd order derivatives and check the
Hessian matrix; method 2 check the flux function at each point against its 8 neighbors,
which is more deterministic.
"""
function find_reconnection_points(ψ; retol=1e-4, method=1)

   indices_x = Matrix{Int64}(undef, 2, 0)
   indices_o = Matrix{Int64}(undef, 2, 0)

   if method == 1
      ∂ψ = gradient(ψ)
      ∂²ψ = gradient(view(∂ψ,1,:,:))
      # fill boundary layers
      ∂²ψ[:,2,:] = ∂²ψ[:,1,:]
      ∂²ψ[:,end-1,:] = ∂²ψ[:,end,:]
      ∂²ψ[:,:,2] = ∂²ψ[:,:,1]
      ∂²ψ[:,:,end-1] = ∂²ψ[:,:,end]
      # 2nd derivatives
      ∂²ψ∂x² = ∂²ψ[1,:,:]
      ∂²ψ∂x∂y = ∂²ψ[2,:,:] # == ∂²ψ∂y∂x = ∂²ψ[1,:,:]

      ∂²ψ = gradient(view(∂ψ,2,:,:))
      # fill boundary layers
      ∂²ψ[:,2,:] = ∂²ψ[:,1,:]
      ∂²ψ[:,end-1,:] = ∂²ψ[:,end,:]
      ∂²ψ[:,:,2] = ∂²ψ[:,:,1]
      ∂²ψ[:,:,end-1] = ∂²ψ[:,:,end]

      ∂²ψ∂y² = ∂²ψ[2,:,:]

      ∂ψmag² = [∂ψ[1,i,j]^2 + ∂ψ[2,i,j]^2 for j in axes(∂ψ,3), i in axes(∂ψ,2)]
      ∂ψmean = mean(∂ψmag²)

      for j in axes(∂ψ,3)[4:end-3], i in axes(∂ψ,2)[4:end-3]
         ∂ψ[1,i,j]^2 + ∂ψ[2,i,j]^2 > retol*∂ψmean && continue
         # Hessian matrix det(H) < 0 => saddle point (X-point)
         if ∂²ψ∂x²[i,j] * ∂²ψ∂y²[i,j] - ∂²ψ∂x∂y[i,j]^2 < 0
            indices_x = hcat(indices_x, [i,j])
         else # O-point
            indices_o = hcat(indices_o, [i,j])
         end
      end
   elseif method == 2
      for j in axes(ψ,2)[4:end-3], i in axes(ψ,1)[4:end-3]
         if issaddle(@views ψ[i-1:i+1,j-1:j+1])
            indices_x = hcat(indices_x, [i,j])
         elseif isextrema(@views ψ[i-1:i+1,j-1:j+1])
            indices_o = hcat(indices_o, [i,j])
         end
      end
   end

   return indices_x, indices_o
end

"Check if the center point in a 3x3 matrix `ψ` is a saddle point."
function issaddle(ψ)
   nmax, nmin = 0, 0
   minmax1 = extrema(@views ψ[:,2])
   minmax2 = extrema(@views ψ[2,:])
   minmax3 = extrema([ψ[1,1],ψ[2,2],ψ[3,3]])
   minmax4 = extrema([ψ[1,3],ψ[2,2],ψ[3,1]])
   if ψ[2,2] == minmax1[1]
      nmin += 1
   elseif ψ[2,2] == minmax1[2]
      nmax += 1
   end
   if ψ[2,2] == minmax2[1]
      nmin += 1
   elseif ψ[2,2] == minmax2[2]
      nmax += 1
   end
   if ψ[2,2] == minmax3[1]
      nmin += 1
   elseif ψ[2,2] == minmax3[2]
      nmax += 1
   end
   if ψ[2,2] == minmax4[1]
      nmin += 1
   elseif ψ[2,2] == minmax4[2]
      nmax += 1
   end
   if nmin > 0 && nmax > 0
      return true
   else
      return false
   end
end

"Check if the center point in a 3x3 matrix `ψ` is an extrema point."
function isextrema(ψ)
   if ψ[2,2] == maximum(ψ) || ψ[2,2] == minimum(ψ)
      return true
   else
      return false
   end
end


text:
module Quadtrees

export position,
    Boundary,
    Rect,
    Circle,
    Quadtree,
    subdivide!,
    isleaf,
    query,
    query!,
    search,
    search!

import Base: push!, insert!, append!, length, in, +, -, *, /

const Point = Tuple{Real, Real}

"""
        distance²((x₁, y₁), (x₂, y₂)) -> Real

Returns the square distance between two points.

```math
d²(p₁, p₂) = (x₂ - x₁)² + (y₂ - y₁)²
```
"""
distance²((x₁, y₁)::Point, (x₂, y₂)::Point) = (x₂ - x₁) ^ 2 + (y₂ - y₁) ^ 2
"""
        distance((x₁, y₁), (x₂, y₂)) -> Real

Returns the distance between two points.

```math
distance(p₁, p₂) = \\sqrt{(x₂ - x₁)² + (y₂ - y₁)²}
```
"""
distance(P₁::Point, P₂::Point) = √distance²(P₁, P₂)

"""
        position(value::T) -> Tuple{Real, Real}

This function is used for the [in](@ref) method for determine
the intersections between the [Boundary](@ref) objects and points.

This function is not mented to be used directly, but to be extended
for the objects `T` that your [quadtree](@ref) carries.

This function **must** return an tuple of two [Real](@ref) numbers.

For example,

```julia
mutble struct MyType
    pos::Vector{Real}
    # other fields...
end

Quadtrees.position(obj::MyType) = (obj.pos[1], obj.pos[2])
```
"""
function position(::T) where {T}
    throw(ErrorException("You need to implement the Quadtrees.position method for the type $T"))
end

position((x, y)::Point) = (x, y)

function _position(value::T) where {T}
    pos = position(value)
    @assert typeof(pos) <: Point "The return value of the position method must be
    a tuple of two Real numbers, but got type $(typeof(pos)) instead."
    return pos;
end

abstract type Boundary end
"""
Represents a 2D rectangle boundary.
"""
struct Rect <: Boundary
    x::Real
    y::Real
    w::Real
    h::Real
end

+(r::Rect, (x, y)::Point) = Rect(r.x + x, r.y + y, r.w, r.h)
-(r::Rect, (x, y)::Point) = Rect(r.x - x, r.y - y, r.w, r.h)
*(r::Rect, n::Real) = Rect(r.x, r.y, n * r.w, n * r.h)
/(r::Rect, n::Real) = Rect(r.x, r.y, r.w / n, r.h / n)

in(item::T, R::Rect) where {T} = begin
    x, y = _position(item)
    rx, ry = R.w / 2, R.h / 2

    x ≥ R.x - rx &&
    x ≤ R.x + rx &&
    y ≥ R.y - ry &&
    y ≤ R.y + ry
end

function in(A::Rect, B::Rect)
    Arx, Ary = A.w / 2, A.h / 2
    Brx, Bry = B.w / 2, B.h / 2
    return !(
        A.x - Arx > B.x + Brx ||
        A.x + Arx < B.x - Brx ||
        A.y - Ary > B.y + Bry ||
        A.y + Ary < B.y - Bry
    )
end

"""
        vertices(rect) -> Vector

Returns the vertices of a rectangle.
"""
function vertices(R::Rect)
    x, y = R.x, R.y
    rx, ry = R.w / 2, R.h / 2
    return Point[
        # Top-Left corner
        (x - rx, y - ry),
        # Top-Right corner
        (x + rx, y - ry),
        # Bottom-Right corner
        (x + rx, y + ry),
        # Botton-Left corner
        (x - rx, y + ry)
    ]
end

"""
Represents a 2D circular boundary.
"""
struct Circle <: Boundary
    x::Real
    y::Real
    r::Real
end

+(c::Circle, (x, y)::Point) = Circle(c.x + x, c.y + y, c.r)
-(c::Circle, (x, y)::Point) = Circle(c.x - x, c.y - y, c.r)
*(c::Circle, n::Real) = Circle(c.x, c.y, n * c.r)
/(c::Circle, n::Real) = Circle(c.x, c.y, c.r / n)

Base.convert(::Type{Rect}, C::Circle) = Rect(C.x, C.y, C.r, C.r)

in(item::T, C::Circle) where {T} = begin
    x, y = _position(item)
    d² = distance²((x, y), (C.x, C.y))
    return d² ≤ C.r ^ 2
end

in(A::Circle, B::Circle) = begin
    d² = distance²((B.x, B.y), (A.x, A.y))
    return d² < A.r + B.r
end

in(self::Circle, range::Rect) = begin
    dx = abs(range.x - self.x)
    dy = abs(range.y - self.y)

    r = self.r
    w, h = range.w/2, range.h/2
    edges = (dx - w)^2 + (dy - h)^2

    if dx > r + w || dy > r + h return false end
    if dx <= w || dy <= h return true end;

    return edges <= r^2
end

in(R::Rect, C::Circle) = C in R

mutable struct Quadtree{T}
    datatype::Type{T}
    capacity::Int
    bounds::Rect
    items::Vector{T}
    size::Int
    sections::Vector{Quadtree{T}}
    Quadtree{T}(init::Function) where {T} = (self = new{T}(); init(self); self)
end

function Quadtree{T}(bounds::Rect, capacity::Int) where {T}
    return Quadtree{T}() do self
        self.datatype = T
        self.capacity = capacity
        self.bounds = bounds
        self.size = 0
        self.items = Vector{T}()
        self.sections = Vector{Quadtree{T}}()
    end
end

function Quadtree{T}((x, y)::Point, width::Real, height::Real, capacity::Int) where {T}
    return Quadtree{T}(Rect(x, y, width, height), capacity)
end

function Quadtree{T}(init::Function, range::AbstractRange, bounds::Rect, capacity::Int) where {T}
    qt = Quadtree{T}(bounds, capacity)

    foreach(range) do i
        push!(qt, init(i))
    end

    return qt
end

"""
        isleaf(qt::Quadtree) -> Bool

Return true if `qt` is leaf (has no children quadtree).
"""
isleaf(qt::Quadtree) = isempty(qt.sections)

"""
        insert!(qt::Quadtree{T}, item::T) -> Bool

Insert an `item` in the quadtree `qt`.

Returns true if the item was inserted successfully, false otherwise.
"""
function insert!(qt::Quadtree{T}, item::T)::Bool where {T}
    # do not add if its not
    # in the current bounds
    !(item in qt.bounds) && return false;

    if length(qt.items) < qt.capacity
        push!(qt.items, item)
        qt.size += 1
        return true;
    else
        isleaf(qt) && subdivide!(qt)
        for section in qt.sections
            if insert!(section, item)
                qt.size += 1
                return true;
            end
        end
    end
end

"""
        subdivide!(qt::Quadtree)

Subdivides the quadtree `qt` in four sections.
"""
function subdivide!(qt::Quadtree{T}) where {T}
    w, h = qt.bounds.w, qt.bounds.h
    P = vertices(qt.bounds / 2)
    qt.sections = Quadtree{T}.(P, w/2, h/2, qt.capacity)
end

"""
        push!(qt, items...)

Inserts all `items` on `qt`.
"""
push!(qt::Quadtree{T}, items::T...) where {T} = begin
    for item in items
        insert!(qt, item)
    end
end

append!(qt::Quadtree{T}, collection::Vector{T}...) where {T} = begin
    for items in collection
        push!(qt, items...)
    end
end

"""
        length(qt::Quadtree) -> Int

Returns the number of items in `qt`.
"""
length(qt::Quadtree) = qt.size
"""
        size(qt::Quadtree) -> Int

Returns the number of items in `qt`.
"""
size(qt::Quadtree) = qt.size

"""
        query(qt::Quadtree{T}, bounds::Boundary) -> Vector{T}

Returns all items from `qt` that are in the `bounds` boundary or an empty array,
otherwise.
"""
function query(qt::Quadtree{T}, bounds::Boundary)::Vector{T} where {T}
    result = Vector{T}()

    if bounds in qt.bounds
        append!(result, filter(item -> item in bounds, qt.items))
        if !isleaf(qt)
            foreach(qt.sections) do section
                append!(result, query(section, bounds))
            end
        end
    end

    return result;
end
"""
        query!(array, quadtree, bounds)

In place version of the [query](@ref) method.
"""
query!(array::Vector{T}, qt::Quadtree{T}, range::Boundary) where {T} = empty!(array) && append!(array, query(qt, range))

"""
        search(quadtree, x, y, radius) -> Vector

Returns all items from the `quadtree` that falls on a circular boundary,
where `(x, y)` is the position of the search and `r` is the radius.
"""
function search(qt::Quadtree, x::Real, y::Real, radius::Real)
    return query(qt, Circle(x, y, radius))
end

"""
        search!(array, quadtree, x, y, radius)

In place version of the [search](@ref) method.
"""
function search!(array::Vector{T}, qt::Quadtree{T}, x::Real, y::Real, radius::Real) where {T}
    query!(array, qt, Circle(x, y, radius))
end

end # module


text:
"""
    AbstractDimStack

Abstract supertype for dimensional stacks.

These have multiple layers of data, but share dimensions.
"""
abstract type AbstractDimStack{L} end

data(s::AbstractDimStack) = s.data
dims(s::AbstractDimStack) = s.dims
refdims(s::AbstractDimStack) = s.refdims
layerdims(s::AbstractDimStack) = s.layerdims
layerdims(s::AbstractDimStack, key::Symbol) = dims(s, layerdims(s)[key])
metadata(s::AbstractDimStack) = s.metadata
layermetadata(s::AbstractDimStack) = s.layermetadata
layermetadata(s::AbstractDimStack, key::Symbol) = layermetadata(s)[key]
metadata(s::AbstractDimStack, key::Symbol) = layermetadata(s)[key]

layerdims(A::AbstractDimArray) = basedims(A)

@inline Base.keys(s::AbstractDimStack) = keys(data(s))
Base.haskey(s::AbstractDimStack, k) = k in keys(s)
Base.values(s::AbstractDimStack) = values(layers(s))
Base.first(s::AbstractDimStack) = s[first(keys(s))]
Base.last(s::AbstractDimStack) = s[last(keys(s))]
# Only compare data and dim - metadata and refdims can be different
Base.:(==)(s1::AbstractDimStack, s2::AbstractDimStack) =
    data(s1) == data(s2) && dims(s1) == dims(s2) && layerdims(s1) == layerdims(s2)
Base.length(s::AbstractDimStack) = length(keys(s))
Base.size(s::AbstractDimStack) = map(length, dims(s))
Base.size(A::AbstractDimStack, dims::DimOrDimType) = size(A, dimnum(A, dims))
Base.size(A::AbstractDimStack, dims::Integer) = size(A)[dims]
Base.axes(s::AbstractDimStack) = map(first ∘ axes, dims(s))
Base.axes(A::AbstractDimStack, dims::DimOrDimType) = axes(A, dimnum(A, dims))
Base.axes(A::AbstractDimStack, dims::Integer) = axes(A)[dims]
Base.iterate(s::AbstractDimStack, args...) = iterate(layers(s), args...)

function rebuild(
    s::AbstractDimStack, data, dims=dims(s), refdims=refdims(s),
    layerdims=layerdims(s), metadata=metadata(s), layermetadata=layermetadata(s)
)
    basetypeof(s)(data, dims, refdims, layerdims, metadata, layermetadata)
end

function rebuildsliced(f::Function, s::AbstractDimStack, layers, I)
    layerdims = map(basedims, layers)
    dims, refdims = slicedims(f, s, I)
    rebuild(s; data=map(parent, layers), dims=dims, refdims=refdims, layerdims=layerdims)
end

"""
    rebuild_from_arrays(s::AbstractDimStack, das::NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractDimArray}}}; kw...)

Rebuild an `AbstractDimStack` from a `NamedTuple` of `AbstractDimArray`
and an existing stack.

# Keywords

Keywords are simply the fields of the stack object:

- `data`
- `dims`
- `refdims`
- `metadata`
- `layerdims`
- `layermetadata`
"""
function rebuild_from_arrays(
    s::AbstractDimStack, das::NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractDimArray}}};
    refdims=refdims(s),
    metadata=DD.metadata(s),
    data=map(parent, das),
    dims=DD.combinedims(das...),
    layerdims=map(DD.basedims, das),
    layermetadata=map(DD.metadata, das),
)
    rebuild(s; data, dims, refdims, layerdims, metadata, layermetadata)
end

function layers(s::AbstractDimStack{<:NamedTuple{Keys}}) where Keys
    NamedTuple{Keys}(map(K -> s[K], Keys))
end

Adapt.adapt_structure(to, s::AbstractDimStack) = map(A -> Adapt.adapt(to, A), s)

# Dipatch on Tuple of Dimension, and map
for func in (:index, :mode, :metadata, :sampling, :span, :bounds, :locus, :order)
    @eval ($func)(s::AbstractDimStack, args...) = ($func)(dims(s), args...)
end

"""
    Base.copy!(dst::AbstractArray, src::DiskGeoStack, key::Key)

Copy the stack layer `key` to `dst`, which can be any `AbstractArray`.

## Example

Copy the `:humidity` layer from `stack` to `array`.

```julia
copy!(array, stack, :humidity)
```
"""
Base.copy!(dst::AbstractArray, src::AbstractDimStack, key) = copy!(dst, src[key])

"""
    Base.copy!(dst::AbstractDimStack, src::AbstractDimStack, [keys=keys(dst)])

Copy all or a subset of layers from one stack to another.

## Example

Copy just the `:sea_surface_temp` and `:humidity` layers from `src` to `dst`.

```julia
copy!(dst::AbstractDimStack, src::AbstractDimStack, keys=(:sea_surface_temp, :humidity))
```
"""
function Base.copy!(dst::AbstractDimStack, src::AbstractDimStack, keys=keys(dst))
    # Check all keys first so we don't copy anything if there is any error
    for key in keys
        key in Base.keys(dst) || throw(ArgumentError("key $key not found in dest keys"))
        key in Base.keys(src) || throw(ArgumentError("key $key not found in source keys"))
    end
    for key in keys
        copy!(dst[key], src[key])
    end
end

"""
    Base.map(f, s::AbstractDimStack)

Apply functrion `f` to each layer of the stack `s`, and rebuild it.

If `f` returns `DimArray`s the result will be another `DimStack`.
Other values will be returned in a `NamedTuple`.
"""
Base.map(f, s::AbstractDimStack...) = _maybestack(s[1], map(f, map(layers, s)...))

_maybestack(s::AbstractDimStack, x::NamedTuple) = x
function _maybestack(
    s::AbstractDimStack, das::NamedTuple{K,<:Tuple{Vararg{<:AbstractDimArray}}}
) where K
    rebuild_from_arrays(s, das; refdims=DD.refdims(s), metadata=DD.metadata(s))
end

"""
    Base.cat(stacks::AbstractGeoStack...; [keys=keys(stacks[1])], dims)

Concatenate all or a subset of layers for all passed in stacks.

# Keywords

- `keys`: `Tuple` of `Symbol` for the stack keys to concatenate.
- `dims`: Dimension of child array to concatenate on.

# Example

Concatenate the :sea_surface_temp and :humidity layers in the time dimension:

```julia
cat(stacks...; keys=(:sea_surface_temp, :humidity), dims=Ti)
```
"""
function Base.cat(s1::AbstractDimStack, stacks::AbstractDimStack...; keys=keys(s1), dims)
    vals = Tuple(cat((s[k] for s in (s1, stacks...))...; dims=dims) for k in keys)
    DD.basetypeof(s1)(NamedTuple{keys}(vals))
end

# Array methods

# Methods with no arguments that return a DimStack
for (mod, fnames) in
    (:Base => (:inv, :adjoint, :transpose), :LinearAlgebra => (:Transpose,))
    for fname in fnames
        @eval ($mod.$fname)(s::AbstractDimStack) = map(A -> ($mod.$fname)(A), s)
    end
end

# Methods with an argument that return a DimStack
for fname in (:rotl90, :rotr90, :rot180, :PermutedDimsArray, :permutedims)
    @eval (Base.$fname)(s::AbstractDimStack, args...) =
        map(A -> (Base.$fname)(A, args...), s)
end

# Methods with keyword arguments that return a DimStack
for (mod, fnames) in
    (:Base => (:sum, :prod, :maximum, :minimum, :extrema, :dropdims),
     :Statistics => (:mean, :median, :std, :var))
    for fname in fnames
        @eval function ($mod.$fname)(s::AbstractDimStack; dims=:, kw...)
            map(s) do A
                # Ignore dims not found in layer
                if dims isa Union{Colon,Int}
                    ($mod.$fname)(A; dims, kw...)
                else
                    ldims = commondims(DD.dims(A), dims)
                    # With no matching dims we do nothing
                    ldims == () ? A : ($mod.$fname)(A; dims=ldims, kw...)
                end
            end
        end
    end
end
for fname in (:cor, :cov)
    @eval function (Statistics.$fname)(s::AbstractDimStack; dims=1, kw...)
        map(s) do A
            if dims isa Int
                (Statistics.$fname)(A; dims, kw...)
            else
                ldims = only(commondims(DD.dims(A), dims))
                ldims == () ? A : (Statistics.$fname)(A; dims=ldims, kw...)
            end
        end
    end
end

# Methods that take a function
for (mod, fnames) in (:Base => (:reduce, :sum, :prod, :maximum, :minimum, :extrema),
                      :Statistics => (:mean,))
    for fname in fnames
        _fname = Symbol(:_, fname)
        @eval function ($mod.$fname)(f::Function, s::AbstractDimStack; dims=Colon())
            map(A -> ($mod.$fname)(f, A; dims=dims), s)
        end
                # ($_fname)(f, s, dims)
            # Colon returns a NamedTuple
            # ($_fname)(f::Function, s::AbstractDimStack, dims::Colon) =
                # map(A -> ($mod.$fname)(f, A), data(s))
            # Otherwise maybe return a DimStack
            # function ($_fname)(f::Function, s::AbstractDimStack, dims) =
            # end
        # end
    end
end

"""
    DimStack <: AbstractDimStack

    DimStack(data::AbstractDimArray...)
    DimStack(data::Tuple{Vararg{<:AbstractDimArray}})
    DimStack(data::NamedTuple{Keys,Vararg{<:AbstractDimArray}})
    DimStack(data::NamedTuple, dims::DimTuple; metadata=NoMetadata())

DimStack holds multiple objects sharing some dimensions, in a `NamedTuple`.
Indexing operates as for [`AbstractDimArray`](@ref), except it occurs for all
data layers of the stack simulataneously. Layer objects can hold values of any type.

DimStack can be constructed from multiple `AbstractDimArray` or a `NamedTuple`
of `AbstractArray` and a matching `dims` `Tuple`. If `AbstractDimArray`s have
the same name they will be given the name `:layer1`, substitiuting the
layer number for `1`.

`getindex` with `Int` or `Dimension`s or `Selector`s that resolve to `Int` will
return a `NamedTuple` of values from each layer in the stack. This has very good
performace, and usually takes less time than the sum of indexing each array
separately.

Indexing with a `Vector` or `Colon` will return another `DimStack` where
all data layers have been sliced.  `setindex!` must pass a `Tuple` or `NamedTuple` maching
the layers.

Most `Base` and `Statistics` methods that apply to `AbstractArray` can be used on
all layers of the stack simulataneously. The result is a `DimStack`, or
a `NamedTuple` if methods like `mean` are used without `dims` arguments, and
return a single non-array value.

## Example

```jldoctest
julia> using DimensionalData

julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0];

julia> dimz = (X([:a, :b]), Y(10.0:10.0:30.0))

  X: Symbol[a, b],
  Y: 10.0:10.0:30.0

julia> da1 = DimArray(1A, dimz; name=:one);

julia> da2 = DimArray(2A, dimz; name=:two);

julia> da3 = DimArray(3A, dimz; name=:three);

julia> s = DimStack(da1, da2, da3);

julia> s[At(:b), At(10.0)]
(one = 4.0, two = 8.0, three = 12.0)

julia> s[X(At(:a))] isa DimStack
true
```

"""
struct DimStack{L,D<:Tuple,R<:Tuple,LD<:NamedTuple,M,LM<:NamedTuple} <: AbstractDimStack{L}
    data::L
    dims::D
    refdims::R
    layerdims::LD
    metadata::M
    layermetadata::LM
end
DimStack(das::AbstractDimArray...; kwargs...) = DimStack(das; kwargs...)
function DimStack(das::Tuple{Vararg{<:AbstractDimArray}}; kwargs...)
    DimStack(NamedTuple{uniquekeys(das)}(das); kwargs...)
end
function DimStack(das::NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractDimArray}}};
    data=map(parent, das), dims=combinedims(das...), layerdims=map(basedims, das),
    refdims=(), metadata=NoMetadata(), layermetadata=map(DD.metadata, das)
)
    DimStack(data, dims, refdims, layerdims, metadata, layermetadata)
end
# Same sized arrays
function DimStack(data::NamedTuple, dims::Tuple;
    refdims=(), metadata=NoMetadata(), layermetadata=map(_ -> NoMetadata(), data)
)
    all(map(d -> axes(d) == axes(first(data)), data)) || _stack_size_mismatch()
    layerdims = map(_ -> basedims(dims), data)
    DimStack(data, formatdims(first(data), dims), refdims, layerdims, metadata, layermetadata)
end

@noinline _stack_size_mismatch() = throw(ArgumentError("Arrays must have identical axes. For mixed dimensions, use DimArrays`"))


text:
using OffsetArrays
using .SetProb

# ============================================================================
function rp1(parms::CParam,
             maxm::Int, meqn::Int, mwaves::Int, maux::Int, mbc::Int, mx::Int,
             ql::OffsetArray{Float64,2}, qr::OffsetArray{Float64,2},
             auxl::OffsetArray{Float64}, auxr::OffsetArray{Float64},
             wave::OffsetArray{Float64,3}, s::OffsetArray{Float64,2},
             amdq::OffsetArray{Float64,2}, apdq::OffsetArray{Float64,2})
# ============================================================================

# subroutine rp1(maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)

# Riemann solver for the acoustics equations in 1d,

# On input, ql contains the state vector at the left edge of each cell
#           qr contains the state vector at the right edge of each cell

# On output, wave contains the waves,
#            s the speeds,
#
#            amdq = A^- Delta q,
#            apdq = A^+ Delta q,
#                   the decomposition of the flux difference
#                       f(qr(i-1)) - f(ql(i))
#                   into leftgoing and rightgoing parts respectively.
#

# Note that the i'th Riemann problem has left state qr(i-1,:)
#                                    and right state ql(i,:)
# From the basic clawpack routines, this routine is called with ql = qr

#    dimension wave(meqn, mwaves, 1-mbc:maxm+mbc)
#    dimension    s(mwaves,1-mbc:maxm+mbc)
#    dimension   ql(meqn, 1-mbc:maxm+mbc)
#    dimension   qr(meqn, 1-mbc:maxm+mbc)
#    dimension apdq(meqn, 1-mbc:maxm+mbc)
#    dimension amdq(meqn, 1-mbc:maxm+mbc)

#     local arrays
#     ------------
#    dimension delta(2)
delta = Vector{Float64}(undef, 2)

# density, bulk modulus, and sound speed, and impedence of medium:
# (should be set in setprob.f)
#    common /cparam/ rho,bulk,cc,zz

    # split the jump in q at each interface into waves

    # find a1 and a2, the coefficients of the 2 eigenvectors:

    for i = 2-mbc : mx+mbc
        delta[1] = ql[1,i] - qr[1,i-1]
        delta[2] = ql[2,i] - qr[2,i-1]
        a1 = (-delta[1] + parms.zz*delta[2]) / (2.0*parms.zz)
        a2 = ( delta[1] + parms.zz*delta[2]) / (2.0*parms.zz)

        # Compute the waves.

        wave[1,1,i] = -a1*parms.zz
        wave[2,1,i] = a1
        s[1,i] = -parms.cc

        wave[1,2,i] = a2*parms.zz
        wave[2,2,i] = a2
        s[2,i] = parms.cc

    end

    # compute the leftgoing and rightgoing flux differences:
    # Note s(1,i) < 0   and   s(2,i) > 0.

    for m = 1 : meqn
        for i = 2-mbc : mx+mbc
            amdq[m,i] = s[1,i]*wave[m,1,i]
            apdq[m,i] = s[2,i]*wave[m,2,i]
        end
    end

end #subroutine rp1


text:
module BenchCumsum

using BenchmarkTools
using ThreadedScans

const CACHE = Ref{Any}()

function setup(n = 2^23)
    CACHE[] = rand(-2^5:2^5, n)
    T = typeof(CACHE[])

    suite = BenchmarkGroup()
    suite["base"] = @benchmarkable(cumsum!(ys, xs), setup = begin
        xs = CACHE[]::$T
        ys = similar(xs)
    end)
    suite["simple!"] = @benchmarkable(ThreadedScans.simple!(+, xs), setup = begin
        xs = copy(CACHE[]::$T)
    end)
    suite["dac!"] = @benchmarkable(ThreadedScans.dac!(+, xs), setup = begin
        xs = copy(CACHE[]::$T)
    end)
    suite["linear!"] = @benchmarkable(ThreadedScans.linear!(+, xs), setup = begin
        xs = copy(CACHE[]::$T)
    end)
    suite["partitioned_hillis_steele!"] = @benchmarkable(
        ThreadedScans.partitioned_hillis_steele!(+, xs),
        setup = begin
            xs = copy(CACHE[]::$T)
        end
    )

    return suite
end

function clear()
    CACHE[] = nothing
end

end  # module


text:
export Ray

struct Ray{T}
    orig::Point3{T}
    dir::Vec3{T}
end
at(r::Ray, t) = r.orig + t * r.dir


text:
"""
For a prime number p and a multiple q, the wheel index encodes the prime number index of
p and q modulo 30, which can be encoded to a single number from 1 ... 64. This allows us
to jump immediately into the correct loop at the correct offset.
"""
create_jump(wheel_index, i) = :($wheel_index === $i && @goto $(Symbol(:x, i)))

create_label(wheel_index) = :(@label $(Symbol(:x, wheel_index)))

wheel_mask(prime_mod_30)::UInt8 = ~(0x01 << (to_idx(prime_mod_30) - 1))

"""
For any prime number `p` we compute its prime number index modulo 30 (here `wheel`) and we
generate the loop that crosses of the next 8 multiples that, modulo 30, are
p * {1, 7, 11, 13, 17, 19, 23, 29}.
"""
function unrolled_loop(p_idx)
    p = ps[p_idx]

    # First push the stopping criterion
    unrolled_loop_body = Any[:(byte_idx > unrolled_max && break)]

    # Cross off the 8 next multiples
    for q in ps
        div, rem = divrem(p * q, 30)
        bit = wheel_mask(rem)
        push!(unrolled_loop_body, :(xs[byte_idx + increment * $(q - 1) + $div] &= $bit))
    end

    # Increment the byte index to where the next / 9th multiple is located
    push!(unrolled_loop_body, :(byte_idx += increment * 30 + $p))

    quote
        while true
            $(unrolled_loop_body...)
        end
    end
end

"""
The fan-in / fan-out phase that crosses off one multiple and then checks bounds; this is
before and after the unrolled loop starts and finishes respectively.
"""
function single_loop_item_not_unrolled(p_idx, q_idx, save_on_exit = true)
    # Our prime number modulo 30
    p = ps[p_idx]

    ps_next = (1, 7, 11, 13, 17, 19, 23, 29, 31)

    # Label name
    jump_idx = 8 * (p_idx - 1) + q_idx

    # Current and next multiplier modulo 30
    q_curr, q_next = ps_next[q_idx], ps_next[q_idx + 1]

    # Get the bit mask for crossing off p * q_curr
    div_curr, rem_curr = divrem(p * q_curr, 30)
    bit = wheel_mask(rem_curr)

    # Compute the increments for the byte index for the next multiple
    incr_bytes = p * q_next ÷ 30 - div_curr
    incr_multiple = q_next - q_curr

    quote
        # Todo: this generates an extra jump, maybe conditional moves are possible?
        if byte_idx > n_bytes

            # For a segmented sieve we store where we exit the loop, since that is the
            # entrypoint in the next loop; to avoid modulo computation to find the offset
            $(save_on_exit ? :(last_idx = $jump_idx) : nothing)
            @goto out
        end

        # Cross off the multiple
        xs[byte_idx] &= $bit

        # Increment the byte index to where the next multiple is located
        byte_idx += increment * $incr_multiple + $incr_bytes
    end
end

"""
Full loop generates a potentially unrolled loop for a particular wheel
that may or may not save the exit point.
"""
function full_loop_for_wheel(wheel, unroll = true, save_on_exit = true)
    loop_statements = []

    for i = 1 : 8
        push!(loop_statements, create_label(8 * (wheel - 1) + i))
        unroll && i == 1 && push!(loop_statements, unrolled_loop(wheel))
        push!(loop_statements, single_loop_item_not_unrolled(wheel, i, save_on_exit))
    end

    quote
        while true
            $(loop_statements...)
        end
    end
end

"""
Generates a sieving loop that crosses off multiples of a given prime number.

    @sieve_loop :unroll :save_on_exit
    @sieve_loop
"""
macro sieve_loop(options...)
    unroll, save_on_exit = :(:unroll) in options, :(:save_on_exit) in options

    # When crossing off p * q where `p` is the siever prime and `q` the current multiplier
    # we have that p and q are {1, 7, 11, 13, 17, 19, 23, 29} mod 30.
    # For each of these 8 possibilities for `p` we create a loop, and per loop we
    # create 8 entrypoints to jump into. The first entrypoint is the unrolled loop for
    # whenever we can remove 8 multiples at the same time when all 8 fit in the interval
    # between byte_start:byte_next_start-1. Otherwise we can only remove one multiple at
    # a time. With 8 loops and 8 entrypoints per loop we have 64 different labels, numbered
    # x1 ... x64.

    # As an example, take p = 7 as a prime number and q = 23 as the first multiplier, and
    # assume our number line starts at 1 (so byte 1 represents 1:30, byte 2 represent 31:60).
    # We have to cross off 7 * 23 = 161 first, which has byte index 6. Our prime number `p`
    # is in the 2nd spoke of the wheel and q is in the 7th spoke. This means we have to jump
    # to the 7th label in the 2nd loop; that is label 8 * (2 - 1) + 7 = 15. There we cross
    # off the multiple (since 161 % 30 = 11 is the 3rd spoke, we "and" the byte with 0b11011111)
    # Then we move to 7 * 29 (increment the byte index accordingly), cross it off as well.
    # And now we enter the unrolled loop where 7 * {31, 37, ..., 59} are crossed off, then
    # 7 * {61, 67, ..., 89} etc. Lastly we reach the end of the sieving interval, we cross
    # off the remaining multiples one by one, until the byte index is passed the end.
    # When that is the case, we save at which multiple / label we exited, so we can jump
    # there without computation when the next interval of the number line is sieved.

    esc(quote
        $(unroll ? :(unrolled_max = n_bytes - increment * 28 - 28) : nothing)

        # Create jumps inside loops
        $([create_jump(:wheel_idx, i) for i = 1 : 64]...)

        # # Create loops
        $([full_loop_for_wheel(wheel, unroll, save_on_exit) for wheel in 1 : 8]...)

        # Point of exit
        @label out
    end)
end


text:
module DebbyPackerTests
using DebbyPacker
using FactCheck: facts, context, @fact, greater_than, not, @fact_throws, anyof, anything

import Base: ismatch
ismatch(regex::String) = x -> ismatch(Regex(regex), x)
ismatch(regex::Regex) = x -> ismatch(regex, x)
contains(item) = x -> item ∈ x

#= include("GitSource.jl") =#
include("debianize.jl")

end

max_stars_repo_name: UCL/DebbyPacker

max_issues_repo_name: UCL/DebbyPacker

max_forks_repo_name: UCL/DebbyPacker


text:
using LinearAlgebra
using Base: require_one_based_indexing

"""
    systemSolver(A::AbstractMatrix, b::Vector, x_initial::Vector) -> Vector

    Compute the orthogonal projection of a given point to the
solution set of a system of linear equations. Also a direct
method for solving systems of linear equations. The output
is either the projection or inconsistency of the system.

# Arguments
 - `A::AbstractMatrix{Real}`: coefficients of linear equations
 - `b::AbstractVector{Real}`: intended results of equations
 - `x0::AbstractVector{Real}`: initial point for projection
 - `β::Real`: set in order to obtain comparable summands in norm

# References
> Ján Plesník, Finding the orthogonal projection of a point onto
> an affine subspace, Linear Algebra and its Applications, Volume 422,
> Issues 2–3, 2007, Pages 455-470
"""
function systemSolver(A::AbstractMatrix{<:Real},    # in: equations matrix (m x n) rows = coefficients
                    b::AbstractVector{<:Real},      # in: equations results
                    x0::AbstractVector{<:Real},     # in: initial point
                    β::Real=1.0)                    # in: optional normalizer

    require_one_based_indexing(A, b, x0)
    m, n = size(A)
    if m != length(b) || n != length(x0)
        throw(DimensionMismatch("dimensions of `A`, $m * $n, does not match one of the lengths of `x0`, $(length(x0)), or `b`, $(length(b))"))
    end
    for i in 1:m
        iszero(A[i, :]) && throw(ArgumentError("`A` rows must be nonzero"))
    end
    iszero(β) && throw(ArgumentError("`β` must be nonzero"))

    if m == 0
        return x0
    end

    k = 1
    p = 1
    V = hcat(A[1, :])
    x = x0 + (b[1] - dot(A[1, :], x0))/(dot(A[1, :], V[:, 1])) * V[:, 1]

    denominators = zeros(Real, m-1)

    while k != m
        k = k + 1
        p = p + 1

        denominators[k-1] = dot((A[k, :] - A[1, :]), V[:, k-1])
        if denominators[k-1] == 0
            A[k, :] = - A[k, :]
            b[k] = - b[k]
            denominators[k-1] = dot((A[k, :] - A[1, :]), V[:, k-1])
        end

        w = x0 + β * A[k, :]
        for j in 2:k
            w = w + (b[j]-b[1]-dot(A[j, :]-A[1, :], w))/denominators[j-1] * V[:, j-1]
        end

        u = x + (b[k] - dot(A[k, :], x))/denominators[k-1] * V[:, k-1]
        if u != w
            V = hcat(V, w - u)
            x = u + (b[1] - dot(A[1, :], u))/dot(A[1, :], V[:, k]) * V[:, k]
        elseif dot(A[k, :], x) == b[k]
            if k != m
                A = A[setdiff(1:end, k), :]
                b = b[setdiff(1:end, k)]
                k = k - 1
                m = m - 1
            end
        else
            throw(DomainError("the ($p)th equation is inconsistent"))
        end
    end
    return x
end

max_stars_repo_name: amir-ni/ALA-extra

max_issues_repo_name: amir-ni/ALA-extra

max_forks_repo_name: amir-ni/ALA-extra


text:
module SubPkgB

println("EVALUATING SubPkgB!")

fnB() = println("SubPkgB.fnB()")

function __init__()
	println("Initializing SubPkgB!")
end

end # module


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

JLLWrappers.@generate_main_file_header("DecFP")
JLLWrappers.@generate_main_file("DecFP", UUID("47200ebd-12ce-5be5-abb7-8e082af23329"))
end  # module DecFP_jll


text:
@testset "multivariate" begin
    struct CholeskyMvNormal{TL,Tm} <: Distributions.AbstractMvNormal
        m::Tm
        L::TL
    end
    Distributions.mean(p::CholeskyMvNormal) = p.m
    Distributions.cov(p::CholeskyMvNormal) = p.L * p.L'
    Distributions.rand(p::CholeskyMvNormal, n::Int) = p.m .+ p.L * randn(length(p), n)
    Distributions.length(p::CholeskyMvNormal) = length(p.m)
    function Distributions.logpdf(p::CholeskyMvNormal, x::AbstractVector)
        return -0.5 * (length(p) * log(2π) + 2 * logdet(p.L) + sum(abs2, p.L \ (x .- p.m)))
    end
    @testset "AbstractMvNormal" begin
        n_dim = 2
        X1 = cholesky(Matrix(0.5 * I(n_dim))).L
        X2 = cholesky(Matrix(0.3 * I(n_dim))).L
        p = CholeskyMvNormal(zeros(n_dim), X1)
        q = CholeskyMvNormal(ones(n_dim), X2)
        @test KL(p, q) > 0
        @test KL(p, q) ≈ KL(p, q, 100_000) atol = 0.1
    end
    @testset "MvNormal" begin
        n_dim = 2
        p = MvNormal(zeros(n_dim), Matrix(0.5 * I(n_dim)))
        q = MvNormal(ones(n_dim), Matrix(0.3 * I(n_dim)))
        @test KL(p, q) > 0
        @test KL(p, q) ≈ KL(p, q, 100_000) atol = 0.1
    end
end


text:
struct Run
    lap::Vector{Integer}
    date::Vector{DateTime}
    lon::Vector{AbstractFloat}
    lat::Vector{AbstractFloat}
    ele::Vector{AbstractFloat}
    hr::Vector{Integer}
    function Run(date, lon, lat,
        ele=Vector{AbstractFloat}[], lap=Vector{Integer}[], hr=Vector{Integer}[])
        n = length(lon)
        if length(lat) != n
            error("'lon' and 'lat' differ in length")
        end
        if length(date) != n
            error("length of 'data' is not equal to the length of 'lon' and 'lat'")
        end
        if length(lap) == 0
            lap = one(n)
        elseif length(lap) != n
            error("length of 'lap' is not equal to the length of 'lon' and 'lat'")
        end
        if length(ele) != n & length(ele) != 0
            error("length of 'ele' is not zero or equal to the length of 'lon' and 'lat'")
        end
        if length(hr) != n & length(hr) != 0
            error("length of 'hr' is not zero or equal to the length of 'lon' and 'lat'")
        end
        new(lap, date, lon, lat, ele, hr)
    end
end

# Methods
Base.length(run::Run) = length(run.date)

## Indexing
Base.endof(run::Run) = length(run)
function Base.getindex(run::Run, i::Number)
    i = convert(Int, i)
    1 <= i <= length(run) || throw(BoundsError(run, i))
    lap = [run.lap[i]]
    date = [run.date[i]]
    lon = [run.lon[i]]
    lat = [run.lat[i]]
    if length(run.ele) == 0
        ele = run.ele
    else
        ele = [run.ele[i]]
    end
    if length(run.hr) == 0
        hr = run.hr
    else
        hr = [run.hr[i]]
    end
    Run(lap, date, lon, lat, ele, hr)
end

function Base.getindex(run::Run, i)
    1 <= minimum(i) <= length(run) || throw(BoundsError(run, i))
    1 <= maximum(i) <= length(run) || throw(BoundsError(run, i))
    lap = [run.lap[j] for j in i]
    date = [run.date[j] for j in i]
    lon = [run.lon[j] for j in i]
    lat = [run.lat[j] for j in i]
    if length(run.ele) == 0
        ele = run.ele
    else
        ele = [run.ele[j] for j in i]
    end
    if length(run.hr) == 0
        hr = run.hr
    else
        hr = [run.hr[j] for j in i]
    end
    Run(lap, date, lon, lat, ele, hr)
end

## Iteration
Base.start(run::Run) = 1
Base.next(run::Run, state) = (run[1], state + 1)
Base.done(run::Run, state) = state > length(run)


text:
using DensestCommonSubgraph

for i=1:10
  d_true = 0.0;
  G = get_dcs_test_instance(15);
  (s_true, d_true) = dcs_parallel(G, 2);
  p_dcs_lp, t_dcs, X_dcs, Y_dcs = dcs_lp(G);
  p_charikar_lp, t_c, X_c, Y_c = charikar_lp(dcs_intersection_graph(G));
  p_avg_lp, t_avg, X_avg, Y_avg = charikar_lp(dcs_avg_graph(G));
  println("opt: $d_true DCS: $(p_dcs_lp.optval) CHARIKAR: $(p_charikar_lp.optval) AVG: $(p_avg_lp.optval)")
  println("opt: $d_true DCS: $(p_dcs_lp.optval) CHARIKAR: $(p_charikar_lp.optval)")
end


text:
@testset "unit_value" begin
    @testset "EnumParameter" begin
        p = EnumParameter([StringParameter("value", "i0"),
                           FloatParameter(1.3, 56.3, 3.3, "i1"),
                           IntegerParameter(1, 40, 3, "i3")],
                           "test")
        c = p.value
        u = unit_value(p)
        v = unit_value!(p, u)
        @test unit_value!(p, u) == c
        @test 1 <= v <= length(p.values)
    end
    @testset "IntegerParameter" begin
        p = IntegerParameter(1, 40, 23, "i0")
        c = p.value
        u = unit_value(p)
        v = unit_value!(p, u)
        @test unit_value!(p, u) == c
        @test p.min <= v <= p.max
    end
    @testset "FloatParameter" begin
        p = FloatParameter(1.0, 40.0, 23.0, "i0")
        c = p.value
        u = unit_value(p)
        v = unit_value!(p, u)
        # Inexact conversion.
        @test p.min <= v <= p.max
        @test unit_value!(p, u) == c
    end
end


text:
export MtlArgumentEncoder

const MTLArgumentEncoder = Ptr{MtArgumentEncoder}

mutable struct MtlArgumentEncoder
    handle::MTLArgumentEncoder
end

function MtlArgumentEncoder(fun::MtlFunction, entry::Integer)
    handle = mtNewArgumentEncoderWithBufferIndexFromFunction(fun, entry-1)
    obj = MtlArgumentEncoder(handle)
    finalizer(unsafe_destroy!, obj)
    return obj
end

function unsafe_destroy!(enc::MtlArgumentEncoder)
    mtRelease(enc.handle)
end

Base.unsafe_convert(::Type{MTLArgumentEncoder}, enc::MtlArgumentEncoder) = enc.handle

Base.:(==)(a::MtlArgumentEncoder, b::MtlArgumentEncoder) = a.handle == b.handle
Base.hash(fun::MtlArgumentEncoder, h::UInt) = hash(mod.handle, h)

## properties

Base.propertynames(o::MtlArgumentEncoder) = (
    # identification
    #=:device, :label,=#
    # creation
    :encodedLength,
    # alignment
    :alignment,
)

function Base.getproperty(o::MtlArgumentEncoder, f::Symbol)
    if f === :encodedLength
        mtArgumentEncoderLength(o)
    elseif f === :alignment
        mtArgumentEncoderAlignment(o)
    else
        getfield(o, f)
    end
end

Base.sizeof(a::MtlArgumentEncoder) = Int(mtArgumentEncoderLength(a))

## operations

function assign_argument_buffer!(enc::MtlArgumentEncoder, buf::MtlBuffer, offset::Integer=1)
    mtArgumentEncoderSetArgumentBufferWithOffset(enc, buf, offset-1)
end

function assign_argument_buffer!(enc::MtlArgumentEncoder, buf::MtlBuffer, offset::Integer, element::Integer)
    mtArgumentEncoderSetArgumentBufferWithOffsetForElement(enc, buf, offset-1, element)
end

set_buffer!(enc::MtlArgumentEncoder, buf::MtlBuffer, offset::Integer, index::Integer) =
    mtArgumentEncoderSetBufferOffsetAtIndex(enc, buf, offset, index-1)
set_buffers!(enc::MtlArgumentEncoder, bufs::Vector{<:MtlBuffer},
             offsets::Vector{Int}, indices::UnitRange{Int}) =
    mtArgumentSetBuffersOffsetsWithRange(enc, handle_array(bufs), offsets, indices .- 1)

function set_bytes!(enc::MtlArgumentEncoder, src, length::Integer, index::Integer)
    dst = Base.bitcast(typeof(ptr), mtArgumentEncoderConstantDataAtIndex(enc, index-1))
    Base.unsafe_copyto!(dst, src, length)
    return
end

function set_field!(enc::MtlArgumentEncoder, val::Number, index::Integer)
    dst = Base.bitcast(Ptr{typeof(val)}, mtArgumentEncoderConstantDataAtIndex(enc, index-1))
    Base.unsafe_store!(dst, val, 1)
    return
end

function set_field!(enc::MtlArgumentEncoder, val::NTuple{N,T}, index::Integer) where {N,T}
    dst = Base.bitcast(Ptr{typeof(val)}, mtArgumentEncoderConstantDataAtIndex(enc, index-1))
    Base.unsafe_store!(dst, val, 1)
    return
end


text:
using DataFrames
# Test evolve_function.jl in Evolve.jl for the specific purpose of reducing the number of gates
#    in a 2-component goal where circuits has already been evolved to map to each single-component goal.

# Goal is the 2-component goal
function test_evolve_function( goal::Goal, p::Parameters, maxsteps::Int64 )
  @assert p.numouputs==2
  p1 = Parameters( div(p.numinputs,2), 1, p.numlevelsback, p.numinteriors, p.numlevelsback )
  maxtries = 10
  tries = 1
  done = :false
  while !done && tries < maxtries
    ec = random_chromosome(p1)
    (c1,steps,maxsteps) = neutral_evolution( ec, [goal[1]], maxsteps )
    if steps < maxsteps
      done = :true
    end
  end
  if !done
    error("unable to evolve circuit to map to goal ",[goal[1]])
  end
  tries = 1
  done = :false
  while !done && tries < maxtries
    ec = random_chromosome(p1)
    (c2,steps,maxsteps) = neutral_evolution( ec, [goal[2]], maxsteps )
    if steps < maxsteps
      done = :true
    end
  end
  if !done
    error("unable to evolve circuit to map to goal ",[goal[1]])
  end
end

# Combines two 1-output circuits c1 and c2 with the same paramterers and number of gates into a 2-output circuit
# The output values of the result is the 2-output goal whose components are the output values of c1 and of c2.
function combine_circuits( c1::Chromosome, c2::Chromosome )
  # Two functions used only within this function
  inp1(i,numinputs) = i <= numinputs ? i : 2*(i-numinputs)+numinputs-1
  inp2(i,numinputs) = i <= numinputs ? i : 2*(i-numinputs)+numinputs
  @assert c1.params.numinputs == c2.params.numinputs
  @assert c1.params.numoutputs == 1
  @assert c2.params.numoutputs == 1
  if length(c1.interiors) < length(c2.interiors)
    println("  c1ints: ",length(c1.interiors),"  c2ints: ",length(c2.interiors))
    ctemp = c1
    c1 =c2
    c2 = ctemp
  end
  #println("c1ints: ",length(c1.interiors),"  c2ints: ",length(c2.interiors))
  p = Parameters( c1.params.numinputs, 2, c1.params.numinteriors+c2.params.numinteriors,
    c1.params.numlevelsback+c2.params.numlevelsback )
  inputs = c1.inputs
  interiors = InteriorNode[]
  new_numinteriors = 2*length(c2.interiors)
  #println("new_numints: ",new_numinteriors)
  #print_parameters(p)
  for i = 1:length(c2.interiors)
    int_node1 = deepcopy(c1.interiors[i])
    input1 = inp1(c1.interiors[i].inputs[1],p.numinputs)
    input2 = inp1(c1.interiors[i].inputs[2],p.numinputs)
    int_node1.inputs = [input1,input2]
    #println("i: ",i,"  int_node1: ",int_node1)
    push!(interiors,int_node1)
    int_node2 = deepcopy(c2.interiors[i])
    input1 = inp2(c2.interiors[i].inputs[1],p.numinputs)
    input2 = inp2(c2.interiors[i].inputs[2],p.numinputs)
    int_node2.inputs = [input1,input2]
    #println("i: ",i,"  int_node2: ",int_node2)
    push!(interiors,int_node2)
  end
  #println("interiors: ",interiors)
  outputs = [OutputNode(p.numinputs+new_numinteriors-1),OutputNode(p.numinputs+new_numinteriors)]
  Chromosome( p, c1.inputs, interiors, outputs, 0.0, 0.0 )
end

function test_combine_circuits( p::Parameters )
  c1 = random_chromosome(p); print(output_values(c1),"  "); print_circuit(c1)
  c2 = random_chromosome(p); print(output_values(c2),"  "); print_circuit(c2)
  c = combine_circuits(c1,c2); print(output_values(c),"  "); print_circuit(c)
end

# Evaluate the performance of function evolve_reduce_numactive() which is in evolve_funct.jl.
# For each pair of chromosomes
#   Creates two chromosomes using parameter settings p.
#   Combines them into one chromosome with one more input
#   For each trial run evolve_reduce_numactive() and record average and minimum numactive.
function test_reduce_numactive( p::Parameters, numpairs::Int64, numtrials::Int64, maxsteps::Int64;
      num_mutations::Int64=1, csvfile::String="" )
  df = DataFrame()
  df.combined_goal = Goal[]
  df.numinputs = Int64[]
  df.numgates = Int64[]
  df.numlevelsback = Int64[]
  df.numpairs = Int64[]
  df.numtrials = Int64[]
  df.maxsteps = Int64[]
  df.nactive0 = Float64[]
  df.nactive1_mean = Float64[]
  df.nactive1_std = Float64[]
  if num_mutations >= 2
    df.nactive2_mean = Float64[]
    df.nactive2_std = Float64[]
  end
  if num_mutations == 3
    df.nactive3_mean = Float64[]
    df.nactive3_std = Float64[]
  end
  for npair = 1:numpairs
    c1 = random_chromosome(p); print(output_values(c1),"  "); print_circuit(c1)
    c2 = random_chromosome(p); print(output_values(c2),"  "); print_circuit(c2)
    c = combine_circuits(c1,c2)
    rowlist = [ output_values(c), p.numinputs, p.numinteriors, p.numlevelsback, numpairs, numtrials, maxsteps ]
    #@assert [output_values(c1)[1],output_values(c1)[1]] == output_values(c)
    println("output_values(c): ",output_values(c))
    println("goal: ",[output_values(c1)[1],output_values(c2)[1]])
    #nactive_lists = map(i->reduce_nactive_helper(c,maxsteps,num_mutations=num_mutations), collect(1:numtrials ) )
    nactive_lists = pmap(i->reduce_nactive_helper(c,maxsteps,num_mutations=num_mutations), collect(1:numtrials ) )
    nactive0 = vcat( [ na[1] for na in nactive_lists]... )
    nactive1 = vcat( [ na[2] for na in nactive_lists]... )
    if num_mutations == 1
      row_results = [mean(nactive0),mean(nactive1),std(nactive1)]
    end
    if num_mutations >= 2
      nactive2 = vcat( [ na[3] for na in nactive_lists]... )
      row_results = [mean(nactive0),mean(nactive1),std(nactive1),mean(nactive2),std(nactive2)]
    end
    if num_mutations >= 3
      nactive3 = vcat( [ na[4] for na in nactive_lists]... )
      row_results = [mean(nactive0),mean(nactive1),std(nactive1),mean(nactive2),std(nactive2),mean(nactive3),std(nactive3)]
    end
    row = deepcopy(rowlist)
    println("length(row): ",length(row),"  length(row_results): ",length(row_results),"  size(df): ",size(df))
    #println("npair: ",npair,"  mean(nactive1): ",mean(nactive1),"  mean(nactive2): ",mean(nactive2),"  mean(nactive3): ",mean(nactive3))
    push!(df,vcat(row,row_results))
  end
  hostname = chomp(open("/etc/hostname") do f read(f,String) end)
  if length(csvfile) > 0
    open( csvfile, "w" ) do f
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      print_parameters(f,p,comment=true)
      println(f,"# numpairs: ", numpairs)
      println(f,"# numtrials: ", numtrials)
      println(f,"# maxsteps: ", maxsteps)
      CSV.write( f, df, append=true, writeheader=true )
    end
  end
  df
end

# nactive0 is number_active_gates() after evolution to find circuit that maps to goal
#    but before running evolve_reduce_numactive()
function reduce_nactive_helper( c::Chromosome, maxsteps::Int64; num_mutations::Int64=1 )
  #println("num_mutations: ",num_mutations)
  goal = output_values( c )
  (ec,steps) = neutral_evolution( random_chromosome(c.params), goal, maxsteps )
  nactive0 = number_active_gates(ec)
  println("nactive0: ",nactive0)
  nactive1 = Int64[]
  nc = deepcopy(c)
  (ec,cur_nactive) = evolve_reduce_numactive( nc, maxsteps, num_mutations=1 )
  @assert cur_nactive == number_active_gates(ec)
  push!(nactive1,number_active_gates(ec))
  if num_mutations == 1
    return(nactive0, nactive1)
  end
  nactive2 = Int64[]
  nc = deepcopy(c)
  (ec,cur_nactive) = evolve_reduce_numactive( nc, maxsteps, num_mutations=2 )
  push!(nactive2,number_active_gates(ec))
  if num_mutations == 2
    return(nactive0, nactive1, nactive2)
  end
  nactive3 = Int64[]
  nc = deepcopy(c)
  (ec,cur_nactive) = evolve_reduce_numactive( nc, maxsteps, num_mutations=3 )
  push!(nactive3,number_active_gates(ec))
  #println("trial: ",trial,"  mean(nactive1): ",mean(nactive1),"  mean(nactive2): ",mean(nactive2),"  mean(nactive2): ",mean(nactive2))
  (nactive0, nactive1, nactive2, nactive3)
end

# folder_path is relative to data folder
function compare_dataframe( folder_path::String, suffix::String )
  rdf = read_dataframe("$(folder_path)/reduce_numactive_$(folder_path)_$(suffix).csv")
  kdf = read_dataframe("$(folder_path)/k_complexity$(folder_path)$(suffix).csv")
  cdf = DataFrame()
  cdf.goal = kdf.goal
  cdf.kng = kdf.num_gates
  cdf.rng = rdf.nactive1_mean
  cdf.k_minus_r_diff = cdf.kng .- cdf.rng
  cdf
end


text:
using MatrixNetworks
#include("../src/diffusions.jl")

using Base.Test

function _normout(P)
    n = size(P,1)
    colsums = sum(P,1)
    (pi,pj,pv) = findnz(P)
    P = sparse(pi,pj,pv./colsums[pj],n,n)
end

function prlinsys(A::SparseMatrixCSC{Float64,Int64},alpha::Float64,v::Vector{Float64})
    P = full(_normout(A'))
    n = size(A,1)
    z =  (eye(n) - alpha*P) \ v
    z = z /sum(z)
    return z
end

function pagerank_test()

    A = load_matrix_network("celegans")
    x = pagerank(A, 0.85)
    y = pagerank(A, 0.85, 1e.-2)

    n = 10
    A = speye(10)
    v = rand(n)
    vn = v/sum(v)
    x = personalized_pagerank(A,0.85,v)
    @test norm(x-vn,1) <= n*eps(Float64)

    v = 5;
    x = personalized_pagerank(A,0.85,v)
    xtrue = zeros(n)
    xtrue[5] = 1.
    @test norm(x -xtrue,1) <= n*eps(Float64)

    x = seeded_pagerank(A,0.85,Set([5]))

    x = personalized_pagerank(A,0.85,Set([5]))
    @test norm(x -xtrue,1) <= n*eps(Float64)

    x = personalized_pagerank(A,0.85,Dict{Int64,Float64}(5 => 6.))
    @test norm(x -xtrue,1) <= n*eps(Float64)

    x = personalized_pagerank(A,0.85,Set(collect(1:8)))
    @test norm(x - sparsevec(collect(1:8),1./8,n),1) <= n*eps(Float64)

    x = personalized_pagerank(A,0.85,Dict{Int64,Float64}(zip(1:8, ones(8))))
    @test norm(x - sparsevec(collect(1:8),1./8,n),1) <= n*eps(Float64)

    @test_throws ArgumentError personalized_pagerank(A,0.85,sparsevec(Dict{Int64,Float64}(zip(1:8, ones(8)))))

    x = personalized_pagerank(A,0.85,sparsevec(Dict{Int64,Float64}(zip(4:10, ones(7)))))
    @test norm(x - sparsevec(collect(4:10),1./7,n),1) <= n*eps(Float64)

    x = personalized_pagerank(A, 0.85, sparsevec([5], [2.], 10))
    @test norm(x -xtrue,1) <= n*eps(Float64)

    A = speye(Int64,n)
    x = personalized_pagerank(A,0.85,5)
    @test norm(x -xtrue,1) <= n*eps(Float64)

    n = 25
    A = spones(sprand(n,n,2/n))
    x = pagerank(A,0.85)
    xtrue = prlinsys(A,0.85,ones(n))
    @test norm(x -xtrue,1) <= n*eps(Float64)

    x = pagerank(MatrixNetwork(A),0.85)
    @test norm(x -xtrue,1) <= n*eps(Float64)

    n = 10
    A = spdiagm(ones(n-1),-1,n,n)'
    v = 1./n
    tol = 1e-8
    x = pagerank(A,0.85)
    z = zeros(n)
    z[1] = 1.
    z[2] = 1.85
    for i=3:n
        z[i] = z[i-1] + 0.85*(z[i-1]-z[i-2])
    end
    z = z/sum_kbn(z)
    @test norm(x-z,1) <= n*tol

    return true
end

function pagerank_error_test()

    @test_throws DimensionMismatch pagerank(spzeros(6,5),0.85)

    @test_throws ArgumentError pagerank(spzeros(5,5),-0.85)
    @test_throws ArgumentError pagerank(spzeros(5,5),1.0)
    @test_throws ArgumentError pagerank(spzeros(5,5),1.01)

    @test_throws ArgumentError personalized_pagerank(spzeros(5,5), 0.85, rand(11))
    @test_throws DomainError personalized_pagerank(spzeros(5,5), 0.85, -ones(5))
    @test_throws DomainError personalized_pagerank(spzeros(5,5), 0.85, -sparsevec(ones(5)))
end

function pagerank_alg_test()
n = 10
P = speye(n)
x = zeros(n)
y = zeros(n)
v = rand(n)
v = v/sum(v)
tol = 1e-8
maxiter = 1000
iterfunc = MatrixNetworks._noiterfunc

x = pagerank_power!(x,y,P,0.85,v,tol,maxiter,iterfunc)
@test norm(x-v,1) <= n*tol

P = spdiagm(ones(n-1),-1,n,n)
v[:] = 0.
v[1] = 1.
x = pagerank_power!(x,y,P,0.85,v,tol,maxiter,iterfunc)
z = zeros(n)
z[1] = 1.
for i=2:n
    z[i] = 0.85*z[i-1]
end
z = z/sum_kbn(z)
@test norm(x-z,1) <= n*tol

tol = 1e-12
x = pagerank_power!(x,y,P,0.85,v,tol,maxiter,iterfunc)
@test norm(x-z,1) <= n*tol

tol = 1e-12
x = pagerank_power!(x,y,P,0.85,sparsevec(v),tol,maxiter,iterfunc)
@test norm(x-z,1) <= n*tol

tol = 1e-15
x = pagerank_power!(x,y,P,0.85,v,tol,maxiter,iterfunc)
@test norm(x-z,1) <= n*tol

tol = 1e-3
x = pagerank_power!(x,y,P,0.85,v,tol,maxiter,(iter,x) -> @show iter, norm(x,1))
@test norm(x-z,1) <= n*tol

v = 1./n
tol = 1e-8
x = pagerank_power!(x,y,P,0.85,v,tol,maxiter,iterfunc)
z[1] = 1.
z[2] = 1.85
for i=3:n
    z[i] = z[i-1] + 0.85*(z[i-1]-z[i-2])
end
z = z/sum_kbn(z)
@test norm(x-z,1) <= n*tol

return true
end

function pagerank_perf_test()
    maxiter=500

    n = 1_000
    x = zeros(n)
    y = zeros(n)
    z = zeros(n)
    tol = 1e-8
    v = 1./n
    A = spones(sprand(n,n,25/n))
    P = _normout(A')

    x = pagerank_power!(x,y,P,0.85,v,tol,maxiter,MatrixNetworks._noiterfunc)
    P2 = MatrixNetworks._create_stochastic_mult(A)
    y = pagerank_power!(y,z,P2,0.85,v,tol,maxiter,MatrixNetworks._noiterfunc)
    P3 = MatrixNetworks._create_stochastic_mult(MatrixNetwork(A))
    y = pagerank_power!(y,z,P3,0.85,v,tol,maxiter,MatrixNetworks._noiterfunc)

    n = 10_000
    x = zeros(n)
    y = zeros(n)
    z = zeros(n)
    v = 1./n
    tol = 1e-8
    A = spones(sprand(n,n,25/n))
    P = _normout(A')

    dt = @elapsed x = pagerank_power!(x,y,P,0.85,v,tol,maxiter,MatrixNetworks._noiterfunc)
    P2 = MatrixNetworks._create_stochastic_mult(A)
    dt2 = @elapsed y = pagerank_power!(y,z,P2,0.85,v,tol,maxiter,MatrixNetworks._noiterfunc)

    @test dt2 <= 2*dt
    @test norm(x-y) <= n*eps(Float64)

    # now test with matrix networks
    P3 = MatrixNetworks._create_stochastic_mult(MatrixNetwork(A))

    y = zeros(n)
    z = zeros(n)
    dt2 = @elapsed y = pagerank_power!(y,z,P3,0.85,v,tol,maxiter,MatrixNetworks._noiterfunc)

    @test dt2 <= 2*dt
    @test norm(x-y) <= n*eps(Float64)

    dt = @elapsed x = pagerank_power!(x,y,P,0.85,v,eps(Float64),1000,MatrixNetworks._noiterfunc)
    dt2 = @elapsed y = pagerank(A,0.85)
    @test norm(x-y) <= n*eps(Float64)
    @test dt2 <= 2*dt

end

function shkexpm(A,t,v)
    P = full(_normout(A'))
    return expm(-t*(eye(size(A,1))-P))*v
end

function heatkernel_test()
    seeded_stochastic_heat_kernel(speye(5),2.,1)
    seeded_stochastic_heat_kernel(speye(5),2.,Set([1]))
    seeded_stochastic_heat_kernel(speye(5),2.,Dict{Int,Float64}(1 => 1.))
    seeded_stochastic_heat_kernel(speye(5),2.,sparsevec([1],[1.],5))
    seeded_stochastic_heat_kernel(speye(5),2.,sparse([1],[1],[1.], 5, 1))
    v = zeros(5)
    v[1] = 1.
    seeded_stochastic_heat_kernel(speye(5),2.,v)
    seeded_stochastic_heat_kernel(speye(5),2.,1./5.)

    n = 5
    v = rand(5)
    vn = v/sum(v)
    x = seeded_stochastic_heat_kernel(speye(5),2.,vn)
    @test norm(x-vn,1) <= n*eps(Float64)

    n = 25
    A = spones(sprand(n,n,2/n))
    v = ones(n)./n
    x = seeded_stochastic_heat_kernel(A,2.,v)
    xtrue = shkexpm(A,2.,v)
    @test norm(x -xtrue,1) <= n*eps(Float64)

    x = seeded_stochastic_heat_kernel(MatrixNetwork(A),2.,v)
    @test norm(x -xtrue,1) <= n*eps(Float64)

    n = 10
    A = spdiagm(ones(n-1),-1,n,n)'
    v = zeros(n)
    v[1] = 1.
    t = 2.
    x = seeded_stochastic_heat_kernel(A,t,v)
    z = zeros(n)
    z[1] = exp(-t)
    for i=2:n
        z[i] = z[i-1]*t./(i-1.)
    end
    @test norm(x-z,1) <= n*eps(Float64)

    @test abs(seeded_stochastic_heat_kernel(spzeros(1,1),2.,1)[1] - exp(-2.)) <= 10*eps(Float64)
    @test abs(seeded_stochastic_heat_kernel(spzeros(1,1),5.,1)[1] - exp(-5.)) <= 10*eps(Float64)
    @test abs(seeded_stochastic_heat_kernel(speye(1,1),2.,1)[1] - 1.) <= 10*eps(Float64)
end

function diffusions_test()

heatkernel_test()

pagerank_error_test()
pagerank_test()
pagerank_alg_test()
pagerank_perf_test()

return true

end


text:
using MOCNeutronTransport
@testset "Triangle_2D" begin
    for F in [Float32, Float64]
        @testset "Constructors" begin
            p₁ = Point_2D(F, 0)
            p₂ = Point_2D(F, 1)
            p₃ = Point_2D(F, 1, 1)
            tri = Triangle_2D(SVector(p₁, p₂, p₃))
            @test tri.points == SVector(p₁, p₂, p₃)

            # single constructor
            tri = Triangle_2D(p₁, p₂, p₃)
            @test tri.points == SVector(p₁, p₂, p₃)
        end

        @testset "Methods" begin
            p₁ = Point_2D(F, 0)
            p₂ = Point_2D(F, 1)
            p₃ = Point_2D(F, 0, 1)
            tri = Triangle_2D(SVector(p₁, p₂, p₃))

            # interpolation
            @test tri(0, 0) ≈ p₁
            @test tri(1, 0) ≈ p₂
            @test tri(0, 1) ≈ p₃
            @test tri(1//2, 1//2) ≈ Point_2D(F, 1//2, 1//2)

            # area
            a = area(tri)
            @test typeof(a) == F
            @test a == F(1//2)

            # in
            p = Point_2D(F, 1//2, 1//10)
            @test p ∈  tri
            p = Point_2D(F, 1//2, -1//10)
            @test p ∉ tri

            # intersect
            # 3 intersections
            l = LineSegment_2D(p₁, Point_2D(F, 1, 1))
            ipoints, points = intersect(l, tri)
            @test ipoints == 3
            @test points[1] ≈ p₁
            @test points[2] ≈ Point_2D(F, 1//2, 1//2)
            @test points[3] ≈ p₁

            # 2 intersections
            l = LineSegment_2D(Point_2D(F, 0, 1//2), Point_2D(F, 1//2, 0))
            ipoints, points = intersect(l, tri)
            @test ipoints == 2
            @test points[1] ≈ Point_2D(F, 1//2, 0)
            @test points[2] ≈ Point_2D(F, 0, 1//2)

            # 0 intersections
            l = LineSegment_2D(Point_2D(F, -1, -1), Point_2D(F, 2, -1))
            ipoints, points = intersect(l, tri)
            @test ipoints == 0
        end
    end
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport