text:
"""
R_multipole = 1.1 is the constant ratio between scattering disks and the maximal
radius of their particles, and thus half the minimal distance between neighboring
particles. While mathematically this can be reduced to `1 + eps()`, that will
increase the necessary `P`.
"""
const R_multipole = 1.1

"""
	AbstractShapeParams

Abstract type which all shape types inherit from.
"""
abstract type AbstractShapeParams end

"""
	ShapeParams(t,ft,dft)

Returns `ShapeParams` object containing the parametrization of a two-dimensional
shape. `t` is a uniform sampling of [0,2π), `ft = [x(t) y(t)]`, and
`dft = [x'(t) y'(t)]`. The field `R` contains the radius of the shape's
scattering disk.

See also: `CircleParams`,`R_multipole`.
"""
struct ShapeParams <: AbstractShapeParams
	t::Array{Float64,1}
	ft::Array{Float64,2}
	dft::Array{Float64,2}
	R::Float64 #radius of multipole disk

	ShapeParams(t,ft,dft) = new(t,ft,dft,R_multipole*maximum(hypot.(ft[:,1],ft[:,2])))
end

"""
	CircleParams(R)

Returns object for a circular shape, containing its radius in the field `R`
(which is also the radius of the scattering disk).

See also: `ShapeParams`,`R_multipole`.
"""
struct CircleParams <: AbstractShapeParams
	R::Float64 #radius of multipole disk = radius of circle
end

#### FMM #####

mutable struct FMMgroup
    #TODO: make immutable and change divideSpace accordingly.
    point_ids::Vector{Int}
    center::Array{Float64,2}
	size::Int

	FMMgroup(point_ids,center) = new(point_ids,center,length(point_ids))
end

import Base.push!
function push!(g::FMMgroup,p)
	push!(g.point_ids,p)
	g.size = length(g.point_ids)
end

mutable struct FMMbuffer
    rhs::Vector{Complex{Float64}}
    pre_agg::Array{Complex{Float64},2}
    trans::Vector{Complex{Float64}}

    FMMbuffer(Ns::Integer,P::Integer,Q::Integer,Ngroups::Integer) =
                                new(zeros(Complex{Float64},Ns*(2*P+1)),
                                    zeros(Complex{Float64},Q,Ngroups),
                                    zeros(Complex{Float64},Q))
end

mutable struct FMMmatrices
	Agg::Array{Array{Complex{Float64},2},1}
	Trans::Array{Array{Complex{Float64},1},1}
	Disagg::Array{Array{Complex{Float64},2},1}
	Znear::SparseMatrixCSC{Complex{Float64},Int}
	groups::Vector{ParticleScattering.FMMgroup}
	P2::Int
	Q::Int
end

"""
	FMMoptions(FMM; nx = 0, dx = 0.0, acc = 0, tol = 0.0, method = "pre")

Constructor for `struct` containing all FMM options. `FMM` decides if FMM is used,
and the following keyword arguments dictate its behavior:

- `nx::Integer`: number of groups in x direction (for division)
- `dx::Real`: group height/width (alternative division)
- `acc::Integer`: accuracy digits for translation truncation, and also for gmres if `tol` is not given
- `tol::Real`: gmres tolerance
- `method::String`: method used: for now can be "pre". Mainly used for development.
"""
mutable struct FMMoptions
    FMM::Bool       #Is FMM used?
    nx::Integer     #number of groups in x direction (for division)
    dx::Real        #group height/width (alternative division)
    acc::Integer    #accuracy digits for translation truncation, and also for gmres if tol is not given
    tol::Real       #gmres tolerance
    method::String  #method used: can be "pre"
	# symmetric::Bool #are agg = disagg points, and thus Disagg[k] = Agg^*[k] ∀k?

    #empty contructor - for not using FMM
    FMMoptions() = new(false, 0, 0.0, 0, 0.0, "")

    #Full FMMoptions constructor with value checking
    function FMMoptions(FMM; nx = 0, dx = 0.0, acc = 0, tol = 0.0, method = "pre")
        FMM == false && return FMMoptions()
        (nx == 0 && dx == 0.0) &&
            error("FMMoptions: either nx or dx must be specified")
        (nx > 0 && dx > 0.0) &&
            error("FMMoptions: either nx or dx must be specified")
        nx < 0 &&
            error("FMMoptions: nx must be greater than 0")
        dx < 0.0 &&
            error("FMMoptions: dx must be greater than 0")
        !in(acc, 1:16) &&
            error("FMMoptions: accuracy digits must be in [1,16]")
        tol < 0.0 &&
            error("FMMoptions: gmres tolerance must be greater than 0.0")
        !in(method, ["pre"]) &&
            error("FMMoptions: method must be \"pre\" or ...")
        tol == 0.0 && (tol = 10^(-Float64(acc)))
        return new(true, nx, dx, acc, tol, method)
    end
end

"""
	OptimBuffer(Ns::Integer, P::Integer, Npoints::Integer, [J::Integer])

Constructor for the `OptimBuffer` type, which stores some of the buffers and
shared variables necessary for optimization. Includes the cylindrical harmonics
coefficient vector `β`, field values at points of interest (`f`), the partial
derivatives `∂β`, and storage for the various right-hand side vectors used while
solving for `∂β`.

If the number of optimization variables `J` is not supplied, it is assumed that
`J` = `Ns`.
"""
mutable struct OptimBuffer
    β::Vector{Complex{Float64}}
    f::Vector{Complex{Float64}}
    ∂β::Array{Complex{Float64},2}
    rhs_grad::Vector{Complex{Float64}}

    OptimBuffer(Ns::Integer, P::Integer, Npoints::Integer, J = Ns) =
                                new(Array{Complex{Float64}}(undef, Ns*(2*P+1)),
                                    Array{Complex{Float64}}(undef, Npoints),
                                    Array{Complex{Float64}}(undef, Ns*(2*P+1), J),
                                    Array{Complex{Float64}}(undef, Ns*(2*P+1)))
end

"""
	ScatteringProblem(shapes, ids, centers, φs)

Constructor for the `ScatteringProblem` type, including particle shape information
for multiple-scattering problems.
"""
struct ScatteringProblem
	shapes::Vector{T} where T <: AbstractShapeParams
	ids::Vector{Int}
	centers::Array{Float64,2}
	φs::Vector{Float64}

	ScatteringProblem(shapes, ids, centers, φs) =
		((length(ids) == size(centers,1)) &&
		(length(ids) == length(φs)) &&
		(maximum(ids) <= length(shapes))) ? new(shapes,ids,centers,φs) :
		error("ScatteringProblem: size mismatch")
end

import Base.size
size(q::ParticleScattering.ScatteringProblem) = length(q.ids)

abstract type Einc end

"""
	PlaneWave(θi)

Constructor for the `PlaneWave` incident field type, where `θi` is the angle between
the wavevector and the x-axis.
"""
struct PlaneWave <: Einc
	θi::Float64
	PlaneWave(θi = 0) = new(θi)
end

"""
	LineSource(x, y)

Constructor for the `LineSource` type, where `(x,y)` is the coordinate of the
current filament.
"""
struct LineSource <: Einc
	x::Float64
	y::Float64
end

struct CurrentSource <: Einc
	#p is a *straight* line for integration simplicity
	p::Array{Float64,2}
	σ::Vector{Complex{Float64}}
	len::Float64
end

"""
	CurrentSource(x1, y1, x2, y2, σ)

Constructor for the `CurrentSource` type, where `(x1,y1)` and `(x2,y2)` denote the
start and end points of the source, and `σ` contains the potential density.
"""
function CurrentSource(x1, y1, x2, y2, σ)
	len = sqrt((x2 - x1)^2 + (y2 - y1)^2)
	t = range(0, stop=1, length=length(σ))
	p = [(x2 - x1)*t + x1		(y2 - y1)*t + y1]
	CurrentSource(p, σ, len)
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.

module TestNDArray

using MXNet
using Statistics
using LinearAlgebra
using Test

using ..Main: rand_dims

################################################################################
# Test Implementations
################################################################################
rand_tensors(dims::NTuple{N,Int}) where {N} = rand_tensors(mx.MX_float, dims)
function rand_tensors(::Type{T}, dims::NTuple{N,Int}) where {N,T}
  tensor = rand(T, dims)
  array  = copy(tensor, mx.cpu())
  return (tensor, array)
end

function test_constructor()
  @info("NDArray::NDArray(x::AbstractArray)")
  function check_absarray(x)
    y = mx.NDArray(x)
    @test ndims(x)  == ndims(y)
    @test eltype(x) == eltype(y)
    @test x[3]      == y[3][]
  end

  check_absarray(1:10)
  check_absarray(1.0:10)

  @info("NDArray::NDArray(Type, AbstractArray)")
  let
    x = mx.NDArray(Float32, [1, 2, 3])
    @test eltype(x) == Float32
    @test copy(x) == [1, 2, 3]
  end
  let
    x = mx.NDArray(Float32, [1.1, 2, 3])
    @test eltype(x) == Float32
    @test copy(x) ≈ [1.1, 2, 3]
  end

  @info "NDArray::NDArray{T,N}(undef, dims...)"
  let
    x = NDArray{Int,2}(undef, 5, 5)
    @test eltype(x) == Int
    @test size(x) == (5, 5)
    @test x.writable

    y = NDArray{Int,2}(undef, 5, 5, writable = false)
    @test !y.writable

    # dimension mismatch
    @test_throws MethodError NDArray{Int,1}(undef, 5, 5)
  end

  @info "NDArray::NDArray{T,N}(undef, dims)"
  let
    x = NDArray{Int,2}(undef, (5, 5))
    @test eltype(x) == Int
    @test size(x) == (5, 5)
    @test x.writable

    y = NDArray{Int,2}(undef, (5, 5), writable = false)
    @test !y.writable

    # dimension mismatch
    @test_throws MethodError NDArray{Int,1}(undef, (5, 5))
  end

  @info "NDArray::NDArray{T}(undef, dims...)"
  let
    x = NDArray{Int}(undef, 5, 5)
    @test eltype(x) == Int
    @test size(x) == (5, 5)
    @test x.writable

    y = NDArray{Int}(undef, 5, 5, writable = false)
    @test !y.writable
  end

  @info "NDArray::NDArray{T}(undef, dims)"
  let
    x = NDArray{Int}(undef, (5, 5))
    @test eltype(x) == Int
    @test size(x) == (5, 5)
    @test x.writable

    y = NDArray{Int}(undef, (5, 5), writable = false)
    @test !y.writable
  end

  @info "NDArray::NDArray(undef, dims...)"
  let
    x = NDArray(undef, 5, 5)
    @test eltype(x) == mx.MX_float
    @test size(x) == (5, 5)
    @test x.writable

    y = NDArray(undef, 5, 5, writable = false)
    @test !y.writable
  end

  @info "NDArray::NDArray(undef, dims)"
  let
    x = NDArray(undef, (5, 5))
    @test eltype(x) == mx.MX_float
    @test size(x) == (5, 5)
    @test x.writable

    y = NDArray(undef, (5, 5), writable = false)
    @test !y.writable
  end
end  # function test_constructor

function test_ones_zeros_like()
  @info("NDArray::Base.zeros")
  let x = mx.rand(1, 3, 2, 4, low = 1, high = 10)
    y = zeros(x)
    @test sum(copy(y)) == 0

    y = mx.zeros(x)
    @test sum(copy(y)) == 0
  end

  @info("NDArray::Base.ones")
  let x = mx.rand(1, 3, 2, 4, low = 1, high = 10)
    y = ones(x)
    @test sum(copy(y)) == 1 * 3 * 2 * 4

    y = mx.ones(x)
    @test sum(copy(y)) == 1 * 3 * 2 * 4
  end
end  # function test_ones_zeros_like

function test_copy()
  dims    = rand_dims()
  tensor  = rand(mx.MX_float, dims)

  @info("NDArray::copy::dims = $dims")

  # copy to NDArray and back
  array   = copy(tensor, mx.cpu())
  tensor2 = copy(array)
  @test tensor ≈ tensor2

  # copy between NDArray
  array2  = copy(array, mx.cpu())
  tensor2 = copy(array2)
  @test tensor ≈ tensor2

  @info("NDArray::copy::AbstractArray")
  let x = copy(1:4, mx.cpu())
    @test eltype(x) == Int
    @test copy(x) == [1, 2, 3, 4]
  end

  let x = copy(1.:4, mx.cpu())
    @test eltype(x) == Float64
    @test copy(x) ≈ [1., 2, 3, 4]
  end

  @info("NDArray::copy!::AbstractArray")
  let
    x = mx.zeros(4)
    copy!(x, 1:4)

    @test eltype(x) == Float32
    @test copy(x) == [1, 2, 3, 4]
  end
end

function test_deepcopy()
  @info("NDArray::deepcopy")

  x = mx.zeros(2, 5)
  y = deepcopy(x)
  x[:] = 42
  @test copy(x) != copy(y)
end

function test_assign()
  dims    = rand_dims()
  tensor  = rand(mx.MX_float, dims)

  @info("NDArray::assign::dims = $dims")

  # Julia Array -> NDArray assignment
  array    = NDArray(undef, size(tensor)...)
  array[:] = tensor
  @test tensor ≈ copy(array)

  array2  = mx.zeros(size(tensor))
  @test zeros(size(tensor)) ≈ copy(array2)

  array3 = mx.zeros(Float16, size(tensor))
  @test zeros(Float16, size(tensor)) ≈ copy(array2)

  # scalar -> NDArray assignment
  scalar    = rand()
  array2[:] = scalar
  @test zeros(size(tensor)) .+ scalar ≈ copy(array2)

  scalar = rand(Float16)
  array2[:] = scalar
  @test zeros(size(tensor)) .+ scalar ≈ copy(array2)

  scalar = rand(Float64)
  array2[:] = scalar
  array3[:] = scalar
  @test zeros(size(tensor)) .+ scalar ≈ copy(array2)
  @test zeros(Float16, size(tensor)) .+ scalar ≈ copy(array3)

  # NDArray -> NDArray assignment
  array[:]  = array2
  @test zeros(size(tensor)) .+ scalar ≈ copy(array)
end

function test_slice()
  array = mx.zeros((2, 4))
  array[2:3] = ones(2, 2)
  @test copy(array) == [0 1 1 0; 0 1 1 0]
  @test copy(mx.slice(array, 2:3)) == [1 1; 1 1]
end

function test_linear_idx()
  @info("NDArray::getindex::linear indexing")
  let A = reshape(1:30, 3, 10)
    x = mx.NDArray(A)

    @test copy(x) == A
    @test copy(x[1])  == [1]
    @test copy(x[2])  == [2]
    @test copy(x[3])  == [3]
    @test copy(x[12]) == [12]
    @test copy(x[13]) == [13]
    @test copy(x[14]) == [14]

    @test_throws BoundsError x[-1]
    @test_throws BoundsError x[0]
    @test_throws BoundsError x[31]
    @test_throws BoundsError x[42]
  end

  let A = reshape(1:24, 3, 2, 4)
    x = mx.NDArray(A)

    @test copy(x) == A
    @test copy(x[1])  == [1]
    @test copy(x[2])  == [2]
    @test copy(x[3])  == [3]
    @test copy(x[11]) == [11]
    @test copy(x[12]) == [12]
    @test copy(x[13]) == [13]
    @test copy(x[14]) == [14]
  end

  @info("NDArray::setindex!::linear indexing")
  let A = reshape(1:24, 3, 2, 4)
    x = mx.NDArray(A)

    @test copy(x) == A

    x[4] = -4
    @test copy(x[4])  == [-4]

    x[11] = -11
    @test copy(x[11]) == [-11]

    x[24] = 42
    @test copy(x[24]) == [42]
  end

  @info("NDArray::setindex!::type convert")
  let
    x = NDArray([1, 2, 3])
    @test eltype(x) == Int
    x[:] = π
    @test copy(x) == [3, 3, 3]
  end
end  # function test_linear_idx

function test_first()
  @info("NDArray::first")
  let A = reshape(1:30, 3, 10)
    x = mx.NDArray(A)

    @test x[]    == 1
    @test x[5][] == 5

    @test first(x)    == 1
    @test first(x[5]) == 5
  end
end  # function test_first

function test_lastindex()
  @info("NDArray::lastindex")
  let A = [1 2; 3 4; 5 6], x = mx.NDArray(A)
    @test lastindex(A) == lastindex(x)
  end
end  # function test_lastindex

function test_cat()
  function check_cat(f, A, B = 2A)
    C = [A B]
    D = [A; B]
    x = NDArray(A)
    y = NDArray(B)
    z = NDArray(C)
    d = NDArray(D)

    if f == :hcat
      @test copy([x y]) == [A B]
      @test copy([x y 3y x]) == [A B 3B A]
      @test copy([z y x]) == [C B A]
    elseif f == :vcat
      @test copy([x; y]) == [A; B]
      @test copy([x; y; 3y; x]) == [A; B; 3B; A]
      @test copy([x; d]) == [A; D]
      @test copy([d; x]) == [D; A]
    else
      @assert false
    end
  end

  let A = [1, 2, 3, 4]
    @info("NDArray::hcat::1D")
    check_cat(:hcat, A)

    @info("NDArray::vcat::1D")
    check_cat(:vcat, A)
  end

  let A = [1 2; 3 4]
    @info("NDArray::hcat::2D")
    check_cat(:hcat, A)

    @info("NDArray::vcat::2D")
    check_cat(:vcat, A)
  end

  let A = rand(4, 3, 2)
    @info("NDArray::hcat::3D")
    check_cat(:hcat, A)

    @info("NDArray::vcat::3D")
    check_cat(:vcat, A)
  end

  let A = rand(4, 3, 2, 2)
    @info("NDArray::hcat::4D")
    check_cat(:hcat, A)

    @info("NDArray::vcat::4D")
    check_cat(:vcat, A)
  end

  let A = [1, 2, 3, 4]
    @info("NDArray::cat::3D/1D")
    check_cat(:vcat, reshape(A, 4, 1, 1), 2A)
  end
end  # function test_cat

function test_plus()
  dims   = rand_dims()
  t1, a1 = rand_tensors(dims)
  t2, a2 = rand_tensors(dims)
  t3, a3 = rand_tensors(dims)

  @info("NDArray::plus::dims = $dims")

  @test t1 .+ t2 ≈ copy(a1 .+ a2)

  @test t1 .+ t2 .+ t3 ≈ copy(a1 .+ a2 .+ a3)

  # test inplace += operation
  a0 = a1               # keep a reference to a1
  @mx.inplace a1 += a2  # perform inplace +=
  @test a0 == a1        # make sure they are still the same object
  @test copy(a0) ≈ copy(a1)
  @test copy(a1) ≈ t1 .+ t2

  # test scalar
  scalar = rand()
  @test t3 .+ scalar       ≈ copy(a3 .+ scalar)
  @test t2 .+ scalar .+ t3 ≈ copy(a2 .+ scalar .+ a3)

  # test small and large scalar
  t4 = zeros(Float32, dims)
  a4 = copy(t4, mx.cpu())
  scalar_small = 1e-8
  scalar_large = 1e8
  @test t4 .+ scalar_small ≈ copy(a4 .+ scalar_small)
  @test t4 .+ scalar_large ≈ copy(a4 .+ scalar_large)

  t5 = zeros(Float64, dims)
  a5 = copy(t5, mx.cpu())
  scalar_small = 1e-8
  scalar_large = 1e8
  @test t5 .+ scalar_small ≈ copy(a5 .+ scalar_small)
  @test t5 .+ scalar_large ≈ copy(a5 .+ scalar_large)

  t6 = zeros(Float16, dims)
  a6 = copy(t6, mx.cpu())
  scalar_small = Float16(1e-5)
  scalar_large = Float16(1e4)
  @test t6 .+ scalar_small ≈ copy(a6 .+ scalar_small)
  @test t6 .+ scalar_large ≈ copy(a6 .+ scalar_large)

  let x = mx.NDArray([1 2; 3 4]), y = mx.NDArray([1 1; 1 1])
    @test copy(42 .+ x) == [43 44; 45 46]
    @test copy(x .+ 42) == [43 44; 45 46]
    @test copy(0 .+ x .+ y .+ 41) == [43 44; 45 46]
  end

  @info("NDArray::plus::scalar::type convert")
  let x = mx.NDArray([1, 2, 3])
    y = x .+ 0.5
    @test copy(y) == copy(x)

    y = x .+ 2.9
    @test copy(y) == [3, 4, 5]
  end

  @info("NDArray::broadcast_add")
  let
    A = [1 2 3;
         4 5 6]
    B = [1,
         2]
    x = NDArray(A)
    y = NDArray(B)

    z = x .+ y
    @test copy(z) == A .+ B

    # TODO
    # @inplace x .+= y
    # @test copy(x) == A .+ B
  end
end

function test_minus()
  dims   = rand_dims()
  t1, a1 = rand_tensors(dims)
  t2, a2 = rand_tensors(dims)

  @info("NDArray::minus::dims = $dims")

  @test t1 .- t2 ≈ copy(a1 .- a2)

  @test -t1 ≈ copy(-a1)

  # make sure the negation is not in-place, so a1 is not changed after previous
  # statement is executed
  @test t1 ≈ copy(a1)

  # test inplace -= operation
  a0 = a1              # keep a reference to a1
  @mx.inplace a1 -= a2 # perform inplace -=
  @test a0 == a1       # make sure they are still the same object
  @test a0.handle == a1.handle
  @test copy(a0) ≈ copy(a1)
  @test copy(a1) ≈ t1 .- t2

  # test scalar
  scalar = rand()
  @test t2 .- scalar ≈ copy(a2 .- scalar)

  # test small and large scalar
  t4 = zeros(Float32, dims)
  a4 = copy(t4, mx.cpu())
  scalar_small = 1e-8
  scalar_large = 1e8
  @test t4 .- scalar_small ≈ copy(a4 .- scalar_small)
  @test t4 .- scalar_large ≈ copy(a4 .- scalar_large)

  t5 = zeros(Float64, dims)
  a5 = copy(t5, mx.cpu())
  scalar_small = 1e-8
  scalar_large = 1e8
  @test t5 .- scalar_small ≈ copy(a5 .- scalar_small)
  @test t5 .- scalar_large ≈ copy(a5 .- scalar_large)

  t6 = zeros(Float16, dims)
  a6 = copy(t6, mx.cpu())
  scalar_small = Float16(1e-5)
  scalar_large = Float16(1e4)
  @test t6 .- scalar_small ≈ copy(a6 .- scalar_small)
  @test t6 .- scalar_large ≈ copy(a6 .- scalar_large)

  @info("NDArray::minus::scalar::type convert")
  let x = mx.NDArray([1, 2, 3])
    @test copy(x .- π) ≈ [-2, -1, 0]
  end

  @info("NDArray::broadcast_minus")
  let
    A = [1 2 3;
         4 5 6]
    B = [1,
         2]
    x = NDArray(A)
    y = NDArray(B)

    z = x .- y
    @test copy(z) == A .- B

    # TODO
    # @inplace x .-= y
    # @test copy(x) == A .- B
  end

  @info("NDArray::scalar::rminus")
  let
    A = [1 2 3;
         4 5 6]
    B = 10 .- A

    x = NDArray(A)
    y = 10 .- x

    @test copy(y) == B
  end
end

function test_mul()
  dims   = rand_dims()
  t1, a1 = rand_tensors(dims)
  t2, a2 = rand_tensors(dims)
  t3, a3 = rand_tensors(dims)

  @info("NDArray::mul::dims = $dims")

  @test t1 .* t2 ≈ copy(a1 .* a2)

  # test inplace .*= operation
  a0 = a1               # keep a reference to a1
  @mx.inplace a1 .*= a2 # perform inplace .*=
  @test a0 == a1        # make sure they are still the same object
  @test a0.handle == a1.handle
  @test copy(a0) ≈ copy(a1)
  @test copy(a1) ≈ t1 .* t2

  # test scalar
  scalar = mx.MX_float(rand())
  @test t3 .* scalar ≈ copy(a3 .* scalar)

  # test small and large scalar
  t4, a4 = rand_tensors(Float32, dims)
  scalar_small = 1e-8
  scalar_large = 1e8
  @test t4 * scalar_small ≈ copy(a4 .* scalar_small)
  @test t4 * scalar_large ≈ copy(a4 .* scalar_large)

  t5, a5 = rand_tensors(Float64, dims)
  scalar_small = 1e-8
  scalar_large = 1e8
  @test t5 * scalar_small ≈ copy(a5 .* scalar_small)
  @test t5 * scalar_large ≈ copy(a5 .* scalar_large)

  t6, a6 = rand_tensors(Float16, dims)
  scalar_small = Float16(1e-5)
  @test t6 * scalar_small ≈ copy(a6 .* scalar_small)

  @info("NDArray::mul::matrix multiplication")
  let x = mx.NDArray([1.  2])
    y = x' * x
    @test copy(y) == [1. 2; 2 4]
  end

  @info("NDArray::mul::elementwise::issue 253")
  let x = mx.NDArray([1.  2])
    y = x .* x
    @test copy(y) == [1. 4.]
  end

  @info("NDArray::mul::scalar::type convert")
  let x = mx.NDArray([1, 2, 3])
    y = x .* π
    @test eltype(x) == Int
    @test copy(y) == [3, 6, 9]
  end

  @info("NDArray::broadcast_mul")
  let
    A = [1 2 3;
         4 5 6]
    B = [1,
         2]
    x = NDArray(A)
    y = NDArray(B)

    z = x .* y
    @test copy(z) == A .* B

    # TODO
    # @inplace x .*= y
    # @test copy(x) == A .* B
  end
end

function test_div()
  dims   = rand_dims()
  t1, a1 = rand_tensors(dims)
  t2, a2 = rand_tensors(dims)

  @info("NDArray::div::dims = $dims")
  t2             .+= 2  # avoid numerical instability
  @mx.inplace a2 .+= 2

  @test t1 ./ t2 ≈ copy(a1 ./ a2)

  # test inplace -= operation
  a0 = a1                # keep a reference to a2
  @mx.inplace a1 ./= a2  # perform inplace ./=
  @test a0 == a1         # make sure they are still the same object
  @test a0.handle == a1.handle
  @test copy(a0) ≈ copy(a1)
  @test copy(a1) ≈ t1 ./ t2

  # test scalar
  scalar = rand() + 2
  @test t2 ./ scalar ≈ copy(a2 ./ scalar)

  # test small and large scalar
  t4, a4 = rand_tensors(Float32, dims)
  scalar_small = 1e-8
  scalar_large = 1e8
  @test t4 ./ scalar_small ≈ copy(a4 ./ scalar_small)
  @test t4 ./ scalar_large ≈ copy(a4 ./ scalar_large)

  t5, a5 = rand_tensors(Float64, dims)
  scalar_small = 1e-8
  scalar_large = 1e8
  @test t5 ./ scalar_small ≈ copy(a5 ./ scalar_small)
  @test t5 ./ scalar_large ≈ copy(a5 ./ scalar_large)

  t6, a6 = rand_tensors(Float16, dims)
  scalar_large = 1e4
  @test t6 ./ scalar_large ≈ copy(a6 ./ scalar_large)

  @info("NDArray::div::scalar::type convert")
  let x = mx.NDArray([1, 2, 3])
    y = x ./ 1.1
    @test eltype(y) == Int
    @test copy(y) == [1, 2, 3]

    y = x ./ 2
    @test eltype(y) == Int  # this differs from julia
    @test copy(y) == [0, 1, 1]

    @test_throws AssertionError x ./ 0.5
  end

  @info("NDArray::broadcast_div")
  let
    A = Float32[1 2 3;
                4 5 6]
    B = Float32[1,
                2]
    x = NDArray(A)
    y = NDArray(B)

    z = x ./ y
    @test copy(z) == A ./ B

    # TODO
    # @inplace x ./= y
    # @test copy(x) == A ./ B
  end
end

function test_rdiv()
  @info("NDArray::rdiv")

  @info("NDArray::rdiv::Inf16")
  let x = 1 ./ mx.zeros(Float16, 4)
    @test copy(x) == [Inf16, Inf16, Inf16, Inf16]
  end

  @info("NDArray::rdiv::Inf32")
  let x = 1 ./ mx.zeros(Float32, 4)
    @test copy(x) == [Inf32, Inf32, Inf32, Inf32]
  end

  @info("NDArray::rdiv::Inf64")
  let x = 1 ./ mx.zeros(Float64, 4)
    @test copy(x) == [Inf64, Inf64, Inf64, Inf64]
  end

  @info("NDArray::rdiv::Int")
  let x = 1 ./ mx.NDArray([1 2; 3 4])
    @test copy(x) == [1 0; 0 0]
  end

  @info("NDArray::rdiv::Float32")
  let x = 1 ./ mx.NDArray(Float32[1 2; 3 4])
    y = 1 ./ Float32[1 2; 3 4]
    @test copy(x) ≈ y
  end

  @info("NDArray::rdiv::type convert")
  let x = mx.NDArray([1, 2, 3])
    y = 5.5 ./ x
    @test eltype(y) == Int  # this differs from julia
    @test copy(y) == [5, 2, 1]
  end
end  # function test_rdiv

function test_mod()
  @info("NDArray::mod")
  A = [1 2; 3 4]
  B = [1 1; 3 3]

  let x = NDArray(A), y = NDArray(B)
    C = A .% B
    D = B .% A

    w = x .% y
    z = y .% x

    @test copy(w) ≈ C
    @test copy(z) ≈ D
  end

  @info("NDArray::mod::scalar")
  let x = NDArray(A)
    C = A .% 2
    y = x .% 2
    @test copy(y) ≈ C
  end

  @info("NDArray::rmod")
  let x = NDArray(A)
    C = 11 .% A
    y = 11 .% x
    @test copy(y) ≈ C
  end

  @info("NDArray::mod_from!")
  let
    x = NDArray(A)
    y = NDArray(B)
    C = A .% B
    mx.mod_from!(x, y)
    @test copy(x) ≈ C
  end

  let
    x = NDArray(A)
    y = NDArray(B)
    C = B .% A
    mx.mod_from!(y, x)

    @test copy(y) ≈ C
  end

  @info("NDArray::mod_from!::scalar")
  let
    x = NDArray(A)
    C = A .% 2
    mx.mod_from!(x, 2)
    @test copy(x) ≈ C
  end

  @info("NDArray::rmod_from!")
  let
    x = NDArray(A)
    C = 11 .% A
    mx.rmod_from!(11, x)
    @test copy(x) ≈ C
  end

  @info("NDArray::mod_from!::writable")
  let
    x = NDArray(A)
    y = NDArray(B)
    x.writable = false
    y.writable = false
    @test_throws AssertionError mx.mod_from!(x, y)
    @test_throws AssertionError mx.mod_from!(y, x)
    @test_throws AssertionError mx.mod_from!(x, 2)
    @test_throws AssertionError mx.rmod_from!(2, x)
  end

  @info("NDArray::mod::inplace")
  let
    x = NDArray(A)
    y = NDArray(B)
    C = A .% B
    @inplace x .%= y
    @test copy(x) ≈ C
  end

  @info("NDArray::broadcast_mod")
  let
    A = [1 2 3;
         4 5 6]
    B = [1,
         2]
    x = NDArray(A)
    y = NDArray(B)

    z = x .% y
    @test copy(z) == A .% B

    # TODO
    # @inplace x .%= y
    # @test copy(x) == A .% B
  end
end  # function test_mod

function test_gd()
  dims   = rand_dims()
  tw, aw = rand_tensors(dims)
  tg, ag = rand_tensors(dims)

  @info("NDArray::gd::dims = $dims")

  lr = rand()
  wd = rand()

  @mx.inplace aw += -lr * (ag + wd * aw)
  tw += -lr * (tg + wd * tw)
  @test copy(aw) ≈ tw
end

function test_saveload()
  n_arrays = 5
  @info("NDArray::saveload::n_arrays = $n_arrays")
  fname = tempname()

  # save and load a single array
  dims   = rand_dims()
  j_array, nd_array = rand_tensors(dims)
  mx.save(fname, nd_array)
  data = mx.load(fname, mx.NDArray)
  @test data isa Vector{<:mx.NDArray}
  @test length(data) == 1
  @test copy(data[1]) ≈ j_array

  # save and load N arrays of different shape
  arrays = [rand_tensors(rand_dims()) for i = 1:n_arrays]
  nd_arrays = mx.NDArray[x[2] for x in arrays]
  mx.save(fname, nd_arrays)
  data = mx.load(fname, mx.NDArray)
  @test data isa Vector{<:mx.NDArray}
  @test length(data) == n_arrays
  for i = 1:n_arrays
    @test copy(data[i]) ≈ arrays[i][1]
  end

  # save and load dictionary of ndarrays
  names = [Symbol("array$i") for i = 1:n_arrays]
  dict = Dict([(n, v) for (n,v) in zip(names, nd_arrays)])
  mx.save(fname, dict)
  data = mx.load(fname, mx.NDArray)
  @test data isa Dict{Symbol,<:mx.NDArray}
  @test length(data) == n_arrays
  for i = 1:n_arrays
    @test copy(data[names[i]]) ≈ arrays[i][1]
  end

  rm(fname)
end

function test_clip()
  dims = rand_dims()
  @info("NDArray::clip::dims = $dims")

  j_array, nd_array = rand_tensors(dims)
  clip_up   = maximum(abs.(j_array)) / 2
  clip_down = 0
  clipped   = clip(nd_array, clip_down, clip_up)

  # make sure the original array is not modified
  @test copy(nd_array) ≈ j_array

  @test all(clip_down .<= copy(clipped) .<= clip_up)

  @info("NDArray::clip!")
  let
    x = NDArray(1.0:20)
    clip!(x, 5, 15)
    @test all(5 .<= copy(x) .<= 15)
  end
end

function test_power()
  @info("NDArray::power")

  @info("NDArray::power::Int::x .^ n")
  let x = mx.NDArray([1 2; 3 4])
    @test eltype(x) == Int
    @test copy(x .^ -1)  == [1 0; 0 0]
    @test copy(x .^ 0)   == [1 1; 1 1]
    @test copy(x .^ 1)   == [1 2; 3 4]
    @test copy(x .^ 1.1) == [1 2; 3 4]
    @test copy(x .^ 2)   == [1 4; 9 16]
    @test copy(x .^ 2.9) == [1 4; 9 16]
    @test copy(x .^ 3)   == [1 8; 27 64]
  end

  @info("NDArray::power::Int::n .^ x")
  let x = mx.NDArray([1 2; 3 4])
    @test eltype(x) == Int
    @test copy(0   .^ x)   == [0 0; 0 0]
    @test copy(1   .^ x)   == [1 1; 1 1]
    @test copy(1.1 .^ x) == [1 1; 1 1]
    @test copy(2   .^ x)   == [2 4; 8 16]
    @test copy(2.9 .^ x) == [2 4; 8 16]
    @test copy(3   .^ x)   == [3 9; 27 81]
  end

  @info("NDArray::power::Int::x .^ y")
  let x = mx.NDArray([1 2; 3 4]), y = mx.NDArray([2 2; 2 2])
    @test eltype(x) == Int
    @test eltype(y) == Int
    @test copy(x .^ y) == [1 4; 9 16]
    @test copy(y .^ x) == [2 4; 8 16]
  end

  @info("NDArray::power::Float32::x .^ n")
  let x = mx.NDArray(Float32[1 2; 3 4]), A = Float32[1 2; 3 4]
    @test eltype(x) == Float32
    @test copy(x .^ 0) == Float32[1 1; 1 1]
    @test copy(x .^ 1) == Float32[1 2; 3 4]
    @test copy(x .^ 2) == Float32[1 4; 9 16]
    @test copy(x .^ 3) == Float32[1 8; 27 64]

    @test copy(x .^ -1)  ≈ A .^ -1
    @test copy(x .^ 1.1) ≈ A .^ 1.1
    @test copy(x .^ 2.9) ≈ A .^ 2.9
  end

  @info("NDArray::power::Float32::n .^ x")
  let x = mx.NDArray(Float32[1 2; 3 4]), A = Float32[1 2; 3 4]
    @test eltype(x) == Float32
    @test copy(0 .^ x) == Float32[0 0; 0 0]
    @test copy(1 .^ x) == Float32[1 1; 1 1]
    @test copy(2 .^ x) == Float32[2 4; 8 16]
    @test copy(3 .^ x) == Float32[3 9; 27 81]

    @test copy(1.1 .^ x) ≈ 1.1 .^ A
    @test copy(2.9 .^ x) ≈ 2.9 .^ A
  end

  @info("NDArray::power::Float32::x .^ y")
  let x = mx.NDArray(Float32[1 2; 3 4]), y = mx.NDArray(Float32[2 2; 2 2])
    @test eltype(x) == Float32
    @test eltype(y) == Float32
    @test copy(x .^ y) == Float32[1 4; 9 16]
    @test copy(y .^ x) == Float32[2 4; 8 16]
  end

  @info("NDArray::power::ℯ .^ x::x .^ ℯ")
  let x = mx.zeros(2, 3), A = [1 1 1; 1 1 1]
    @test copy(ℯ .^ x) ≈ A
  end

  let A = Float32[1 2; 3 4], x = mx.NDArray(A)
    @test copy(ℯ .^ x) ≈ ℯ .^ A
    @test copy(x .^ ℯ) ≈ A .^ ℯ
  end

  @info("NDArray::power::π .^ x::x .^ π")
  let A = Float32[1 2; 3 4], x = mx.NDArray(A)
    @test copy(π .^ x) ≈ π .^ A
    @test copy(x .^ π) ≈ A .^ π
  end

  # TODO: Float64: wait for https://github.com/apache/incubator-mxnet/pull/8012

  @info("NDArray::broadcast_power")
  let
    A = [1 2 3;
         4 5 6]
    B = [1,
         2]
    x = NDArray(A)
    y = NDArray(B)

    z = x.^y
    @test copy(z) == A.^B

    # TODO
    # @inplace x .^= y
    # @test copy(x) == A.^B
  end
end # function test_power

function test_sqrt()
  dims = rand_dims()
  @info("NDArray::sqrt::dims = $dims")

  j_array, nd_array = rand_tensors(dims)
  sqrt_ed = sqrt(nd_array)
  @test copy(sqrt_ed) ≈ sqrt.(j_array)
end

function test_nd_as_jl()
  dims = (2, 3)
  @info("NDArray::nd_as_jl::dims = $dims")

  x = mx.zeros(dims) + 5
  y = mx.ones(dims)
  z = mx.zeros(dims)
  @mx.nd_as_jl ro=x rw=(y, z) begin
    for i = 1:length(z)
      z[i] = x[i]
    end

    z[:, 1] = y[:, 1]
    y .= 0
  end

  @test sum(copy(y)) == 0
  @test sum(copy(z)[:, 1]) == 2
  @test copy(z)[:, 2:end] ≈ copy(x)[:, 2:end]
end

function test_dot()
  dims1 = (2, 3)
  dims2 = (3, 8)
  @info("NDArray::dot")

  x = mx.zeros(dims1)
  y = mx.zeros(dims2)
  z = mx.dot(x, y)
  @test size(z) == (2, 8)

  x = mx.zeros(1, 2)
  y = mx.zeros(1, 2, 3)
  @test_throws mx.MXError dot(x, y)  # dimension mismatch

  @info("NDArray::matrix mul")
  let
    A = [1. 2 3; 4 5 6]
    B = [-1., -2, -3]
    x = NDArray(A)
    y = NDArray(B)
    z = x * y
    @test copy(z) == A * B
    @test size(z) == (2,)
  end

  let
    A = [1. 2 3; 4 5 6]
    B = [-1. -2; -3 -4; -5 -6]
    x = NDArray(A)
    y = NDArray(B)
    z = x * y
    @test copy(z) == A * B
    @test size(z) == (2, 2)
  end
end

function test_eltype()
  @info("NDArray::eltype")
  dims = (3,3)

  x = NDArray(undef, dims)
  @test eltype(x) == mx.DEFAULT_DTYPE

  for TF in instances(mx.TypeFlag)
    T = mx.fromTypeFlag(TF)
    x = NDArray{T}(undef, dims)
    @test eltype(x) == T
  end
end

function test_reshape()
  @info("NDArray::reshape")
  A = rand(2, 3, 4)

  B = reshape(NDArray(A), 4, 3, 2)
  @test size(B) == (4, 3, 2)
  @test copy(B)[3, 1, 1] == A[1, 2, 1]

  C = reshape(NDArray(A), (4, 3, 2))
  @test size(C) == (4, 3, 2)
  @test copy(C)[3, 1, 1] == A[1, 2, 1]

  @info("NDArray::reshape::reverse")
  A = mx.zeros(10, 5, 4)

  B = reshape(A, -1, 0)
  @test size(B) == (40, 5)

  C = reshape(A, -1, 0, reverse=true)
  @test size(C) == (50, 4)
end

function test_expand_dims()
  @info("NDArray::expand_dims")
  let A = [1, 2, 3, 4], x = NDArray(A)
    @test size(x) == (4,)

    y = expand_dims(x, 1)
    @test size(y) == (1, 4)

    y = expand_dims(x, 2)
    @test size(y) == (4, 1)
  end

  let A = [1 2; 3 4; 5 6], x = NDArray(A)
    @test size(x) == (3, 2)

    y = expand_dims(x, 1)
    @test size(y) == (1, 3, 2)

    y = expand_dims(x, 2)
    @test size(y) == (3, 1, 2)

    y = expand_dims(x, 3)
    @test size(y) == (3, 2, 1)
  end
end  # test_expand_dims

function test_sum()
  @info("NDArray::sum")

  let A = reshape(1.0:8, 2, 2, 2), X = mx.NDArray(A)
    @test copy(sum(X))[]              == sum(A)
    @test copy(sum(X, dims = 1))      == sum(A, dims = 1)
    @test copy(sum(X, dims = 2))      == sum(A, dims = 2)
    @test copy(sum(X, dims = 3))      == sum(A, dims = 3)
    @test copy(sum(X, dims = [1, 2])) == sum(A, dims = [1, 2])
    @test copy(sum(X, dims = (1, 2))) == sum(A, dims = (1, 2))
  end
end

function test_mean()
  @info("NDArray::mean")

  let A = reshape(1.0:8, 2, 2, 2), X = mx.NDArray(A)
    @test copy(mean(X))[]              == mean(A)
    @test copy(mean(X, dims = 1))      == mean(A, dims = 1)
    @test copy(mean(X, dims = 2))      == mean(A, dims = 2)
    @test copy(mean(X, dims = 3))      == mean(A, dims = 3)
    @test copy(mean(X, dims = [1, 2])) == mean(A, dims = [1, 2])
    @test copy(mean(X, dims = (1, 2))) == mean(A, dims = (1, 2))
  end
end

function test_maximum()
  @info("NDArray::maximum")

  let A = reshape(1.0:8, 2, 2, 2), X = mx.NDArray(A)
    @test copy(maximum(X))[]              == maximum(A)
    @test copy(maximum(X, dims = 1))      == maximum(A, dims = 1)
    @test copy(maximum(X, dims = 2))      == maximum(A, dims = 2)
    @test copy(maximum(X, dims = 3))      == maximum(A, dims = 3)
    @test copy(maximum(X, dims = [1, 2])) == maximum(A, dims = [1, 2])
    @test copy(maximum(X, dims = (1, 2))) == maximum(A, dims = (1, 2))
  end

  @info("NDArray::broadcast_maximum")
  let
    A = [1 2 3;
         4 5 6]
    B = [1,
         2]
    x = NDArray(A)
    y = NDArray(B)

    z = max.(x, y)
    @test copy(z) == max.(A, B)
  end
end

function test_minimum()
  @info("NDArray::minimum")

  let A = reshape(1.0:8, 2, 2, 2), X = mx.NDArray(A)
    @test copy(minimum(X))[]              == minimum(A)
    @test copy(minimum(X, dims = 1))      == minimum(A, dims = 1)
    @test copy(minimum(X, dims = 2))      == minimum(A, dims = 2)
    @test copy(minimum(X, dims = 3))      == minimum(A, dims = 3)
    @test copy(minimum(X, dims = [1, 2])) == minimum(A, dims = [1, 2])
    @test copy(minimum(X, dims = (1, 2))) == minimum(A, dims = (1, 2))
  end

  @info("NDArray::broadcast_minimum")
  let
    A = [1 2 3;
         4 5 6]
    B = [1,
         2]
    x = NDArray(A)
    y = NDArray(B)

    z = min.(x, y)
    @test copy(z) == min.(A, B)
  end
end

function test_prod()
  @info("NDArray::prod")

  let A = reshape(1.0:8, 2, 2, 2), X = mx.NDArray(A)
    @test copy(prod(X))[]              == prod(A)
    @test copy(prod(X, dims = 1))      == prod(A, dims = 1)
    @test copy(prod(X, dims = 2))      == prod(A, dims = 2)
    @test copy(prod(X, dims = 3))      == prod(A, dims = 3)
    @test copy(prod(X, dims = [1, 2])) == prod(A, dims = [1, 2])
    @test copy(prod(X, dims = (1, 2))) == prod(A, dims = (1, 2))
  end
end

function test_fill()
  @info("NDArray::fill")

  let x = mx.fill(42, 2, 3, 4)
    @test eltype(x) == Int
    @test size(x) == (2, 3, 4)
    @test copy(x) == fill(42, 2, 3, 4)
  end

  let x = mx.fill(Float32(42), 2, 3, 4)
    @test eltype(x) == Float32
    @test size(x) == (2, 3, 4)
    @test copy(x) ≈ fill(Float32(42), 2, 3, 4)
  end

  let x = mx.fill(42, (2, 3, 4))
    @test eltype(x) == Int
    @test size(x) == (2, 3, 4)
    @test copy(x) == fill(42, 2, 3, 4)
  end

  let x = mx.fill(Float32(42), (2, 3, 4))
    @test eltype(x) == Float32
    @test size(x) == (2, 3, 4)
    @test copy(x) ≈ fill(Float32(42), 2, 3, 4)
  end

  @info("NDArray::fill!::arr")
  let x = fill!(mx.zeros(2, 3, 4), 42)
    @test eltype(x) == Float32
    @test size(x) == (2, 3, 4)
    @test copy(x) ≈ fill(Float32(42), 2, 3, 4)
  end
end  # function test_fill

function test_transpose()
  @info("NDArray::transpose::1D")
  let A = rand(Float32, 4), x = NDArray(A)
    @test size(x) == (4,)
    @test size(x') == (1, 4)
  end

  @info("NDArray::transpose::2D")
  let A = rand(Float32, 2, 3), x = mx.NDArray(A)
    @test size(x) == (2, 3)
    @test size(x') == (3, 2)
  end

  @info("NDArray::permutedims")
  let A = collect(Float32, reshape(1.0:24, 2, 3, 4)), x = mx.NDArray(A)
    A′ = permutedims(A, [2, 1, 3])
    x′ = permutedims(x, [2, 1, 3])
    @test size(A′) == size(x′)
    @test A′ == copy(x′)
  end
end

function test_show()
  @info("NDArray::show::REPL")
  let str = sprint(show, MIME"text/plain"(), mx.NDArray([1 2 3 4]))
    @test occursin("1×4", str)
    @test occursin("NDArray", str)
    @test occursin("Int64", str)
    @test occursin("CPU", str)
    @test match(r"1\s+2\s+3\s+4", str) != nothing
  end

  @info("NDArray::show")
  let str = sprint(show, mx.NDArray([1 2 3 4]))
    @test str == "NDArray([1 2 3 4])"
  end

  let str = sprint(show, mx.zeros(4))
    @test str == "NDArray(Float32[0.0, 0.0, 0.0, 0.0])"
  end
end

function test_size()
  @info("NDArray::size")
  let A = [1 2; 3 4; 5 6], x = mx.NDArray(A)
    @test size(A) == size(x)
    dims = (1, 2, 3, 4, 5)
    @test map(d -> size(A, d), dims) == map(d -> size(x, d), dims)
    @inferred map(d -> size(x, d), dims)
  end
end  # function test_size()

function check_trigonometric(f)
  @info("NDArray::$f")
  let A = [.1 .2; .3 .4], x = mx.NDArray(A)
    B = f.(A)
    y = f.(x)
    @test copy(y) ≈ B
  end

  let A = Float32[.1 .2; .3 .4], x = mx.NDArray(A)
    B = f.(A)
    y = f.(x)
    @test copy(y) ≈ B
  end
end  # function check_trigonometric

function test_trigonometric()
  for f ∈ [sin, cos, tan, asin, acos, atan]
    check_trigonometric(f)
  end
end  # function test_trigonometric

function check_hyperbolic(f, A)
  @info("NDArray::$f")
  let x = NDArray(A)
    B = f.(A)
    y = f.(x)
    @test copy(y) ≈ B
  end

  let A = Float32.(A), x = NDArray(A)
    B = f.(A)
    y = f.(x)
    @test copy(y) ≈ B
  end
end  # function check_hyperbolic

function test_hyperbolic()
  for f ∈ [sinh, cosh, tanh, asinh, acosh, atanh]
    A = if f == acosh
      [1.1, 1.2, 1.3, 1.4]
    else
      [.1, .2, .3, .4]
    end
    check_hyperbolic(f, A)
  end
end  # function test_hyperbolic

function test_act_funcs()
  @info("NDArray::σ/sigmoid")
  let
    A = Float32[.1, .2, -.3, -.4]
    B = @. 1 / (1 + ℯ ^ (-A))
    x = NDArray(A)
    y = σ.(x)
    @test copy(y) ≈ B

    z = sigmoid.(x)
    @test copy(z) ≈ B
  end

  @info("NDArray::relu")
  let
    A = [1, 2, -3, -4]
    B = max.(A, 0)
    x = NDArray(A)
    y = relu.(x)
    @test copy(y) ≈ B
  end

  @info("NDArray::softmax::1D")
  let
    A = Float32[1, 2, 3, 4]
    B = exp.(A) ./ sum(exp.(A))
    x = NDArray(A)
    y = softmax.(x)
    @test copy(y) ≈ B
  end

  @info("NDArray::softmax::2D")
  let
    A = Float32[1 2; 3 4]
    B = exp.(A) ./ sum(exp.(A), dims = 1)
    x = NDArray(A)
    y = softmax.(x, 1)
    @test copy(y) ≈ B

    C = exp.(A) ./ sum(exp.(A), dims = 2)
    z = softmax.(x, 2)
    @test copy(z) ≈ C
  end

  @info("NDArray::log_softmax::1D")
  let
    A = Float32[1, 2, 3, 4]
    B = log.(exp.(A) ./ sum(exp.(A)))
    x = NDArray(A)
    y = log_softmax.(x)
    @test copy(y) ≈ B
  end

  @info("NDArray::log_softmax::2D")
  let
    A = Float32[1 2; 3 4]
    B = log.(exp.(A) ./ sum(exp.(A), dims = 1))
    x = NDArray(A)
    y = log_softmax.(x, 1)
    @test copy(y) ≈ B

    C = log.(exp.(A) ./ sum(exp.(A), dims = 2))
    z = log_softmax.(x, 2)
    @test copy(z) ≈ C
  end
end  # function test_act_funcs

macro check_equal(op)
  quote
    A = [1 2 3
         4 5 6]
    B = [1,
         6]
    x = NDArray(A)
    y = NDArray(B)
    a = broadcast($op, x, y)
    @test copy(a) == broadcast($op, A, B)

    C = [3 2 1
         6 5 4]
    z = NDArray(C)
    b = broadcast($op, x, z)
    @test copy(b) == broadcast($op, A, C)
  end
end

function test_equal()
  @info("NDArray::broadcast_equal")
  @check_equal ==

  @info("NDArray::broadcast_not_equal")
  @check_equal !=

  @info("NDArray::broadcast_greater")
  @check_equal >

  @info("NDArray::broadcast_greater_equal")
  @check_equal >=

  @info("NDArray::broadcast_lesser")
  @check_equal <

  @info("NDArray::broadcast_lesser_equal")
  @check_equal <=
end  # function test_equal

function test_broadcast_to()
  @info("NDArray::broadcast_to")
  A = [1 2 3]
  x = NDArray(A)
  @test mx.broadcast_to(x, (1, 3)) |> copy == A
  @test mx.broadcast_to(x, (5, 3)) |> copy == repeat(A, outer = (5, 1))

  @test mx.broadcast_to(x, 1, 3) |> copy == A
  @test mx.broadcast_to(x, 5, 3) |> copy == repeat(A, outer = (5, 1))
end  # function test_broadcast_to

function test_broadcast_axis()
  @info("NDArray::broadcast_axis")
  A = reshape([1, 2, 3], 1, 3, 1)
  x = NDArray(A)

  @test mx.broadcast_axis(x, 1, 4) |> copy == [A; A; A; A]
  @test mx.broadcast_axis(x, 3, 2) |> copy == cat(A, A, dims = 3)

  @info("NDArray::broadcast_axes")
  @test mx.broadcast_axes(x, 1, 4) |> copy == [A; A; A; A]
  @test mx.broadcast_axes(x, 3, 2) |> copy == cat(A, A, dims = 3)
end  # function test_broadcast_axis

function test_hypot()
  @info("NDArray::hypot")
  A = [3 3 3]
  B = [4, 4]
  C = hypot.(A, B)

  x = NDArray(A)
  y = NDArray(B)
  z = hypot.(x, y)

  @test copy(z) == C
end  # function test_hypot

function test_argmax()
  @info "NDArray::argmax"
  let
    A = [1. 5 3;
         4 2 6]
    x = NDArray(A)

    @test copy(argmax(x, dims = 1)) == [2 1 2]
    @test copy(argmax(x, dims = 2)) == reshape([2, 3], :, 1)
  end

  @info "NDArray::argmax::NaN"
  let
    A = [1.  5 3;
         NaN 2 6]
    x = NDArray(A)

    @test copy(argmax(x, dims = 1)) == [1 1 2]
    @test copy(argmax(x, dims = 2)) == reshape([2, 3], :, 1)
  end
end

function test_argmin()
  @info "NDArray::argmin"
  let
    A = [1. 5 3;
         4 2 6]
    x = NDArray(A)

    @test copy(argmin(x, dims = 1)) == [1 2 1]
    @test copy(argmin(x, dims = 2)) == reshape([1, 2], :, 1)
  end

  @info "NDArray::argmin::NaN"
  let
    A = [1.  5 3;
         NaN 2 6]
    x = NDArray(A)

    @test copy(argmin(x, dims = 1)) == [1 2 1]
    @test copy(argmin(x, dims = 2)) == reshape([1, 2], :, 1)
  end
end

################################################################################
# Run tests
################################################################################
@testset "NDArray Test" begin
  test_constructor()
  test_ones_zeros_like()
  test_assign()
  test_copy()
  test_slice()
  test_linear_idx()
  test_first()
  test_lastindex()
  test_cat()
  test_plus()
  test_minus()
  test_mul()
  test_div()
  test_rdiv()
  test_mod()
  test_gd()
  test_saveload()
  test_clip()
  test_power()
  test_sqrt()
  test_eltype()
  test_nd_as_jl()
  test_dot()
  test_reshape()
  test_expand_dims()
  test_sum()
  test_mean()
  test_maximum()
  test_minimum()
  test_prod()
  test_fill()
  test_transpose()
  test_show()
  test_size()
  test_trigonometric()
  test_hyperbolic()
  test_act_funcs()
  test_equal()
  test_broadcast_to()
  test_broadcast_axis()
  test_hypot()
  test_argmax()
  test_argmin()
end

end

ext: jl

max_stars_repo_name: HyperZealot/mxnet

max_issues_repo_name: HyperZealot/mxnet

max_forks_repo_name: HyperZealot/mxnet


text:
module Loom

using HDF5: HDF5

struct Connection
    file::HDF5.File
end

function Base.getproperty(conn::Connection, name::Symbol)
    if name == :ca
        return ColumnAttributes(conn)
    elseif name == :ra
        return RowAttributes(conn)
    elseif name == :matrix
        return Dataset(conn, "/matrix")
    else
        return getfield(conn, name)
    end
end

function Base.show(io::IO, conn::Connection)
    print(io, summary(conn))
    if get(io, :compact, false)
        print(io, "(<filename=", repr(basename(HDF5.filename(conn.file))), ">)")
    else
        m, n = matrixsize(conn)
        print(
            io,
            summary(conn), ':', '\n',
            "  filename: ", basename(HDF5.filename(conn.file)), '\n',
            "  #genes: ", m, '\n',
            "  #cells: ", n, '\n',
            "  column attributes: ", join(colattrs(conn), ", "), '\n',
            "  row attributes: ", join(rowattrs(conn), ", "),
        )
    end
end

function matrixsize(conn::Connection)
    n, m = size(conn.file["/matrix"])
    return m, n
end

ngenes(conn::Connection) = matrixsize(conn)[1]
ncells(conn::Connection) = matrixsize(conn)[2]

colattrs(conn::Connection) = names(conn.file["/col_attrs"])
rowattrs(conn::Connection) = names(conn.file["/row_attrs"])

open(filepath::AbstractString; mode::AbstractString="r") = Connection(HDF5.h5open(filepath, mode))
Base.close(conn::Connection) = close(conn.file)

readmatrix(conn::Connection) = Matrix(read(conn.file["/matrix"])')

# Dataset
# -------

struct Dataset
    conn::Connection
    path::String
end

const Slice = Union{Colon,AbstractRange{<:Integer}}

# open(filepath::AbstractString; mode::AbstractString="r") = Dataset(Connection(HDF5.h5open(filepath, mode)), filepath)

Base.getindex(dataset::Dataset, i::Integer, j::Integer) = dataset.conn.file[dataset.path][j,i]
Base.getindex(dataset::Dataset, i::Integer, j::Slice) = vec(dataset.conn.file[dataset.path][j,i])
Base.getindex(dataset::Dataset, i::Slice, j::Integer) = vec(dataset.conn.file[dataset.path][j,i])
Base.getindex(dataset::Dataset, i::Slice, j::Slice) = Matrix(dataset.conn.file[dataset.path][j,i]')

# Column Attributes
# -----------------

struct ColumnAttributes
    conn::Connection
end

function Base.show(io::IO, ca::ColumnAttributes)
    print(io, summary(ca), ": ", join(colattrs(ca.conn), ", "))
end

Base.getindex(ca::ColumnAttributes, name::AbstractString) = read(ca.conn.file["/col_attrs/$(name)"])

# Row Attributes
# --------------

struct RowAttributes
    conn::Connection
end

function Base.show(io::IO, ra::RowAttributes)
    print(io, summary(ra), ": ", join(rowattrs(ra.conn), ", "))
end

Base.getindex(ra::RowAttributes, name::AbstractString) = read(ra.conn.file["/row_attrs/$(name)"])

end # module


text:
function run_thermal_bubble_regression_test(arch)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100
    Δt = 6

    grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz))
    closure = ConstantIsotropicDiffusivity(ν=4e-2, κ=4e-2)
    model = Model(architecture=arch, grid=grid, closure=closure, coriolis=FPlane(f=1e-4))

    model.tracers.T.data.parent .= 9.85
    model.tracers.S.data.parent .= 35.0

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
    model.tracers.T.data[i1:i2, j1:j2, k1+1:k2+1] .+= 0.01

    # Uncomment to generate regression data.
    # outputs = Dict("v" => model.velocities.v,
    #                "u" => model.velocities.u,
    #                "w" => model.velocities.w,
    #                "T" => model.tracers.T,
    #                "S" => model.tracers.S)
    #
    # nc_writer = NetCDFOutputWriter(model, outputs,
    #                                filename="data/thermal_bubble_regression_10.nc",
    #                                frequency=10)
    # push!(model.output_writers, nc_writer)

    time_step!(model, 10, Δt)

    regression_data_filepath = joinpath(dirname(@__FILE__), "data", "thermal_bubble_regression_10.nc")
    ds = Dataset(regression_data_filepath, "r")

    u = ds["u"][:]
    v = ds["v"][:]
    w = ds["w"][:]
    T = ds["T"][:]
    S = ds["S"][:]

    u = reverse(u; dims=3)
    v = reverse(v; dims=3)
    T = reverse(T; dims=3)
    S = reverse(S; dims=3)

    w = reverse(w[:, :, 2:Nz]; dims=3)
    w = cat(zeros(Nx, Ny), w; dims=3)

    field_names = ["u", "v", "w", "T", "S"]
    fields = [model.velocities.u, model.velocities.v, model.velocities.w, model.tracers.T, model.tracers.S]
    fields_gm = [u, v, w, T, S]
    for (field_name, φ, φ_gm) in zip(field_names, fields, fields_gm)
        φ_min = minimum(Array(interior(φ)) - φ_gm)
        φ_max = maximum(Array(interior(φ)) - φ_gm)
        φ_mean = mean(Array(interior(φ)) - φ_gm)
        φ_abs_mean = mean(abs.(Array(interior(φ)) - φ_gm))
        φ_std = std(Array(interior(φ)) - φ_gm)
        @info(@sprintf("Δ%s: min=%.6g, max=%.6g, mean=%.6g, absmean=%.6g, std=%.6g\n",
                       field_name, φ_min, φ_max, φ_mean, φ_abs_mean, φ_std))
    end

    # Now test that the model state matches the regression output.
    @test all(Array(interior(model.velocities.u)) .≈ u)
    @test all(Array(interior(model.velocities.v)) .≈ v)
    @test all(Array(interior(model.velocities.w))[:, :, 2:Nz] .≈ w[:, :, 2:Nz])
    @test all(Array(interior(model.tracers.T))    .≈ T)
    @test all(Array(interior(model.tracers.S))    .≈ S)
end


text:
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))

datasets = Vector(undef, 6)

X = Data.roessler_chaotic(; N = 1e4)
εs = MathConstants.e .^ (-8:0.2:5)

Hs = genentropy.(Ref(X), εs)

x = log.(εs)
y = - Hs
figure()
plot(x, y)

lrs, slopes = linear_regions(x, y; tol = 0.25)

plot(x, y)
for i in 1:length(lrs)-1
    plot(
        x[lrs[i]:lrs[i+1]], y[lrs[i]:lrs[i+1]];
        marker = "o", ms = 5
    )
end

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
# =========================================================================
#    BASIC STRUCTURES
# =========================================================================

struct vtx
    #IDENTIFIERS
    name::String
    type::String

    #VALS
    x::Float64 #xcoor
    y::Float64 #ycoor
    MAX::Float64 #max inventory level
    MIN::Float64 #min inventory level
    START::Float64 #starting inventory level

    #COSTS
    h::Float64 #inv cost per unit
end

struct veh
    #IDENTIFIERS
    name::String
    type::String

    #CHAR
    cover::Vector{Int64}
    BP::Dict{Int64,Int64}
    Q::Int64

    #COSTS
    vx::Float64
    vl::Float64
    fp::Float64
    fd::Float64
end

struct col
    #quantity
    u::JuMP.Containers.DenseAxisArray
    v::JuMP.Containers.DenseAxisArray
    l::JuMP.Containers.DenseAxisArray

    #decision
    o::JuMP.Containers.DenseAxisArray
    x::JuMP.Containers.DenseAxisArray
end

struct dv
    #three master constraints
    λ::JuMP.Containers.DenseAxisArray
    δ::JuMP.Containers.SparseAxisArray
    ϵ::JuMP.Containers.DenseAxisArray

    #bounding constraints
    ρ::JuMP.Containers.DenseAxisArray
    σ::JuMP.Containers.DenseAxisArray
end

struct β
    i::Int64
    v::Int64
end

struct S
    k::Int64
    t::Int64
    seq::Vector{β}
end

struct bound
    S::S #component bound sequence
    sense::String #≳ or ≲
    κ::Int64 #value of aggregation
end

struct stabilizer
    slCoeff::Float64
    suCoeff::Float64
    slLim::Dict
    suLim::Dict
end

struct node
    #IDENTIFIER
    parent::UUID
    self::UUID

    #Dynamic SET
    bounds::Vector{bound}
    columns::Vector{Dict}

    #SUPPORT
    stab::stabilizer

    #STATUS
    status::Vector{String}
end


text:
# ============= accuracy tests ================
print("transforms: accuracy tests ...\n")

# test against data made by data/make_filter_data.m in Octave
wname = Any["Daubechies","Coiflet","Haar","Symmlet","Battle","Vaidyanathan","Beylkin"];
wnum = Any[[4:2:20],[2:5],[0],[4:10],[1,3,5],[0],[0]];
wvm = Any[[2:1:10],[4:2:10],[0],[4:10],[2,4,6],[0],[0]];

name = "filter"
data = vec(readdlm(joinpath(dirname(@__FILE__), "data", "filter1d_data.txt"),'\t'))
data2 = readdlm(joinpath(dirname(@__FILE__), "data", "filter2d_data.txt"),'\t')

stderr = 1e-9*sqrt(length(data))
stderr2 = 1e-9*sqrt(length(data2))

for i = 1:length(wname)
    for num = 1:length(wnum[i])
        # expected transforms
        ye = vec(readdlm(joinpath(dirname(@__FILE__), "data", string(name,"1d_",wname[i],wnum[i][num],".txt")),'\t'))
        ye2 = readdlm(joinpath(dirname(@__FILE__), "data", string(name,"2d_",wname[i],wnum[i][num],".txt")),'\t')

        wn = wnum[i][num]
        wn!=0 && (wt = waveletfilter(wname[i],wvm[i][num]))
        wn==0 && (wt = waveletfilter(lowercase(wname[i][1:4])))
        # transform data
        y = dwt(data, wt)
        y2 = dwt(data2, wt)

        @test_vecnorm_eq_eps y ye stderr
        @test_vecnorm_eq_eps y2 ye2 stderr2
        # a few bad cases have quite large norms
        if wname[i] != "Battle" && (wname[i] != "Coiflet" && wvm[i][num]==10)
            @test abs(vecnorm(data)-vecnorm(y)) < 1e-9
            @test abs(vecnorm(data2)-vecnorm(y2)) < 1e-9
            @test_vecnorm_eq_eps idwt(y,wt) data stderr*100
            @test_vecnorm_eq_eps idwt(y2,wt) data2 stderr2*100
        end
    end
end

# 1-d and 2-d lifting and filtering comparison
for WT in ("db1","db2")
    wf = waveletfilter(WT)
    wls = GLS(WT)
    n = 64
    x = randn(n)
    stderr = 1e-10*sqrt(n)
    stderr2 = 1e-10*sqrt(n*n)

    x2 = copy(x)
    tmp = zeros(n)

    for L in (nscales(n),0,1,2)
        yf = dwt(x, wf, L)
        yls = dwt(x, wls, L)
        dwt!(x2, wls, L, true, tmp)

        @test_vecnorm_eq_eps yf yls stderr
        @test_vecnorm_eq_eps yf x2 stderr

        ytf = idwt(yf, wf, L)
        ytls = idwt(yls, wls, L)
        dwt!(x2, wls, L, false, tmp)

        @test_vecnorm_eq_eps ytf x stderr
        @test_vecnorm_eq_eps ytls x stderr
        @test_vecnorm_eq_eps x2 x stderr

    end

    x = randn(n,n)
    x2 = copy(x)
    for L in (nscales(n),0,1,2)
        yf = dwt(x, wf, L)
        yls = dwt(x, wls, L)

        @test_vecnorm_eq_eps yf yls stderr2
        @test_approx_eq yf yls

        ytf = idwt(yf, wf, L)
        ytls = idwt(yls, wls, L)

        @test_vecnorm_eq_eps ytf x stderr2
        @test_vecnorm_eq_eps ytls x stderr2

    end
end

# ============= transform functionality ================
print("transforms: transform functionality ...\n")

# column-wise 1-d
wf = waveletfilter("db2")
x = randn(16,2)
y = copy(x)
y[:,1] = dwt(vec(x[:,1]),wf)
y[:,2] = dwt(vec(x[:,2]),wf)
@test_approx_eq dwtc(x,wf) y
@test_approx_eq dwtc(x',wf,nscales(16),1)' y

# column-wise 2-d
n = 16
x = randn(n,n,2)
y = copy(x)
y[:,:,1] = dwt(reshape(x[:,:,1],n,n),wf)
y[:,:,2] = dwt(reshape(x[:,:,2],n,n),wf)
@test_approx_eq dwtc(x,wf) y
x = randn(n,2,n)
y = copy(x)
y[:,1,:] = dwt(reshape(x[:,1,:],n,n),wf)
y[:,2,:] = dwt(reshape(x[:,2,:],n,n),wf)
@test_approx_eq dwtc(x,wf,nscales(16),2) y

# "inplace" for filter
wf = waveletfilter("db2")
x = randn(16)
@test_approx_eq dwt(x,wf,2) dwt!(copy(x),wf,2,true)

# "out of place" for LS
wt = GLS("db2")
x = randn(16)
@test_approx_eq dwt(x,wt,2) dwt!(similar(x),x,wt,2,true)

# ============= types and sizes ================
print("transforms: types and sizes ...\n")

function makedwt(ft::Type, n, wf, L)
	x0 = rand(-5:5, n)
	x = zeros(ft, n)
	copy!(x,x0)
	return (x, dwt(x, wf, L))
end

# 1-d
n = 8; wf = waveletfilter("db2"); L = 2
sett = (n,wf,L)

ft = Float64; x, y = makedwt(ft, sett...)
@test Array{ft,1} == typeof(y) && length(y) == n
ft = Float32; x, y = makedwt(ft, sett...)
@test Array{ft,1} == typeof(y) && length(y) == n
ft = Int64; x, y = makedwt(ft, sett...)
@test Array{typeof(float(x[1])),1} == typeof(y) && length(y) == n
ft = Int32; x, y = makedwt(ft, sett...)
@test Array{typeof(float(x[1])),1} == typeof(y) && length(y) == n
@test_approx_eq dwt(x,wf) dwt(float(x),wf)

# 2-d
sett = ((n,n),wf,L)
ft = Float64; x, y = makedwt(ft, sett...)
@test Array{ft,2} == typeof(y) && length(y) == n*n
ft = Float32; x, y = makedwt(ft, sett...)
@test Array{ft,2} == typeof(y) && length(y) == n*n
ft = Int64; x, y = makedwt(ft, sett...)
@test Array{typeof(float(x[1])),2} == typeof(y) && length(y) == n*n
ft = Int32; x, y = makedwt(ft, sett...)
@test Array{typeof(float(x[1])),2} == typeof(y) && length(y) == n*n
@test_approx_eq dwt(x,wf) dwt(float(x),wf)

# non-Array type
wt = waveletls("db2")
x = randn(16, 16)
xs = sub(copy(x), 1:16, 1:16)
@test_approx_eq dwt(x,wt,2) dwt(xs,wt,2)

# ============= error tests ================
print("transforms: error tests ...\n")

type wunknownt <: DiscreteWavelet end
uwt = wunknownt()
EE = Exception
@test_throws EE dwt(randn(4),uwt)
@test_throws EE dwt(randn(4,4),uwt)
@test_throws EE waveletfilter("db2asdsad")
@test_throws EE waveletfilter("db2", "ppppp")

# ============= WPT ================
print("transforms: WPT ...\n")

wf = waveletfilter("db2")
x = randn(16)

L = 1
wp = wpt(x,wf,L)
dw = dwt(x,wf,L)
@test_approx_eq wp dw
@test_approx_eq iwpt(wp,wf,L) x

L = 2
wp = wpt(x,wf,L)
dw = dwt(x,wf,L)
dw2 = copy(dw)
dw2[9:end] = dwt(dw[9:end],wf,1)
@test_approx_eq dw[1:8] wp[1:8]
@test_approx_eq dw2 wp
@test_approx_eq iwpt(wp,wf,L) x

L = 3
wp = wpt(x,wf,L)
dw = dwt(x,wf,L)
@test_approx_eq dw[1:4] wp[1:4]
@test_approx_eq dwt(dw2[5:8],wf,1) wp[5:8]
@test_approx_eq dwt(dw2[9:12],wf,1) wp[9:12]
@test_approx_eq dwt(dw2[13:16],wf,1) wp[13:16]
@test_approx_eq iwpt(wp,wf,L) x

wl = waveletls("db2")
x = randn(128)
@test_approx_eq iwpt(wpt(x,wf),wf) x
@test_approx_eq iwpt(wpt(x,wl),wl) x
@test_approx_eq wpt(x,wl,1) wpt(x,wf,1)
@test_approx_eq wpt(x,wl,2) wpt(x,wf,2)
@test_approx_eq wpt(x,wl,4) wpt(x,wf,4)
@test_approx_eq wpt(x,wl) wpt(x,wf)

wl = waveletls("db2")
n = 128
x = randn(n)
for L = 0:nscales(n)
    @test_approx_eq wpt(x, wl, maketree(n, L, :dwt)) dwt(x, wl, L)
    @test_approx_eq iwpt(x, wl, maketree(n, L, :dwt)) idwt(x, wl, L)
    @test_approx_eq wpt(x, wf, maketree(n, L, :dwt)) dwt(x, wf, L)
    @test_approx_eq iwpt(x, wf, maketree(n, L, :dwt)) idwt(x, wf, L)
end

# ============= tranform low level functions ================

#...


text: @test false broken=true


text:
abstract type DynamicComponent <: DeviceParameter end

"""
Abstract type for all dynamic injection types
"""
abstract type DynamicInjection <: Device end

supports_services(d::DynamicInjection) = false

function get_dynamic_components(device::T) where {T <: DynamicInjection}
    return (
        getfield(device, x) for
        (x, y) in zip(fieldnames(T), fieldtypes(T)) if y <: DynamicComponent
    )
end


text:
@with_kw struct BoxNode <: Node
	lb :: Vector{Float64}
	ub :: Vector{Float64}
	subdivision_level :: Int64 = 0
	parent :: Union{BoxNode, Nothing} = nothing
	sub_boxes :: Vector{BoxNode} = []

    # Test point data
	database :: D where D<:SomeDataBase = DataBase()
    db_indices :: Vector{Int} = [];

end

# required Tree methods (see "Trees.jl")
children( bn :: BoxNode ) = bn.sub_boxes;
parent( bn :: BoxNode ) = bn.parent;
depth( bn :: BoxNode ) = bn.subdivision_level;

# utilities
lb( bn :: BoxNode ) = bn.lb; # use methods for `lb` & `ub` …
ub( bn :: BoxNode ) = bn.ub; # … in case we would rather want to store center & radius
n_vars( bn :: BoxNode ) = length( bn.lb );
width( bn :: BoxNode ) = ub(bn) .- lb(bn); # maybe raname?
center(bn :: BoxNode ) = lb(bn) .+ .5 .* width( bn );
edge_length( bn :: BoxNode, dim :: Int) = ub(bn)[dim] - lb(bn)[dim];

@doc """
    contains(box, point)

Return true if `point` lies in `box`.
"""
function contains( bn :: BoxNode, point :: Vector{R} where R<:Real)
    all( lb(bn) .<= point ) && all( point .<= ub(bn) )
end

function vertices( bn :: BoxNode )
    LB = lb(bn);
    UB = ub(bn);

    other_ind = combinations( eachindex(UB) )
    return [ [LB,];
        [ let v = copy(LB); v[ind] .= UB[ind]; v end for ind in other_ind ]
    ]
end

function edges( bn :: BoxNode )
    verts = vertices( bn );
    return combinations( verts, 2 )
end

@doc "Retrieve sites indexed by `db_indices` from referenced DataBase `db`."
function sites( bn :: BoxNode )
    return sites( bn.database, bn.db_indices);
end

@doc "Retrieve evaluation results indexed by `db_indices` from referenced DataBase `db`."
function results( bn :: BoxNode; include_nothing = true )
    return results( bn.database, bn.db_indices; include_nothing);
end

"Evaluate all sites referenced by field `db_indices`."
function eval!( bn :: BoxNode, f :: F where F<:Function; reevaluate = true)
    for db_id in bn.db_indices
        eval!( bn.database, f, db_id; reevaluate );
    end
end

@doc """
    init_box_tree( lb, ub, db = nothing )

Initialize the root node of a box Tree with variable boundaries
`lb` and `ub`.
If `db` is a database then it is referenced in the node.
"""
function init_box_tree(
        lb :: Vector{R} where R<:Real, ub :: Vector{R} where R<:Real,
        db :: Union{Nothing, DataBase} = nothing
    )
    box_tree = BoxNode(;
        lb = lb,
		ub = ub,
        database = isnothing(db) ? DataBase() : db,
    );
	return box_tree
end

@doc """
    add_sub_box!( parent, sub_lb, sub_ub )

Add a new box node with variable boundaries `sub_lb` and `sub_ub`
to parent.
Modifies `sub_boxes` field of `parent` and sets the parent reference
of the new box correctly.
"""
function add_sub_box!(
        parent :: BoxNode, sub_lb :: Vector{Float64},
        sub_ub :: Vector{Float64}
    )
	# check if boundaries of sub-box lay within parent box
	if all( lb(parent) .<= sub_lb ) && all( sub_ub .<= ub(parent) )
		push!( parent.sub_boxes,
			BoxNode(;
				lb = sub_lb,
				ub = sub_ub,
				subdivision_level = parent.subdivision_level + 1,
				parent = parent,
				sub_boxes = BoxNode[],
                database = parent.database,
			)
		);
	else
		@error("Variable boundaries of subbox must not exceed parent boundaries.")
	end
end

# recursive subdivison function
# usually not called directly by the user
"Divide each of `bn`s edges by divisor or the entries of divisor
and add sub boxes as children to the tree."
function subdivide!(
        bn :: BoxNode, divisor :: Union{Int, Vector{Int}},
        LB :: Vector{Float64}, UB :: Vector{Float64};
    )

    # recursion starts with `LB` and `UB` empty
	var_level = length(LB) + 1; # recursion level

    # divisor for this level
    sub_divisor = isa( divisor, Int ) ? divisor : divisor[ var_level ]

    width = edge_length( bn, var_level );
	width_fraction = width / sub_divisor
	for s = 1 : sub_divisor
		sub_lb = [LB; lb( bn )[ var_level ] + (s-1) * width_fraction ]
		sub_ub = [UB; lb( bn )[ var_level ] + s * width_fraction ]
		if var_level == n_vars( bn )
			add_sub_box!( bn, sub_lb, sub_ub )
		else
			subdivide!( bn, divisor, sub_lb, sub_ub )
		end
	end
end

function subdivide!(
        bn :: BoxNode;
        method = :cycle
    )
    if method == :cycle
        divisor = fill(1, n_vars(bn) );
        divisor[ ( depth(bn) % n_vars(bn) + 1 ) ] = 2;
    else
        divisor = 2;
    end
    subdivide!( bn, divisor, Float64[], Float64[] );
end

@doc """
    subdivide!(parent, divisor=2)

Recursively subdivide node `parent`.
`parent` is divided by `divisor` along each box dimension.
"""
function subdivide!( parent :: BoxNode, divisor :: Union{Vector{Int}, Int} = 2)
	lb = Float64[];
	ub = Float64[];
	subdivide!( parent, divisor, lb, ub )
end

function generate_sites!( bn :: BoxNode; method ::Symbol = :edges_center, kwargs...  )
    generate_sites!(bn , Val(method); kwargs...)
end

function generate_sites!( bn :: BoxNode, ::Val{:edges_center}; num_edge_points = -1)
    if num_edge_points < 0
        num_edge_points = n_vars(bn)
    end

    # reset old database ids
    empty!(bn.db_indices);

    # get points for unit HyperCube (cached for future reference)
    unit_points = unit_sites_edges_center( n_vars(bn); num_edge_points = num_edge_points)

    scale_fn(p) = p .* width(bn) .+ lb(bn);
    new_points = scale_fn.(unit_points);

    db_indices = push_sites!( bn.database, new_points );
    push!(bn.db_indices,db_indices...);
    nothing;
end

@memoize function unit_sites_edges_center(n_vars ::Int; num_edge_points ::Int = -1)
    if num_edge_points < 0
        num_edge_points = n_vars
    end

    bn = init_box_tree( zeros(n_vars), ones(n_vars) );

    ℓ = collect( 1 : num_edge_points )
    steps = ( 2 .* ℓ .- 1 ) ./ (2*num_edge_points);

    # N * num_edges + 1
    box_edges = edges(bn)
    num_edges = length(edges(bn))
    num_new_points = num_edge_points * num_edges + 1;

    # pre-allocate site array
    new_points = Vector{Vector{Float64}}(undef, num_new_points);

    new_points[end] = center(bn);

    counter = 1;
    for e ∈ box_edges
        v1, v2 = e;
        for λ ∈ steps
            new_points[counter] = λ .* v1 .+ (1-λ) .* v2;
            counter += 1;
        end
    end

    return new_points
end


text:
# An implementation of a FISTA-like method, where the smooth part of the objective function can be strongly convex.

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalCore: Zero
using LinearAlgebra
using Printf

"""
    SFISTAIteration(; <keyword-arguments>)

Iterator implementing the FISTA-like algorithm in [3].

This iterator solves strongly convex composite optimization problems of the form

    minimize f(x) + g(x),

where g is proper closed convex and f is a continuously differentiable function that is `mf`-strongly convex and whose gradient is
`Lf`-Lipschitz continuous.

The scheme is based on Nesterov's accelerated gradient method [1, Eq. (4.9)] and Beck's method for the convex case [2]. Its full
definition is given in [3, Algorithm 2.2.2.], and some analyses of this method are given in [3, 4, 5]. Another perspective is that
it is a special instance of [4, Algorithm 1] in which μh=0.

See also: [`SFISTA`](@ref).

# Arguments
- `x0`: initial point; must be in the domain of g.
- `f=Zero()`: smooth objective term.
- `g=Zero()`: proximable objective term.
- `mf=0` : strong convexity constant of f (see above).
- `Lf` : Lipschitz constant of ∇f (see above).

# References
1. Nesterov, Y. (2013). Gradient methods for minimizing composite functions. Mathematical Programming, 140(1), 125-161.
2. Beck, A., & Teboulle, M. (2009). A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM journal on imaging sciences, 2(1), 183-202.
3. Kong, W. (2021). Accelerated Inexact First-Order Methods for Solving Nonconvex Composite Optimization Problems. arXiv preprint arXiv:2104.09685.
4. Kong, W., Melo, J. G., & Monteiro, R. D. (2021). FISTA and Extensions - Review and New Insights. arXiv preprint arXiv:2107.01267.
5. Florea, M. I. (2018). Constructing Accelerated Algorithms for Large-scale Optimization-Framework, Algorithms, and Applications.
"""
Base.@kwdef struct SFISTAIteration{R,C<:Union{R,Complex{R}},Tx<:AbstractArray{C},Tf,Th}
    x0::Tx
    f::Tf = Zero()
    g::Th = Zero()
    Lf::R
    mf::R = real(eltype(Lf))(0.0)
end

Base.IteratorSize(::Type{<:SFISTAIteration}) = Base.IsInfinite()

Base.@kwdef mutable struct SFISTAState{R, Tx}
    λ::R                                # stepsize.
    yPrev::Tx                           # previous main iterate.
    y:: Tx = zero(yPrev)                # main iterate.
    xPrev::Tx = copy(yPrev)             # previous auxiliary iterate.
    x::Tx  = zero(yPrev)                # auxiliary iterate (see [3]).
    xt::Tx = zero(yPrev)                # prox center used to generate main iterate y.
    τ::R = real(eltype(yPrev))(1.0)     # helper variable (see [3]).
    a::R = real(eltype(yPrev))(0.0)     # helper variable (see [3]).
    APrev::R = real(eltype(yPrev))(1.0) # previous A (helper variable).
    A::R = real(eltype(yPrev))(0.0)     # helper variable (see [3]).
    gradf_xt::Tx = zero(yPrev)          # array containing ∇f(xt).
end

function Base.iterate(
    iter::SFISTAIteration,
    state::SFISTAState = SFISTAState(λ=1/iter.Lf, yPrev=copy(iter.x0))
)
    # Set up helper variables.
    state.τ = state.λ * (1 + iter.mf * state.APrev)
    state.a = (state.τ + sqrt(state.τ ^ 2 + 4 * state.τ * state.APrev)) / 2
    state.A = state.APrev + state.a
    state.xt .= (state.APrev / state.A) .* state.yPrev + (state.a / state.A) .* state.xPrev
    gradient!(state.gradf_xt, iter.f, state.xt)
    λ2 = state.λ / (1 + state.λ * iter.mf)
    # FISTA acceleration steps.
    prox!(state.y, iter.g, state.xt - λ2 * state.gradf_xt, λ2)
    state.x .= state.xPrev .+ (state.a / (1 + state.A * iter.mf)) .* ((state.y .- state.xt) ./ state.λ .+ iter.mf .* (state.y .- state.xPrev))
    # Update state variables.
    state.yPrev .= state.y
    state.xPrev .= state.x
    state.APrev = state.A
    return state, state
end

# Different stopping conditions (sc). Returns the current residual value and whether or not a stopping condition holds.
function check_sc(state::SFISTAState, iter::SFISTAIteration, tol, termination_type)
    if termination_type == "AIPP"
        # AIPP-style termination [4]. The main inclusion is: r ∈ ∂_η(f + h)(y).
        r = (iter.y0 - state.x) / state.A
        η = (norm(iter.y0 - state.y) ^ 2 - norm(state.x - state.y) ^ 2) / (2 * state.A)
        res = (norm(r) ^ 2 + max(η, 0.0)) / max(norm(iter.y0 - state.y + r) ^ 2, 1e-16)
    else
        # Classic (approximate) first-order stationary point [4]. The main inclusion is: r ∈ ∇f(y) + ∂h(y).
        λ2 = state.λ / (1 + state.λ * iter.mf)
        gradf_y, = gradient(iter.f, state.y)
        r = gradf_y - state.gradf_xt + (state.xt - state.y) / λ2
        res = norm(r)
    end
    return res, (res <= tol || res ≈ tol)
end

default_solution(::SFISTAIteration, state::SFISTAState) = state.y

"""
    SFISTA(; <keyword-arguments>)

Constructs the the FISTA-like algorithm in [3].

This algorithm solves strongly convex composite optimization problems of the form

    minimize f(x) + h(x),

where h is proper closed convex and f is a continuously differentiable function that is μ-strongly convex and whose gradient is
Lf-Lipschitz continuous.

The scheme is based on Nesterov's accelerated gradient method [1, Eq. (4.9)] and Beck's method for the convex case [2]. Its full
definition is given in [3, Algorithm 2.2.2.], and some analyses of this method are given in [3, 4, 5]. Another perspective is that
it is a special instance of [4, Algorithm 1] in which μh=0.

The returned object has type `IterativeAlgorithm{SFISTAIteration}`,
and can be called with the problem's arguments to trigger its solution.

See also: [`SFISTAIteration`](@ref), [`IterativeAlgorithm`](@ref).

# Arguments
- `maxit::Int=10_000`: maximum number of iteration
- `tol::1e-6`: tolerance for the default stopping criterion
- `stop::Function`: termination condition, `stop(::T, state)` should return `true` when to stop the iteration
- `solution::Function`: solution mapping, `solution(::T, state)` should return the identified solution
- `verbose::Bool=false`: whether the algorithm state should be displayed
- `freq::Int=100`: every how many iterations to display the algorithm state
- `display::Function`: display function, `display(::Int, ::T, state)` should display a summary of the iteration state
- `kwargs...`: additional keyword arguments to pass on to the `SFISTAIteration` constructor upon call

# References
1. Nesterov, Y. (2013). Gradient methods for minimizing composite functions. Mathematical Programming, 140(1), 125-161.
2. Beck, A., & Teboulle, M. (2009). A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM journal on imaging sciences, 2(1), 183-202.
3. Kong, W. (2021). Accelerated Inexact First-Order Methods for Solving Nonconvex Composite Optimization Problems. arXiv preprint arXiv:2104.09685.
4. Kong, W., Melo, J. G., & Monteiro, R. D. (2021). FISTA and Extensions - Review and New Insights. arXiv preprint arXiv:2107.01267.
5. Florea, M. I. (2018). Constructing Accelerated Algorithms for Large-scale Optimization-Framework, Algorithms, and Applications.
"""
SFISTA(;
    maxit=10_000,
    tol=1e-6,
    termination_type="",
    stop=(iter, state) -> check_sc(state, iter, tol, termination_type)[2],
    solution=default_solution,
    verbose=false,
    freq=100,
    display=(it, iter, state) -> @printf("%5d | %.3e\n", it, check_sc(state, iter, tol, termination_type)[1]),
    kwargs...
) = IterativeAlgorithm(SFISTAIteration; maxit, stop, solution, verbose, freq, display, kwargs...)


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

JLLWrappers.@generate_main_file_header("Mesa")
JLLWrappers.@generate_main_file("Mesa", UUID("78dcde23-ec64-5e07-a917-6fe22bbc0f45"))
end  # module Mesa_jll


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

using Markdown
using InteractiveUtils

# ╔═╡ 62b0cf82-e00c-498b-a6d8-124ea9950ef4
using Manopt, Manifolds, Random, BenchmarkTools

# ╔═╡ d5e636f0-533a-11ec-168a-e118d53012cb
md"""
# Illustration how to use mutating gradient functions

When it comes to time critital operations, a main ingredient in Julia are
mutating functions, i.e. those that compute in place without additional Memory
allocations. In the following the illustrate how to do this with `Manopt.jl`.

Let's start with the same function as in [Get Started: Optimize!](https://manoptjl.org/stable/tutorials/MeanAndMedian.html)
and compute the mean of some points. Just that here we use the sphere ``\mathbb S^{30}``
and `n=800` points.

From the just mentioned example, the implementation looks like
"""

# ╔═╡ 967e5caa-1d27-4c24-a882-6f7126053754
begin
    Random.seed!(42)
    m = 30
    M = Sphere(m)
    n = 800
    σ = π / 8
    x = zeros(Float64, m + 1)
    x[2] = 1.0
    data = [exp(M, x, random_tangent(M, x, Val(:Gaussian), σ)) for i in 1:n]
end;

# ╔═╡ 1436cdc7-1e84-4438-937d-2211856348de
md"""
## Classical definition

The variant from the previous tutorial defines a cost ``F(x)`` and its gradient ``gradF(x)``
"""

# ╔═╡ 6f1b04cf-2a76-4660-916c-3b8b8472d3c0
F(M, x) = sum(1 / (2 * n) * distance.(Ref(M), Ref(x), data) .^ 2)

# ╔═╡ e487f38c-334e-4819-b775-a9a44632d7ff
gradF(M, x) = sum(1 / n * grad_distance.(Ref(M), data, Ref(x)))

# ╔═╡ e4a7b318-65e7-40bb-899d-8bdeeefbf304
md"""
we further set the stopping criterion to be a little more strict, then we obtain
"""

# ╔═╡ 5a2d6d9f-20dc-41a2-bfce-32509a1ef106
begin
    sc = StopWhenGradientNormLess(1e-10)
    x0 = random_point(M)
    m1 = gradient_descent(M, F, gradF, x0; stopping_criterion=sc)
    @benchmark gradient_descent($M, $F, $gradF, $x0; stopping_criterion=$sc)
end

# ╔═╡ 306c7706-f399-4104-9159-c2e6b0bca189
md"""
## Inplace computation of the gradient

We can reduce the memory allocations, by implementing the gradient as a [functor](https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects).
The motivation is twofold: On the one hand, we want to avoid variables from global scope,
for example the manifold `M` or the `data`, to be used within the function
For more complicated cost functions it might also be worth considering to do the same.

Here we store the data (as reference) and one temporary memory in order to avoid
reallocation of memory per `grad_distance` computation. We have
"""

# ╔═╡ 4e464076-d477-4e38-aac6-79639c1cb9b5
begin
    struct grad!{TD,TTMP}
        data::TD
        tmp::TTMP
    end
    function (gradf!::grad!)(M, X, x)
        fill!(X, 0)
        for di in gradf!.data
            grad_distance!(M, gradf!.tmp, di, x)
            X .+= gradf!.tmp
        end
        X ./= length(gradf!.data)
        return X
    end
end

# ╔═╡ 15d4c49c-7178-4be1-8e83-00900e384831
md"""
Then we just have to initialize the gradient and perform our final benchmark.
Note that we also have to interpolate all variables passed to the benchmark with a `$`.
"""

# ╔═╡ 90f40202-69ec-4dba-b1d0-905c916a3642
begin
    gradF2! = grad!(data, similar(data[1]))
    m2 = deepcopy(x0)
    gradient_descent!(
        M, F, gradF2!, m2; evaluation=MutatingEvaluation(), stopping_criterion=sc
    )
    @benchmark gradient_descent!(
        $M, $F, $gradF2!, m2; evaluation=$(MutatingEvaluation()), stopping_criterion=$sc
    ) setup = (m2 = deepcopy($x0))
end

# ╔═╡ 3e01469b-670b-47a6-a379-199b70c59207
md"""Mote that the results `m1`and `m2` are of course (approximately) the same."""

# ╔═╡ 10c7b9e5-024c-4815-be58-c524a7a08d25
distance(M, m1, m2)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
Manopt = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[compat]
BenchmarkTools = "~1.2.0"
Manifolds = "~0.7.1"
Manopt = "~0.3.14"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractFFTs]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.0.1"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "62e51b39331de8911e4a7ff6f5aaf38a5f4cc0ae"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.2.0"

[[deps.ArrayInterface]]
deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "265b06e2b1f6a216e0e8f183d28e4d354eab3220"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.2.1"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
git-tree-sha1 = "61adeb0823084487000600ef8b1c00cc2474cd47"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.2.0"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "f885e7e7c124f8c92650d61b9477b9ac2ee607dd"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.11.1"

[[deps.ChangesOfVariables]]
deps = ["LinearAlgebra", "Test"]
git-tree-sha1 = "9a1d594397670492219635b35a3d830b04730d62"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.1"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "a851fec56cb73cfdf43762999ec72eff5b86882a"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.15.0"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dce3e3fea680869eaa0b774b2e8343e9ff442313"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.40.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.CovarianceEstimation]]
deps = ["LinearAlgebra", "Statistics", "StatsBase"]
git-tree-sha1 = "bc3930158d2be029e90b7c40d1371c4f54fa04db"
uuid = "587fd27a-f159-11e8-2dae-1979310e6154"
version = "0.2.6"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "7f3bec11f4bcd01bc1f507ebce5eadf1b0a78f47"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.34"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.Einsum]]
deps = ["Compat"]
git-tree-sha1 = "4a6b3eee0161c89700b6c1949feae8b851da5494"
uuid = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0"
version = "0.4.1"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "8756f9935b7ccc9064c6eef0bff0ad643df733a3"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.12.7"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Graphs]]
deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "92243c07e786ea3458532e199eb3feee0e7e08eb"
uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
version = "1.4.1"

[[deps.HybridArrays]]
deps = ["LinearAlgebra", "Requires", "StaticArrays"]
git-tree-sha1 = "c1d5b1dcdf2140644e1c6beb9ca09fbed601c241"
uuid = "1baab800-613f-4b0a-84e4-9cd3431bfbb9"
version = "0.4.9"

[[deps.IfElse]]
git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.1"

[[deps.Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.2"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.Kronecker]]
deps = ["LinearAlgebra", "NamedDims", "SparseArrays", "StatsBase"]
git-tree-sha1 = "a51f46415c844dee694cb8b20a3fcbe6dba342c2"
uuid = "2c470bb0-bcc8-11e8-3dad-c9649493f05e"
version = "0.5.1"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "be9eef9f9d78cecb6f262f3c10da151a6c5ab827"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.5"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Manifolds]]
deps = ["Colors", "Distributions", "Einsum", "Graphs", "HybridArrays", "Kronecker", "LinearAlgebra", "ManifoldsBase", "Markdown", "Random", "RecipesBase", "RecursiveArrayTools", "Requires", "SimpleWeightedGraphs", "SpecialFunctions", "StaticArrays", "Statistics", "StatsBase"]
git-tree-sha1 = "3811935fd2549b0f5f9b365d6e7173fcdb98ec9a"
uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
version = "0.7.1"

[[deps.ManifoldsBase]]
deps = ["LinearAlgebra", "Markdown"]
git-tree-sha1 = "77a5949567437d185ee929c405e3c6c0768118ea"
uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
version = "0.12.9"

[[deps.Manopt]]
deps = ["ColorSchemes", "ColorTypes", "Colors", "DataStructures", "Dates", "LinearAlgebra", "ManifoldsBase", "Markdown", "Random", "Requires", "SparseArrays", "StaticArrays", "Statistics", "Test"]
git-tree-sha1 = "e87c818e19f79444fff1c77880cdb43b51974dfe"
uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
version = "0.3.14"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NamedDims]]
deps = ["AbstractFFTs", "ChainRulesCore", "CovarianceEstimation", "LinearAlgebra", "Pkg", "Requires", "Statistics"]
git-tree-sha1 = "4f34e3ff2fa7f2a1c03fb2c4fdd637380d760bbc"
uuid = "356022a1-0364-5f58-8944-0da4b18d706f"
version = "0.2.42"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "ee26b350276c51697c9c2d88a072b339f9f03d73"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.5"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "ae4bbcadb2906ccc085cf52ac286dc1377dceccc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.2"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.Profile]]
deps = ["Printf"]
uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.RecursiveArrayTools]]
deps = ["ArrayInterface", "ChainRulesCore", "DocStringExtensions", "FillArrays", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
git-tree-sha1 = "c944fa4adbb47be43376359811c0a14757bdc8a8"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.20.0"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[deps.SimpleWeightedGraphs]]
deps = ["Graphs", "LinearAlgebra", "Markdown", "SparseArrays", "Test"]
git-tree-sha1 = "a6f404cc44d3d3b28c793ec0eb59af709d827e4e"
uuid = "47aef6b3-ad0c-573a-a1e2-d07658019622"
version = "1.2.1"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "f0bccf98e16759818ffc5d97ac3ebf87eb950150"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.8.1"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "e7bc80dc93f50857a5d1e3c8121495852f407e6a"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.4.0"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StatsAPI]]
git-tree-sha1 = "0f2aa8e32d511f758a2ce49208181f7733a0936a"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.1.0"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "2bb0cb32026a66037360606510fca5984ccc6b75"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.13"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "bedb3e17cc1d94ce0e6e66d3afa47157978ba404"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.14"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.2"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─d5e636f0-533a-11ec-168a-e118d53012cb
# ╠═62b0cf82-e00c-498b-a6d8-124ea9950ef4
# ╠═967e5caa-1d27-4c24-a882-6f7126053754
# ╟─1436cdc7-1e84-4438-937d-2211856348de
# ╠═6f1b04cf-2a76-4660-916c-3b8b8472d3c0
# ╠═e487f38c-334e-4819-b775-a9a44632d7ff
# ╟─e4a7b318-65e7-40bb-899d-8bdeeefbf304
# ╠═5a2d6d9f-20dc-41a2-bfce-32509a1ef106
# ╟─306c7706-f399-4104-9159-c2e6b0bca189
# ╠═4e464076-d477-4e38-aac6-79639c1cb9b5
# ╟─15d4c49c-7178-4be1-8e83-00900e384831
# ╠═90f40202-69ec-4dba-b1d0-905c916a3642
# ╟─3e01469b-670b-47a6-a379-199b70c59207
# ╠═10c7b9e5-024c-4815-be58-c524a7a08d25
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using JEMSS

println("Loading sim.")
sim = initSim(joinpath(@__DIR__, "input/sim_config.xml"))

println("Simulating.")
simulateReps!(sim)

println("Writing stats files.")
writeStatsFiles(sim, sim.reps)


text:
using MultiScaleArrays
using OrdinaryDiffEq, DiffEqBase, Test, StochasticDiffEq

#=
struct Cell{B} <: AbstractMultiScaleArrayLeaf{B}
    values::Vector{B}
end
struct Population{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end
struct Tissue{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end
struct Embryo{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArrayHead{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end
=#

cell1 = Cell([1.0;2.0;3.0])
cell2 = Cell([4.0;5])
population = construct(Population,deepcopy([cell1,cell2]),[11.0,12.0,13.0]) # Make a Population from cells
cell3 = Cell([3.0;2.0;5.0])
cell4 = Cell([4.0;6])
population2 = construct(Population,deepcopy([cell3,cell4]),[8.0,9.0])

@test population[1] == 1.0
@test population[2] == 2.0
@test population[3] == 3.0
@test population[4] == 4.0
@test population[5] == 5.0
@test population[6] == 11.0
@test population[7] == 12.0
@test population[8] == 13.0

tissue1 = construct(Tissue,deepcopy([population,population2])) # Make a Tissue from Populations
tissue2 = construct(Tissue,deepcopy([population2,population]))
embryo = construct(Embryo,deepcopy([tissue1,tissue2]),[31.0,32.0]) # Make an embryo from Tissues

@test tissue1[1] == 1.0
@test tissue1[2] == 2.0
@test tissue1[3] == 3.0
@test tissue1[4] == 4.0
@test tissue1[5] == 5.0
@test tissue1[6] == 11.0
@test tissue1[7] == 12.0
@test tissue1[8] == 13.0

@test tissue1[9]  == 3.0
@test tissue1[10] == 2.0
@test tissue1[11] == 5.0
@test tissue1[12] == 4.0
@test tissue1[13] == 6.0
@test tissue1[14] == 8.0
@test tissue1[15] == 9.0

@test embryo[1] == 1.0
@test embryo[2] == 2.0
@test embryo[3] == 3.0
@test embryo[4] == 4.0
@test embryo[5] == 5.0
@test embryo[6] == 11.0
@test embryo[7] == 12.0
@test embryo[8] == 13.0

@test embryo[9]  == 3.0
@test embryo[10] == 2.0
@test embryo[11] == 5.0
@test embryo[12] == 4.0
@test embryo[13] == 6.0
@test embryo[14] == 8.0
@test embryo[15] == 9.0

@test embryo[16] == 3.0
@test embryo[17] == 2.0
@test embryo[18] == 5.0
@test embryo[19] == 4.0
@test embryo[20] == 6.0
@test embryo[21] == 8.0
@test embryo[22] == 9.0

@test embryo[23] == 1.0
@test embryo[24] == 2.0
@test embryo[25] == 3.0
@test embryo[26] == 4.0
@test embryo[27] == 5.0
@test embryo[28] == 11.0
@test embryo[29] == 12.0
@test embryo[30] == 13.0

@test embryo[31] == 31.0
@test embryo[32] == 32.0


text:
"""
```
newnote(dsn, job_id, username, note)
```

Create a new note in the `notes` table. A note is associated with a particular
job and therefore a `job_id` is required, along with a `username` to indicate
who wrote the note.
"""
function newnote(dsn, job_id, username, note)
    ODBC.execute!(dsn, """
        INSERT INTO notes (job_id, username, dt, note) VALUES
            ($job_id, '$username', '$(now())', '$note');
    """)
end

"""
```
listnotes(dsn; job_id=0, username="")
```

If keyword arguments `job_id` and/or `username` are provided, list the notes
in the `notes` table that are at the intersection of those conditions.
Otherwise, list all notes in the `notes` table.
"""
function listnotes(dsn; job_id=0, username="")
    if job_id != 0 && username != ""
        wherestr = "WHERE "
        if job_id != 0
            wherestr *= "job_id = $job_id"
        end
        if username != ""
            wherestr *= ", username = '$username'"
        end
    else
        wherestr = ""
    end

    ODBC.query(dsn, """
        SELECT job_id, username, dt, note FROM notes $wherestr;
    """)
end


text:
function mj_model(model; timestep=0.01)
    module_dir()
    jm = jlModel(joinpath(module_dir(), "model/$model"))
    MuJoCo.@set!! jm.opt.timestep=timestep
    jd = jlData(jm)
    mjsim = MJSim(jm, jd)
    return jm, jd, mjsim
end

function angular_momentum(jm, jd)
    MuJoCo.mj_subtreeVel(jm, jd)
    p = jd.subtree_angmom[:,1]
    return p
end

function linear_momentum(jm, jd)
    MuJoCo.mj_subtreeVel(jm, jd)
	mass = robot_mass(jm)
    p = mass * jd.subtree_linvel[:,1]
    return p
end

function robot_mass(jm)
	return sum(jm.body_mass[2:end])
end

function energy(jm, jd)
	MuJoCo.mj_energyPos(jm, jd)
	MuJoCo.mj_energyVel(jm, jd)
	return sum(jd.d.energy)
end