text:
function solveMUMPS(A::SparseMatrixCSC{T1}, rhs::AbstractArray{T2}, sym::Int=0,ooc::Int=0,tr::Int=0) where {T1,T2}

	x = zeros(promote_type(T1,T2),size(rhs))

	return solveMUMPS!(A,rhs,x,sym,ooc,tr)
end

function solveMUMPS!(A::SparseMatrixCSC, rhs::AbstractArray, x::Array, sym::Int=0,ooc::Int=0,tr::Int=0)

# factorization
factor = factorMUMPS(A,sym,ooc)

# solve system
x = applyMUMPS!(factor,rhs,x,tr)

# free memory
destroyMUMPS(factor)
return x
end

function factorMUMPS(A::SparseMatrixCSC{ComplexF64},sym::Int=0,ooc::Int=0)
    # Generate LU-factorization of complex matrix A, 'sym' can be: 0=unsymmetric, 1=symm. pos def, 2=general symmetric
    if size(A,1) != size(A,2)
	error("factorMUMPS: Matrix must be square!")
    end
    n  = size(A,1);
    mumpsstat = [0];
    facTime = @elapsed begin
    p  = ccall( (:factor_mumps_cmplx_, MUMPSlibPath),
    		 Int64, ( Ref{Int64}, Ref{Int64}, Ref{Int64}, Ref{ComplexF64}, Ref{Int64}, Ref{Int64}, Ref{Int64}),
                 n, sym, ooc,  convert(Ref{ComplexF64}, A.nzval), A.rowval, A.colptr, mumpsstat);
    checkMUMPSerror(mumpsstat);
	end
    factor = MUMPSfactorization(p,myid(),n,A.nzval[1],facTime);
    return factor
end

function factorMUMPS(A::SparseMatrixCSC{Float64},sym::Int=0,ooc::Int=0)
    # Generate LU-factorization of a real matrix A, 'sym' can be: 0=unsymmetric, 1=symm. pos def, 2=general symmetric
    if size(A,1) != size(A,2)
		error("factorMUMPS: Matrix must be square!")
    end
    n  = size(A,1);
    mumpsstat = [0];
    facTime = @elapsed begin
    p  = ccall( (:factor_mumps_, MUMPSlibPath),
 	       Int64, ( Ref{Int64}, Ref{Int64}, Ref{Int64}, Ref{Float64}, Ref{Int64}, Ref{Int64}, Ref{Int64}),
               n, sym, ooc,  A.nzval, A.rowval, A.colptr, mumpsstat);
    checkMUMPSerror(mumpsstat);
	end
    factor = MUMPSfactorization(p,myid(),n,A.nzval[1],facTime);
    return factor
end

function checkMUMPSerror(mumpsstat)
	# check if MUMPS reported error

	if mumpsstat[1]==-10
	     error("MUMPS: Numerically singular matrix.");
	elseif mumpsstat[1]==-13
	     error("MUMPS: memory allocation error")
	elseif mumpsstat[1]==-40
	     error("MUMPS: matrix is not positive definite")
	elseif mumpsstat[1]==-90
	     error("MUMPS: Error in out-of-core management.Probably there is not enough disk space.")
	elseif mumpsstat[1]<0
	     error( @sprintf("MUMPS: error --> %d <--. Please refer to Ch. 7 of MUMPS User's guide!",mumpsstat[1]))
	end
end

function applyMUMPS(factor::MUMPSfactorization{T1},rhs::AbstractArray{T2,N},x::Array=zeros(promote_type(T1,T2),size(rhs)),tr::Int=0) where {T1,T2,N}

	id1 = myid(); id2 = factor.worker
	if id1 != id2
		@warn "Worker $id1 has no access to MUMPS factorization stored on $id2. Trying to remotecall!"
		return remotecall_fetch(applyMUMPS,factor.worker,factor,rhs,x,tr)::Array{promote_type(T1,T2),N}
	end

	if size(rhs,1) != factor.n;
		error("applyMUMPS: wrong size of rhs, size(A)=$(factor.n), size(rhs)=$(size(rhs,1)) x $(size(rhs,2)).");
	end

	if size(x)!=size(rhs);
		if isempty(x)
			x = zeros(promote_type(T1,T2),size(rhs))
		else
			error("applyMUMPS: wrong size of x, size(A)=$(factor.n), size(rhs)=$(size(rhs)), size(x)=$(size(x)) provided");
		end
	end

	return applyMUMPS!(factor,rhs,x,tr)::Array{promote_type(T1,T2),N}
end

function applyMUMPS!(factor::MUMPSfactorization{Float64},rhs::Array{Float64},
                      x::Array{Float64},tr::Int=0)

	 nrhs = size(rhs,2)
	ptr = factor.ptr
	ccall( (:solve_mumps_, MUMPSlibPath),
				Int64, (Ref{Int64}, Ref{Int64}, Ref{Float64}, Ref{Float64}, Ref{Int64} ),
							ptr,      nrhs, 	 	rhs,      x, tr)
	return x
end

function applyMUMPS!(factor::MUMPSfactorization{Float64},rhs::SparseMatrixCSC{Float64}, x::Array{Float64,2},tr::Int=0)
	nrhs = size(rhs,2)
	nzrhs = nnz(rhs)
	ptr = factor.ptr
	ccall( (:solve_mumps_sparse_rhs_, MUMPSlibPath),
              Nothing, (Ref{Int64}, Ref{Int64}, Ref{Int64}, Ref{Float64}, Ref{Int64}, Ref{Int64},
	      	  Ref{Float64}, Ref{Int64} ),
		      ptr, nzrhs, nrhs, rhs.nzval, rhs.rowval, rhs.colptr, x, tr)

	return x
end

function applyMUMPS!(factor::MUMPSfactorization{Float64},rhs::SparseVector{Float64}, x::Array{Float64,1},tr::Int=0)
	nrhs = 1
	nzrhs = nnz(rhs)
	colptr = [1; nzrhs+1]
	ptr = factor.ptr
	ccall( (:solve_mumps_sparse_rhs_, MUMPSlibPath),
              Nothing, (Ref{Int64}, Ref{Int64}, Ref{Int64}, Ref{Float64}, Ref{Int64}, Ref{Int64},
	      	  Ref{Float64}, Ref{Int64} ),
		      ptr, nzrhs, nrhs, rhs.nzval, rhs.nzind, colptr, x, tr)

	return x
end

function applyMUMPS!(factor::MUMPSfactorization{ComplexF64},rhs::Array{ComplexF64},
                      x::Array{ComplexF64},tr::Int=0)
	n     = size(rhs,1)
	nrhs  = size(rhs,2)
	ptr = factor.ptr
	ccall( (:solve_mumps_cmplx_, MUMPSlibPath),
				Int64, (Ref{Int64}, Ref{Int64}, Ref{ComplexF64}, Ref{ComplexF64}, Ref{Int64} ),
				ptr,        nrhs,   convert(Ref{ComplexF64}, rhs),   convert(Ref{ComplexF64}, x),   tr)
	return x
end

function applyMUMPS!(factor::MUMPSfactorization{ComplexF64},rhs::SparseMatrixCSC{ComplexF64},
                      x::Array{ComplexF64,2},tr::Int=0)
nrhs = size(rhs,2)
nzrhs = nnz(rhs)
ptr = factor.ptr
ccall( (:solve_mumps_cmplx_sparse_rhs_, MUMPSlibPath),
			Nothing, (Ref{Int64}, Ref{Int64}, Ref{Int64}, Ref{ComplexF64}, Ref{Int64},
			Ref{Int64}, Ref{ComplexF64}, Ref{Int64} ),
			ptr, nzrhs,nrhs, convert(Ref{ComplexF64}, pointer(rhs.nzval)), rhs.rowval,
			rhs.colptr, convert(Ref{ComplexF64}, pointer(x)), tr)
return x
end

function applyMUMPS!(factor::MUMPSfactorization{ComplexF64},rhs::SparseVector{ComplexF64},
                      x::Array{ComplexF64,1},tr::Int=0)
nrhs = 1
nzrhs = nnz(rhs)
colptr = [1; nzrhs + 1]
ptr = factor.ptr
ccall( (:solve_mumps_cmplx_sparse_rhs_, MUMPSlibPath),
			Nothing, (Ref{Int64}, Ref{Int64}, Ref{Int64}, Ref{ComplexF64}, Ref{Int64},
			Ref{Int64}, Ref{ComplexF64}, Ref{Int64} ),
			ptr, nzrhs,nrhs, convert(Ref{ComplexF64}, pointer(rhs.nzval)), rhs.nzind,
			colptr, convert(Ref{ComplexF64}, pointer(x)), tr)
return x
end

function destroyMUMPS(factor::MUMPSfactorization{Float64})
	 #  free memory
	id1 = myid(); id2 = factor.worker;
	if myid() != factor.worker
		remotecall_fetch(destroyMUMPS,id2,factor)
		return
	end
	ccall( (:destroy_mumps_, MUMPSlibPath),
			Int64, (Ref{Int64}, ), factor.ptr )
	factor.ptr = -1
	factor.n    = -1
	factor.time = -1.0

end

function destroyMUMPS(factor::MUMPSfactorization{ComplexF64})
	 #  free memory
	id1 = myid(); id2 = factor.worker;
	if myid() != factor.worker
		remotecall_fetch(destroyMUMPS,id2,factor)
		return
	end
	 ccall( (:destroy_mumps_cmplx_, MUMPSlibPath),
	 		Int64, (Ref{Int64}, ), factor.ptr )
	factor.ptr = -1
	factor.n    = -1
	factor.time = -1.0

end


text:
#using PyPlot
import Colors: RGB
using PyCall
@pyimport matplotlib.patches as mpatches
@pyimport matplotlib.collections as mcol
@pyimport matplotlib.lines as mlines

type Circle <: Grob
    c::Point
    r::Float64
    style::Style
    # by default it seems svg has stroke=nothing, fill=NC"black", not sure if that is good
    function Circle(p, r; style=Style(stroke=nothing, fill=NC"black"))
        new(p, r, style)
    end
end

function render(go::Circle)
    mpatches.Circle(go.c, go.r; style2kw(go.style)...)
end

type PolyLine <: Grob
    points::Array{Point}
    style::Style
    function PolyLine(points; style=Style())
        new(points, style)
    end
end

function render(go::PolyLine)
    #TODO: likely super slow to have to reshape every time. Think of general way
    # looks like I could do Line2D(go.points...; style2kw_l(go.style)...) the only issue
    # I see with this is that that could be a stupid number of arguments being passed
    n = length(go.points)
    x = zeros(n)
    y = zeros(n)
    for (i, p) in enumerate(go.points)
        x[i] = p.x
        y[i] = p.y
    end
    mlines.Line2D(x, y; style2kw_l(go.style)...)
end

type Polygon <: Grob
    points::Array{Point}
    style::Style
    function Polygon(points; style=Style(fill=NC"black"))
        new(points, style)
    end
end

function render(go::Polygon)
    mpatches.Polygon(go.points; style2kw(go.style)...)
end

ext: jl

max_stars_repo_name: gabrielgellner/OdinSon

max_issues_repo_name: gabrielgellner/OdinSon

max_forks_repo_name: gabrielgellner/OdinSon


text:
#!/usr/bin/env julia

using SegmentTrees
import SegmentTrees: Slice, InternalNode, LeafNode, Segment, SegmentBTree
using Test
using Random: seed!, shuffle!

@testset "All tests" begin
    include("creation_tests.jl")
    include("insertion_deletion_tests.jl")
    include("intersection_tests.jl")
    include("search_tests.jl")
    include("iteration.jl")
    include("low_level_tests.jl")
end


text:
@testset "pcl_square_sum" begin
    using Random; Random.seed!(233)
    y = placeholder(rand(10))
    loss = sum((y-rand(10))^2)
    g = tf.convert_to_tensor(gradients(loss, y))
    function test_f(x0)
        run(sess, g, y=>x0), pcl_square_sum(length(y))
    end
    err1, err2 = test_jacobian(test_f, rand(10); showfig=false)
    @test all(err1.>err2)
end

@testset "pl_hessian" begin
    using Random; Random.seed!(233)
    x = placeholder(rand(10))
    z = (x[1:5]^2 + x[6:end]^3) * sum(x)
    y = [reshape(sum(z[1:3]), (-1,));z]
    loss = sum((y-rand(6))^2)
    H, dW = pcl_hessian(y, x, loss)

    g = gradients(loss, x)
    function test_f(x0)
        y0, H0 = run(sess, [g, H], feed_dict = Dict(
            x=>x0,
            dW => pcl_square_sum(length(y))
        ))
    end
    err1, err2 = test_hessian(test_f, rand(length(x)); showfig = false)
    @test all(err1.>err2)
end

@testset "pcl_linear_op" begin
    using Random; Random.seed!(233)
    x = placeholder(rand(10))
    A = rand(6, 10)
    y = A * x + rand(6)
    loss = sum((y-rand(6))^2)
    J = jacobian(y, x)
    g = gradients(loss, x)
    function test_f(x0)
        gval, Jval = run(sess, [g, J], x=>x0)
        H = pcl_linear_op(Array(Jval'), pcl_square_sum(6))
        return gval, H
    end
    err1, err2 = test_hessian(test_f, rand(length(x)); showfig = false)
    @test all(err1.>err2)
end

@testset "pcl_sparse_solve" begin
    using Random; Random.seed!(233)
    A = sprand(10,10,0.6)
    II, JJ, VV = findnz(A)
    x = placeholder(rand(length(VV)))
    temp = SparseTensor(II,JJ,x,10,10)
    indices = run(sess, temp.o.indices)
    B = RawSparseTensor(constant(indices), x, 10, 10)

    rhs = rand(10)
    sol = B\rhs
    loss = sum(sol^2)
    g = gradients(loss, sol)
    G = gradients(loss, x)

    indices = run(sess, B.o.indices) .+ 1
    function test_f(x0)
        s, gval, Gval = run(sess, [sol, g, G], x=>x0)
        H = pcl_sparse_solve(indices, x0,
            s,
            pcl_square_sum(length(sol)), gval)
        return Gval, H
    end
    err1, err2 = test_hessian(test_f, rand(length(x)); showfig = false)
    @test all(err1.>err2)
end

@testset "pcl_compress" begin
    using Random; Random.seed!(233)
    if !Sys.iswindows()
        indices = [
            1 1
            1 1
            1 1
            2 2
            3 3
        ]
        pl = placeholder(rand(5))
        A = RawSparseTensor(constant(indices)-1, pl, 3, 3)
        B = compress(A)
        loss = sum(B.o.values^2)
        g = gradients(loss, pl)
        function test_f(x0)
            G = run(sess, g, pl=>x0)
            W = pcl_square_sum(3)
            J = pcl_compress(indices)
            H = pcl_linear_op(J, W)
            G, H
        end

        x0 = rand(length(pl))
        err1, err2 = test_hessian(test_f, x0; showfig = false)
        @test all(err1.>err2)
    end
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

import Base: +

"""
    +(args...)
    .+(args...)

Summation. Multiple arguments of either scalar or `NDArray` could be
added together. Note at least the first or second argument needs to be an
`NDArray` to avoid ambiguity of built-in summation.
"""
+(x::NDArray)             = x
+(x::NDArray, y::NDArray) = _plus(x, y)
+(x::NDArray, y::Real)    = _plus_scalar(x, scalar = y)
+(y::Real,    x::NDArray) = _plus_scalar(x, scalar = y)

broadcasted(::typeof(+), x::NDArray{T,N}, y::NDArray{T,M}) where {T,N,M} =
  _broadcast_add(x, y)

"""
    sub_from!(dst::NDArray, args::NDArrayOrReal...)

Subtract a bunch of arguments from `dst`. Inplace updating.
"""
function sub_from!(dst::NDArray, arg::NDArrayOrReal)
  @assert dst.writable
  if isa(arg, Real)
    _minus_scalar(dst, scalar = arg, out = dst)
  else
    _minus!(dst, arg)
  end
  dst
end

import Base: -

"""
    -(x::NDArray)
    -(x, y)
    .-(x, y)

Subtraction `x - y`, of scalar types or `NDArray`.
Or create the negative of `x`.
"""
-(x::NDArray)             = _mul_scalar(x, scalar = -one(eltype(x)))
-(x::NDArray, y::NDArray) = _minus(x, y)
-(x::NDArray, y::Real)    = _minus_scalar(x, scalar = y)
-(y::Real, x::NDArray)    = _rminus_scalar(x, scalar = y)

broadcasted(::typeof(-), x::NDArray{T,N}, y::NDArray{T,M}) where {T,N,M} =
  _broadcast_minus(x, y)

"""
    mul_to!(dst::NDArray, arg::NDArrayOrReal)

Elementwise multiplication into `dst` of either a scalar or an `NDArray` of the same shape.
Inplace updating.
"""
function mul_to!(dst::NDArray, arg::NDArrayOrReal)
  @assert dst.writable
  if isa(arg, Real)
    _mul_scalar(dst, scalar = arg, out = dst)
  else
    _mul(dst, arg, out = dst)
  end
  dst
end

import Base: *

"""
    .*(x, y)

Elementwise multiplication for `NDArray`.
"""
*(x::NDArray, y::Real)  = _mul_scalar(x, scalar = y)
*(y::Real, x::NDArray)  = _mul_scalar(x, scalar = y)

broadcasted(::typeof(*), x::NDArray{T,N}, y::NDArray{T,N}) where {T,N} =
  _mul(x, y)
broadcasted(::typeof(*), x::NDArray{T,N}, y::NDArray{T,M}) where {T,N,M} =
  _broadcast_mul(x, y)

"""
    *(A::NDArray, B::NDArray)

Matrix/tensor multiplication.
"""
*(x::NDArray{T}, y::NDArray{T}) where T = x ⋅ y

LinearAlgebra.adjoint(x::NDArray{T,1}) where T = transpose(x)
LinearAlgebra.adjoint(x::NDArray{T,2}) where T = transpose(x)

"""
    div_from!(dst::NDArray, arg::NDArrayOrReal)

Elementwise divide a scalar or an `NDArray` of the same shape from `dst`. Inplace updating.
"""
function div_from!(dst::NDArray, arg::NDArrayOrReal)
  @assert dst.writable
  if isa(arg, Real)
    _div_scalar(dst, scalar = arg, out = dst)
  else
    _div(dst, arg, out = dst)
  end
  dst
end

function div_from!(dst::NDArray{T}, arg::Real) where {T<:Integer}
  @assert dst.writable
  @assert(round(T, arg) != zero(T), "Integer divided by zero")
  _div_scalar(dst, scalar = arg, out = dst)
  dst
end

"""
    rdiv_from!(x:: Real, y::NDArray)

Elementwise divide a scalar by an `NDArray`. Inplace updating.
"""
function rdiv_from!(x::Real, y::NDArray)
  @assert y.writable
  _rdiv_scalar(y, scalar = x, out = y)
  y
end

import Base: /

"""
    ./(x::NDArray, y::NDArray)
    ./(x::NDArray, y::Real)
    ./(x::Real, y::NDArray)

* Elementwise dividing an `NDArray` by a scalar or another `NDArray`
of the same shape.

* Elementwise divide a scalar by an `NDArray`.

* Matrix division (solving linear systems) is not implemented yet.
"""
/(x::NDArray, y::Real) = _div_scalar(x, scalar = y)

broadcasted(::typeof(/), y::Real, x::NDArray) = _rdiv_scalar(x, scalar = y)
broadcasted(::typeof(/), x::NDArray{T,N}, y::NDArray{T,N}) where {T,N} =
  _div(x, y)
broadcasted(::typeof(/), x::NDArray{T,N}, y::NDArray{T,M}) where {T,N,M} =
  _broadcast_div(x, y)

function broadcasted(::typeof(/), x::NDArray{T}, y::Real) where {T<:Integer}
  @assert(round(T, y) != zero(T), "Integer divided by zero")
  _div_scalar(x, scalar = y)
end

"""
    mod_from!(x::NDArray, y::NDArray)
    mod_from!(x::NDArray, y::Real)

Elementwise modulo for `NDArray`.
Inplace updating.
"""
mod_from!(x::NDArray, y::NDArray) = _mod!(x, y)
mod_from!(x::NDArray, y::Real)    = _mod_scalar!(x, y)

"""
    rmod_from!(y::Real, x::NDArray)

Elementwise modulo for `NDArray`.
Inplace updating.
"""
rmod_from!(y::Real, x::NDArray) = _rmod_scalar!(x, y)

import Base: %

"""
    .%(x::NDArray, y::NDArray)
    .%(x::NDArray, y::Real)
    .%(x::Real, y::NDArray)

Elementwise modulo for `NDArray`.
"""
%(x::NDArray, y::Real) = _mod_scalar(x, y)

broadcasted(::typeof(%), y::Real, x::NDArray) = _rmod_scalar(x, y)
broadcasted(::typeof(%), x::NDArray{T,N}, y::NDArray{T,N}) where {T,N} =
  _mod(x, y)
broadcasted(::typeof(%), x::NDArray{T,N}, y::NDArray{T,M}) where {T,N,M} =
  _broadcast_mod(x, y)

# document of `.^` is merged into SymbolicNode's

broadcasted(::typeof(Base.literal_pow), ::typeof(^), x::NDArray, ::Val{s}) where {s} =
  _power_scalar(x, scalar = s)
broadcasted(::typeof(^), x::NDArray, s::Real) = _power_scalar(x,  scalar = s)
broadcasted(::typeof(^), s::Real, x::NDArray) = _rpower_scalar(x, scalar = s)

broadcasted(::typeof(^), ::Irrational{:ℯ}, x::NDArray) = exp(x)
broadcasted(::typeof(^), x::NDArray, s::Irrational)    = _power_scalar(x, scalar = s)
broadcasted(::typeof(^), s::Irrational, x::NDArray)    = _rpower_scalar(x, scalar = s)

broadcasted(::typeof(^), x::NDArray{T,N}, y::NDArray{T,N}) where {T,N} =
  _power(x, y)
broadcasted(::typeof(^), x::NDArray{T,N}, y::NDArray{T,M}) where {T,N,M} =
  _broadcast_power(x, y)

_nddoc[:clip] = _nddoc[:clip!] =
"""
    clip(x::NDArray, min, max)
    clip!(x::NDArray, min, max)

Clips (limits) the values in `NDArray`.
Given an interval, values outside the interval are clipped to the interval edges.
Clipping `x` between `min` and `x` would be:

```julia
clip(x, min_, max_) = max(min(x, max_), min_))
```

```jldoctest
julia> x = NDArray(1:9);

julia> mx.clip(x, 2, 8)'
1×9 mx.NDArray{Int64,2} @ CPU0:
 2  2  3  4  5  6  7  8  8
```

The storage type of clip output depends on storage types of inputs and the
`min`, `max` parameter values:

- clip(default) = default
- clip(row_sparse, min <= 0, max >= 0) = row_sparse
- clip(csr, min <= 0, max >= 0) = csr
- clip(row_sparse, min < 0, max < 0) = default
- clip(row_sparse, min > 0, max > 0) = default
- clip(csr, min < 0, max < 0) = csr
- clip(csr, min > 0, max > 0) = csr
"""
@_remap clip(x::NDArray, min::Real, max::Real) clip(x; a_min = min, a_max = max)
@_remap clip!(x::NDArray, min::Real, max::Real) clip(x; a_min = min, a_max = max)

################################################################################
# remapping to solving type unstablility
################################################################################

@_remap _plus(x::NDArray, y::NDArray)  _plus(x, y)
@_remap _plus!(x::NDArray, y::NDArray) _plus(x, y)

@_remap _minus(x::NDArray, y::NDArray)  _minus(x, y)
@_remap _minus!(x::NDArray, y::NDArray) _minus(x, y)

@_remap _mod(x::NDArray, y::NDArray)  _mod(x, y)
@_remap _mod!(x::NDArray, y::NDArray) _mod(x, y)

@_remap _mod_scalar(x::NDArray, y::Real)  _mod_scalar(x; scalar = y)
@_remap _mod_scalar!(x::NDArray, y::Real) _mod_scalar(x; scalar = y)

@_remap _rmod_scalar(x::NDArray, y::Real)  _rmod_scalar(x; scalar = y)
@_remap _rmod_scalar!(x::NDArray, y::Real) _rmod_scalar(x; scalar = y)

@_remap _broadcast_add(x::NDArray, y::NDArray)  broadcast_add(x, y)
@_remap _broadcast_add!(x::NDArray, y::NDArray) broadcast_add(x, y)

@_remap _broadcast_minus(x::NDArray, y::NDArray)  broadcast_minus(x, y)
@_remap _broadcast_minus!(x::NDArray, y::NDArray) broadcast_minus(x, y)

@_remap _broadcast_mul(x::NDArray, y::NDArray)  broadcast_mul(x, y)
@_remap _broadcast_mul!(x::NDArray, y::NDArray) broadcast_mul(x, y)

@_remap _broadcast_div(x::NDArray, y::NDArray)  broadcast_div(x, y)
@_remap _broadcast_div!(x::NDArray, y::NDArray) broadcast_div(x, y)

@_remap _broadcast_mod(x::NDArray, y::NDArray)  broadcast_mod(x, y)
@_remap _broadcast_mod!(x::NDArray, y::NDArray) broadcast_mod(x, y)

@_remap _broadcast_power(x::NDArray, y::NDArray)  broadcast_power(x, y)
@_remap _broadcast_power!(x::NDArray, y::NDArray) broadcast_power(x, y)

max_stars_repo_name: HyperZealot/mxnet

max_issues_repo_name: HyperZealot/mxnet

max_forks_repo_name: HyperZealot/mxnet


text:
using MinimalPerfectHash
using Test
using Random

@testset "MinimalPerfectHash.jl" begin

    @testset "EmptyConstructor" begin
        chd = MinimalPerfectHash.CHD()
        @test typeof(chd) == MinimalPerfectHash.CHD{Any, Any}
        @test length(chd) == 0
        for (_, _) in chd
            @test false
        end
        chd = MinimalPerfectHash.CHD( () )
        @test typeof(chd) == MinimalPerfectHash.CHD{Any, Any}
        @test length(chd) == 0
        for (_, _) in chd
            @test false
        end

    end

    let
        chd = MinimalPerfectHash.CHD(1=>"Bye")
        @test typeof(chd) == MinimalPerfectHash.CHD{Int, String}
        chd = MinimalPerfectHash.CHD("Hello" => "World")
        @test typeof(chd) == MinimalPerfectHash.CHD{String, String}
        MinimalPerfectHash.CHD{UInt, Int}(0x1 => 2)
        MinimalPerfectHash.CHD(UInt8(1) => 2, UInt16(2) => 3, UInt(3) => 10)
        chd2 = copy(chd)
        @test chd2 !== chd
        @test chd2 == chd
    end

    @testset "Exceptions" begin
        @test_throws ArgumentError MinimalPerfectHash.CHD(1=>2, 1=>3)
        @test_throws ArgumentError MinimalPerfectHash.CHD([1,2,3])
    end

    for K in [Int, UInt8, UInt64, Float64]
        for V in [Int, Float64]
            for n in [0, 1, 2, 3, 4, 5, 8, 9, 27, 32, 33, 64, 65, 1024, 1025, 65536, 65537]
                ks = collect(Set(Random.rand(K, n)))
                vs = rand(V, length(ks))
                chd = MinimalPerfectHash.CHD((k, v) for (k, v) in zip(ks, vs))
                chd = MinimalPerfectHash.CHD{K, V}((k, v) for (k, v) in zip(ks, vs))

                @test all(haskey(chd, k) for k in ks)
                @test all((k => v) in chd for (k, v) in zip(ks, vs))
                @test all(chd[k] == v for (k, v) in zip(ks, vs))
                @test length(chd) == length(ks)
                @test all(get(chd, k, nothing) == v for (k, v) in zip(ks, vs))
                @test all(getkey(chd, k, nothing) == k for k in ks)
                @test all(k in keys(chd) for k in ks)
                if n == 0
                    @test isempty(chd)
                else
                    @test sum(1 for (_, _) in chd) == length(ks)
                end
            end
        end
    end

    let K=Int, V=Int, n=(1<<3)
        ks = sort(unique(Random.rand(K, n)))
        vs = rand(V, length(ks))
        chd = MinimalPerfectHash.CHD((k, v) for (k, v) in zip(ks, vs))
        chd = MinimalPerfectHash.CHD{K, V}((k, v) for (k, v) in zip(ks, vs))

        @test all(haskey(chd, k) for k in ks)
        @test all((k => v) in chd for (k, v) in zip(ks, vs))
        @test all(chd[k] == v for (k, v) in zip(ks, vs))
        @test length(chd) == length(ks)
        @test all(get(chd, k, nothing) == v for (k, v) in zip(ks, vs))
        @test all(getkey(chd, k, nothing) == k for k in ks)
        @test all(k in keys(chd) for k in ks)
        if n == 0
            @test isempty(chd)
        else
            @test sum(1 for (_, _) in chd) == length(ks)
        end
    end

end

max_stars_repo_name: kyungminlee/MinimalPerfectHashing

max_issues_repo_name: kyungminlee/MinimalPerfectHashing

max_forks_repo_name: kyungminlee/MinimalPerfectHashing


text:
## Use this functions to make creating VirtualArrays easier.

function virtual_cat{T}(expanded_dim::Int, parents::AbstractArray{T}...)
    dim = max(expanded_dim, largest_dim(parents...))
    return VirtualArray{T, dim}(expanded_dim, parents...)
end

function virtual_cat(expanded_dim::Int)
    return virtual_cat(expanded_dim, [])
end

function virtual_vcat{T}(parents::AbstractArray{T}...)
    return virtual_cat(1, parents...)
end

function virtual_hcat{T}(parents::AbstractArray{T}...)
    return virtual_cat(2, parents...)
end

function largest_dim(parents::AbstractArray...)
    largest = 0
    for parent in parents
        largest = max(largest, length(size(parent)))
    end
    return largest
end


text:
using Test

@testset "GraphicalLinearAlgebra" begin
  include("GLA.jl")
  include("StructuredGLA.jl")
end


text:
function reverse(mg::MetaDiGraph{T,U}) where T <: Integer where U <: Real
    rg = reverse(mg.graph)
    map = nv(mg):-1:1
    rvprops = Dict{T,PropDict}()
    reprops = Dict{SimpleEdge{T},PropDict}()
    rgprops = mg.gprops
    rweightfield = mg.weightfield
    rdefaultweight = mg.defaultweight
    rindices = mg.indices

    for v in keys(mg.vprops)
        rvprops[map[v]] = mg.vprops[v]
    end

    for e in keys(mg.eprops)
        reprops[reverse(e)] = mg.eprops[e]
    end

    rmg = MetaDiGraph(rg,
        rvprops,
        reprops,
        rgprops,
        rweightfield,
        rdefaultweight,
        MetaDict(),
        Set{Symbol}())
    for p in mg.indices
        set_indexing_prop!(rmg, p)
    end

    return rmg
end


text:
using Base.Test
using MathProgBase
using MathProgBase.SolverInterface

function quadprogtest(solver=MathProgBase.defaultQPsolver)
    println("Testing quadprog with solver ", string(typeof(solver)))

    sol = quadprog([0., 0., 0.],[2. 1. 0.; 1. 2. 1.; 0. 1. 2.],[1. 2. 3.; 1. 1. 0.],'>',[4., 1.],-Inf,Inf,solver)
    @test sol.status == :Optimal
    @test_approx_eq_eps sol.objval 130/70 1e-6
    @test_approx_eq_eps norm(sol.sol[1:3] - [0.5714285714285715,0.4285714285714285,0.8571428571428572]) 0.0 1e-6

    let
    m = LinearQuadraticModel(solver)
    loadproblem!(m, [1. 2. 3.; 1. 1. 0.],[-Inf,-Inf,-Inf],[Inf,Inf,Inf],[0.,0.,0.],[4., 1.],[Inf,Inf], :Min)

    setquadobj!(m,diagm([10.0,10.0,10.0]))
    rows = [1, 2, 2, 2, 3, 3, 3]
    cols = [1, 1, 1, 2, 2, 3, 3]
    vals = Float64[2, 0.5, 0.5, 2, 1, 1, 1]
    setquadobj!(m,rows,cols,vals)
    optimize!(m)
    stat = status(m)
    @test stat == :Optimal
    @test_approx_eq_eps getobjval(m) 130/70 1e-6
    @test_approx_eq_eps norm(getsolution(m) - [0.5714285714285715,0.4285714285714285,0.8571428571428572]) 0.0 1e-6
    end

    let
    m = LinearQuadraticModel(solver)
    loadproblem!(m, [-1. 1.; 1. 1.], [0.,0.], [Inf,Inf], [1.,1.], [0.,0.], [Inf,Inf], :Max)
    addquadconstr!(m, [2], [1.], [1], [1], [1.], '<', 2)
    optimize!(m)
    stat = status(m)
    @test stat == :Optimal
    @test_approx_eq_eps getobjval(m) 2.25 1e-6
    @test_approx_eq_eps norm(getsolution(m) - [0.5,1.75]) 0.0 1e-3
    end
    println("Done")
end

function qpdualtest(solver=MathProgBase.defaultQPsolver)
    println("Testing QP duals with solver ", string(typeof(solver)))
    # max x
    # s.t. x^2 <= 2
    m = LinearQuadraticModel(solver)
    loadproblem!(m, Array(Float64,0,1), [-Inf], [Inf], [1.0], Float64[], Float64[], :Max)
    addquadconstr!(m, [], [], [1], [1], [1.0], '<', 2.0)
    optimize!(m)
    stat = status(m)

    @test numlinconstr(m) == 0
    @test numquadconstr(m) == 1
    @test numconstr(m) == 1
    @test stat == :Optimal
    @test_approx_eq_eps getobjval(m) sqrt(2) 1e-6
    @test_approx_eq_eps getsolution(m)[1] sqrt(2) 1e-6
    @test_approx_eq_eps getquadconstrduals(m)[1] 0.5/sqrt(2) 1e-6

    # min -x
    # s.t. x^2 <= 2
    m = LinearQuadraticModel(solver)
    loadproblem!(m, Array(Float64,0,1), [-Inf], [Inf], [-1.0], Float64[], Float64[], :Min)
    addquadconstr!(m, [], [], [1], [1], [1.0], '<', 2.0)
    optimize!(m)
    stat = status(m)

    @test numlinconstr(m) == 0
    @test numquadconstr(m) == 1
    @test numconstr(m) == 1
    @test stat == :Optimal
    @test_approx_eq_eps getobjval(m) -sqrt(2) 1e-6
    @test_approx_eq_eps getsolution(m)[1] sqrt(2) 1e-6
    @test_approx_eq_eps getquadconstrduals(m)[1] -0.5/sqrt(2) 1e-6

    println("Done")
end

function socptest(solver=MathProgBase.defaultQPsolver)
    println("Testing SOCP interface with solver ", string(typeof(solver)))

    # min t
    # s.t. x + y >= 1
    #      x^2 + y^2 <= t^2
    #      t >= 0
    m = LinearQuadraticModel(solver)
    loadproblem!(m, [ 1.0 1.0 0.0 ], [-Inf,-Inf,0.0], [Inf,Inf,Inf], [0.0,0.0,1.0], [1.0],[Inf], :Min)
    addquadconstr!(m, [], [], [1,2,3], [1,2,3], [1.0,1.0,-1.0],'<',0.0)
    optimize!(m)
    stat = status(m)

    @test stat == :Optimal
    @test_approx_eq_eps getobjval(m) sqrt(1/2) 1e-6
    @test_approx_eq_eps norm(getsolution(m) - [0.5,0.5,sqrt(1/2)]) 0.0 1e-3
    println("Done")

end


text:
module GADM

using DataDeps
using ArchGDAL
using Logging
using Tables
using GeoInterface

function __init__()
    # make sure geometries are always downloaded
    # without user interaction from DataDeps.jl
    ENV["DATADEPS_ALWAYS_ACCEPT"] = true
end

"""
    isvalidcode(str)

Tells whether or not `str` is a valid
ISO 3166 Alpha 3 country code. Valid
code examples are "IND", "USA", "BRA".
"""
isvalidcode(str) = match(r"\b[A-Z]{3}\b", str) !== nothing

"""
    download(country)

Downloads data for `country` using DataDeps.jl and returns path.
"""
function download(country)
    ID = "GADM_$country"
    try
        # if data is already on disk
        # we just return the path
        @datadep_str ID
    catch KeyError
        # otherwise we register the data
        # and download using DataDeps.jl
        register(DataDep(ID,
            "Geographic data for country $country provided by the https://gadm.org project.",
            "https://biogeo.ucdavis.edu/data/gadm3.6/gpkg/gadm36_$(country)_gpkg.zip",
            post_fetch_method=DataDeps.unpack))
        @datadep_str ID
    end
end

"""
    dataread(path)

Read data in `path` returned by [`download`](@ref).
"""
function dataread(path)
    files = readdir(path; join=true)

    isgpkg(f) = last(splitext(f)) == ".gpkg"
    gpkg = files[findfirst(isgpkg, files)]

    data = ArchGDAL.read(gpkg)

    !isnothing(data) ? data : throw(ArgumentError("failed to read data from disk"))
end

"""
    getdataset(country)

Downloads and extracts dataset of the given country code
"""
function getdataset(country)
    isvalidcode(country) || throw(ArgumentError("please provide standard ISO 3 country codes"))
    data = country |> download |> dataread
end

"""
    getlayer(data, level)

Get layer of the desired `level` from the `data`.
"""
function getlayer(data, level)
    nlayers = ArchGDAL.nlayer(data)
    for l = 0:nlayers - 1
        layer = ArchGDAL.getlayer(data, l)
        lname = ArchGDAL.getname(layer)
        llevel = last(split(lname, "_"))
        string(level) == llevel && return layer
    end
    throw(ArgumentError("asked for level $(level), valid levels are 0-$(nlayers - 1)"))
end

"""
    get(country, subregions...; children=false)
Returns a Tables.jl columntable for the requested region
Geometry of the region(s) can be accessed with the key `geom`
The geometries are GeoInterface compliant Polygons/MultiPolygons.

1. country: ISO 3166 Alpha 3 country code
2. subregions: Full official names in hierarchial order (provinces, districts, etc.)
3. children: When true, function returns two columntables -> parent, children.
Eg. when children is set true when querying just the country,
second return parameter are the states/provinces.

## Examples

```julia
# columntable of size 1, data of India's borders
data = get("IND")
# parent -> state data, children -> table of all districts inside Uttar Pradesh
parent, children = get("IND", "Uttar Pradesh"; children=true)
```
"""
function get(country, subregions...; children=false)
    data = getdataset(country)
    nlayers = ArchGDAL.nlayer(data)

    function filterlayer(layer, key, value, all=false)
        filtered = []
        for row in layer
            index = ArchGDAL.findfieldindex(row, Symbol(key))
            field = ArchGDAL.getfield(row, index)
            if all || occursin(lowercase(value), lowercase(field))
                push!(filtered, row)
            end
        end
        return Tables.columntable(filtered)
    end

    # p -> parent, is the requested region
    plevel = length(subregions)
    plevel >= nlayers && throw(ArgumentError("more subregions provided than actual"))
    pname = isempty(subregions) ? "" : last(subregions)
    player = getlayer(data, plevel)
    p = filterlayer(player, "NAME_$(plevel)", pname, iszero(plevel))
    isempty(p) && throw(ArgumentError("could not find required region"))

    !children && return p

    # c -> children, is the region 1 level lower than p
    clevel = plevel + 1
    clevel == nlayers && return (p, Tables.rowtable([]))
    clayer = getlayer(data, clevel)
    c = filterlayer(clayer, "NAME_$(plevel)", pname, iszero(plevel))

    return p, c
end

"""
    coordinates(country, subregions...)

Returns a deep array of coordinates for the requested region.
Input: ISO 3166 Alpha 3 Country Code, and further full official names of subdivisions

## Example:

```julia
c = coordinates("VAT") # Returns a deep array of Vatican city

1-element Array{Array{Array{Float64,1},1},1}:
 [[12.455550193786678, 41.90755081176758], ..., [12.454191207885799, 41.90721130371094], [12.455550193786678, 41.90755081176758]]

c = coordinates("IND") # returns an array of coordinates of India's boundary

763-element Array{Array{Array{Array{Float64,1},1},1},1}:
 [[[93.78772736, 6.85264015], [93.7884903, 6.85257101], ...]]

 ```
"""
function coordinates(country, subregions...)
    p = get(country, subregions...)
    GeoInterface.coordinates(p.geom[1])
end

end


text:
#
# This file contains all structre types and functions necessary for building
# the CellList and CellListPair structures.
#

"""

$(TYPEDEF)

Internal function or structure - interface may change.

# Extended help

$(TYPEDFIELDS)

Copies particle coordinates and associated index, to build contiguous particle lists
in memory when building the cell lists. This strategy duplicates the particle coordinates
data, but is probably worth the effort.

"""
struct ParticleWithIndex{N,T}
    index::Int
    real::Bool
    coordinates::SVector{N,T}
end
Base.zero(::Type{ParticleWithIndex{N,T}}) where {N,T} =
    ParticleWithIndex{N,T}(0,false,zeros(SVector{N,T}))

"""

$(TYPEDEF)

Internal function or structure - interface may change.

# Extended help

$(TYPEDFIELDS)

Structure to define the number of batches used in the parallel splitting of the calculations
of the cell list construction and of the `map_pairwise` computation. It is initialized with
a standard heuristic that returns at most the number of threads, but may return a smaller
number if the system is small. The two parameters can be tunned for optimal performance of each
step of the calculation (cell list construction and mapping of interactions). The construction
of the cell lists require a larger number of particles for threading to be effective, Thus
by default the system size that allows multi-threading is greater for this part of the calculation.

"""
struct NumberOfBatches
    build_cell_lists::Int
    map_computation::Int
end
NumberOfBatches(t::Tuple{Int,Int}) = NumberOfBatches(t[1],t[2])
Base.zero(::Type{NumberOfBatches}) = NumberOfBatches(0,0)
Base.iszero(x::NumberOfBatches) = (iszero(x.build_cell_lists) && iszero(x.map_computation))
function Base.show(io::IO,::MIME"text/plain",nbatches::NumberOfBatches)
    println(io,typeof(nbatches))
    println(io,"  Number of batches for cell list construction: $(nbatches.build_cell_lists)")
    print(io,"  Number of batches for function mapping: $(nbatches.map_computation)")
end

"""

$(TYPEDEF)

Internal function or structure - interface may change.

# Extended help

$(TYPEDFIELDS)

This structure contains the cell linear index and the information
about if this cell is in the border of the box (such that its
neighboring cells need to be wrapped)

"""
Base.@kwdef struct Cell{N,T}
    linear_index::Int = 0
    cartesian_index::CartesianIndex{N} = CartesianIndex{N}(ntuple(i->0,N))
    center::SVector{N,T} = zeros(SVector{N,T})
    contains_real::Bool = false
    n_particles::Int = 0
    particles::Vector{ParticleWithIndex{N,T}} = Vector{ParticleWithIndex{N,T}}(undef,0)
end
function Cell{N,T}(cartesian_index::CartesianIndex,box::Box; sizehint::Int=0) where {N,T}
    return Cell{N,T}(
        linear_index=cell_linear_index(box.nc,cartesian_index),
        cartesian_index=cartesian_index,
        center=cell_center(cartesian_index,box),
        particles = Vector{ParticleWithIndex{N,T}}(undef,sizehint)
    )
end

"""

$(TYPEDEF)

Internal function or structure - interface may change.

# Extended help

$(TYPEDFIELDS)

Auxiliary structure to contain projected particles. Types of
scalars are chosen such that with a `SVector{3,Float64}` the
complete struct has 32bytes.

"""
Base.@kwdef struct ProjectedParticle{N,T}
    index::Int
    xproj::T
    coordinates::SVector{N,T}
end

"""

$(TYPEDEF)

Internal function or structure - interface may change.

# Extended help

$(TYPEDFIELDS)

Structure that contains the cell lists information.

"""
Base.@kwdef struct CellList{N,T}
    " Number of real particles. "
    n_real_particles::Int = 0
    " Number of cells. "
    number_of_cells::Int = 0
    " *mutable* number of particles in the computing box. "
    n_particles::Int = 0
    " *mutable* number of cells with real particles. "
    n_cells_with_real_particles::Int = 0
    " *mutable* number of cells with particles, real or images. "
    n_cells_with_particles::Int = 0
    " Auxiliary array that contains the indexes in list of the cells with particles, real or images. "
    cell_indices::Vector{Int} = zeros(Int,number_of_cells)
    " Auxiliary array that contains the indexes in the cells with real particles. "
    cell_indices_real::Vector{Int} = zeros(Int,0)
    " Vector containing cell lists of cells with particles. "
    cells::Vector{Cell{N,T}} = Cell{N,T}[]
    " Number of batches for the parallel calculations. "
    nbatches::NumberOfBatches = zero(NumberOfBatches)
    " Auxiliar array to store projected particles. "
    projected_particles::Vector{Vector{ProjectedParticle{N,T}}} =
        Vector{Vector{ProjectedParticle{N,T}}}(undef,0)
end
function Base.show(io::IO,::MIME"text/plain",cl::CellList)
    println(io,typeof(cl))
    println(io,"  $(cl.n_real_particles) real particles.")
    println(io,"  $(cl.n_cells_with_real_particles) cells with real particles.")
    print(io,"  $(cl.n_particles) particles in computing box, including images.")
end

"""

$(TYPEDEF)

Internal function or structure - interface may change.

# Extended help

$(TYPEDFIELDS)

Structure that will cointain the cell lists of two independent sets of
particles for cross-computation of interactions

"""
@with_kw struct CellListPair{V,N,T}
    ref::V
    target::CellList{N,T}
    swap::Bool
end
function Base.show(io::IO,::MIME"text/plain",cl::CellListPair)
    print(io,typeof(cl),"\n")
    print(io,"   $(length(cl.ref)) particles in the reference vector.\n")
    print(io,"   $(cl.target.n_cells_with_real_particles) cells with real particles of target vector.")
end

"""

```
set_number_of_batches!(cl,nbatches::Tuple{Int,Int}=(0,0))
```

Internal function or structure - interface may change.

# Extended help

Functions that set the default number of batches for the construction of the cell lists,
and mapping computations. This is of course heuristic, and may not be the best choice for
every problem. See the parameter `nbatches` of the construction of the cell lists for
tunning this.

"""
function set_number_of_batches!(cl::CellList{N,T},nbatches::Tuple{Int,Int}=(0,0)) where {N,T}
    nbatches = NumberOfBatches(nbatches)
    if nbatches.build_cell_lists < 1
        n1 = _nbatches_build_cell_lists(cl.n_real_particles)
    else
        n1 = nbatches.build_cell_lists
    end
    if nbatches.map_computation < 1
        n2 = _nbatches_map_computation(cl.n_real_particles)
    else
        n2 = nbatches.map_computation
    end
    nbatches = NumberOfBatches(n1,n2)
    @set! cl.nbatches = nbatches
    for _ in 1:cl.nbatches.map_computation
        push!(cl.projected_particles,Vector{ProjectedParticle{N,T}}(undef,0))
    end
    return cl
end
# Heuristic choices for the number of batches, for an atomic system
_nbatches_build_cell_lists(n::Int) = min(n,min(8,nthreads()))
_nbatches_map_computation(n::Int) = min(n,min(floor(Int,2^(log10(n)+1)),nthreads()))

function set_number_of_batches!(cl::CellListPair{N,T},nbatches::Tuple{Int,Int}=(0,0)) where {N,T}
    nbatches = NumberOfBatches(nbatches)
    if nbatches.build_cell_lists < 1
        n1 = _nbatches_build_cell_lists(cl.target.n_real_particles)
    else
        n1 = nbatches.build_cell_lists
    end
    if nbatches.map_computation < 1
        n2 = _nbatches_map_computation(length(cl.ref))
    else
        n2 = nbatches.map_computation
    end
    nbatches = NumberOfBatches(n1,n2)
    target = cl.target
    @set! target.nbatches = nbatches
    @set! cl.target = target
    return cl
end

"""

```
nbatches(cl)
```

Returns the number of batches for parallel processing that will be used in the pairwise function mappings associated to cell list `cl`.
It returns the `cl.nbatches.map_computation` value. This function is important because it must be used to set the number of copies
of custom preallocated output arrays.

A second argument can be provided, which may be `:map` or `:build`, in which case the function returns either the number of batches used
for pairwise mapping or for the construction of the cell lists. Since this second value is internal and does not affect the interface,
it can be usually ignored.

## Example

```julia-repl
julia> x = rand(3,1000); box = Box([1,1,1],0.1);

julia> cl = CellList(x,box,nbatches=(2,16));

julia> nbatches(cl)
16

julia> nbatches(cl,:map)
16

julia> nbatches(cl,:build)
2
```

"""
nbatches(cl::CellList) = cl.nbatches.map_computation
function nbatches(cl::CellList,s::Symbol)
    s == :map_computation || s == :map && return cl.nbatches.map_computation
    s == :build_cell_lists || s == :build && return cl.nbatches.build_cell_lists
end
nbatches(cl::CellListPair) = nbatches(cl.target)
nbatches(cl::CellListPair,s::Symbol) = nbatches(cl.target,s)

"""

$(TYPEDEF)

Internal function or structure - interface may change.

# Extended help

$(TYPEDFIELDS)

Auxiliary structure to carry threaded lists and ranges of particles to
be considered by each thread on parallel construction.

"""
@with_kw struct AuxThreaded{N,T}
    particles_per_batch::Int
    idxs::Vector{UnitRange{Int}} = Vector{UnitRange{Int}}(undef,0)
    lists::Vector{CellList{N,T}} = Vector{CellList{N,T}}(undef,0)
end
function Base.show(io::IO,::MIME"text/plain",aux::AuxThreaded)
    println(io,typeof(aux))
    print(io," Auxiliary arrays for nbatches = ", length(aux.lists))
end

"""

```
AuxThreaded(cl::CellList{N,T}) where {N,T}
```

Constructor for the `AuxThreaded` type, to be passed to `UpdateCellList!` for in-place
update of cell lists.

## Example
```julia-repl
julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(3) for _ in 1:100_000 ];

julia> cl = CellList(x,box);

julia> aux = CellListMap.AuxThreaded(cl)
CellListMap.AuxThreaded{3, Float64}
 Auxiliary arrays for nthreads = 8

julia> cl = UpdateCellList!(x,box,cl,aux)
CellList{3, Float64}
  100000 real particles.
  31190 cells with real particles.
  1134378 particles in computing box, including images.

```
"""
function AuxThreaded(cl::CellList{N,T};particles_per_batch=10_000) where {N,T}
    nbatches = cl.nbatches.build_cell_lists
    aux = AuxThreaded{N,T}(
        particles_per_batch=particles_per_batch,
        idxs=Vector{UnitRange{Int}}(undef,nbatches),
        lists=Vector{CellList{N,T}}(undef,nbatches)
    )
    # If the calculation is not parallel, no need to initialize this
    nbatches == 1 && return aux
    for ibatch in 1:nbatches
        cl_batch = CellList{N,T}(
            n_real_particles=particles_per_batch, # this is reset before filling, in UpdateCellList!
            number_of_cells=cl.number_of_cells,
        )
        aux.lists[ibatch] = cl_batch
    end
    # Indices of the atoms that will be added by each thread
    nrem = cl.n_real_particles%nbatches
    nperthread = (cl.n_real_particles-nrem)÷nbatches
    first = 1
    for ibatch in 1:nbatches
        nx = nperthread
        if ibatch <= nrem
            nx += 1
        end
        aux.idxs[ibatch] = first:(first-1)+nx
        first += nx
    end
    return aux
end

"""

```
AuxThreaded(cl::CellListPair{N,T}) where {N,T}
```

Constructor for the `AuxThreaded` type for lists of disjoint particle sets,
to be passed to `UpdateCellList!` for in-place update of cell lists.

## Example
```julia-repl
julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(3) for i in 1:50_000 ];

julia> y = [ 250*rand(3) for i in 1:10_000 ];

julia> cl = CellList(x,y,box);

julia> aux = CellListMap.AuxThreaded(cl)
CellListMap.AuxThreaded{3, Float64}
 Auxiliary arrays for nthreads = 8

julia> cl = UpdateCellList!(x,box,cl,aux)
CellList{3, Float64}
  100000 real particles.
  31190 cells with real particles.
  1134378 particles in computing box, including images.

```
"""
AuxThreaded(cl_pair::CellListPair;particles_per_batch=10_000) =
    AuxThreaded(cl_pair.target,particles_per_batch=particles_per_batch)

"""

```
CellList(
    x::AbstractVector{AbstractVector},
    box::Box{UnitCellType,N,T};
    parallel::Bool=true,
    nbatches::Tuple{Int,Int}=(0,0)
) where {UnitCellType,N,T}
```

Function that will initialize a `CellList` structure from scracth, given a vector
or particle coordinates (a vector of vectors, typically of static vectors)
and a `Box`, which contain the size ofthe system, cutoff, etc. Except for small
systems, the number of parallel batches is equal to the number of threads, but it can
be tunned for optimal performance in some cases.

## Example

```julia-repl
julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:100000 ];

julia> cl = CellList(x,box)
CellList{3, Float64}
  100000 real particles.
  15600 cells with real particles.
  126276 particles in computing box, including images.

```

"""
function CellList(
    x::AbstractVector{<:AbstractVector},
    box::Box{UnitCellType,N,T};
    parallel::Bool=true,
    nbatches::Tuple{Int,Int} = (0,0)
) where {UnitCellType,N,T}
    cl = CellList{N,T}(
        n_real_particles=length(x),
        number_of_cells=prod(box.nc),
    )
    cl = set_number_of_batches!(cl,nbatches)
    return UpdateCellList!(x,box,cl,parallel=parallel)
end

"""

```
function CellList(x::AbstractMatrix, box::Box{UnitCellType,N,T}; kargs...) where {UnitCellType,N,T}
```

Reinterprets the matrix `x` as vectors of static vectors and calls the
equivalent function with the reinterprted input. The first dimension of the
matrix must be the dimension of the points (`2` or `3`).

"""
function CellList(x::AbstractMatrix, box::Box{UnitCellType,N,T}; kargs...) where {UnitCellType,N,T}
    @assert size(x,1) == N "First dimension of input matrix must be $N"
    x_re = reinterpret(reshape, SVector{N,eltype(x)}, x)
    return CellList(x_re, box; kargs...)
end

"""

```
reset!(cl::CellList{N,T},box,n_real_particles) where{N,T}
```

Internal function or structure - interface may change.

# Extended help

Resets a cell list, by setting everything to zero, but retaining
the allocated `particles` and `projected_particles` vectors.

"""
function reset!(cl::CellList{N,T},box,n_real_particles) where{N,T}
    new_number_of_cells = prod(box.nc)
    if new_number_of_cells > cl.number_of_cells
        resize!(cl.cell_indices,new_number_of_cells)
    end
    for i in 1:length(cl.cells)
        cl.cells[i] = Cell{N,T}(particles=cl.cells[i].particles)
    end
    @. cl.cell_indices = 0
    @. cl.cell_indices_real = 0
    cl = CellList{N,T}(
        n_real_particles = n_real_particles,
        n_particles = 0,
        number_of_cells = new_number_of_cells,
        n_cells_with_real_particles = 0,
        n_cells_with_particles = 0,
        cell_indices = cl.cell_indices,
        cell_indices_real = cl.cell_indices_real,
        cells=cl.cells,
        nbatches=cl.nbatches,
        projected_particles = cl.projected_particles
    )
    return cl
end

"""

```
CellList(
    x::AbstractVector{<:AbstractVector},
    y::AbstractVector{<:AbstractVector},
    box::Box{UnitCellType,N,T};
    parallel::Bool=true,
    nbatches::Tuple{Int,Int}=(0,0),
    autoswap::Bool=true
) where {UnitCellType,N,T}
```

Function that will initialize a `CellListPair` structure from scracth, given two vectors
of particle coordinates and a `Box`, which contain the size of the system, cutoff, etc.
By default, the cell list will be constructed for smallest vector, but this is not always
the optimal choice. Using `autoswap=false` the cell list is constructed for the second (`y`)

## Example

```julia-repl
julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:1000 ];

julia> y = [ 250*rand(SVector{3,Float64}) for i in 1:10000 ];

julia> cl = CellList(x,y,box)
CellListMap.CellListPair{Vector{SVector{3, Float64}}, 3, Float64}
   10000 particles in the reference vector.
   961 cells with real particles of target vector.

julia> cl = CellList(x,y,box,autoswap=false)
CellListMap.CellListPair{Vector{SVector{3, Float64}}, 3, Float64}
   1000 particles in the reference vector.
   7389 cells with real particles of target vector.

```

"""
function CellList(
    x::AbstractVector{<:AbstractVector},
    y::AbstractVector{<:AbstractVector},
    box::Box{UnitCellType,N,T};
    parallel::Bool=true,
    nbatches::Tuple{Int,Int} = (0,0),
    autoswap=true
) where {UnitCellType,N,T}
    if !autoswap || length(x) >= length(y)
        ref = [ SVector{N,T}(ntuple(i->el[i],N)...) for el in x ]
        target = CellList(y,box,parallel=parallel)
        swap = false
    else
        ref = [ SVector{N,T}(ntuple(i->el[i],N)...) for el in y ]
        target = CellList(x,box,parallel=parallel)
        swap = true
    end
    cl_pair = CellListPair(ref=ref,target=target,swap=swap)
    cl_pair = set_number_of_batches!(cl_pair,nbatches)
    return cl_pair
end

"""

```
CellList(x::AbstractMatrix, y::AbstractMatrix, box::Box{UnitCellType,N,T}; kargs...) where {UnitCellType,N,T}
```

Reinterprets the matrices `x` and `y` as vectors of static vectors and calls the
equivalent function with the reinterprted input. The first dimension of the
matrices must be the dimension of the points (`2` or `3`).

"""
function CellList(x::AbstractMatrix, y::AbstractMatrix, box::Box{UnitCellType,N,T}; kargs...) where {UnitCellType,N,T}
    @assert size(x,1) == N "First dimension of input matrix must be $N"
    @assert size(y,1) == N "First dimension of input matrix must be $N"
    x_re = reinterpret(reshape, SVector{N,eltype(x)}, x)
    y_re = reinterpret(reshape, SVector{N,eltype(y)}, y)
    CellList(x_re,y_re, box; kargs...)
end

"""

```
UpdateCellList!(
    x::AbstractVector{<:AbstractVector},
    box::Box,
    cl:CellList,
    parallel=true
)
```

Function that will update a previously allocated `CellList` structure, given new
updated particle positions. This function will allocate new threaded auxiliary
arrays in parallel calculations. To preallocate these auxiliary arrays, use
the `UpdateCellList!(x,box,cl,aux)` method instead.

## Example

```julia-repl
julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:1000 ];

julia> cl = CellList(x,box);

julia> box = Box([260,260,260],10);

julia> x = [ 260*rand(SVector{3,Float64}) for i in 1:1000 ];

julia> cl = UpdateCellList!(x,box,cl); # update lists

```

"""
function UpdateCellList!(
    x::AbstractVector{<:AbstractVector},
    box::Box,
    cl::CellList;
    parallel::Bool=true
)
    aux = AuxThreaded(cl)
    return UpdateCellList!(x,box,cl,aux,parallel=parallel)
end

"""

```
function UpdateCellList!(
    x::AbstractMatrix,
    box::Box,
    cl::CellList{N,T};
    parallel::Bool=true
) where {N,T}
```

Reinterprets the matrix `x` as vectors of static vectors and calls the
equivalent function with the reinterprted input. The first dimension of the
matrix must be the dimension of the points (`2` or `3`).

"""
function UpdateCellList!(
    x::AbstractMatrix,
    box::Box,
    cl::CellList{N,T};
    parallel::Bool=true
) where {N,T}
    @assert size(x,1) == N "First dimension of input matrix must be $N"
    x_re = reinterpret(reshape, SVector{N,eltype(x)}, x)
    return UpdateCellList!(x_re,box,cl,parallel=parallel)
end

"""

```
function UpdateCellList!(
    x::AbstractVector{<:AbstractVector},
    box::Box,
    cl::CellList{N,T},
    aux::AuxThreaded{N,T};
    parallel::Bool=true
) where {N,T}
```

Function that updates the cell list `cl` new coordinates `x` and possibly a new
box `box`, and receives a preallocated `aux` structure of auxiliary vectors for
threaded cell list construction. Given a preallocated `aux` vector, allocations in
this function should be minimal, only associated with the spawning threads, or
to expansion of the cell lists if the number of cells or number of particles
increased.

## Example

```julia-repl
julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:100000 ];

julia> cl = CellList(x,box);

julia> aux = CellListMap.AuxThreaded(cl)
CellListMap.AuxThreaded{3, Float64}
 Auxiliary arrays for nthreads = 8

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:100000 ];

julia> UpdateCellList!(x,box,cl,aux)
CellList{3, Float64}
  100000 real particles.
  15599 cells with real particles.
  125699 particles in computing box, including images.

```

To illustrate the expected ammount of allocations, which are a consequence
of thread spawning only:

```julia-repl
julia> using BenchmarkTools

julia> @btime UpdateCellList!(\$x,\$box,\$cl,\$aux)
  16.384 ms (41 allocations: 3.88 KiB)
CellList{3, Float64}
  100000 real particles.
  15599 cells with real particles.
  125699 particles in computing box, including images.

julia> @btime UpdateCellList!(\$x,\$box,\$cl,\$aux,parallel=false)
  20.882 ms (0 allocations: 0 bytes)
CellList{3, Float64}
  100000 real particles.
  15603 cells with real particles.
  125896 particles in computing box, including images.

```

"""
function UpdateCellList!(
    x::AbstractVector{<:AbstractVector},
    box::Box,
    cl::CellList{N,T},
    aux::AuxThreaded{N,T};
    parallel::Bool=true
) where {N,T}

    # Add particles to cell list
    nbatches = cl.nbatches.build_cell_lists
    if !parallel || nbatches == 1
        cl = reset!(cl,box,length(x))
        cl = add_particles!(x,box,0,cl)
    else
        # Reset cell list
        cl = reset!(cl,box,0)
        @sync for ibatch in 1:nbatches
            Threads.@spawn begin
                prange = aux.idxs[ibatch]
                aux.lists[ibatch] = reset!(aux.lists[ibatch],box,length(prange))
                xt = @view(x[prange])
                aux.lists[ibatch] = add_particles!(xt,box,prange[begin]-1,aux.lists[ibatch])
            end
        end
        for ibatch in 1:nbatches
            cl = merge_cell_lists!(cl,aux.lists[ibatch])
        end
    end

    # allocate, or update the auxiliary projected_particles arrays
    maxnp = 0
    for i in 1:cl.n_cells_with_particles
        maxnp = max(maxnp,cl.cells[i].n_particles)
    end
    for i in 1:length(cl.projected_particles)
        if maxnp > length(cl.projected_particles[i])
            resize!(cl.projected_particles[i],maxnp)
        end
    end

    return cl
end

"""

```
function UpdateCellList!(
    x::AbstractMatrix,
    box::Box,
    cl::CellList{N,T},
    aux::AuxThreaded{N,T};
    parallel::Bool=true
) where {N,T}
```

Reinterprets the matrix `x` as vectors of static vectors and calls the
equivalent function with the reinterprted input. The first dimension of the
matrix must be the dimension of the points (`2` or `3`).

"""
function UpdateCellList!(
    x::AbstractMatrix,
    box::Box,
    cl::CellList{N,T},
    aux::AuxThreaded{N,T};
    parallel::Bool=true
) where {N,T}
    @assert size(x,1) == N "First dimension of input matrix must be $N"
    x_re = reinterpret(reshape, SVector{N,eltype(x)}, x)
    return UpdateCellList!(x_re,box,cl,aux,parallel=parallel)
end

"""

```
add_particles!(x,box,ishift,cl::CellList{N,T}) where {N,T}
```

Internal function or structure - interface may change.

# Extended help

Add all particles in vector `x` to the cell list `cl`. `ishift` is the shift in particle
index, meaning that particle `i` of vector `x` corresponds to the particle with original
index `i+ishift`. The shift is used to construct cell lists from fractions of the original
set of particles in parallel list construction.

"""
function add_particles!(x,box,ishift,cl::CellList{N,T}) where {N,T}
    for ip in eachindex(x)
        xp = x[ip]
        # This converts the coordinates to static arrays, if necessary, and possibly
        # removes decorations (like units) from the quantities
        p = SVector{N,T}(ntuple(i->xp[i],N))
        p = wrap_to_first(p,box)
        cl = add_particle_to_celllist!(ishift+ip,p,box,cl) # add real particle
        cl = replicate_particle!(ishift+ip,p,box,cl) # add virtual particles to border cells
    end
    return cl
end

"""

```
copydata!(cell1::Cell,cell2::Cell)
```

Internal function or structure - interface may change.

# Extended help

Copies the data from `cell2` to `cell1`, meaning that particles are
copied element-wise from `cell2` to `cell1`, with the `particles` array
of `cell1` being resized (increased) if necessary.

"""
function copydata!(cell1::Cell,cell2::Cell)
    @set! cell1.linear_index = cell2.linear_index
    @set! cell1.cartesian_index = cell2.cartesian_index
    @set! cell1.center = cell2.center
    @set! cell1.n_particles = cell2.n_particles
    @set! cell1.contains_real = cell2.contains_real
    for ip in 1:cell2.n_particles
        p = cell2.particles[ip]
        if ip > length(cell1.particles)
            push!(cell1.particles,p)
        else
            cell1.particles[ip] = p
        end
    end
    return cell1
end

"""
```
append_particles!(cell1::Cell,cell2::Cell)
```

Internal function or structure - interface may change.

# Extended help

Add the particles of `cell2` to `cell1`, updating the cell data and, if necessary,
resizing (increasing) the `particles` array of `cell1`

"""
function append_particles!(cell1::Cell,cell2::Cell)
    if cell2.contains_real
        @set! cell1.contains_real = true
    end
    n_particles_old = cell1.n_particles
    @set! cell1.n_particles += cell2.n_particles
    if cell1.n_particles > length(cell1.particles)
        resize!(cell1.particles,cell1.n_particles)
    end
    for ip in 1:cell2.n_particles
        cell1.particles[n_particles_old+ip] = cell2.particles[ip]
    end
    return cell1
end

"""

```
merge_cell_lists!(cl::CellList,aux::CellList)
```

Internal function or structure - interface may change.

# Extended help

Merges an auxiliary `aux` cell list to `cl`, and returns the modified `cl`. Used to
merge cell lists computed in parallel threads.

"""
function merge_cell_lists!(cl::CellList,aux::CellList)
    # One should never get here if the lists do not share the same # computing box
    if cl.number_of_cells != aux.number_of_cells
        error("cell lists must have the same number of cells to be merged.")
    end
    @set! cl.n_particles += aux.n_particles
    @set! cl.n_real_particles += aux.n_real_particles
    for icell in 1:aux.n_cells_with_particles
        aux_cell = aux.cells[icell]
        linear_index = aux_cell.linear_index
        cell_index = cl.cell_indices[linear_index]
        # If cell was yet not initialized in merge, push it to the list
        if cell_index == 0
            @set! cl.n_cells_with_particles += 1
            cell_index = cl.n_cells_with_particles
            cl.cell_indices[linear_index] = cell_index
            if cell_index > length(cl.cells)
                push!(cl.cells,deepcopy(aux_cell))
            else
                cl.cells[cell_index] = copydata!(cl.cells[cell_index],aux_cell)
            end
            if aux_cell.contains_real
                @set! cl.n_cells_with_real_particles += 1
                if cl.n_cells_with_real_particles > length(cl.cell_indices_real)
                    push!(cl.cell_indices_real,cell_index)
                else
                    cl.cell_indices_real[cl.n_cells_with_real_particles] = cell_index
                end
            end
        # Append particles to initialized cells
        else
            # If the previous cell didn't contain real particles, but the current one
            # does, update the list information
            if !cl.cells[cell_index].contains_real && aux_cell.contains_real
                @set! cl.n_cells_with_real_particles += 1
                if cl.n_cells_with_real_particles > length(cl.cell_indices_real)
                    push!(cl.cell_indices_real,cell_index)
                else
                    cl.cell_indices_real[cl.n_cells_with_real_particles] = cell_index
                end
            end
            cl.cells[cell_index] = append_particles!(cl.cells[cell_index],aux_cell)
        end
    end
    return cl
end

"""

```
add_particle_to_celllist!(
    ip,
    x::SVector{N,T},
    box,
    cl::CellList{N,T};
    real_particle::Bool=true
) where {N,T}
```

Internal function or structure - interface may change.

# Extended help

Adds one particle to the cell lists, updating all necessary arrays.

"""
function add_particle_to_celllist!(
    ip,
    x::SVector{N,T},
    box,
    cl::CellList{N,T};
    real_particle::Bool=true
) where {N,T}
    @unpack n_cells_with_real_particles,
            n_cells_with_particles,
            cell_indices,
            cell_indices_real,
            cells = cl

    # Increase the counter of number of particles of this list
    @set! cl.n_particles += 1
    # Cell of this particle
    cartesian_index = particle_cell(x,box)
    linear_index = cell_linear_index(box.nc,cartesian_index)

    # Check if this is the first particle of this cell, if it is,
    # initialize a new cell, or reset a previously allocated one
    cell_index = cell_indices[linear_index]

    if cell_index == 0
        n_cells_with_particles += 1
        cell_index = n_cells_with_particles
        cell_indices[linear_index] = cell_index
        if cell_index > length(cells)
            particles_sizehint = cl.n_real_particles ÷ prod(box.nc)
            push!(cells,Cell{N,T}(cartesian_index,box,sizehint=particles_sizehint))
        else
            cell = cells[cell_index]
            @set! cell.linear_index = linear_index
            @set! cell.cartesian_index = cartesian_index
            @set! cell.center = cell_center(cell.cartesian_index,box)
            @set! cell.contains_real = false
            @set! cell.n_particles = 0
            cells[cell_index] = cell
        end
    end

    # Increase particle counter for this cell
    cell = cells[cell_index]
    @set! cell.n_particles += 1

    #
    # Cells with real particles are annotated to be run over
    #
    if real_particle && (!cell.contains_real)
        @set! cell.contains_real = true
        n_cells_with_real_particles += 1
        if n_cells_with_real_particles > length(cell_indices_real)
            push!(cell_indices_real,cell_index)
        else
            cell_indices_real[n_cells_with_real_particles] = cell_index
        end
    end

    #
    # Add particle to cell list
    #
    p = ParticleWithIndex(ip,real_particle,x)
    if cell.n_particles > length(cell.particles)
        push!(cell.particles,p)
    else
        cell.particles[cell.n_particles] = p
    end

    #
    # Update (imutable) cell in list
    #
    @set! cl.n_cells_with_particles = n_cells_with_particles
    @set! cl.n_cells_with_real_particles = n_cells_with_real_particles
    cells[cell_index] = cell

    return cl
end

"""

```
UpdateCellList!(
  x::AbstractVector{<:AbstractVector},
  y::AbstractVector{<:AbstractVector},
  box::Box,
  cl:CellListPair,
  parallel=true
)
```

Function that will update a previously allocated `CellListPair` structure, given
new updated particle positions, for example. This method will allocate new
`aux` threaded auxiliary arrays. For a non-allocating version, see the
`UpdateCellList!(x,y,box,cl,aux)` method.

```julia-repl
julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:1000 ];

julia> y = [ 250*rand(SVector{3,Float64}) for i in 1:10000 ];

julia> cl = CellList(x,y,box);

julia> cl = UpdateCellList!(x,y,box,cl); # update lists

```

"""
function UpdateCellList!(
    x::AbstractVector{<:AbstractVector},
    y::AbstractVector{<:AbstractVector},
    box::Box,
    cl_pair::CellListPair;
    parallel::Bool=true
)
    aux = AuxThreaded(cl)
    return UpdateCellList!(x,y,box,cl_pair,aux,parallel=parallel)
end

"""

```
function UpdateCellList!(
    x::AbstractMatrix,
    y::AbstractMatrix,
    box::Box{UnitCellType,N,T},
    cl_pair::CellListPair;
    parallel::Bool=true
) where {UnitCellType,N,T}
```

Reinterprets the matrices `x` and `y` as vectors of static vectors and calls the
equivalent function with the reinterprted input. The first dimension of the
matrices must be the dimension of the points (`2` or `3`).

"""
function UpdateCellList!(
    x::AbstractMatrix,
    y::AbstractMatrix,
    box::Box{UnitCellType,N,T},
    cl_pair::CellListPair;
    parallel::Bool=true
) where {UnitCellType,N,T}
    @assert size(x,1) == N "First dimension of input matrix must be $N"
    @assert size(y,1) == N "First dimension of input matrix must be $N"
    x_re = reinterpret(reshape, SVector{N,eltype(x)}, x)
    y_re = reinterpret(reshape, SVector{N,eltype(y)}, y)
    return UpdateCellList!(x_re,y_re,box,cl_pair,parallel=parallel)
end

"""

```
function UpdateCellList!(
    x::AbstractVector{<:AbstractVector},
    y::AbstractVector{<:AbstractVector},
    box::Box,
    cl_pair::CellListPair,
    aux::AuxThreaded;
    parallel::Bool=true
)
```

This function will update the `cl_pair` structure that contains the cell lists
for disjoint sets of particles. It receives the preallocated `aux` structure to
avoid reallocating auxiliary arrays necessary for the threaded construct of the
lists.

## Example

```julia-repl
julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(3) for i in 1:50_000 ];

julia> y = [ 250*rand(3) for i in 1:10_000 ];

julia> cl = CellList(x,y,box)
CellListMap.CellListPair{Vector{SVector{3, Float64}}, 3, Float64}
   50000 particles in the reference vector.
   7381 cells with real particles of target vector.

julia> aux = CellListMap.AuxThreaded(cl)
CellListMap.AuxThreaded{3, Float64}
 Auxiliary arrays for nthreads = 8

julia> x = [ 250*rand(3) for i in 1:50_000 ];

julia> y = [ 250*rand(3) for i in 1:10_000 ];

julia> cl = UpdateCellList!(x,y,box,cl,aux)
CellList{3, Float64}
  10000 real particles.
  7358 cells with real particles.
  12591 particles in computing box, including images.

```
To illustrate the expected ammount of allocations, which are a consequence
of thread spawning only:

```julia-repl
julia> using BenchmarkTools

julia> @btime UpdateCellList!(\$x,\$y,\$box,\$cl,\$aux)
  715.661 μs (41 allocations: 3.88 KiB)
CellListMap.CellListPair{Vector{SVector{3, Float64}}, 3, Float64}
   50000 particles in the reference vector.
   7414 cells with real particles of target vector.

julia> @btime UpdateCellList!(\$x,\$y,\$box,\$cl,\$aux,parallel=false)
   13.042 ms (0 allocations: 0 bytes)
 CellListMap.CellListPair{Vector{SVector{3, Float64}}, 3, Float64}
    50000 particles in the reference vector.
    15031 cells with real particles of target vector.

```

"""
function UpdateCellList!(
    x::AbstractVector{<:AbstractVector},
    y::AbstractVector{<:AbstractVector},
    box::Box,
    cl_pair::CellListPair,
    aux::AuxThreaded;
    parallel::Bool=true
)
    if !cl_pair.swap
        target = UpdateCellList!(x,box,cl_pair.target,aux,parallel=parallel)
    else
        target = UpdateCellList!(y,box,cl_pair.target,aux,parallel=parallel)
    end
    cl_pair = CellListPair(ref=cl_pair.ref,target=target,swap=cl_pair.swap)
    return cl_pair
end

"""

```
function UpdateCellList!(
    x::AbstractMatrix,
    y::AbstractMatrix,
    box::Box,
    cl_pair::CellListPair,
    aux::AuxThreaded;
    parallel::Bool=true
) where {UnitCellType,N,T}
```

Reinterprets the matrices `x` and `y` as vectors of static vectors and calls the
equivalent function with the reinterprted input. The first dimension of the
matrices must be the dimension of the points (`2` or `3`).

"""
function UpdateCellList!(
    x::AbstractMatrix,
    y::AbstractMatrix,
    box::Box{UnitCellType,N,T},
    cl_pair::CellListPair{N,T},
    aux::AuxThreaded{N,T};
    parallel::Bool=true
) where {UnitCellType,N,T}
    @assert size(x,1) == N "First dimension of input matrix must be $N"
    @assert size(y,1) == N "First dimension of input matrix must be $N"
    x_re = reinterpret(reshape, SVector{N,eltype(x)}, x)
    y_re = reinterpret(reshape, SVector{N,eltype(y)}, y)
    return UpdateCellList!(x_re,y_re,box,cl_pair,aux,parallel=parallel)
end

"""

```
particles_per_cell(cl)
```

Returns the average number of real particles per computing cell.

"""
particles_per_cell(cl::CellList) = cl.n_real_particles / cl.number_of_cells
particles_per_cell(cl::CellListPair) = particles_per_cell(cl.target)


text:
module Bonsai

using HTTP, JSON3, StructTypes, Dates, LRUCache, FilePaths
import HTTP: Request, Response, Stream

export ResponseCodes

include("ResponseCodes.jl")

include("mime_types.jl")
include("http.jl")
include("exceptions.jl")
include("json_schema.jl")
include("handlers.jl")
include("io.jl")
include("cancel_token.jl")
include("router.jl")
include("middleware.jl")
include("openapi.jl")
include("app.jl")
include("server.jl")
include("docs.jl")

end

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: braket
using AWS.Compat
using AWS.UUIDs

"""
    CancelQuantumTask()

Cancels the specified task.

# Required Parameters
- `clientToken`: The client token associated with the request.
- `quantumTaskArn`: The ARN of the task to cancel.

"""
cancel_quantum_task(clientToken, quantumTaskArn; aws_config::AWSConfig=global_aws_config()) = braket("PUT", "/quantum-task/$(quantumTaskArn)/cancel", Dict{String, Any}("clientToken"=>clientToken); aws_config=aws_config)
cancel_quantum_task(clientToken, quantumTaskArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("PUT", "/quantum-task/$(quantumTaskArn)/cancel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientToken"=>clientToken), args)); aws_config=aws_config)

"""
    CreateQuantumTask()

Creates a quantum task.

# Required Parameters
- `action`: The action associated with the task.
- `clientToken`: The client token associated with the request.
- `deviceArn`: The ARN of the device to run the task on.
- `outputS3Bucket`: The S3 bucket to store task result files in.
- `outputS3KeyPrefix`: The key prefix for the location in the S3 bucket to store task results in.
- `shots`: The number of shots to use for the task.

# Optional Parameters
- `deviceParameters`: The parameters for the device to run the task on.
- `tags`: Tags to be added to the quantum task you're creating.
"""
create_quantum_task(action, clientToken, deviceArn, outputS3Bucket, outputS3KeyPrefix, shots; aws_config::AWSConfig=global_aws_config()) = braket("POST", "/quantum-task", Dict{String, Any}("action"=>action, "clientToken"=>clientToken, "deviceArn"=>deviceArn, "outputS3Bucket"=>outputS3Bucket, "outputS3KeyPrefix"=>outputS3KeyPrefix, "shots"=>shots); aws_config=aws_config)
create_quantum_task(action, clientToken, deviceArn, outputS3Bucket, outputS3KeyPrefix, shots, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("POST", "/quantum-task", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("action"=>action, "clientToken"=>clientToken, "deviceArn"=>deviceArn, "outputS3Bucket"=>outputS3Bucket, "outputS3KeyPrefix"=>outputS3KeyPrefix, "shots"=>shots), args)); aws_config=aws_config)

"""
    GetDevice()

Retrieves the devices available in Amazon Braket.

# Required Parameters
- `deviceArn`: The ARN of the device to retrieve.

"""
get_device(deviceArn; aws_config::AWSConfig=global_aws_config()) = braket("GET", "/device/$(deviceArn)"; aws_config=aws_config)
get_device(deviceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("GET", "/device/$(deviceArn)", args; aws_config=aws_config)

"""
    GetQuantumTask()

Retrieves the specified quantum task.

# Required Parameters
- `quantumTaskArn`: the ARN of the task to retrieve.

"""
get_quantum_task(quantumTaskArn; aws_config::AWSConfig=global_aws_config()) = braket("GET", "/quantum-task/$(quantumTaskArn)"; aws_config=aws_config)
get_quantum_task(quantumTaskArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("GET", "/quantum-task/$(quantumTaskArn)", args; aws_config=aws_config)

"""
    ListTagsForResource()

Shows the tags associated with this resource.

# Required Parameters
- `resourceArn`: Specify the resourceArn for the resource whose tags to display.

"""
list_tags_for_resource(resourceArn; aws_config::AWSConfig=global_aws_config()) = braket("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("GET", "/tags/$(resourceArn)", args; aws_config=aws_config)

"""
    SearchDevices()

Searches for devices using the specified filters.

# Required Parameters
- `filters`: The filter values to use to search for a device.

# Optional Parameters
- `maxResults`: The maximum number of results to return in the response.
- `nextToken`: A token used for pagination of results returned in the response. Use the token returned from the previous request continue results where the previous request ended.
"""
search_devices(filters; aws_config::AWSConfig=global_aws_config()) = braket("POST", "/devices", Dict{String, Any}("filters"=>filters); aws_config=aws_config)
search_devices(filters, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("POST", "/devices", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("filters"=>filters), args)); aws_config=aws_config)

"""
    SearchQuantumTasks()

Searches for tasks that match the specified filter values.

# Required Parameters
- `filters`: Array of SearchQuantumTasksFilter objects.

# Optional Parameters
- `maxResults`: Maximum number of results to return in the response.
- `nextToken`: A token used for pagination of results returned in the response. Use the token returned from the previous request continue results where the previous request ended.
"""
search_quantum_tasks(filters; aws_config::AWSConfig=global_aws_config()) = braket("POST", "/quantum-tasks", Dict{String, Any}("filters"=>filters); aws_config=aws_config)
search_quantum_tasks(filters, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("POST", "/quantum-tasks", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("filters"=>filters), args)); aws_config=aws_config)

"""
    TagResource()

Add a tag to the specified resource.

# Required Parameters
- `resourceArn`: Specify the resourceArn of the resource to which a tag will be added.
- `tags`: Specify the tags to add to the resource.

"""
tag_resource(resourceArn, tags; aws_config::AWSConfig=global_aws_config()) = braket("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

Remove tags from a resource.

# Required Parameters
- `resourceArn`: Specify the resourceArn for the resource from which to remove the tags.
- `tagKeys`: pecify the keys for the tags to remove from the resource.

"""
untag_resource(resourceArn, tagKeys; aws_config::AWSConfig=global_aws_config()) = braket("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)


text:
"""
    circ_xcorr2(u,v; norm=true, center=true, phat=false)

Perform circular cross-correlation using the FFT algorithm, shifting `v`
relative to `u`, so delaying `u` will shift the result to the right. Index 1
corresponds to zero lag.

Another interpretation is that index of the peak of this cross correlation
(minus 1) represents the delay of `u` with respect to `v`, or alternatively the
position of `v` within `u`.

If `phat=true`, this implements GCC-PHAT, where the amplitudes are set to one
and only the phases are used.
"""
function circ_xcorr2(u,v; norm=true, center=true, phat=false)
    L = length(u)
    @assert length(v) == L
    if center
        u .-= mean(u)
        v .-= mean(v)
    end
    p = plan_rfft(u)
    U = p*u
    V = p*v
    if phat
        @. U = U * conj(V) / (abs(U) * abs(V))
    else
        @. U *= conj(V)
    end
    result = irfft(U, L)

    if norm
        result ./= sqrt(sum(x->x^2, u)*sum(x->x^2, v))
    end

    result
end

"""
    xcorr2(u,v;
           lagbounds=(-length(v)+1, length(u)-1),
           unbiased=true,
           norm=true,
           center=true,
           phat=false,
           plantype=FFTW.ESTIMATE)

Perform linear cross-correlation using the FFT algorithm, shifting `v` relative
to `u`, so delaying `u` will shift the result to the right. Index
`-lagbounds[1]+1` (default `length(v)`) corresponds to zero lag.

If `unbiased` is true (the default) then this will compensate for lost energy
due to lags where the signals don't entirely overlap. This is generally appropriate
when the target signals being correlated occupy a substantial fraction of `u` and `v`.
It also helps keep the noise power constant for the uncorrelated lags. Note that there
may be situations where this is not appropriate - more investigation needs to be done.

If `phat=true`, this implements GCC-PHAT, where the amplitudes are set to one
and only the phases are used.

`plantype` can be set to one of:
    - `FFTW.ESTIMATE` (the default)
    - `FFTW.MEASURE`
    - `FFTW.PATIENT`
    - `FFTW.EXHAUSTIVE`

See the FFTW documentation for more information on these settings. The plan information
is saved internally within FFTW on each session, so you only pay the cost of planning
the first time `xcorr2` is called.

!!! WARNING: Currently (Feb 2020) I get very long pauses and sometimes hangs with any
plan other than `FFTW.ESTIMATE`
"""
function xcorr2(u,v;
                lagbounds=(-length(v)+1, length(u)-1),
                unbiased=true,
                norm=true,
                center=true,
                phat=false,
                plantype=FFTW.ESTIMATE)
    # TODO make sure it works for zero-length vectors
    # TODO: handle minlag == maxlag
    # TODO: really `unbiased` is a form of normalization, where we normalize
    # each lag by the number of positions that overlap. This also suggests yet
    # another normalization approach where we normalize each lag by geometric
    # mean (similar to the current `norm=true` but we only include the parts
    # that overlap)
    su = length(u)
    sv = length(v)
    minlag, maxlag = lagbounds
    minlag <= maxlag || throw(ArgumentError(
        "maximum lag must be greater or equal to minimum lag"))
    # if `u` and `v` are different lengths and the lags don't cover the whole
    # signal then we might not need all of the signals. E.g. if `u` is longer
    # than `v` by more than `minlag` then the end of the `u` is never used,
    # so we can chop it off.
    su, sv = if su > sv
        clamp(sv + maxlag, sv, su), sv
    elseif su < sv
        su, clamp(su - minlag, su, sv)
    else
        su, sv
    end
    nlags = maxlag - minlag + 1
    nfft = nextfastfft(max(su - minlag, sv + maxlag))
    upad = zeros(float(eltype(u)), nfft)
    vpad = zeros(float(eltype(v)), nfft)

    copyto!(upad, 1, u, 1, su)
    # pre-shift `v` to correspond to the `minlag` position
    if minlag < 0
        copyto!(vpad, 1, v, -minlag+1, sv+minlag)
        copyto!(vpad, nfft+minlag+1, v, 1, -minlag)
    else
        copyto!(vpad, minlag+1, v, 1, sv)
    end

    if center
        umean = mean(u)
        upad[-minlag+1:sv] .-= umean
        upad[1:-minlag] .-= umean
        vpad[1:sv] .-= mean(v)
    end
    # we don't normalize if we're using GCC-PHAT because we set all the amplitudes
    # to 1
    normval = if norm && !phat
        sqrt(sum(abs2, upad)*sum(abs2, vpad))
    else
        one(eltype(upad))
    end

    p = plan_rfft(upad; flags=plantype | FFTW.DESTROY_INPUT)
    uspec = p * upad
    vspec = p * vpad
    magthresh = sqrt(eps())
    for i in eachindex(uspec)
        @inbounds uspec[i] *= conj(vspec[i])
        if phat
            @inbounds mag = abs(uspec[i])
            @inbounds mag > magthresh && (uspec[i] /= mag)
        end
    end

    result = resize!(inv(p) * uspec, nlags)

    if unbiased
        # maximum overlap when lag == 0
        maxoverlap = min(su, sv)
        for (i, lag) in enumerate(UnitRange(lagbounds...))
            overlap = if lag < 0
                max(0, maxoverlap+lag)
            elseif lag < su - sv
                maxoverlap
            else
                max(0, maxoverlap - lag + su - sv)
            end
            if overlap > 0
                # this power is empirical - I haven't throught through whether
                # there's a theoretical reason, or whether it is signal-dependent.
                # this is probably another good reason to normalize each lag based
                # on the actual energy in the overlap instead of this estimate
                result[i] *= (maxoverlap / overlap)^0.5
            end
        end
    end

    if norm
        result ./= normval
    end
    result
end

"""
    xc_align(x, y, lagbounds; phat=true)

Returns `(vx, vy, os)` with views into `x` and `y` that are aligned to maximize the
cross-correlation, as well as delay (in samples) of `y` relative to `x`. The views will
be of equal length.

`lagbounds` should be a `Tuple` with the minimum and maximum shift of `y` relative to
`x`.
`phat` enables the GCC-PHAT algorithm for sharper cross-correlations
"""
function xc_align(x, y, lagbounds; phat=true)
    xc = xcorr2(y, x; phat=phat, lagbounds=lagbounds)
    offset = findmax(abs2.(xc))[2] + lagbounds[1] - 1
    vx, vy = if offset >= 0
        L = min(length(x), length(y)-offset)
        @views x[1:L], y[(1:L).+offset]
    else
        L = min(length(x)+offset, length(y))
        @views x[(1:L).-offset], y[1:L]
    end

    vx, vy, offset
end

max_stars_repo_name: ssfrr/SFRDSP

max_issues_repo_name: ssfrr/SFRDSP

max_forks_repo_name: ssfrr/SFRDSP


text:
using QuDiffEq, Yao
using Test, Random, LinearAlgebra

@testset "HHL_tests" begin
    include("HHL_tests.jl")
end

@testset "PhaseEstimation_tests" begin
   include("PhaseEstimation_tests.jl")
end

@testset "QuLDE_test" begin
    include("QuLDE_tests.jl")
end

@testset "QuNLDE_test" begin
    include("QuNLDE_tests.jl")
end

@testset "QuDiffHHL_tests" begin
    include("QuDiffHHL_tests.jl")
end

@testset "TaylorTrunc_tests" begin
    include("TaylorTrunc_tests.jl")
end


text:
# util.jl
# Utility code for QuickAPI.jl

"""
@trycatch expr

Convenience macro which wraps an expression in a `try/catch` block.
If the expression throws an error, it will be printed with the
`@error` macro.
# Examples
julia> @trycatch 1 + 1
2
julia> @trycatch 1 + "dog"
┌ Error: MethodError(+, (1, "dog"), 0x0000000000006cbc)
└ @ Main REPL[13]:6
"""
macro trycatch(expr)
    quote
        try
            $expr
        catch err
            @error err
        end
    end
end

function iscallable(s::Symbol)
    Bool(length(methods(eval(s))))
end

macro whatisit(e)
    if typeof(e) == Symbol
        if iscallable(e)
            println("$e looks like an already-defined function to me!")
            return
        else
            println("$e isn't a function!")
        end
    else if typeof(e) == Expr
        if
    else
        println("$e isn't a Symbol or an Expr! It appears to be a $(typeof($e))")

    end
end

macro getargmap(e::Expr)
    @assert e.head == :->
    args = [a.args for a in e.args[1].args]
    argmap = Dict([(a => eval(b)) for (a,b) in args])
    quote
        $([ "$k::String" for k in keys(argmap) ])
    end
end

macro fillintheblanks(args::Dict{String,Any},fn::Function)
    fmt_args = join(["$a=$v" for (a,v) = $args]
    quote quote
        $fn($(join(
            "$a=$v" for (a,v) = $args
        )))
    end end
end

function fn(a::Int = 1, b::String = "2")
    @show (a, b)
end

args = Dict("a"=>1,"b"=>"test")

@fillintheblanks args fn

macro passiton(fn,args)
    quote
        fn($args...)
    end
end