text:
using Distributions, Statistics, PyPlot, Random
rng = MersenneTwister(01012019);

#This defines the tree structure.
mutable struct Tree
    name::String                        # name of the tree
    parent::Vector{Int64}               # parents of nodes in the tree
    children::Vector{Vector{Int64}}     # successor nodes of each parent
    state::Matrix{Float64}              # states of nodes in the tree
    probability::Matrix{Float64}        # probability to go from one node to another

    """
    	child(parent::Vector{Int64})

    Returns a vector of successors of each of the parent nodes.
    """
    function child(parent::Vector{Int64})
        allchildren = Vector{Vector{Int64}}([])
        for node in unique(parent)
            push!(allchildren, [i for i = 1 : length(parent) if parent[i] == node])
        end
        return allchildren
    end
    Tree(name::String,parent::Vector{Int64},
    children::Vector{Vector{Int64}},
    state::Matrix{Float64}, probability::Matrix{Float64})=new(name, parent, child(parent), state, probability)

    """
    	Tree(identifier::Int64)

    Returns some examples of predefined trees.
    - These are (0, 302, 303, 304, 305, 306, 307, 402, 404, 405).
    - You can call any of the above tree and plot to see the properties of the tree.
    """
    function Tree(identifier::Int64)
        self= new()
        if identifier==0
            self.name = "Empty Tree"
            self.parent = zeros(Int64 , 0)
            self.children = child(self.parent)
            self.state = zeros(Float64, 0)'
            self.probability = ones(Float64, 0)'
        elseif identifier == 302
            self.name = "Tree 1x2x2"
            self.parent =[0,1,1,2,2,3,3]
            self.children = child(self.parent)
            self.state = [2.0 2.1 1.9 4.0 1.0 3.0]'
            self.probability = [1.0 0.5 0.5 0.5 0.5 0.5 0.5]'
        elseif identifier == 303
            self.name = "Tree 1x1x4"
            self.parent = [0,1,2,2,2,2]
            self.children = child(self.parent)
            self.state = [3.0 3.0 6.0 4.0 2.0 0.0]'
            self.probability = [1.0 1.0 0.25 0.25 0.25 0.25]'
        elseif identifier == 304
            self.name = "Tree 1x4x1x1"
            self.parent = [0,1,2,0,4,5,0,7,8,0,10,11]
            self.children = child(self.parent)
            self.state = [0.1 2.1 3.0 0.1 1.9 1.0 0.0 -2.9 -1.0 -0.1 -3.1 -4.0]'
            self.probability = [0.14 1.0 1.0 0.06 1.0 1.0 0.48 1.0 1.0 0.32 1.0 1.0]'
        elseif identifier == 305
            self.name = "Tree 1x1x4"
            self.parent = [0,1,2,2,2,2]
            self.children = child(self.parent)
            self.state = [0.0 10.0 28.0 22.0 21.0 20.0]'
            self.probability = [1.0 1.0 0.25 0.25 0.25 0.25]'
        elseif identifier == 306
            self.name = "Tree 1x2x2"
            self.parent = [0,1,1,2,2,3,3]
            self.children = child(self.parent)
            self.state = [0.0 10.0 10.0 28.0 22.0 21.0 20.0]'
            self.probability = [1.0 0.5 0.5 0.5 0.5 0.5 0.5]'
        elseif identifier == 307
            self.name = "Tree 1x4x1"
            self.parent = [0,1,1,1,1,2,3,4,5]
            self.children = child(self.parent)
            self.state = [0.0 10.0 10.0 10.0 10.0 28.0 22.0 21.0 20.0]'
            self.probability = [1.0 0.25 0.25 0.25 0.25 1.0 1.0 1.0 1.0]'
        elseif identifier == 401
            self.name = "Tree 1x1x2x2"
            self.parent = [0,1,2,2,3,3,4,4]
            self.children = child(self.parent)
            self.state = [10.0 10.0 8.0 12.0 9.0 6.0 10.0 13.0]'
            self.probability = [1.0 1.0 0.66 0.34 0.24 0.76 0.46 0.54]'
        elseif identifier == 402
            self.name = "Tree 1x2x2x2"
            self.parent = [0,1,1,2,2,3,3,4,4,5,5,6,6,7,7]
            self.children = child(self.parent)
            self.state = [10.0 12.0 8.0 15.0 11.0 9.0 5.0 18.0 16.0 13.0 11.0 10.0 7.0 6.0 3.0]'
            self.probability = [1.0 0.8 0.7 0.3 0.2 0.8 0.4 0.6 0.2 0.5 0.5 0.4 0.6 0.7 0.3]'
        elseif identifier == 4022
            self.name = "Tree 1x2x2x2"
            self.parent = [0,1,1,2,2,3,3,4,4,5,5,6,6,7,7]
            self.children = child(self.parent)
            self.state = [10.0 0.0; 12.0 1.0; 8.0 -1.0; 15.0 2.0; 11.0 1.0; 9.0 -0.5; 5.0 -2.0; 18.0 3.0; 16.0 1.8; 13.0 0.9; 11.0 0.2; 10.0 0.0; 7.0 -1.2; 6.0 -2.0; 3.0 -3.2]
            self.probability = [1.0 0.8 0.7 0.3 0.2 0.8 0.4 0.6 0.2 0.5 0.5 0.4 0.6 0.7 0.3]'
        elseif identifier ==404
            self.name = "Tree 1x2x2x2"
            self.parent = [0,1,1,2,2,3,3,4,4,5,5,6,6,7,7]
            self.children = child(self.parent)
            self.state = (0.2+0.6744).*[0.0 1.0 -1.0 2.0 0.1 0.0 -2.0 3.0 1.1 0.9 -1.1 1.2 -1.2 -0.8 -3.2]'
            self.probability = [1.0 0.3 0.7 0.2 0.8 0.1 0.9 0.5 0.5 0.6 0.4 0.4 0.6 0.3 0.7]'
        elseif identifier == 405
            self.name = "Tree 5"
            self.parent = [0,1,2,3,3,2,6,6,6,1,10,10,11,12,12,12]
            self.children = child(self.parent)
            self.state = [7.0 12.0 14.0 15.0 13.0 9.0 11.0 10.0 8.0 4.0 6.0 2.0 5.0 0.0 1.0 3.0]'
            self.probability = [1.0 0.7 0.4 0.7 0.3 0.6 0.2 0.3 0.5 0.3 0.2 0.8 1.0 0.3 0.5 0.2]'
        end
        return self
    end

    """
    	Tree(bstructure::Vector{Int64}, dimension=Int64[])

    Returns a tree according to the specified branching vector and the dimension.

    Args:
    - bstructure - the branching structure of the scenario tree.
    - dimension - describes how many values (states) in each node.

    The branching vector must start with 1 for the root node.
    """
    function Tree(bstructure::Vector{Int64}, dimension=Int64[])
        self = new()
        if isempty(dimension)
            dimension = 1
        end
        leaves = bstructure
        for stage = 1 : length(bstructure)
            if stage == 1
                leaves = 1
                self.name = "Tree $(bstructure[1])"
                self.parent = zeros(Int64, bstructure[1])
                self.children = child(self.parent)
                self.state = randn(rng, bstructure[1], dimension)
                self.probability = ones(bstructure[1], 1)
            else
                leaves = leaves * bstructure[stage-1]
                newleaves = vec(ones(Int64, bstructure[stage]) .* transpose(length(self.parent) .+ (1 .- leaves : 0)))
                self.parent =  vcat(self.parent, newleaves)
                self.children = child(self.parent)
                self.state = vcat(self.state, randn(length(newleaves), dimension))
                self.name = "$(self.name)x$(bstructure[stage])"
                tmp = rand(rng, Uniform(0.3, 1.0), bstructure[stage], leaves)
                tmp = tmp ./ (transpose(ones(1, bstructure[stage])) .* sum(tmp, dims = 1))
                self.probability = vcat(self.probability, vec(tmp))
            end
        end
        return self
    end
end

"""
	stage(trr::Tree, node=Int64[])

Returns the stage of each node in the tree.

Args:
- trr - an instance of a Tree.
- node - the number of node in the scenario tree you want to know its stage.
"""
function stage(trr::Tree, node = Int64[])
    if isempty(node)
        node = 1 : length(trr.parent)
    elseif isa(node, Int64)
        node = Int64[node]
    end
    #stage = zeros(length(collect(node))) #you can also use this
    stage = zero(node)
    for i = 1 : length(node)
        pred = node[i]
        while pred > 0 && trr.parent[pred] > 0
            pred = trr.parent[pred]
            stage[i] += 1
        end
    end
    return stage
end

"""
	height(trr::Tree)

Returns the height of the tree which is just the maximum number of the stages of each node.

Args:
- trr - an instance of a Tree.
"""
function height(trr::Tree)
    return maximum(stage(trr))
end

"""
	leaves(trr::Tree,node=Int64[])

Returns the leaf nodes, their indexes and the conditional probabilities.

Args:
- trr - an instance of a Tree.
- node - a node in the tree you want to its children.
"""
function leaves(trr::Tree, node = Int64[])
    nodes = 1 : length(trr.parent)
    leaves = setdiff(nodes, trr.parent)
    #leaves are all nodes which are not parents
    #setdiff(A,B) finds all the elements that belongs to A and are not in B
    omegas = 1 : length(leaves)
    if !isempty(node) && isa(node,Int64)
        node = Int64[node]
    end
    if !isempty(node) && (0 ∉ node)
        omegas = Set{Int64}()
        nodes = leaves
        while any(j !=0 for j ∈ nodes)
            omegas = union(omegas, (ind for (ind, j) ∈ enumerate(nodes) if j ∈ node))
            nodes = Int64[trr.parent[max(0 , j)] for j ∈ nodes]
        end
        omegas = collect(omegas)
    end
    leaves = Int64[leaves[j] for j ∈ omegas]
    prob = ones(Float64 , length(leaves))
    nodes = leaves
    while any(j !=0 for j ∈ nodes)
        prob = prob .* trr.probability[[j for j in nodes]]
        nodes = Int64[trr.parent[max(0 , j)] for j ∈ nodes]
    end
    return leaves, omegas, prob
end

"""
	nodes(trr::Tree,t=Int64[])

Returns the nodes in the tree, generally the range of the nodes in the tree.

Args:
- trr - an instance of a Tree.
- t  - stage in the tree.

Example : nodes(trr,2) - gives all nodes at stage 2.
"""
function nodes(trr::Tree, t = Int64[])
    nodes = 1 : length(trr.parent)
    if isempty(t) #if stage t is not given, return all nodes of the tree
        return nodes
    else # else return nodes at the given stage t
        stg = stage(trr)
        return Int64[i for i in nodes if stg[i] == t]
    end
end

"""
	root(trr::Tree,nodes=Int64[])
Returns the root of the tree if the node is not specified.

Args:
- trr - an instance of Tree.
- nodes - node in the tree you want to know the sequence from the root.

If `nodes` is not specified, it returns the root of the tree.

If `nodes` is specified, it returns a sequence of nodes from the root to the specified node.
"""
function root(trr::Tree, nodes = Int64[])
    if isempty(nodes)
        nodes = trr.children[1]
    elseif isa(nodes , Int64)
        nodes = Int64[nodes]
    end
    root = 1 : length(trr.parent)
    for i in nodes
        iroot = Vector{Int64}([])
        tmp = i
        while tmp > 0
            push!(iroot , tmp)
            tmp = trr.parent[tmp]
        end
        root = Int64[i for i in reverse(iroot) if i in root]
    end
    return root
end

"""
	part_tree(trr::Tree)

Returns a vector of trees in d-dimension.

Args:
- trr - an instance of Tree.
"""
function part_tree(trr::Tree)
    trees = Tree[]
    for col = 1:size(trr.state , 2)
        subtree = Tree("Tree of state $col", trr.parent, trr.children, hcat(trr.state[:,col]), trr.probability)
        push!(trees,subtree)
    end
    return trees
end

"""
	build_probabilities!(trr::Tree,probabilities::Array{Float64,2})

Returns the probabilities of the nodes without probabilities
if the array of probabilities is less than the length of parents
in the stochastic approximation procedure.
"""
function build_probabilities!(trr::Tree, probabilities::Array{Float64,2})
    leaf, omegas, probaLeaf = leaves(trr)
    if length(probabilities) == length(nodes(trr))
        trr.probability .= probabilities
    else
        probabilities .= max.(0.0,probabilities)
        j = leaf; i = trr.parent[j]; j .= j[i .> 0]; i .= i[i .> 0]
        trr.probability .= zeros(length(trr.state[:,1]), 1)
        trr.probability[j] = probabilities
        while !isempty(i)
            for k = 1:length(i)
                trr.probability[i[k]] = trr.probability[i[k]] + trr.probability[j[k]]
            end
            trr.probability[j] .= trr.probability[j] ./ trr.probability[i]
            trr.probability[isnan.(trr.probability)] .= 0.0
            j = unique(i); i = trr.parent[j]; j = j[i .> 0] ; i = i[i .> 0]
        end
    end
    return trr.probability
end

"""
	tree_plot(trr::Tree,fig=1)

Returns the plot of the input tree annotated with density of
probabilities of reaching the leaf nodes in the tree.
Args:
- trr - A scenario tree.
- fig - Specifies the size of the image you want to be returned, default = 1.

Using the matplotlib version, fig usually takes a tuple. For example, fig = (8,6) produces an image with length = 8 and width = 6.
The user is not obliged to give this specifications as already there is a default value.
"""
function tree_plot(trr::Tree, fig= 1)
    if !isempty(fig)
        figure(fig)
    end
    #suptitle(trr.name, fontsize = 14)
    trs = subplot2grid((1,4), (0,0), colspan=3)
    title("states")
    stg = stage(trr)
    xticks(1 : height(trr) + 1)         # Set the ticks on the x-axis
    xlabel("stage,time",fontsize=12)
    trs.spines["top"].set_visible(false)         # remove the line of the box at the top
    trs.spines["right"].set_visible(false)		 # remove the line of the box at the right
    for i = 1 : length(trr.parent)
        if stg[i] > 0
            trs.plot([stg[i],stg[i]+1],[trr.state[trr.parent[i]],trr.state[i]],linewidth=1.5)
        end
    end
    prs = subplot2grid((1,4), (0,3))
    title("probabilities")
    prs.spines["top"].set_visible(false)
    prs.spines["left"].set_visible(false)
    prs.spines["right"].set_visible(false)

    (Yi,_,probYi) = leaves(trr)
    Yi = [trr.state[i] for i in Yi]
    nY = length(Yi)
    h = 1.05 * std(Yi) / (nY^0.2) + 1e-3         #Silverman rule of thumb
    #trs.set_ylim(minimum(Yi)-h, maximum(Yi)+h)
    #prs.set_ylim(minimum(Yi)-h, maximum(Yi)+h)
    yticks(())                                                                                #remove the ticks on probability plot
    t = LinRange(minimum(Yi)-h, maximum(Yi)+h, 100)  #100 points on probability plot
    density = zero(t)
    for (i, ti) in enumerate(t)
        for (j, xj) in enumerate(Yi)
            tmp = (xj - ti) / h
            density[i] += probYi[j]* 35/32 * max(1.0 -tmp^2, 0.0)^3 /h #triweight kernel
        end
    end
    plot(density, t)
    prs.fill_betweenx(t, 0 , density)
end

"""
	plot_hd(newtree::Tree)

Returns a plots of trees in higher dimensions.
"""
function plot_hd(newtree::Tree)
    fig = figure(figsize = (10,6))
    stg = stage(newtree)
    for rw = 1:size(newtree.state,2)
      ax = subplot(1,size(newtree.state,2),rw)
      ax.spines["top"].set_visible(false)                                                  #remove the box top
      ax.spines["right"].set_visible(false)                                               #remove the box right
      for i in range(1,stop = length(newtree.parent))
          if stg[i] > 0
              ax.plot([stg[i], stg[i]+1],[newtree.state[:,rw][newtree.parent[i]], newtree.state[:,rw][i]])
          end
          xlabel("stage, time")
          ylabel("states")
      end
      xticks(unique(stg))
    end
end


text:
# test attn:
#
using Knet, NNHelferlein

depth, mb, t = 3, 4, 6

h_enc = randn(Float32, depth, mb, t)
h_t = randn(Float32, depth, mb)

function test_attn(attn)
    a = attn(depth, depth)
    c,α = a(h_t, h_enc)

    return size(c) == (3,4)
end

function test_attnDot()
    a = AttnDot()
    c,α = a(h_t, h_enc)

    return size(c) == (3,4)
end

function test_attnLocation()
    a = AttnLocation(t, depth)
    c,α = a(h_t, h_enc)

    return size(c) == (3,4)
end

function test_attnInFeed()
    a = AttnInFeed(t, depth, depth)
    c,α = a(h_t, init0(depth), h_enc)

    return size(c) == (3,4)
end

# transformer tests:
#
function test_dpa()
    kqv = rand(Float32, 3,8,10)
    c,a = dot_prod_attn(kqv, kqv, kqv)

    return size(c) == (3,8,10) && size(a) == (8,8,10)
end

function test_masks()
    seqs = rand(1:16, 4,6)
    el = Embed(16,8)
    seqs_e = el(seqs)

    pl = PositionalEncoding()
    pos_enc = pl(seqs)   # asser 4x6

    peek_ah = mk_peek_ahead_mask(seqs)
    padd = mk_padding_mask(seqs)

    return size(pos_enc) == (4,6) &&
           size(peek_ah) == (4,4) &&
           size(padd) == (4,1,1,6)
end

function test_dotp_attn()

    function separate_heads(x, n)
        depth, seq, mb = size(x)
        mh_depth = depth ÷ n
        x = reshape(x, mh_depth, n, :, mb)     # split depth in 2 separate dims for heads
        return permutedims(x, (1,3,2,4))       # bring seq-len back to 2nd dim
    end

    seqs = rand(1:16, 4,6)
    el = Embed(16,8)
    emb = el(seqs)
    heads = separate_heads(emb, 2)
    padd = mk_padding_mask(seqs)

    dpa = dot_prod_attn(heads, heads, heads, mask=padd)

    return size(dpa[1]) == (4,4,2,6)
end

ext: jl


text:
module AttachConstraintsTests

using Test
using FillArrays
using LinearAlgebra
using Gridap.Helpers
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Arrays
using Gridap.Integration
using Gridap.CellData
using Gridap.TensorValues

using Gridap.Fields: MockField, MockBasis, OtherMockBasis

ncells = 10
ndofs = 3
ndofs_c = 4
cellmat = [ rand(ndofs,ndofs) for cell in 1:ncells]
cellvec = [ rand(ndofs) for cell in 1:ncells]
cellmatvec = pair_arrays(cellmat,cellvec)
cellconstr = [ rand(ndofs_c,ndofs) for cell in 1:ncells]

a = attach_constraints_rows(cellvec,cellconstr)
r = map( (vec,constr) -> constr*vec ,cellvec,cellconstr)
test_array(a,r)

a = attach_constraints_rows(cellmat,cellconstr)
r = map( (mat,constr) -> constr*mat ,cellmat,cellconstr)
test_array(a,r)

a = attach_constraints_rows(cellmatvec,cellconstr)
r = map( (mat,vec,constr) -> (constr*mat,constr*vec),cellmat,cellvec,cellconstr)
test_array(a,r)

a = attach_constraints_cols(cellmat,cellconstr)
r = map( (mat,constr) -> mat*transpose(constr) ,cellmat,cellconstr)
test_array(a,r)

a = attach_constraints_cols(cellmatvec,cellconstr)
r = map( (mat,vec,constr) -> (mat*transpose(constr),vec),cellmat,cellvec,cellconstr)
test_array(a,r)

#cache = array_cache(a)
#using BenchmarkTools
#@btime getindex!($cache,$a,2)

ndofs = 4

z = 2.0
v = VectorValue(3.0,1.5)
w = VectorValue(3.4,3.5)
a = MockBasis{2}(v,ndofs)
b = MockBasis{2}(w,ndofs)
c = fill(1.0,ndofs)
f = OtherMockBasis{2}(ndofs)
g = MockField{2}(v)

l = 10
zl = [ z for  i in 1:l]
cl = fill(c,l)
fl = Fill(f,l)
ϕl = lincomb(fl,cl)
gl = fill(g,l)
al = Fill(a,l)
bl = fill(b,l)

cell_axes = Fill((Base.OneTo(ndofs),),l)
gf = GenericCellField(gl,ϕl,Val(true))
af = GenericCellField(al,ϕl,Val(true),cell_axes,Val((:,)))
bf = GenericCellField(bl,ϕl,Val(true),cell_axes,Val((:,)))
zf = convert_to_cell_field(zl,ϕl)
df = af*zf
dft = trialize_cell_basis(df)

degree = 3
quad = CellQuadrature(degree,[QUAD,],ones(Int,l))

ndofs_c = 5
cellconstr = [ rand(ndofs_c,ndofs) for i in 1:l]
cellvals = [ rand(ndofs) for i in 1:l]
cellvec = integrate(bf⋅v,ϕl,quad)
cellmat = integrate(∇(bf)⊙∇(dft),ϕl,quad)
cellmatvec = pair_arrays(cellmat,cellvec)
cellmatvec = attach_dirichlet(cellmatvec,cellvals)
cellmatvec = attach_constraints_rows(cellmatvec,cellconstr)
cellmatvec = attach_constraints_cols(cellmatvec,cellconstr)

@test size(cellmatvec[1][1]) == (ndofs_c,ndofs_c)
@test size(cellmatvec[1][2]) == (ndofs_c,)

cellconstr = identity_constraints(cell_axes)
test_array(cellconstr,Fill(Matrix(I,ndofs,ndofs),l))
@test isa(cellconstr,Fill)

# Test at skeleton
aS = merge_cell_fields_at_skeleton(af,af)
dSt = merge_cell_fields_at_skeleton(dft,dft)

axesL = Fill((Base.OneTo(ndofs),),l)
axesR = axesL
cellvalsL = [ rand(ndofs) for i in 1:l]
cellvalsR = cellvalsL
cellvals = merge_cell_dofs_at_skeleton(cellvalsL,cellvalsR,axesL,axesR)
cellconstrL = [ rand(ndofs_c,ndofs) for i in 1:l]
cellconstrR = cellconstrL
axesL_rows = Fill((Base.OneTo(ndofs_c),),l)
axesR_rows = axesL_rows
axesL_cols = axesL
axesR_cols = axesR
cellconstr = merge_cell_constraints_at_skeleton(cellconstrL,cellconstrR,axesL_rows,axesR_rows,axesL_cols,axesR_cols)

cellvec = integrate(jump(aS⋅v),ϕl,quad)
cellmat = integrate( jump(aS⋅v)*(w⋅dSt.⁻),ϕl,quad)
cellmatvec = pair_arrays(cellmat,cellvec)
cellmatvec = attach_dirichlet(cellmatvec,cellvals)
cellmatvec = attach_constraints_rows(cellmatvec,cellconstr)
cellmatvec = attach_constraints_cols(cellmatvec,cellconstr)
test_array(cellmatvec,collect(cellmatvec))

cellconstr = identity_constraints(cell_axes)
cellconstr = merge_cell_constraints_at_skeleton(cellconstr,cellconstr,cell_axes,cell_axes,cell_axes,cell_axes)
@test isa(cellconstr,VectorOfBlockArrayCoo)

##a = cellmatvec
#a = cellconstr
#cache = array_cache(a)
#using BenchmarkTools
#@btime getindex!($cache,$a,2)

end # module


text:
@testset "device" begin

    @test device(rand(2)) == Val(:cpu)
    @test device(Dense(2, 3)) == Val(:cpu)
    @test device(Conv((2, 2), 1 => 16, relu)) == Val(:cpu)
    @test device(Chain(x -> x .^ 2, Dense(2, 3))) == Val(:cpu)

    if CUDA.functional()
        @test device(rand(2) |> gpu) isa CuDevice
        @test device(Dense(2, 3) |> gpu) isa CuDevice
        @test device(Conv((2, 2), 1 => 16, relu) |> gpu) isa CuDevice
        @test device(Chain(x -> x .^ 2, Dense(2, 3)) |> gpu) isa CuDevice
    end

end


text:
for (elty, relty, ext) in ((:Float32, :Float32, :s),
                           (:Float64, :Float64, :d))
    for (matA, matB, sym) in ((:Matrix, :Matrix, "_"),
                              (:DistMatrix, :DistMatrix, "Dist_"),
                              (:SparseMatrix, :Matrix, "Sparse_"),
                              (:DistSparseMatrix, :DistMultiVec, "DistSparse_"))
        @eval begin
            function leastSquares!(A::$matA{$elty}, B::$matB{$elty}, X::$matB{$elty};
                orientation::Orientation = NORMAL)
                ElError(ccall(($(string("ElLeastSquares", sym, ext)), libEl), Cuint,
                    (Cuint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
                    orientation, A.obj, B.obj, X.obj))
                return X
            end
        end
    end
end

function leastSquares(A::DistMatrix{T}, B::DistMatrix{T};
    orientation::Orientation = NORMAL) where {T}
    X = DistMatrix(T, MC, MR, A.g)
    return leastSquares!(A, B, X, orientation = orientation)
end
function leastSquares(A::DistSparseMatrix{T}, B::DistMultiVec{T};
    orientation::Orientation = NORMAL) where {T}
    X = DistMultiVec(T, comm(A))
    return leastSquares!(A, B, X, orientation = orientation)
end

for (elty, relty, ext) in ((:Float32, :Float32, :s),
                           (:Float64, :Float64, :d))
    for (matA, matB, sym) in ((:(Elemental.Matrix), :(Elemental.Matrix), "_"),
                              (:DistMatrix, :DistMatrix, "Dist_"),
                              (:(Elemental.SparseMatrix), :(Elemental.Matrix), "Sparse_"),
                              (:DistSparseMatrix, :DistMultiVec, "DistSparse_"))
        @eval begin
            function bpdn!(A::$matA{$elty}, B::$matB{$elty}, lambda::$elty, X::$matB{$elty})
                ElError(ccall(($(string("ElBPDNX", sym, ext)), libEl), Cuint,
                    (Ptr{Cvoid}, Ptr{Cvoid}, $elty, Ptr{Cvoid}, Cuint),
                    A.obj, B.obj, lambda, X.obj, 0))
                return X
            end
        end
    end
end

function bpdn(A::DistMatrix{T}, B::DistMatrix{T}, lambda::T) where {T}
    X = DistMatrix(T, MC, MR, A.g)
    return bpdn!(A, B, lambda, X)
end
function bpdn(A::DistSparseMatrix{T}, B::DistMultiVec{T}, lambda) where {T}
    X = DistMultiVec(T, comm(A))
    return bpdn!(A, B, lambda, X)
end

@enum RidgeAlg RIDGE_CHOLESKY RIDGE_QR RIDGE_SVD

for (elty, relty, ext) in ((:Float32, :Float32, :s),
                           (:Float64, :Float64, :d))
    for (matA, matB, sym) in ((:(Elemental.Matrix), :(Elemental.Matrix), "_"),
                              (:DistMatrix, :DistMatrix, "Dist_"),
                              (:(Elemental.SparseMatrix), :(Elemental.Matrix), "Sparse_"),
                              (:DistSparseMatrix, :DistMultiVec, "DistSparse_"))
        @eval begin
            function ridge!(A::$matA{$elty}, B::$matB{$elty}, gamma::$elty, X::$matB{$elty};
                orientation::Orientation = NORMAL, alg::RidgeAlg = RIDGE_CHOLESKY)
                ElError(ccall(($(string("ElRidge", sym, ext)), libEl), Cuint,
                    (Cuint, Ptr{Cvoid}, Ptr{Cvoid}, $elty, Ptr{Cvoid}, Cuint),
                    orientation, A.obj, B.obj, gamma, X.obj, alg))
                return X
            end
        end
    end
end

function ridge(A::DistMatrix{T}, B::DistMatrix{T}, gamma::T; ka...) where {T}
    X = DistMatrix(T, MC, MR, A.g)
    return ridge!(A, B, gamma, X; ka...)
end
function ridge(A::DistSparseMatrix{T}, B::DistMultiVec{T}, gamma::T; ka...) where {T}
    X = DistMultiVec(T, comm(A))
    return ridge!(A, B, gamma, X; ka...)
end


text:
##
# VOO Routines

struct VOOActionGenerator{D1 <: Distribution, D2 <: Distribution, R <: AbstractRNG} <: Function
	exploration_prob::Float64
	exploration_sampler::D1
    voronoi_sampler::D2
    action_rng::R
    acceptance_radius::Float64
    early_halt::Bool
    halt_count::Int64
end

# If only using acceptance radius
function VOOActionGenerator(exploration_prob::Float64,
                            exploration_sampler::Distribution,
                            voronoi_sampler::Distribution,
                            action_rng::AbstractRNG,
                            acceptance_radius::Float64)
    early_halt = false
    halt_count = 1
    return VOOActionGenerator(exploration_prob, exploration_sampler, voronoi_sampler, action_rng, acceptance_radius, early_halt, halt_count)
end

# If only using halt count
function VOOActionGenerator(exploration_prob::Float64,
                            exploration_sampler::Distribution,
                            voronoi_sampler::Distribution,
                            early_halt::Bool,
                            halt_count::Int64)
    acceptance_radius = 0.0
    return VOOActionGenerator(exploration_prob, exploration_sampler, voronoi_sampler, action_rng, acceptance_radius, early_halt, halt_count)
end

function voo_sample(f::VOOActionGenerator, a_list, q_list, a_space)
	# determine explore vs. exploit
	w = rand(f.exploration_sampler)

	if w <= f.exploration_prob || length(a_list) <= 1
		# uniform sample from action space
		a_new = rand(f.action_rng, a_space)
	else
	    # find the best value
	    a_new = voronoi_sample_centered(f, a_list, q_list, a_space)
	end
    return a_new
end

function voo_sample(f::VOOActionGenerator, problem::Union{MDP,POMDP}, b, h::DPWStateNode)
    # extract actions
    t = h.tree
    a_inds = t.children[h.index]
    a_list = t.a_labels[a_inds]
    a_space = actions(problem, b)

    # determine explore vs. exploit
    w = rand(f.exploration_sampler)

    if w <= f.exploration_prob || length(a_list) <= 1
        # uniform sample from action space
        a_new = rand(f.action_rng, a_space)
    else
        # extract q values
        q_list = [t.q[a] for a in a_inds]

        # find the best value
        a_new = voronoi_sample_centered(f, a_list, q_list, a_space)
    end
    return a_new
end

function voo_sample(f::VOOActionGenerator, problem::Union{MDP,POMDP}, b, h::POWTreeObsNode)
	# extract actions
	t = h.tree
    a_inds = t.tried[h.node]
    a_list = t.a_labels[a_inds]
    a_space = actions(problem, b)

    # determine explore vs. exploit
	w = rand(f.exploration_sampler)

	if w <= f.exploration_prob || length(a_list) <= 1
		# uniform sample from action space
		a_new = rand(f.action_rng, a_space)
	else
	    # extract q values
		q_list = [t.v[a] for a in a_inds]

	    # find the best value
	    a_new = voronoi_sample_centered(f, a_list, q_list, a_space)
	end
    return a_new
end

function voronoi_sample_centered(f::VOOActionGenerator, a_list, q_list, a_space)
	best_a = a_list[argmax(q_list)]
	a_new = best_a
    closest = false

    a_closest = best_a
    dist_closest = Inf

    # iterate until it lies within the best Voronoi cell
    n = 1
    while .!closest
    	closest = true
        a_new = action_sample_centered(best_a, f.voronoi_sampler, a_space)
        dist_to_best = action_dist(a_new, best_a)

        if dist_to_best < dist_closest
            a_closest = a_new
            dist_closest = dist_to_best
        end

        # if early halting, return the closest one to the center
        if f.early_halt && n >= f.halt_count
            return a_closest
        end

        # if auto acceptance radius, return if the sample is within the radius
        if dist_to_best < f.acceptance_radius
        	return a_new
        end

        # iterate over all actions to see if it is closer to another Voronoi cell
        for a in a_list
        	if action_dist(a_new, a) < dist_to_best
        		closest = false
        	end
        end
        n += 1
    end
    return a_new
end

function (f::VOOActionGenerator)(pomdp, b, h)
    a_new = voo_sample(f, pomdp, b, h)
    return a_new
end

# Next action function that uses the rollout policy as the first action
abstract type AbstractRolloutFirst end

struct RolloutFirst{R <: AbstractRNG, P <: Policy} <: AbstractRolloutFirst
    rng::R
    rollout_policy::P
end

function MCTS.next_action(gen::RolloutFirst, p::POMDP, b, h)
    if isroot(h) && n_children(h) < 1
        return POMDPs.action(gen.rollout_policy, b)
    else
        return rand(gen.rng, actions(p))
    end
end

MCTS.next_action(gen::RolloutFirst, p::GenerativeBeliefMDP, b, h) = next_action(gen, p.pomdp, b, h)

struct VOORolloutFirst{V <: VOOActionGenerator, P <: Policy} <: AbstractRolloutFirst
    voo::V
    rollout_policy::P
end

function MCTS.next_action(gen::VOORolloutFirst, p::POMDP, b, h)
    if isroot(h) && n_children(h) < 1
        return POMDPs.action(gen.rollout_policy, b)
    else
        return voo_sample(gen.voo, p, b, h)
    end
end


text:
"""
    $(TYPEDEF)
Simple in-memory data store with a specified data type and a specified key type.
"""
struct InMemoryDataStore{T, E <: AbstractBigDataEntry} <: AbstractDataStore{T}
    data::Dict{Symbol, T}
    entries::Dict{Symbol, E}
end

"""
    $(SIGNATURES)
Create an in-memory store using a specific data type.
"""
function InMemoryDataStore{T, E}() where {T, E <: AbstractBigDataEntry}
    return InMemoryDataStore{T, E}(Dict{Symbol, T}(), Dict{Symbol, E}())
end

"""
    $(SIGNATURES)
Create an in-memory store using binary data (UInt8) as a type.
"""
function InMemoryDataStore()
    return InMemoryDataStore{Vector{UInt8}, GeneralBigDataEntry}()
end


text:
abstract type AbstractPortfolioFormulation end
abstract type AbstractMeanVariance <: AbstractPortfolioFormulation end
abstract type AbstractSampleBased <: AbstractPortfolioFormulation end


text:
using Pkg
Pkg.activate(".")
include("docs/make.jl")


text:
# Solve the Landau-Zener problem i ψ' = H(t) ψ, with H(t) = [t 1;1 -t]

using Test
using StaticArrays, LinearAlgebra
using OrdinaryDiffEq, DiffEqBase

H(t) = -im*(@SMatrix [t 1;1 -t])

fun(ψ,p,t) = oftype(ψ, H(t)*ψ)
fun_inplace(dψ,ψ,p,t) = (dψ .= H(t)*ψ)

T = 0.1
tspan = (0,T)
explicit = [Midpoint, RK4, DP5, Tsit5, Vern7]
implicit = [ImplicitEuler, Trapezoid, Kvaerno3, Rosenbrock23]

@testset "Complex Tests on Explicit Methods. alg=$alg" for alg in explicit
  for f in (fun, fun_inplace)
    ψ0 = [1.0 + 0.0im; 0.0]
    prob = ODEProblem(f, ψ0, (-T, T))
    sol = solve(prob, alg())
    @test norm(sol(T)) ≈ 1 atol=1e-2
  end
  ψ0 = @SArray [1.0 + 0.0im; 0.0]
  prob = ODEProblem(fun, ψ0, (-T, T))
  sol = solve(prob, alg())
  @test norm(sol(T)) ≈ 1 atol=1e-2
end

@testset "Complex Tests on Implicit Autodiff Methods. alg=$alg" for alg in implicit
  @test_broken begin
    for f in (fun, fun_inplace)
      ψ0 = [1.0 + 0.0im; 0.0]
      prob = ODEProblem(f, ψ0, (-T, T))
      sol = solve(prob, alg())
      @test norm(sol(T)) ≈ 1 atol=1e-2
    end
    ψ0 = @SArray [1.0 + 0.0im; 0.0]
    prob = ODEProblem(fun, ψ0, (-T, T))
    sol = solve(prob, alg())
    @test norm(sol(T)) ≈ 1 atol=1e-2
  end
end

@testset "Complex Tests on Implicit Finite Diff Methods. alg=$alg" for alg in implicit
  ψ0 = [1.0 + 0.0im; 0.0]
  prob = ODEProblem(fun_inplace, ψ0, (-T, T))
  sol = solve(prob, alg(autodiff=false))
  @test norm(sol(T)) ≈ 1 atol=1e-2
end

@testset "Complex Tests on Implicit Finite Diff Out-of-place Methods. alg=$alg" for alg in implicit
  ψ0 = [1.0 + 0.0im; 0.0]
  prob = ODEProblem(fun, ψ0, (-T, T))
  sol = solve(prob, alg(autodiff=false))
  @test norm(sol(T)) ≈ 1 atol=1e-2
end

@testset "Complex Tests on Implicit Finite Diff Out-of-place Methods SArray. alg=$alg" for alg in implicit
  ψ0 = @SArray [1.0 + 0.0im; 0.0]
  prob = ODEProblem(fun, ψ0, (-T, T))
  sol = solve(prob, alg(autodiff=false))
  @test norm(sol(T)) ≈ 1 atol=1e-2
end


text:
using DimensionalData, Test, Unitful, OffsetArrays, SparseArrays, Dates, Random, ArrayInterface
using DimensionalData: layerdims

using DimensionalData.LookupArrays, DimensionalData.Dimensions

a = [1 2; 3 4]
a2 = [1 2 3 4
      3 4 5 6
      4 5 6 7]
xmeta = Metadata(:meta => "X")
ymeta = Metadata(:meta => "Y")
ameta = Metadata(:meta => "da")
dimz = (X(Sampled(143.0:2.0:145.0; order=ForwardOrdered(), metadata=xmeta)),
        Y(Sampled(-38.0:2.0:-36.0; order=ForwardOrdered(), metadata=ymeta)))
dimz2 = (Dim{:row}(10:10:30), Dim{:column}(-20:10:10))
refdimz = (Ti(1:1),)
da = @test_nowarn DimArray(a, dimz; refdims=refdimz, name=:test, metadata=ameta)
val(dims(da, 1)) |> typeof
da2 = DimArray(a2, dimz2; refdims=refdimz, name=:test2)

@testset "size and axes" begin
    row_dims = (1, Dim{:row}(), Dim{:row}, :row, dimz2[1])
    for dim in row_dims
        @test size(da2, dim) == 3
        @test axes(da2, dim) == 1:3
        @test firstindex(da2, dim) == 1
        @test lastindex(da2, dim) == 3
    end
    @test size(da2, :column) == 4
    @test axes(da2, :column) == 1:4
    @test size(da2) == (3, 4)
    @test axes(da2) == (1:3, 1:4)
    @test firstindex(da2) == 1
    @test lastindex(da2) == 12
    @inferred axes(da2, Dim{:column})
    @test IndexStyle(da) == IndexLinear()
end

@testset "interface methods" begin
    lx = Sampled(143.0:2.0:145.0, ForwardOrdered(), Regular(2.0), Points(), xmeta)
    ly = Sampled(-38.0:2.0:-36.0, ForwardOrdered(), Regular(2.0), Points(), ymeta)
    @test dims(da) == (X(lx), Y(ly))
    @test dims(da, X) == X(lx)
    @test refdims(da) == refdimz
    @test name(da) == :test
    @test metadata(da) == ameta
    @test lookup(da) == (lx, ly)
    @test order(da) == (ForwardOrdered(), ForwardOrdered())
    @test sampling(da) == (Points(), Points())
    @test span(da) == (Regular(2.0), Regular(2.0))
    @test bounds(da) == ((143.0, 145.0), (-38.0, -36.0))
    @test locus(da) == (Center(), Center())
    @test layerdims(da) == (X(), Y())
    @test index(da, Y) == LinRange(-38.0, -36.0, 2)
end

@testset "copy and friends" begin
    dac = copy(da2)
    @test dac == da2
    @test dims(dac) == dims(da2)
    @test refdims(dac) == refdims(da2) == (Ti(1:1),)
    @test name(dac) == name(da2) == :test2
    @test metadata(dac) == metadata(da2)
    dadc = deepcopy(da2)
    @test dadc == da2
    @test dims(dadc) == dims(da2)
    @test refdims(dadc) == refdims(da2) == (Ti(1:1),)
    @test name(dadc) == name(da2) == :test2
    @test metadata(dadc) == metadata(da2)

    o = one(da)
    @test o == [1 0; 0 1]
    @test dims(o) == dims(da)

    ou = oneunit(da)
    @test ou == [1 0; 0 1]
    @test dims(ou) == dims(da)

    z = zero(da)
    @test z == [0 0; 0 0]
    @test dims(z) == dims(da)

    @test Array(da) == [1 2; 3 4]
    @test Array(da) isa Array{Int,2}
    @test collect(da) == [1 2; 3 4]
    @test collect(da) isa Array{Int,2}
    @test vec(da) == [1, 3, 2, 4]
    @test vec(da) isa Array{Int,1}

    # This should do nothing
    A = read(da2)
    @test A === da2
end

@testset "similar" begin
    da_sim = similar(da2)
    @test eltype(da_sim) == eltype(da2)
    @test size(da_sim) == size(da2)
    @test dims(da_sim) == dims(da2)
    @test refdims(da_sim) == refdims(da2)

    da_float = similar(da2, Float64)
    @test eltype(da_float) == Float64
    @test size(da_float) == size(da2)
    @test dims(da_float) == dims(da2)
    @test refdims(da_float) == refdims(da2)

    # Changing the axis size removes dims.
    # TODO we can keep dims, but with NoLookup?
    da_size_float = similar(da2, Float64, (10, 10))
    @test eltype(da_size_float) == Float64
    @test size(da_size_float) == (10, 10)
    @test typeof(da_size_float) <: Array{Float64,2}
    da_size_float_splat = similar(da2, Float64, 10, 10)
    @test size(da_size_float_splat) == (10, 10)
    @test typeof(da_size_float_splat)  <: Array{Float64,2}

    sda = DimArray(sprand(Float64, 10, 10, 0.5), (X, Y))
    sparse_size_int = similar(sda, Int64, (5, 5))
    @test eltype(sparse_size_int) == Int64 != eltype(sda)
    @test size(sparse_size_int) == (5, 5)
    @test sparse_size_int isa SparseMatrixCSC

    @test dims(da_float) == dims(da2)
    @test refdims(da_float) == refdims(da2)
end

@testset "replace" begin
    dar = replace(da, 1 => 99)
    dar isa DimArray
    @test dar == [99 2; 3 4]

    dar = copy(da)
    replace!(dar, 4 => 99)
    @test dar == [1 2; 3 99]
end

@testset "broadcast" begin
    da = DimArray(ones(Int, 5, 2, 4), (Y(10:2:18), Ti(10:11), X(1:4)))
    dab = da .* 2.0
    @test dab == fill(2.0, 5, 2, 4)
    @test eltype(dab) <: Float64
    @test dims(dab) ==
        (Y(Sampled(10:2:18, ForwardOrdered(), Regular(2), Points(), NoMetadata())),
         Ti(Sampled(10:11, ForwardOrdered(), Regular(1), Points(), NoMetadata())),
         X(Sampled(1:4, ForwardOrdered(), Regular(1), Points(), NoMetadata())))
    dab = da .+ fill(10, 5, 2, 4)
    @test dab == fill(11, 5, 2, 4)
    @test eltype(dab) <: Int
end

@testset "eachindex" begin
    # Should have linear index
    da = DimArray(ones(5, 2, 4), (Y(10:2:18), Ti(10:11), X(1:4)))
    @test eachindex(da) == eachindex(parent(da))
    # Should have cartesian index
    sda = DimArray(sprand(10, 10, .1), (Y(1:10), X(1:10)))
    @test eachindex(sda) == eachindex(parent(sda))
end

@testset "convert" begin
    ac = convert(Array, da2)
    @test ac isa Array{Int,2}
    @test ac == a2
end

if VERSION > v"1.1-"
    @testset "copy!" begin
        dimz = dims(da2)
        A = zero(a2)
        sp = sprand(Int, 4, 0.5)
        db = DimArray(deepcopy(A), dimz)
        dc = DimArray(deepcopy(A), dimz)

        copy!(A, da2)
        @test A == parent(da2)
        copy!(db, da2)
        @test parent(db) == parent(da2)
        copy!(dc, a2)
        @test parent(db) == a2
        # Sparse vector has its own method for ambiguity
        copy!(sp, da2[1, :])
        @test sp == parent(da2[1, :])

        @testset "vector copy! (ambiguity fix)" begin
            v = zeros(3)
            dv = DimArray(zeros(3), X)
            copy!(v, DimArray([1.0, 2.0, 3.0], X))
            @test v == [1.0, 2.0, 3.0]
            copy!(dv, DimArray([9.9, 9.9, 9.9], X))
            @test dv == [9.9, 9.9, 9.9]
            copy!(dv, [5.0, 5.0, 5.0])
            @test dv == [5.0, 5.0, 5.0]
        end

    end
end

@testset "copyto!" begin
    A = zero(a2)
    da = DimArray(ones(size(A)), dims(da2))
    copyto!(A, da)
    @test all(A .== 1)
    copyto!(da, 1, zeros(5, 5), 1, 12)
    @test all(da .== 0)
    x = reshape(10:10:40, 1, 4)
    copyto!(da, CartesianIndices(view(da, 1:1, 1:4)), x, CartesianIndices(x))
    @test da[1, 1:4] == 10:10:40
    copyto!(A, CartesianIndices(view(da, 1:1, 1:4)), DimArray(x, (X, Y)), CartesianIndices(x))
    @test A[1, 1:4] == 10:10:40
end

@testset "constructor" begin
    da = DimArray(; data=rand(5, 4), dims=(X, Y))
    @test_throws DimensionMismatch DimArray(1:5, X(1:6))
    @test_throws DimensionMismatch DimArray(1:5, (X(1:5), Y(1:2)))
    da_reconstructed = DimArray(da)
    @test da == da_reconstructed
    @test dims(da) == dims(da_reconstructed)
end

@testset "fill constructor" begin
    da = fill(5.0, X(4), Y(40.0:10.0:80.0))
    @test parent(da) == fill(5.0, (4, 5))
    @test dims(da) == (
         X(NoLookup(Base.OneTo(4))),
         Y(Sampled(40.0:10.0:80.0, ForwardOrdered(), Regular(10.0), Points(), NoMetadata()))
    )
    @test_throws ArgumentError fill(5.0, (X(:e), Y(8)))
end

@testset "ones, zeros, trues, falses constructors" begin
    da = zeros(X(4), Y(40.0:10.0:80.0))
    @test eltype(da) <: Float64
    @test all(==(0), da)
    @test size(da) == (4, 5)
    @test dims(da) == (
         X(NoLookup(Base.OneTo(4))),
         Y(Sampled(40.0:10.0:80.0, ForwardOrdered(), Regular(10.0), Points(), NoMetadata()))
    )
    ti = Ti(Date(2001):Year(1):Date(2004))
    da = ones(Int32, ti)
    @test size(da) == (4,)
    @test eltype(da) <: Int32
    @test all(==(1), da)
    @test dims(da) == (Ti(Sampled(Date(2001):Year(1):Date(2004), ForwardOrdered(), Regular(Year(1)), Points(), NoMetadata())),)
    da = trues(ti)
    @test size(da) == (4,)
    @test eltype(da) <: Bool
    @test all(==(true), da)
    @test dims(da) == (Ti(Sampled(Date(2001):Year(1):Date(2004), ForwardOrdered(), Regular(Year(1)), Points(), NoMetadata())),)
    da = falses(ti)
    @test size(da) == (4,)
    @test eltype(da) <: Bool
    @test all(==(false), da)
    @test dims(da) == (Ti(Sampled(Date(2001):Year(1):Date(2004), ForwardOrdered(), Regular(Year(1)), Points(), NoMetadata())),)
end

@testset "rand constructors" begin
    da = rand(1:10, X(8), Y(11:20))
    @test size(da) == (8, 10)
    @test eltype(da) <: Int
    @test dims(da) == (
         X(NoLookup(Base.OneTo(8))),
         Y(Sampled(11:20, ForwardOrdered(), Regular(1), Points(), NoMetadata()))
    )
    da = rand(X([:a, :b]), Y(3))
    @test size(da) == (2, 3)
    @test eltype(da) <: Float64
    da = rand(Bool, X([:a, :b]), Y(3))
    @test size(da) == (2, 3)
    @test eltype(da) <: Bool
    da = rand(MersenneTwister(), Float32, X([:a, :b]), Y(3))
    @test size(da) == (2, 3)
    @test eltype(da) <: Float32
end

@testset "OffsetArray" begin
    oa = OffsetArray(a2, -1:1, 5:8)
    @testset "Regular dimensions don't work: axes must match" begin
        dimz = (X(100:100:300), Y([:a, :b, :c, :d]))
        @test_throws DimensionMismatch DimArray(oa, dimz)
    end
    odimz = (X(OffsetArray(100:100:300, -1:1)), Y(OffsetArray([:a, :b, :c, :d], 5:8)))
    oda = DimArray(oa, odimz)
    @testset "Indexing and selectors work with offsets" begin
        @test axes(oda) == (-1:1, 5:8)
        @test oda[-1, 5] == oa[-1, 5] == 1
        @test oda[Near(105), At(:a)] == oa[-1, 5] == 1
        @test oda[Between(100, 250), At(:a)] == oa[-1:0, 5] == [1, 3]
    end
    @testset "Subsetting reverts to a regular array and dims" begin
        @test axes(oda[0:1, 7:8]) == (1:2, 1:2)
        @test axes.(dims(oda[0:1, 7:8])) == ((1:2,), (1:2,))
    end
end

@testset "ArrayInterface" begin
    @test ArrayInterface.parent_type(da) == Matrix{Int}
end


text:
abstract type AbstractDrift end

get(drift::AbstractDrift) = drift.μ

@with_kw struct ConstDrift <: AbstractDrift
    μ::Real
end

ConstDrift() = ConstDrift(1)

get(drift::ConstDrift) = [drift.μ]

@with_kw struct VarDrift <: AbstractDrift
    μ::Array{Float64,1}
end


text:
module IntegrationTests

using Test

@testset "Quadratures" begin include("QuadraturesTests.jl") end

@testset "TensorProductQuadratures" begin include("TensorProductQuadraturesTests.jl") end

@testset "DuffyQuadratures" begin include("DuffyQuadraturesTests.jl") end

end # module


text:
using ImageMagick
using TestImages
using RegisterQD
using LinearAlgebra
using ImageMetadata
import AxisArrays
using AxisArrays: AxisArray, Axis
using Unitful: μm, mm, cm, km, s

@testset "default_minwidth_rot" begin
    img = rand(3, 10)
    ci = CartesianIndices(img)
    θ = RegisterQD.default_minrot(ci)
    @test θ ≈ 0.01 rtol=0.1
    θ = RegisterQD.default_minrot(ci, [1 0; 0 2])
    @test θ ≈ 0.005 rtol=0.1
    θ = RegisterQD.default_minrot(ci, [3 0; 0 1])
    @test θ ≈ 0.007 rtol=0.1
    θ = RegisterQD.default_minrot(ci, [10 0; 0 1])
    @test θ ≈ 0.01/3 rtol=0.1
    img = rand(3, 10, 5)
    ci = CartesianIndices(img)
    θ = RegisterQD.default_minrot(ci)
    @test θ ≈ 0.1/sqrt(3^2 + 10^2 + 5^2) rtol=1e-3
end

@testset "getSD" begin
    #test that getSD deals with arbitrary dimensions
    A2 = rand(10,10)
    @test getSD(A2) == I

    A3 = rand(10,10,10)
    @test getSD(A3) == I

    A5 = rand(10,10,10,10,10)
    @test getSD(A5) == I

    Ax3 = AxisArray(A3, 1:1:10, 1:2:20, 1:3:30)
    @test getSD(Ax3) == Diagonal([1.0,2.0,3.0])

    #test that getSD works with test images
    mri = testimage("mri-stack.tif")
    @test getSD(mri) == Diagonal([1.0, 1.0, 5.0])

    #test that getSD deals with images with arbitrary space directions
    skewed = ImageMeta(rand(10,10,10))
    skewmatrix = rand(3,3)
    skewed.spacedirections = (Tuple(skewmatrix[1,:]), Tuple(skewmatrix[2,:]), Tuple(skewmatrix[3,:]))
    getSD(skewed) == skewmatrix

    #test that getSD can reconcile units of different magnitudes
    badsampling = AxisArray(rand(10,10,10), (:x,:y,:z), (1mm, 2km, 3.4cm))
    badsampling = ImageMeta(badsampling)
    @test getSD(badsampling) == Diagonal([1, 2e6, 34])

    #test that getSD ignores the time-axis
    timedarray = AxisArray(rand(10,10,10), (:x, :y, :time), (1μm, 1μm, 1s))
    @test size(getSD(timedarray)) == (2,2)
end

#TODO add a testset for other support functions
#rotations
#arrayscale
#pscale
#restrict


text:
@testset "Symmetry Elements" begin
    for t in allofem
        @test Symmetry.pg_to_symels(t) == eval(Meta.parse(t*"s"))
    end
end


text:
"""
# module cyk

- Julia version:
- Author: bramb
- Date: 2019-02-20

# Examples

```jldoctest
julia>using HyperGraphTools.CYK
julia>g = Grammar()
julia>add_rule!(g, Rule("S",("NP","VP"))) # S -> NP VP
julia>add_rule!(g, Rule("VP",("VP","PP")))
julia>add_rule!(g, Rule("VP",("V","NP")))
julia>add_rule!(g, Rule("VP","eats"))
julia>add_rule!(g, Rule("PP",("P","NP")))
julia>add_rule!(g, Rule("NP",("DET","N")))
julia>add_rule!(g, Rule("NP","she"))
julia>add_rule!(g, Rule("V","eats"))
julia>add_rule!(g, Rule("P","with"))
julia>add_rule!(g, Rule("N","fish"))
julia>add_rule!(g, Rule("N","fork"))
julia>add_rule!(g, Rule("DET","a"))

julia>validate(["she", "eats", "a", "fish", "with", "a", "fork"],g)
true

julia>validate(["with", "a", "fork", "eats", "she", "a", "fish"],g)
false
```
"""
module CYK

export
    Rule,
    Grammar,
    add_rule!,
    validate

DEBUG = false

struct Rule
    lhs::String
    rhs::Union{String,Tuple{String,String}} # either a single terminal or a tuple of 2 nonterminals
end
# struct Rule{X}
#     lhs::X
#     rhs::Union{X,Tuple{X,X}}
# end

struct Grammar
  nonTerminals::Array{String,1}
  rules::Array{Rule,1}
  ruleIndex::Dict{Int,Array{Int,1}}

  function Grammar()
      new([],[],Dict{Int64,Array{Int64,1}}())
  end
end
# struct Grammar{X}
#   nonTerminals::Array{X,1}
#   rules::Array{Rule{X},1}
#   lhs_index::Dict{Int,Array{Int,1}}
#
# #   function Grammar()
# #       new{X}([],[],Dict{Int,Array{Int}})
# #   end
# end

_is_nonterminal(s::String) = all(c->isuppercase(c), s)

function non_terminals(r::Rule)
    ntSet = []
    push!(ntSet,r.lhs) # by definition, the lhs is a nonterminal

    potential_nonterminals = []
    if (typeof(r.rhs) == String)
        push!(potential_nonterminals,r.rhs)
    elseif (typeof(r.rhs) == Tuple{String,String})
        push!(potential_nonterminals, r.rhs[1])
        push!(potential_nonterminals, r.rhs[2])
    end
    for e in potential_nonterminals
        if _is_nonterminal(e)
            push!(ntSet,e)
        end
    end
    return collect(ntSet)
end

function add_rule!(g::Grammar, r::Rule)
    for nt in non_terminals(r)
        if !(nt in g.nonTerminals)
            push!(g.nonTerminals,nt)
        end
    end

    push!(g.rules,r)

    lhs_index = indexin([r.lhs],g.nonTerminals)[1]
    if (!haskey(g.ruleIndex,lhs_index))
        g.ruleIndex[lhs_index] = []
    end
    push!(g.ruleIndex[lhs_index],length(g.rules));
end

function rhs_matches_terminal(rhs, terminal)
    rhs == terminal
end

function lhs_index(grammar, terminal)
    index = []
    for rule in grammar.rules
        if rhs_matches_terminal(rule.rhs,terminal)
            push!(index,indexin([rule.lhs],grammar.nonTerminals)[1])
        end
    end
    if (length(index) == 0)
        println(stderr, "ERROR: No rule found for terminal $terminal")
    end
    return index
end

function index_triples(grammar)
    triples = []
    for a in keys(grammar.ruleIndex)
        for ri in grammar.ruleIndex[a]
            rhs = grammar.rules[ri].rhs
            if (typeof(rhs) == Tuple{String,String})
                b = indexin([rhs[1]],grammar.nonTerminals)[1]
                c = indexin([rhs[2]],grammar.nonTerminals)[1]
                push!(triples,(a,b,c))
            end
        end
    end
    return sort(triples)
end

function validate(tokens, grammar)
    # grammar rules need to be in chomsky normal form:
    # https://en.wikipedia.org/wiki/Chomsky_normal_form
    # TODO: check the grammar rules are in CNF

    # algorithm based on the pseudocode from https://en.wikipedia.org/wiki/CYK_algorithm

    # let the input be a string I consisting of n characters: a1 ... an.
    n = length(tokens)
    DEBUG && println("n=$n")

    # let the grammar contain r nonterminal symbols R1 ... Rr, with start symbol R1.
    r = length(grammar.nonTerminals)
    DEBUG && println("r=$r")

    # let P[n,n,r] be an array of booleans. Initialize all elements of P to false.
    P = Array{Bool,3}[]
    P = [false for x in 1:n, y in 1:n, z in 1:r]

    # for each s = 1 to n
    for s in 1:n
    #   for each unit production Rv → as
        for v in lhs_index(grammar,tokens[s])
    #     set P[1,s,v] = true
            P[1,s,v] = true
            DEBUG && println("P[1,$s,$v] = true: token $s:$(tokens[s]) handled by rule $(grammar.nonTerminals[v]) -> ...")
         end
    end

    triples = index_triples(grammar)
    # for each l = 2 to n -- Length of span
    for length in 2:n
    #   for each s = 1 to n-l+1 -- Start of span
        DEBUG && println("length=$length")
        for s in 1:(n-length+1)
            DEBUG && println("  start=$s: \"$(join(tokens[s:(s+length-1)]))\"")
    #       for each p = 1 to l-1 -- Partition of span
            for p in 1:(length-1)
                DEBUG && println("    partition=$p: \"$(join(tokens[s:s+p-1]))\" + \"$(join(tokens[(s+p):s+length-1]))\"")
    #           for each production Ra  → Rb Rc
                for (a,b,c) in triples
    #         if P[p,s,b] and P[l-p,s+p,c] then set P[l,s,a] = true
                    DEBUG && println("      (a,b,c)=($a,$b,$c)=$(grammar.nonTerminals[a])->$(grammar.nonTerminals[b]) $(grammar.nonTerminals[c]), P[p,s,b]=P[$p,$s($(tokens[s])),$b($(grammar.nonTerminals[b]))]=$(P[p,s,b]), P[length-p,s+p,c]=P[$(length-p),$(s+p)($(tokens[s+p])),$c($(grammar.nonTerminals[c]))]=$(P[length-p,s+p,c])")
                    if (P[p,s,b] && P[length-p,s+p,c])
                        P[length,s,a] = true
                        DEBUG && println("      P[$length,$s($(tokens[s])),$a($(grammar.nonTerminals[a]))]=>true")
                    end
                end
            end
        end
    end
    # if P[n,1,1] is true then
    #   I is member of language
    # else
    #   I is not member of language
    return P[n,1,1]
end

function show(r::Rule)
    println("$(r.lhs) -> $(r.rhs)")
end

end

max_stars_repo_name: brambg/HyperGraphTools

max_issues_repo_name: brambg/HyperGraphTools

max_forks_repo_name: brambg/HyperGraphTools


text:
using Statistics: mean

struct LinearModel
    X::Array
    Y::Vector
    beta::Vector #FIMXE: of Real
end

# OLS estimation using (\)
# https://github.com/giob1994/Alistair.jl/blob/3a11c19150169695581b46e4d1895f0641a4c29d/src/linregress.jl#L38-L41
function ols(X::Array, Y::Array)::Array
    return (X' * X) \ (X' * Y)
end

"""Fit ordinary linear regression for observations."""
function ols(sample::Sample; intercept::Bool=false)::LinearModel
    X = (intercept ? add_intercept : identity)(sample.X)
    beta_hat = ols(X, sample.Y)
    return LinearModel(X, sample.Y, beta_hat)
end

# another interface
function ols(x, y; intercept::Bool=false)::LinearModel
    return ols(Sample(x,y), intercept=intercept)
end

"""Return fitted dependent variable Y."""
function yhat(lm::LinearModel)
    return lm.X * lm.beta
end

sum_of_squares(x::Vector)::Real = sum(x .^ 2)

"""Residual sum of squares."""
rss(lm::LinearModel) = sum_of_squares(yhat(lm) - lm.Y)

"""Total sum of squares for Y - mean.
   Also equals var(Y)*n.
"""
tss(lm::LinearModel) = sum_of_squares(lm.Y .- mean(lm.Y)) #

"""R2 (unadjusted) = 1-(RSS/TSS)"""
r2(lm::LinearModel) = 1 - rss(lm)/tss(lm)

betas(lm::LinearModel) = lm.beta

# display LinearModel
function has_intercept(lm::LinearModel)
    return all(lm.X[:,1] .== 1.0)
end

function equation(lm::LinearModel):: String
    equation(has_intercept(lm), lm.beta)
end

function pretty_round(beta, precision::Int=3)
    map(x -> round(x, digits=precision), beta)
end

"""Make string like 3 + 0.5*X1 - 0.1*X2"""
function equation(has_intercept, betas)
    fmt = b -> (b >= 0) ? (" + $b") : (" - $(abs(b))")
    vars = ["*X$i" for i in 1:length(betas)]
    if has_intercept
        vars = pushfirst!(vars[1:end-1], "")
    end
    gen = zip(map(fmt, pretty_round(betas)), vars)
    return "Y =" * join(["$k$x" for (k, x) in gen], "")
end

function desc(lm::LinearModel)::String
    by_line(args...) = join(args, "\n")
    quack_ = has_intercept(lm) ? " " : " no "
    return by_line("Linear model with$(quack_)intercept: $(equation(lm))",
                   "Coefficients: $(lm.beta)",
                   "Sorry I do not know p-values yet.",
                   "R-squared: $(pretty_round(r2(lm)))")
end

show(lm::LinearModel) = println(desc(lm))


text:
module Taarruz
using Knet
using LinearAlgebra

_etype = gpu() >= 0 ? Float32 : Float64
_atype = gpu() >= 0 ? KnetArray{_etype} : Array{_etype}

"""
    FloatArray
A KnetArray or Array. Its elements should be Float32 or Float64.
"""
const FloatArray = Union{KnetArray{F}, Array{F}} where F <: AbstractFloat

include("attacks.jl"); export FGSM, IterativeFGSM, MomentumIterativeFGSM
include("lenet.jl"); export Lenet

end # module


text:
using Plots
using DataFrames
using JSON
using LightOSM
using DataStructures

results = JSON.parsefile(joinpath(@__DIR__, "benchmark_results.json"))

time_units_mapping = Dict(
    "s" => 1,
    "ms" => 1e-3,
    "μs" => 1e-9
)

memory_units_mapping = Dict(
    "MiB" => 1,
    "KiB" => 1e-3,
    "GiB" => 1e+3,
)

df_dict_display = DefaultDict(Dict)
df_dict_time = DefaultDict(Dict)
df_dict_memory = DefaultDict(Dict)

for (algorithm, r) in results
    n_paths = collect(keys(r))
    n_paths_numeric = parse.(Float64, n_paths)

    df_dict_display[algorithm] = DefaultDict(Any[], "No. Paths" => n_paths)
    df_dict_time[algorithm] = DefaultDict(Number[], "No. Paths" => n_paths_numeric)
    df_dict_memory[algorithm] = DefaultDict(Number[], "No. Paths" => n_paths_numeric)

    for (n, packages) in r
        for (p, items) in packages
            mean_time = items["mean time"]
            gc_match = match(r"\((.*?)\)", mean_time)
            gc = LightOSM.remove_non_numeric(gc_match[1]) # %
            mean_time = replace(mean_time, gc_match.match => "")
            mean_time_units = strip(LightOSM.remove_numeric(mean_time))
            mean_time_units_factor = time_units_mapping[mean_time_units]
            mean_time = LightOSM.remove_non_numeric(mean_time) * mean_time_units_factor # s

            allocs = items["allocs estimate"]
            allocs = LightOSM.remove_non_numeric(allocs)

            memory = items["memory estimate"]
            memory_units = strip(LightOSM.remove_numeric(memory))
            memory_units_factor = memory_units_mapping[memory_units]
            memory = LightOSM.remove_non_numeric(memory) * memory_units_factor # MiB

            metrics = "$(round(mean_time, digits=4))s ($(round(memory, digits=2)) MiB)"
            push!(df_dict_display[algorithm][p], metrics)
            push!(df_dict_time[algorithm][p], mean_time)
            push!(df_dict_memory[algorithm][p], memory)

        end
    end
end

# Summary Tables
dijkstra_results = DataFrame(;[Symbol(k) => v for (k, v) in df_dict_display["dijkstra"]]...)
sort!(dijkstra_results, "No. Paths")
dijkstra_results = dijkstra_results[:, [Symbol("No. Paths"), :losm, :losm_precompute, :osmx, :lg]]
astar_results = DataFrame(;[Symbol(k) => v for (k, v) in df_dict_display["astar"]]...)
sort!(astar_results, "No. Paths")
astar_results = astar_results[:, [Symbol("No. Paths"), :losm, :losm_precompute, :osmx, :lg]]

# Dijsktra Time Comparison
dijkstra_df_time = DataFrame(;[Symbol(k) => v for (k, v) in df_dict_time["dijkstra"]]...)
sort!(dijkstra_df_time, "No. Paths")

dijkstra_time_x = dijkstra_df_time["No. Paths"]
dijkstra_time_y = Matrix(select!(dijkstra_df_time, Not(Symbol("No. Paths"))))
dijkstra_time_labels = permutedims(names(dijkstra_df_time))
plot_dijkstra_time() = plot(dijkstra_time_x, dijkstra_time_y, title="Dijkstra Shortest Path - Runtime (s) vs. No. Paths", label=dijkstra_time_labels, ylabel="Runtime (s)", xlabel="No. Paths", lw=3, fmt=:png)

# Dijsktra Memory Comparison
dijkstra_df_memory = DataFrame(;[Symbol(k) => v for (k, v) in df_dict_memory["dijkstra"]]...)
sort!(dijkstra_df_memory, "No. Paths")

dijkstra_memory_x = dijkstra_df_memory["No. Paths"]
dijkstra_memory_y = Matrix(select!(dijkstra_df_memory, Not(Symbol("No. Paths"))))
dijkstra_memory_labels = permutedims(names(dijkstra_df_memory))
plot_dijkstra_memory() = plot(dijkstra_memory_x, dijkstra_memory_y, title="Dijkstra Shortest Path - Memory (MiB) vs. No. Paths", label=dijkstra_memory_labels, ylabel="Memory (MiB)", xlabel="No. Paths", lw=3, fmt=:png)

# A* Time Comparison
astar_df_time = DataFrame(;[Symbol(k) => v for (k, v) in df_dict_time["astar"]]...)
sort!(astar_df_time, "No. Paths")

astar_time_x = astar_df_time["No. Paths"]
astar_time_y = Matrix(select!(astar_df_time, Not(Symbol("No. Paths"))))
astar_time_labels = permutedims(names(astar_df_time))
plot_astar_time() = plot(astar_time_x, astar_time_y, title="A* Shortest Path - Runtime (s) vs. No. Paths", label=astar_time_labels, ylabel="Runtime (s)", xlabel="No. Paths", lw=3, fmt=:png)

# A* Memory Comparison
astar_df_memory = DataFrame(;[Symbol(k) => v for (k, v) in df_dict_memory["astar"]]...)
sort!(astar_df_memory, "No. Paths")

astar_memory_x = astar_df_memory["No. Paths"]
astar_memory_y = Matrix(select!(astar_df_memory, Not(Symbol("No. Paths"))))
astar_memory_labels = permutedims(names(astar_df_memory))
plot_astar_memory() = plot(astar_memory_x, astar_memory_y, title="A* Shortest Path - Memory (MiB) vs. No. Paths", label=astar_memory_labels, ylabel="Memory (MiB)", xlabel="No. Paths", lw=3, fmt=:png)


text:
# utilities

function elMax!{T<:Real}(x::Array{T},y::T)
  x[:] = max.(x,y)
end

function elMin!{T<:Real}(x::Array{T},y::T)
  x[:] = min.(x,y)
end

function resetArray!{T<:Real}(A::Array{T})
  A[:] = zero(eltype(A))
end

function logSumExpWeightsNorm(log_w::Vector{T} where T <: AbstractFloat)
  # function uses log-sum-exp trick to compute normalized weights of a
  # vector of numbers stored as logs, avoiding underflow
  max_entry = maximum(log_w)

  shifted_weights = log_w - max_entry
  return normalized_weights = exp.(shifted_weights)./sum(exp.(shifted_weights))

end

function logSumExpWeightsNorm!(p::Vector{T},log_w::Vector{T}) where T <: AbstractFloat
  # function uses log-sum-exp trick to compute normalized weights of a
  # vector of numbers stored as logs, avoiding underflow
  max_entry = maximum(log_w)

  shifted_weights = log_w - max_entry
  p[:] = exp.(shifted_weights)./sum(exp.(shifted_weights))

end

function logSumExpWeights(log_w::Vector{T} where T <: AbstractFloat)
  # function uses log-sum-exp trick to compute UNnormalized weights of a
  # vector of numbers stored as logs, avoiding underflow
  max_entry = maximum(log_w)

  shifted_weights = exp.(log_w .- max_entry)

end

function logSumExpWeights!(lse_w::Vector{T},log_w::Vector{T}) where T <: AbstractFloat
  # function uses log-sum-exp trick to compute UNnormalized weights of a
  # vector of numbers stored as logs, avoiding underflow
  max_entry = maximum(log_w)

  lse_w[:] = exp.(log_w .- max_entry)

end

function OneHot(n::Int64,k::Int64)

  onehotvector = spzeros(n)
  onehotvector[k] = 1
  return onehotvector

end

function stratifiedResample(weights::Vector{Float64},n_resample::Int64)::Vector{Int64}
  sample_idx = zeros(Int64,n_resample)
  stratifiedResample!(sample_idx,weights,n_resample)
  return sample_idx
end

function stratifiedResample!(sample_idx::Vector{Int64},weights::Vector{Float64},n_resample::Int64)
"""
 Implements stratified resampling as described in Appendix B of
 Fearnhead and Clifford (2003), JRSSB 65(4) 887--899

 weights is a vector of weights; n_resample is the number of resampled
 indices to return
"""
  norm_weights = cumsum(weights./sum(weights))::Vector{Float64} # normalize
  buckets = collect(0.0:(1.0/n_resample):1.0)[1:n_resample]
  locs = rand(n_resample)./n_resample .+ buckets

  sample_idx[:] = [ findfirst( locs[i] .<= norm_weights ) for i in 1:n_resample ]
  # return sample_idx
end

function systematicResample(weights::Vector{Float64},n_resample::Int64)::Vector{Int64}
  sample_idx = zeros(Int64,n_resample)
  systematicResample!(sample_idx,weights,n_resample)
  return sample_idx
end

function systematicResample!(sample_idx::Vector{Int64},weights::Vector{Float64},n_resample::Int64)

  norm_weights = cumsum(weights./sum(weights))
  buckets = collect(0.0:(1.0/n_resample):1.0)[1:n_resample]
  U = rand()*1.0/n_resample
  locs = (buckets .+ U)

  sample_idx[:] = [ findfirst( locs[i] .<= norm_weights ) for i in 1:n_resample ]
  # return sample_idx
end

function tally_ints(Z::Vector{Int},K::Int)
    """
    counts occurrences in `Z` of integers 1 to `K`
    - `Z`: vector of integers
    - `K`: maximum value to count occurences in `Z`
    """
    ret = zeros(Int,K)
    n = size(Z,1)
    idx_all = 1:n
    idx_j = trues(n)
    for j in 1:K
      for i in idx_all[idx_j]
        if Z[i]==j
          ret[j] += 1
          idx_j[i] = false
        end
      end
    end
    return ret
end

function adj2edgelist(A::T where T<:Union{Array{Int64,2},SparseMatrixCSC{Int64,Int64},SparseMatrixCSC{Float64,Int64}})::Array{Int64,2}
  adj2edgelist(A,false)
end

function adj2edgelist(A::T where T<:Union{Array{Int64,2},SparseMatrixCSC{Int64,Int64},SparseMatrixCSC{Float64,Int64}},multigraph::Bool)::Array{Int64,2}
"""
  Convert binary adjacency matrix A to edge list
"""

  M,N,V = findnz(tril(A))
  if multigraph
    return hcat(M,N),V
  else
    return hcat(M,N)
  end
end

function edgelist2adj(edgelist::Array{Int64,2})::SparseMatrixCSC{Int64,Int64}
"""
  Convert edge list to binary adjacency matrix
"""

  edgelist2adj(edgelist,Int64)
end

function edgelist2adj(edgelist::Array{Int64,2},value_type::DataType)
"""
  Convert edge list to binary adjacency matrix, with non-zero values of type value_type
"""
  nv = maximum(edgelist)
  A = sparse([edgelist[:,1];edgelist[:,2]],[edgelist[:,2];edgelist[:,1]],ones(value_type,2*size(edgelist,1)),nv,nv)
end

function edgelist2adj!(A::Union{Array{Float64,2},SparseMatrixCSC{Int64,Int64},SparseMatrixCSC{Float64,Int64}},edge_list::Array{Int64,2})
"""
  Convert edge list to binary adjacency matrix, with non-zero values of type value_type
"""
  et = eltype(A)
  for i=1:size(edge_list,1)
    A[edge_list[i,1],edge_list[i,2]] = one(et)
    A[edge_list[i,2],edge_list[i,1]] = one(et)
  end

end

function normalizedLaplacian(edgelist::Array{Int64,2})::SparseMatrixCSC{Int64,Int64}
"""
  Construct (sparse) normalized Laplacian matrix
"""
  degrees = getDegrees(edgelist)
  nv = maximum(edgelist)
  L = sparse(1*I,nv,nv) - sparse(Diagonal( degrees.^(-1/2) )) * edgelist2adj(edgelist,Float64) * sparse(Diagonal( degrees.^(-1/2) ))
end

function denseNormalizedLaplacian!(L::Array{Float64,2},adj::SparseMatrixCSC{Float64,Int64},degrees::Array{Float64,1},nv::Int64)
"""
  Construct a symmetric (dense) normalized Laplacian matrix
"""
  # degrees = getDegrees(edgelist)
  # nv = maximum(edgelist)
  L[1:nv,1:nv] = Array(sparse(1*I,nv,nv) - sparse(Diagonal( degrees.^(-1/2) )) * adj * sparse(Diagonal( degrees.^(-1/2) )))
end

function denseNormalizedLaplacian!(L::Array{Float64,2},adj::SparseMatrixCSC{Int,Int64},degrees::Array{Float64,1},nv::Int64)
"""
  Construct a symmetric (dense) normalized Laplacian matrix
"""
  # degrees = getDegrees(edgelist)
  # nv = maximum(edgelist)
  L[1:nv,1:nv] = Array(sparse(1*I,nv,nv) - sparse(Diagonal( degrees.^(-1/2) )) * adj * sparse(Diagonal( degrees.^(-1/2) )))
end

function denseNormalizedLaplacian!(L::Array{Float64,2},edge_list::Array{Int64,2},degrees::Array{Float64,1},nv::Int64)
"""
  Construct a symmetric (dense) normalized Laplacian matrix. L should be all zeros.
"""
  # degrees = getDegrees(edgelist)
  # nv = maximum(edgelist)

  # set diagonal
  for n = 1:nv
    L[n,n] = 1.0
  end

  rt_deg = degrees.^(-1/2)
  for m = 1:size(edge_list,1)
    L[edge_list[m,1],edge_list[m,2]] = -(rt_deg[edge_list[m,1]]*rt_deg[edge_list[m,2]])
    L[edge_list[m,2],edge_list[m,1]] = -(rt_deg[edge_list[m,1]]*rt_deg[edge_list[m,2]])
  end

  # L[1:nv,1:nv] = Array(sparse(1*I,nv,nv) - sparse(Diagonal( degrees.^(-1/2) )) * adj * sparse(Diagonal( degrees.^(-1/2) )))
end

function getDegrees(edgelist::Array{Int64,2})::Array{Float64,1}
"""
  Compute degrees of edge list
"""
  return convert(Array{Float64,1},tally_ints(edgelist[:],maximum(edgelist)))
end

function generateEigenSystem!(L::Array{Float64,2},nv::Int64)::Tuple{Array{Float64,1},Array{Float64,2}}
  LL = Symmetric(L[1:nv,1:nv])
  eigvec_tr = zeros(Float64,nv,nv)
  # THIS OVERWRITES L!
  esys_val,esys_vec = LAPACK.syevr!('V','A',LL.uplo,LL.data,-0.0,0.0,0,0,-1.0)
  transpose!(eigvec_tr,esys_vec) # transposed to optimize computation

  return esys_val,eigvec_tr
end

function generateEigenSystem(L::Array{Float64,2},nv::Int64)::Tuple{Array{Float64,1},Array{Float64,2}}
  Lcopy = copy(L)
  LL = Symmetric(Lcopy[1:nv,1:nv])
  eigvec_tr = zeros(Float64,nv,nv)
  # THIS OVERWRITES Lcopy!
  esys_val,esys_vec = LAPACK.syevr!('V','A',LL.uplo,LL.data,-0.0,0.0,0,0,-1.0)
  transpose!(eigvec_tr,esys_vec) # transposed to optimize computation

  return esys_val,eigvec_tr
end

function updateEigenSystem!(L::Array{Float64,2},nv::Int64,p_state::ParticleState)

  denseNormalizedLaplacian!(L,p_state.vertex_unmap[p_state.edge_list],p_state.degrees,nv)
  LL = Symmetric(L[1:nv,1:nv])
  # THIS OVERWRITES L!
  esys_val,esys_vec = LAPACK.syevr!('V','A',LL.uplo,LL.data,-0.0,0.0,0,0,-1.0)
  # esys = eigfact!(Symmetric(L[1:nv,1:nv]))
  esys_val[1] = zero(Float64) # for numerical stability
  esys_val[end] = min(esys_val[end],Float64(2.0)) # for numerical stability

  for i=1:nv
    p_state.eig_vals[i] = esys_val[i]
    for j=1:nv
      # THIS TRANSPOSES THE EIGENVECTOR MATRIX TO OPTIMIZE COMPUTATION
      # ROWS OF p_state.eig_vecs CORRESPOND TO EIGENVECTORS OF L
      p_state.eig_vecs[j,i] = esys_vec[i,j]
    end
  end
  p_state.has_eigensystem[:] = true

end

function saveSamples!(alpha_samples::Array{Float64,1},lambda_samples::Array{Float64,1},
                      particle_trajectory_samples::Array{Int64,2},edge_sequence_samples::Array{Int64,2},
                      s_state::SamplerState,particle_container::Array{Array{ParticleState,1},1},
                      p_idx::Int64,n_sample::Int64)

  alpha_samples[n_sample] = s_state.α[1]
  lambda_samples[n_sample] = s_state.λ[1]
  particle_trajectory_samples[n_sample,:] = s_state.particle_path[:]
  edge_sequence_samples[n_sample,:] = particle_container[p_idx][end].edge_idx_list[:]

end


text:
@safetestset "isbtcaddress" begin
    using Validator

    @test Validator.isbtcaddress("1MUz4VMYui5qY1mxUiG8BQ1Luv6tqkvaiL") === true
    @test Validator.isbtcaddress("3J98t1WpEZ73CNmQviecrnyiWrnqRhWNLy") === true
    @test Validator.isbtcaddress("bc1qar0srrr7xfkvy5l643lydnw9re59gtzzwf5mdq") === true
    @test Validator.isbtcaddress("14qViLJfdGaP4EeHnDyJbEGQysnCpwk3gd") === true
    @test Validator.isbtcaddress("35bSzXvRKLpHsHMrzb82f617cV4Srnt7hS") === true
    @test Validator.isbtcaddress("17VZNX1SN5NtKa8UQFxwQbFeFc3iqRYhemt") === true
    @test Validator.isbtcaddress("bc1qw508d6qejxtdg4y5r3zarvary0c5xw7kv8f3t4") === true

    @test Validator.isbtcaddress("4J98t1WpEZ73CNmQviecrnyiWrnqRhWNLy") === false
    @test Validator.isbtcaddress("0x56F0B8A998425c53c75C4A303D4eF987533c5597") === false
    @test Validator.isbtcaddress("pp8skudq3x5hzw8ew7vzsw8tn4k8wxsqsv0lt0mf3g") === false
    @test Validator.isbtcaddress("17VZNX1SN5NlKa8UQFxwQbFeFc3iqRYhem") === false
    @test Validator.isbtcaddress("BC1QW508D6QEJXTDG4Y5R3ZARVAYR0C5XW7KV8F3T4") === false
end


text:
module ValidateRCall_ape
using Compat.Test
using Compat: @warn
using Phylo
using RCall

# Create a temporary directory to work in
libdir = mktempdir();

global skipR = false
if !rcopy(R"require(ape)")
    rcall(Symbol(".libPaths"), libdir);
    reval("install.packages(\"ape\", lib=\"$libdir\", " *
          "repos=\"http://cran.r-project.org\")");
    global skipR = !rcopy(R"require(ape, lib.loc=c(\"$libdir\", .libPaths()))") &&
        !mustCrossvalidate;
    skipR && @warn "ape R package not installed and would not install, " *
        "skipping R crossvalidation"
end

if !skipR
    # Run tree comparisons on increasing numbers of tips
    @testset "RCall - testing Phylo vs ape" begin
        @testset "Testing with R rtree($i)" for i in 5:5:50
            rt = rcall(:rtree, i)
            jt = rcopy(NamedTree, rt)
            jl = Set(getleafnames(jt))
            rl = Set(rcopy(rcall(Symbol("[["), rt, "tip.label")))
            @test jl == rl
            rt2 = RObject(jt)
            @test rcopy(rcall(Symbol("all.equal"), rt, rt2))
            @rput jt
            reval("jt2=jt");
            @rget jt2
            jt3=jt2;
            @rput jt3
            @rput rt
            @test rcopy(rcall(Symbol("all.equal"), rt, jt3))
        end

        @testset "Testing with julia rand(Nonultrametric($i))" for i in 5:5:50
            jt = rand(Nonultrametric(i))
            rt = RObject(jt)
            @test Set(getleafnames(jt)) ==
                Set(rcopy(rcall(Symbol("[["), rt, "tip.label")))
            jt2 = rcopy(NamedTree, rt)
            @test getleafnames(jt) == getleafnames(jt2)
            @test rcopy(rcall(Symbol("all.equal"), rt, RObject(jt2)))
        end

        @testset "Testing reading in newick trees from disk" begin
            jtree = open(parsenewick, Phylo.path("H1N1.newick"))
            rtree = rcall(Symbol("read.tree"), Phylo.path("H1N1.newick"))
            @test rcopy(rcall(Symbol("all.equal"), jtree, rtree))
        end

        @testset "Testing reading in nexus trees from disk" begin
            jts = open(parsenexus, Phylo.path("H1N1.trees"))
            rtree1 = R"read.nexus($(Phylo.path(\"H1N1.trees\")))$TREE1"
            jtree1 = jts["TREE1"]
            @test rcopy(rcall(Symbol("all.equal"), jtree1, rtree1))
            @test "H1N1_A_MIYAGI_3_2000" ∈ nodenameiter(jtree1)
            @test collect(keys(first(treeinfoiter(jts)))) == ["lnP"]
        end
    end
end
rm(libdir, force=true, recursive=true);

end


text:
function fixargs(args::Vector{Any}, stub::Any)
    n = length(args)
    data = Array(Any, n)
    na = BitArray(n)
    for i in 1:n
        if args[i] == :NA
            data[i] = stub
            na[i] = true
        else
            data[i] = args[i]
            na[i] = false
        end
    end
    return data, na
end

# We assume that data has at least one "value" that isn't NA
function findstub_vector(ex::Expr)
    if ex.args[1] != :NA
        return ex.args[1]
    end
    for i in 2:length(ex.args)
        if ex.args[i] != :NA
            return ex.args[i]
        end
    end
    return NA
end

# We assume that data has at least one "value" that isn't NA
function findstub_matrix(ex::Expr)
    if ex.args[1].args[1] != :NA
        return ex.args[1].args[1]
    end
    nrows = length(ex.args)
    for row in 1:nrows
        subex = ex.args[row]
        for i in 1:length(subex.args)
            if subex.args[i] != :NA
                return subex.args[i]
            end
        end
    end
    return NA
end

function parsevector(ex::Expr)
    if ex.head in (:ref, :typed_hcat, :typed_vcat)
        data, na = fixargs(ex.args[2:end], :(zero($(ex.args[1]))))
        return Expr(ex.head, ex.args[1], data...), ex.head == :typed_hcat ? na' : na
    else
        stub = findstub_vector(ex)
        data, na = fixargs(ex.args, stub)
        if ex.head == :hcat
            na = na'
        end

        if isequal(stub, NA)
            return Expr(ex.head == :hcat ? (:typed_hcat) : (:ref), Any, data...), na
        else
            return Expr(ex.head, data...), na
        end
    end
end

function parsematrix(ex::Expr)
    if ex.head == :typed_vcat
        stub = :(zero($(ex.args[1])))
        rows = 2:length(ex.args)
    else
        stub = findstub_matrix(ex)
        rows = 1:length(ex.args)
    end

    nrows = length(rows)
    datarows = Array(Expr, nrows)
    narows = Array(Expr, nrows)
    for irow in 1:nrows
        data, na = fixargs(ex.args[rows[irow]].args, stub)
        datarows[irow] = Expr(:row, data...)
        narows[irow] = Expr(:row, na...)
    end
    if ex.head == :typed_vcat
        return Expr(:typed_vcat, ex.args[1], datarows...), Expr(:vcat, narows...)
    elseif isequal(stub, NA)
        return Expr(:typed_vcat, Any, datarows...), Expr(:vcat, narows...)
    else
        return Expr(:vcat, datarows...), Expr(:vcat, narows...)
    end
end

function parsedata(ex::Expr)
    if length(ex.args) == 0
        return :([]), Expr(:call, :Array, :Bool, 0)
    end
    if ex.head == :typed_vcat || (isa(ex.args[1], Expr) && ex.args[1].head == :row)
        return parsematrix(ex)
    else
        return parsevector(ex)
    end
end

macro data(ex)
    if !(ex.head in (:vect, :vcat, :hcat, :ref, :typed_vcat, :typed_hcat))
        return quote
            tmp = $(esc(ex))
            DataArray(tmp, bitbroadcast(x->isequal(x, NA), tmp))
        end
    end
    dataexpr, naexpr = parsedata(ex)
    return Expr(:call, :DataArray, esc(dataexpr), esc(naexpr))
end

macro pdata(ex)
    if !(ex.head in (:vect, :vcat, :hcat, :ref, :typed_vcat, :typed_hcat))
        return quote
            tmp = $(esc(ex))
            PooledDataArray(tmp, bitbroadcast(x->isequal(x, NA), tmp))
        end
    end
    dataexpr, naexpr = parsedata(ex)
    return Expr(:call, :PooledDataArray, esc(dataexpr), esc(naexpr))
end


text:
"""
    tabylnSQ

Container for the squared spherical harmonics prefactor |Y_l^n|^2.
ATTENTION, harmonic index starts at l = 0.
We pre-compute all the values of |Y_l^n (pi/2,pi/2)|^2 rather having to load any special functions librairies in order to save time.

Use `tabylnSQ!()` to fill it in.
"""
const tabylnSQ = zeros(Float64,lmax+1,lmax+1)

"""
    tabylnSQ!()

Fills the array (of size lmax*lmax) tabylnSQ with the values of the |Y_l^n|^2, the squared spherical harmonics prefactors.
tabylnSQ is indexed by two integers [l,n], with l,n going from 1 to lmax+1 since arrays begin at 1.
tabylnSQ[l+1,n+1] contains the value `|Y_l^n|^2 = |Yln(pi/2,pi/2)|^2`.
Those are computed by recursion.
"""
function tabylnSQ!()
    # Making sure that the array is initially filled with zeros
    for n=1:(lmax+1), l=1:(lmax+1) # ATTENTION, the size of the array is (lmax+1) in both directions
        tabylnSQ[l,n] = 0.0
    end
    #####
    # Computing all the squared elements by recurrence !! ATTENTION, this corresponds to the squared elements
    tabylnSQ[1,1] = 1.0/(4.0*pi) # |Yln(pi/2,pi/2)|^(2) for (l,n)=(0,0) !! HARMONIC INDEX STARTS AT \ell=0
    # First loop to fill in the diagonal elements
    for l=0:(lmax-1) #!! ATTENTION TO THE BOUNDS OF THE LOOP
        tabylnSQ[(l+1)+1,(l+1)+1] = ((2.0*l+3.0)/(2.0*l+2.0))*tabylnSQ[l+1,l+1]
    end
    # Second loop to fill in the remaining terms
    for n=0:lmax # We perform the loop at fixed m
        for l=n:2:(lmax-2) # At fixed n, we scan the \ell by step of 2 !! ATTENTION TO THE BOUNDS OF THE LOOP
            tabylnSQ[(l+2)+1,n+1] = ((2.0*l+5.0)/(2.0*l+1.0))*(((l-n+1.0)*(l+n+1.0))/((l-n+2.0)*(l+n+2.0)))*tabylnSQ[l+1,n+1]
        end
    end
end
##################################################
tabylnSQ!() # Initialising tabylnSQ with the pre-computed coefficients
##################################################
"""
    ylnSQ(l,n)

Wrapped function that returns the squared spherical harmonics prefactor |Y_l^n (pi/2,0)|^2.
Caution: this is the SQUARED quantity.
"""
function ylnSQ(l::Int64,n::Int64)
    return tabylnSQ[l+1,abs(n)+1] # Using the symmetry |Yl,n(pi/2,pi/2)|^2 = |Yl,-n(pi/2,pi/2)|^2
end

##################################################
# Definitions of useful quantities
##################################################

"""
    rg

Schwartzschild (gravitational) radius.
Defined as `(G*mBH)/(cvel^(2))`, where `G` is Newton's constant, `mBH` the mass of the SMBH and `cvel` the light celerity.
"""
const rg = (G*mBH)/(cvel^(2))

"""
    Jc(a)

Angular momentum of a circular orbit with radius a.
"""
function Jc(a::Float64)
    return sqrt(G*mBH*a)
end

"""
    jlc(a)

Reduced angular momentum of the loss-cone for a semi-major axis a.
For j<jlc(a), there are no allowed orbits
"""
function jlc(a::Float64)
    return 4.0*sqrt(rg/a)
end

"""
    ecc_to_j(ecc::Float64)

Converts eccentricity to normalized angular momentum.
"""
function ecc_to_j(ecc::Float64)
    return sqrt(1.0-ecc^2)
end

"""
    uncertainty_ecc_to_j(ecc::Float64,de::Float64)

Converts uncertainty in eccentricity to uncertainty in normalized angular momentum.
"""
function uncertainty_ecc_to_j(ecc::Float64,de::Float64)
    j = ecc_to_j(ecc)
    dj = ecc/j * de
    return dj
end

"""
    nuKep(a)

Fast Keplerian frequency of an orbit with semi-major axis a.
"""
function nuKep(a::Float64)
    return sqrt((G*mBH)/(a^(3)))
end

"""
    nuGR(a,j)

Schwarzschild precession frequency for an orbit of semi-major axis a and reduced angular momentum J/Jc = j.
"""
function nuGR(a::Float64,j::Float64)
    return 3.0*(rg/a)*nuKep(a)/(j^(2))
end

"""
    dnuGRdj(a,j)

j-partial gradient of the Schwarzschild precession frequency for an orbit of semi-major axis a and reduced angular momentum J/Jc = j.
"""
function dnuGRdj(a::Float64,j::Float64)
    return -6.0*(rg/a)*nuKep(a)/(j^(3))
end

"""
    d2nuGRdj2(a,j)

jj-second-order gradient of the Schwarzschild precession frequency for an orbit of semi-major axis a and reduced angular momentum J/Jc = j.
"""
function d2nuGRdj2(a::Float64,j::Float64)
    return 18.0*(rg/a)*nuKep(a)/(j^(4))
end

"""
    LegendreP(n,x)

Wrapped function to evalue the Legendre functions with non-integer n at argument x.
This uses the Hypergeometric function expression `Pn(x) = ₂F₁(-n,n+1.0,1.0,(1.0-x)*0.5) `.
"""
function LegendreP(n::Float64,x::Float64)
    return _₂F₁(-n,n+1.0,1.0,(1.0-x)*0.5)
end

##################################################
# Computation of actual physically relevant functions
##################################################

"""
    nuMbar(M_below_a,a)

SMA-dependent part of themass precession frequency.
Defined as `nuKep*M(<a)/mBH`, where `M(<a)` is the mass enclosed within semi-major axis a and `mBH` the mass of the SMBH.
"""
function nuMbar(M_below_a::Float64,a::Float64)
    return nuKep(a)*(M_below_a)/(mBH) # Output
end

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

"""
    nbjInt

Number of j used for the interpolation of the functions j->hM(j) and its derivatives.
Set as 1000 by default.
"""
const nbjInt = 1000

"""
    jmaxhM

For j>jmaxhM, we use a linear approximation for j->hM(j) and its derivatives.
Set as 0.99 by default.
"""
const jmaxhM = 0.99

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

"""
    _hM(gamma,j,p1,p2)

j-dependent part of the mass precession frequency.
We should take p1 = LegendreP[1-gamma,1/j] and p2 = LegendreP[2-gamma,1/j].
Depends on the cusp index gamma and the reduced angular momentum j.

# Remarks
- We use a quadratic approximation for j>jmaxhM
- This function is not continuous since we use a linear approximation for j>jmaxhM. By default, jmaxhM=0.99.
"""
function _hM(gamma::Float64,j::Float64,p1::Float64,p2::Float64)
    if (j > jmaxhM) # We use the linear approximation
        return (1.0/2.0)*(-3.0+gamma)-
               (1.0/8.0)*(-12.0+gamma+4.0*gamma^(2)-gamma^(3))*(1.0-j)+
               (1.0/96.0)*(-3.0+gamma)*(-1.0+gamma)*gamma*(-16.0+gamma*(3.0+gamma))*(1.0-j)^(2) # Output
    else # We can use the generic expression
        return (j^(4.0-gamma))/(1.0-j^(2))*(p1 - p2/(j)) # Output
    end
end

"""
    _dhMdj(gamma,j,p1,p2)

j-gradient of the j-dependent part of the mass precession frequency.
We should take p1 = LegendreP[1-gamma,1/j] and p2 = LegendreP[2-gamma,1/j].
Depends on the cusp index gamma and the reduced angular momentum j.

# Remarks
- We use a quadratic approximation for j>jmaxhM
- This function is not continuous since we use a linear approximation for j>jmaxhM. By default, jmaxhM=0.99.
"""
function _dhMdj(gamma::Float64,j::Float64,p1::Float64,p2::Float64)
    if (j > jmaxhM) # We use the linear approximation
        return (1.0/8.0)*(-12.0+gamma+4.0*gamma^(2)-gamma^(3))-
               (1.0/48.0)*(-3.0+gamma)*(-1.0+gamma)*gamma*(-16.0+gamma*(3.0+gamma))*(1.0-j)-
               (1.0/768)*((-3.0+gamma)*(-1.0+gamma)*gamma*(72.0+gamma*(-94+gamma*(-13.0+gamma*(10.0+gamma)))))*(1.0-j)^(2) # Output
    else # We can use the generic expression
        return (j^(2.0-gamma))/((1.0-j^(2))^(2))*(j*(2.0+(2.0-gamma)*(1.0-j^(2)))*p1 - (1.0+j^(2))*p2) # Output
    end
end

"""
    _d2hMdj2(gamma,j,p1,p2)

Second-order j-gradient of the j-dependent part of the mass precession frequency.
We should take p1 = LegendreP[1-gamma,1/j] and p2 = LegendreP[2-gamma,1/j].
Depends on the cusp index gamma and the reduced angular momentum j.

# Remarks
- We use a quadratic approximation for j>jmaxhM
- This function is not continuous since we use a linear approximation for j>jmaxhM. By default, jmaxhM=0.99.
"""
function _d2hMdj2(gamma::Float64,j::Float64,p1::Float64,p2::Float64)
    if (j > jmaxhM) # We use the linear approximation
        return (1.0/48.0)*(-3.0+gamma)*(-1.0+gamma)*gamma*(-16.0+gamma*(3.0+gamma))+
               (1.0/384.0)*(-3.0+gamma)*(-1.0+gamma)*gamma*(72.0+gamma*(-94.0+gamma*(-13.0+gamma*(10.0+gamma))))*(1.0-j)+
               (1.0/7680.0)*(-3.0+gamma)*(-1.0+gamma)*gamma*(1.0+gamma)*(2.0+gamma)*(360.0+gamma*(-358.0+gamma*(19.0+gamma*(18.0+gamma))))*(1.0-j)^(2) # Output
    else # We can use the generic expression
        return (j^(2.0-gamma))/((1.0-j^(2))^(3))*((2.0*(3.0-gamma)^(2)-3.0*j^(2)*(4.0-gamma)*(1.0-gamma)+j^(4)*(2.0-gamma)*(1.0-gamma))*p1+
                                                  j*(-12.0+j^(2)*(4.0-gamma)*(1.0-gamma)+(5.0-gamma)*gamma)*p2) # Output
    end
end

"""
    gethMInt(iBath)

Function that returns an interpolation function for the function j->hM(j) for the bath population indexed by iBath in the source file.
hM(j) is the j-dependent part of the mass precession frequency.
"""
function gethMInt(iBath::Int64)
    jminInt = jlc(rh) # Minimum value of j for the interpolation, fixed by the influence radius
    jmaxInt = 1.0 # Maximum value of j for the interpolation
    #####
    rangejInt = range(jminInt,length=nbjInt,jmaxInt) # Range in j over which the interpolation is performed
    tabjInt = collect(rangejInt) # Table of the values of j for which the functions must be pre-computed
    tabhMInt = zeros(Float64,nbjInt) # Container for the values of hM(j)
    ####
    gamma = tabgammaBath[iBath] # Cusp index of the current bath population
    #####
    for indj=1:nbjInt # Loop over the j of the grid
        jloc = tabjInt[indj] # Current value of j
        p1 = LegendreP(1.0-gamma,1.0/(jloc)) # Current value of LegendreP[1-gamma,1/j]
        p2 = LegendreP(2.0-gamma,1.0/(jloc)) # Current value of LegendreP[2-gamma,1/j]
        tabhMInt[indj] = _hM(gamma,jloc,p1,p2) # Pre-computing the value of hM(j)
    end
    #####
    inthM = Interpolations.scale(interpolate(tabhMInt, BSpline(Cubic(Line(OnGrid())))),rangejInt) # Constructing the interpolation function for j->hM(j)
    #####
    return inthM # Returning the interpolation function
end

"""
    getdhMdjInt(iBath)

Function that returns an interpolation function for the function j->dhMdj(j) for the bath population indexed by iBath in the source file.
dhMdj(j) is the j-gradient of the j-dependent part of the mass precession frequency.
"""
function getdhMdjInt(iBath::Int64)
    jminInt = jlc(rh) # Minimum value of j for the interpolation, fixed by the influence radius
    jmaxInt = 1.0 # Maximum value of j for the interpolation
    #####
    rangejInt = range(jminInt,length=nbjInt,jmaxInt) # Range in j over which the interpolation is performed
    tabjInt = collect(rangejInt) # Table of the values of j for which the functions must be pre-computed
    tabdhMdjInt = zeros(Float64,nbjInt) # Container for the values of dhMdj(j)
    ####
    gamma = tabgammaBath[iBath] # Cusp index of the current bath population
    #####
    for indj=1:nbjInt # Loop over the j of the grid
        jloc = tabjInt[indj] # Current value of j
        p1 = LegendreP(1.0-gamma,1.0/(jloc)) # Current value of LegendreP[1-gamma,1/j]
        p2 = LegendreP(2.0-gamma,1.0/(jloc)) # Current value of LegendreP[2-gamma,1/j]
        tabdhMdjInt[indj] = _dhMdj(gamma,jloc,p1,p2) # Pre-computing the value of hM(j)
    end
    #####
    intdhMdj = Interpolations.scale(interpolate(tabdhMdjInt, BSpline(Cubic(Line(OnGrid())))),rangejInt) # Constructing the interpolation function for j->hM(j)
    #####
    return intdhMdj # Returning the interpolation function
end

"""
    getd2hMdj2Int(iBath)

Function that returns an interpolation function for the function j->d2hMdj2(j) for the bath population indexed by iBath in the source file.
d2hMdj2(j) is the jj-second-order gradient of the j-dependent part of the mass precession frequency.
"""
function getd2hMdj2Int(iBath::Int64)
    jminInt = jlc(rh) # Minimum value of j for the interpolation, fixed by the influence radius
    jmaxInt = 1.0 # Maximum value of j for the interpolation
    #####
    rangejInt = range(jminInt,length=nbjInt,jmaxInt) # Range in j over which the interpolation is performed
    tabjInt = collect(rangejInt) # Table of the values of j for which the functions must be pre-computed
    tabd2hMdj2Int = zeros(Float64,nbjInt) # Container for the values of dhMdj(j)
    ####
    gamma = tabgammaBath[iBath] # Cusp index of the current bath population
    #####
    for indj=1:nbjInt # Loop over the j of the grid
        jloc = tabjInt[indj] # Current value of j
        p1 = LegendreP(1.0-gamma,1.0/(jloc)) # Current value of LegendreP[1-gamma,1/j]
        p2 = LegendreP(2.0-gamma,1.0/(jloc)) # Current value of LegendreP[2-gamma,1/j]
        tabd2hMdj2Int[indj] = _d2hMdj2(gamma,jloc,p1,p2) # Pre-computing the value of hM(j)
    end
    #####
    intd2hMdj2 = Interpolations.scale(interpolate(tabd2hMdj2Int, BSpline(Cubic(Line(OnGrid())))),rangejInt) # Constructing the interpolation function for j->hM(j)
    #####
    return intd2hMdj2 # Returning the interpolation function
end

##################################################
# Preparing containers for the interpolation function of j->hM(j) and j->dhMdj(j)
# See https://discourse.julialang.org/t/array-of-functions-is-there-a-way-to-avoid-allocations-performance-penalty/24471

"""
    tabhMInt

Static container for the interpolation function of j->hM(j).
Set equal to `SVector{nbBath}([gethMInt(iBath) for iBath=1:nbBath])`.
See `https://discourse.julialang.org/t/array-of-functions-is-there-a-way-to-avoid-allocations-performance-penalty/24471`.
"""
const tabhMInt = SVector{nbBath}([gethMInt(iBath) for iBath=1:nbBath])

"""
    tabdhMIntdj

Static container for the interpolation function of j->dhMdj(j).
Set equal to `SVector{nbBath}([getdhMdjInt(iBath) for iBath=1:nbBath])`.
See `https://discourse.julialang.org/t/array-of-functions-is-there-a-way-to-avoid-allocations-performance-penalty/24471`.
"""
const tabdhMdjInt = SVector{nbBath}([getdhMdjInt(iBath) for iBath=1:nbBath])

"""
    tabd2hMIntdj2

Static container for the interpolation function of j->d2hMdj2(j).
Set equal to `SVector{nbBath}([getd2hMdj2Int(iBath) for iBath=1:nbBath])`.
See `https://discourse.julialang.org/t/array-of-functions-is-there-a-way-to-avoid-allocations-performance-penalty/24471`.
"""
const tabd2hMdj2Int = SVector{nbBath}([getd2hMdj2Int(iBath) for iBath=1:nbBath])

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

"""
    nuMBath(iBath,a,j)

Mass precession frequency for an orbit with semi-major axis a and reduced angular momentum j, for the family indexed by iBath in the source file.
Uses pre-computed interpolations functions stacked into arrays, as well as an @inline function for the evaluation to be fast.
"""
@inline function nuMBath(iBath::Int64,a::Float64,j::Float64)
    Menclosed = MenclosedBath(iBath,a) # Enclosed mass M(<a) of the current bath population
    #####
    return nuMbar(Menclosed,a)*tabhMInt[iBath](j) # Returning the mass precession frequency !! ATTENTION, we use the interpolation function for j->hM(j)
end

"""
    dnuMBathdj(iBath,a,j)

j-gradient of the mass precession frequency for an orbit with semi-major axis a and reduced angular momentum j, for the family indexed by iBath in the source file.
Uses pre-computed interpolations functions stacked into arrays, as well as an @inline function for the evaluation to be fast.
"""
@inline function dnuMBathdj(iBath::Int64,a::Float64,j::Float64)
    Menclosed = MenclosedBath(iBath,a) # Enclosed mass M(<a) of the current bath population
    #####
    return nuMbar(Menclosed,a)*tabdhMdjInt[iBath](j) # Returning the mass precession frequency !! ATTENTION, we use the interpolation function for j->dhMdj(j)
end

"""
    d2nuMBathdj2(iBath,a,j)

jj-second-ordre gradient of the mass precession frequency for an orbit with semi-major axis a and reduced angular momentum j, for the family indexed by iBath in the source file.
Uses pre-computed interpolations functions stacked into arrays, as well as an @inline function for the evaluation to be fast.
"""
@inline function d2nuMBathdj2(iBath::Int64,a::Float64,j::Float64)
    Menclosed = MenclosedBath(iBath,a) # Enclosed mass M(<a) of the current bath population
    #####
    return nuMbar(Menclosed,a)*tabd2hMdj2Int[iBath](j) # Returning the mass precession frequency !! ATTENTION, we use the interpolation function for j->dhMdj(j)
end

"""
    nuMtot(a,j)

Total mass precession frequency for an orbit with semi-major axis a and reduced angular momentum j.
Uses pre-computed interpolations functions stacked into arrays, as well as an @inline function for the evaluation to be fast.
"""
@inline function nuMtot(a::Float64,j::Float64)
    res = 0.0 # Initialising the result
    #####
    for iBath=1:nbBath # Loop over the bath particles
        res += nuMBath(iBath,a,j) # Contribution from the ith bath population
    end
    #####
    return res # Output
end

"""
    dnuMtotdj(a,j)

j-gradient of the total mass precession frequency for an orbit with semi-major axis a and reduced angular momentum j.
Uses pre-computed interpolations functions stacked into arrays, as well as an @inline function for the evaluation to be fast.
"""
@inline function dnuMtotdj(a::Float64,j::Float64)
    res = 0.0 # Initialising the result
    #####
    for iBath=1:nbBath # Loop over the bath particles
        res += dnuMBathdj(iBath,a,j) # Contribution from the ith bath population
    end
    #####
    return res # Output
end

"""
    d2nuMtotdj2(a,j)

jj-second-order gradient of the total mass precession frequency for an orbit with semi-major axis a and reduced angular momentum j.
Uses pre-computed interpolations functions stacked into arrays, as well as an @inline function for the evaluation to be fast.
"""
@inline function d2nuMtotdj2(a::Float64,j::Float64)
    res = 0.0 # Initialising the result
    #####
    for iBath=1:nbBath # Loop over the bath particles
        res += d2nuMBathdj2(iBath,a,j) # Contribution from the ith bath population
    end
    #####
    return res # Output
end

"""
    nup(a,j)

Total precession frequency for an orbit with semi-major axis a and reduced angular momentum j.
Uses pre-computed interpolations functions stacked into arrays, as well as an @inline function for the evaluation to be fast.
"""
@inline function nup(a::Float64,j::Float64)
    return nuGR(a,j) + nuMtot(a,j) # Returning the value of the total precession frequency
end

"""
    dnupdj(a,j)

j-gradient of the total precession frequency for an orbit with semi-major axis a and reduced angular momentum j.
Uses pre-computed interpolations functions stacked into arrays, as well as an @inline function for the evaluation to be fast.
"""
@inline function dnupdj(a::Float64,j::Float64)
    return dnuGRdj(a,j) + dnuMtotdj(a,j) # Returning the gradient of the total precession frequency
end

"""
    d2nupdj2(a,j)

jj-second-order gradient of the total precession frequency for an orbit with semi-major axis a and reduced angular momentum j.
Uses pre-computed interpolations functions stacked into arrays, as well as an @inline function for the evaluation to be fast.
"""
@inline function d2nupdj2(a::Float64,j::Float64)
    return d2nuGRdj2(a,j) + d2nuMtotdj2(a,j) # Returning the second gradient of the total precession frequency
end


text:
#####
##### OneDimensionalTimeSeries
#####

"""
    struct OneDimensionalTimeSeries{FT, F, G, C, D, UU, VV, TT, SS}

A time series of horizontally-averaged observational or LES data
gridded as Oceananigans fields.
"""
struct OneDimensionalTimeSeries{P, F, G, T, I, R, N} <: AbstractTimeSeries
    file_path::P
    fields::F
    grid::G
    times::T
    time_range::R
    name::N
end

const OneDimensionalTimeSeriesBatch = Vector{<:OneDimensionalTimeSeries}

function OneDimensionalTimeSeries(file_path, field_names=(:u, :v, :b, :e), time_range=nothing, Nz=64)

    file = jldopen(file_path, "r")

    αg = file["buoyancy/equation_of_state/α"] * file["buoyancy/gravitational_acceleration"]

    name = get_parameter(datapath, "parameters", "name", "")

    iters = get_iterations(datapath)

    simulation_grid = grid_type(datapath)
    u = [ new_field(XFaceField, simulation_grid, get_data("u", datapath, iter)) for iter in iters ]
    v = [ new_field(YFaceField, simulation_grid, get_data("v", datapath, iter)) for iter in iters ]
    b = [ new_field(CenterField, simulation_grid, get_data("T", datapath, iter) .* αg) for iter in iters ]
    e = [ new_field(CenterField, simulation_grid, get_data("e", datapath, iter)) for iter in iters ]

    times = get_times(datapath)

    model_grid = OneDimensionalEnsembleGrid(datapath; size=(1,1,Nz))

    # deserialize as much as possible into metadata
    meta_data = Dict(group => file[group] for group in filter(n -> n != "timeseries", keys(file)))

    fields = NamedTuple(field_name => eval(field_name) for field_name in field_names)

    time_range = isnothing(time_range) ? range(times[1], stop=times[end]) : time_range

    observation = OneDimensionalTimeSeries(file_path, fields, simulation_grid, times, time_range, name)

    return OneDimensionalTimeSeries(observation, model_grid)
end

"""
    OneDimensionalTimeSeries(observation::OneDimensionalTimeSeries, grid)

Returns `observation::OneDimensionalTimeSeries` interpolated to `grid`.
"""
function OneDimensionalTimeSeries(data::OneDimensionalTimeSeries, grid::AbstractGrid)

    U = [ XFaceField(grid) for t in data.times ]
    V = [ YFaceField(grid) for t in data.times ]
    B = [ CenterField(grid) for t in data.times ]
    E = [ CenterField(grid) for t in data.times ]

    for i = 1:length(data.t)
        set!(U[i], data.fields.u[i])
        set!(V[i], data.fields.v[i])
        set!(B[i], data.fields.b[i])
        set!(E[i], data.fields.e[i])
    end

    return OneDimensionalTimeSeries(data.file_path,
                      data.fields,
                      grid,
                      data.times,
                      data.time_range,
                      data.name)
end


text:
# Prefer to parse as integer, parse as float if not
function parseNumber(n::String)
    try
        return parse(Int, n)

    catch ArgumentError
        return parse(Float64, n)
    end
end

# Base.isnumber doesn't cover floats
function isNumber(n::String)
    try
        parseNumber(n)

        return true

    catch ArgumentError
        return false
    end
end

function hasExt(fn::String, ext::String, ignorecase::Bool=Sys.iswindows())
    path, fnext = splitext(fn)

    if ignorecase
        return lowercase(fnext) == lowercase(ext)

    else
        return fnext == ext
    end
end

function humanizeVariableName(s::String)
    text = replace(s, "_" => " ")
    text = replace(text,  "/" => " ")

    prevUpper = false
    nextLower = false
    len = length(s)

    res = Char[]

    for (i, c) in enumerate(text)
        thisUpper = isuppercase(c)
        nextLower = i < len && islowercase(text[i + 1])

        if thisUpper && !prevUpper
            push!(res, ' ')
        end

        if thisUpper && nextLower && res[end] != ' '
            push!(res, ' ')
        end

        prevUpper = thisUpper

        push!(res, c)
    end

    return titlecase(strip(join(res)))
end

function useIfApplicable(f, args...)
    if applicable(f, args...)
        f(args...)

        return true
    end

    return false
end

lerp(n::Number, m::Number, a::Number) = n + (n - m) * a