text:
using EffectiveWaves, Test, LinearAlgebra

spatial_dim = 3
medium = Acoustic(spatial_dim; ρ=1.2, c=1.5)

# ms = MultipleScattering # just in case Circle and Sphere conflict with definitions from other packages.

s1 = Specie(
    Acoustic(spatial_dim; ρ=10.2, c=10.1), Sphere(0.004);
    # Acoustic(spatial_dim; ρ=10.2, c=10.1), Sphere(0.4);
    volume_fraction=0.2
    # volume_fraction=0.01
);
s2 = Specie(
    Acoustic(spatial_dim; ρ=5.2, c=4.1), Sphere(0.003);
    # Acoustic(spatial_dim; ρ=5.2, c=4.1), Sphere(0.3);
    volume_fraction=0.1
);
species = [s1,s2]
species = [s1]

ω = 1e-1
# ω = 0.2
tol = 1e-7
basis_order = 1
# basis_order = 2

θ = 0.0
psource = PlaneSource(medium, [sin(θ),0.0,cos(θ)]);
source = plane_source(medium; direction = [sin(θ),0.0,cos(θ)])
material = Material(Sphere(4.0),species);

# x = rand(3)
# field(source,x,ω) - field(psource,x,ω)

R = outer_radius(material.shape)
# basis_field_order = estimate_regular_basisorder(typeof(source.medium), R * k_eff )
# basis_field_order = 3
basis_field_order = 3

k = ω / medium.c
ko = k * medium.c / s1.particle.medium.c

k_phi = wavenumber_low_volumefraction(ω, medium, species;
    basis_order = basis_order)

# AP_kps0 =  wavenumbers_bisection_robust(ω, medium, species;
#     bisection_mesh_points = 30,
#     basis_order = basis_order,
#     num_wavenumbers = 3, tol = tol,
#     symmetry = PlanarAzimuthalSymmetry())

opts = Dict(
    :tol => tol, :num_wavenumbers => 6,
    :mesh_size => 2.0, :mesh_points => 20,
    :basis_order => basis_order, :basis_field_order => basis_field_order
);
AP_kps = wavenumbers(ω, medium, species; symmetry = PlanarAzimuthalSymmetry(), opts...)

# wavemodes = WaveModes(ω, source, material; opts...)

k_eff = AP_kps[1]
# k_eff = A_kps[1]
(imag(k_eff) < -tol) && (k_eff = - k_eff)

eff_medium = effective_medium(medium, species)
k_low = ω/eff_medium.c;
eff_medium.c
ω / k_eff
ω / k_phi

Symmetry(source, material)
Symmetry(psource, material)

k_eff = AP_kps[2]
k_eff = AP_kps[1]

A_wave = WaveMode(ω, k_eff, psource, material;
    basis_order = basis_order, basis_field_order = basis_field_order);

scat_azi = material_scattering_coefficients(A_wave);

R_wave = WaveMode(ω, k_eff, source, material;
    basis_order = basis_order, basis_field_order = basis_field_order);

scat_reg = material_scattering_coefficients(R_wave);

norm(scat_azi - scat_reg) / norm(scat_azi)

Linc = basis_field_order + basis_order

effective_sphere = Particle(eff_medium, material.shape)

medium3 = Acoustic(spatial_dim; ρ=0.2, c=0.1)
medium3 = eff_medium
effective_sphere = Particle(medium3, material.shape)

Tmat = MultipleScattering.t_matrix(effective_sphere, medium, ω, Linc)
t_diag = [Tmat[i,i] for i in axes(Tmat,1)]

# medium2 = Acoustic(2; ρ=0.2, c=0.1)
# effective_sphere = Particle(medium2, Circle(outer_radius(material.shape)))
# outer_medium2 = Acoustic(2; ρ=1.0, c=1.0)
# Tmat = MultipleScattering.t_matrix(effective_sphere, outer_medium2, ω, Linc)

p = effective_sphere
outer_medium = medium
basis_order = Linc

n_to_l = [l for l = 0:Linc for m = -l:l];

inds_azi = findall([
    m == 0
for l = 0:Linc for m = -l:l]);

source_coefficients = regular_spherical_coefficients(source)(Linc,zeros(3),ω);

t_diag = [Tmat[i,i] for i in axes(Tmat,1)];

norm(scat_azi - t_diag[n_to_l .+ 1] .* source_coefficients) / norm(scat_azi)

MA = eigensystem(ω, psource, material;
    basis_order = basis_order,
    basis_field_order = basis_field_order
)

MM = eigensystem(ω, source, material;
    basis_order = basis_order,
    basis_field_order = basis_field_order
)

det(MM(k_eff))
det(MM(k_low))
det(MM(k_phi))

det(MA(k_eff))
det(MA(k_low))
det(MA(k_phi))

#NOTE: MM(k_eff) ≈ MM_svd.U * diagm(0 => MM_svd.S) * MM_svd.Vt
MM_svd = svd(MM(k_eff))
MM_svd.S
inds = findall(MM_svd.S .< 1e-4)
# 9 -> 48
eigvectors = MM_svd.V[:,inds]

MM_svd = svd(MA(k_eff))
MM_svd.S
inds = findall(MM_svd.S .< 1e-4)
# 3 -> 10
eigvectors = MM_svd.V[:,inds]

eigvectors = eigenvectors(ω, k_eff, source, material;
        basis_order = basis_order,
        basis_field_order = basis_field_order
)

eigvectors = eigenvectors(ω, k_eff, psource, material;
        basis_order = basis_order,
        basis_field_order = basis_field_order
)

sum(eigvectors[i] * α[i[2]] for i in CartesianIndices(eigvectors), dims = 2)

T = Float64

norm(MM(k_eff) * eigvectors[:,1] - MM_svd.S[inds[1]] .* eigvectors[:,1]) / norm(eigvectors[:,1])

norm(MM(k_eff) * eigvectors[:,end] - MM_svd.S[inds[end]] .* eigvectors[:,end]) / norm(eigvectors[:,end])

L = basis_order
L1 = basis_field_order
L2 = L1
# or is it this:
# L2 = L + L1

# l,dl <= L
# l1 <= L1
# l2 <= L1 + L
# l3 <= 2L # assuming L<= L1

# k1 = 1e-8
# k2 = 1e-8 + 1e-6im
# kernelN3D.(0:5, k1, k2) ./ (k2^2 - k1^2)
#
#
# outer_medium = medium
# p = s1.particle
# m = 0
#
# ak = outer_radius(p)*ω/outer_medium.c
#
# q = impedance(p.medium)/impedance(outer_medium) # Impedance ratio
# γ = outer_medium.c / p.medium.c #speed ratio
# numer = q * diffsbesselj(m, ak) * sbesselj(m, γ * ak) - sbesselj(m, ak)*diffsbesselj(m, γ * ak)
# denom = q * diffshankelh1(m, ak) * sbesselj(m, γ * ak) - shankelh1(m, ak)*diffsbesselj(m, γ * ak)
#
# - numer / denom
# t_matrix(s1.particle, medium, ω, 0)
#
#
# -(SphericalBesselJ[0, z]/(2 z)) +
#  1/2 (SphericalBesselJ[-1, z] - SphericalBesselJ[1, z])
#
# m=0
# (- sbesselj(m,ak) + (ak) * (sbesselj(m-1,ak) - sbesselj(m+1,ak))) / (2 * (ak))

max_stars_repo_name: arturgower/EffectiveWaves


text:
import LowRankModels: objective
import FirstOrderOptimization: LowRankOperator

export objective

## Objective evaluation
# we're not going to bother checking whether W is psd or not
# when evaluating the objective; in the course of the prisma
# algo this makes no difference
function objective(gfrm::GFRM, W::Array{Float64,2};
                   yidxs=get_yidxs(gfrm.losses),
                   include_regularization=true)
    # W is the symmetric parameter; U is the upper right block
    m,n = size(gfrm.A)
    UW = W[1:m, m+1:end]
    err = 0.0
    for j=1:n
        for i in gfrm.observed_examples[j]
            err += evaluate(gfrm.losses[j], UW[i,yidxs[j]], gfrm.A[i,j])
        end
    end
    if include_regularization
        err += evaluate(gfrm.r, W)
    end
    return err
end
function objective(gfrm::GFRM; kwargs...)
    objective(gfrm::GFRM, gfrm.W; kwargs...)
end

function objective(gfrm::GFRM, W::LowRankOperator;
                   yidxs=get_yidxs(gfrm.losses),
                   include_regularization=true)
    m,n = size(gfrm.A)
    err = 0.0
    for j=1:n
        for i in gfrm.observed_examples[j]
            err += evaluate(gfrm.losses[j], W[i,yidxs[j]], gfrm.A[i,j])
        end
    end
    if include_regularization
        err += evaluate(gfrm.r, W)
    end
    return err
end

# objective(gfrm::GFRM, args...; kwargs...) = objective(gfrm, gfrm.W, args...; kwargs...)

ext: jl


text:
# # Working with Strings
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__notebooks/06_Strings.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__notebooks/06_Strings.ipynb)
#
# ## Regex
#
# You can also use [regular expressions](https://docs.julialang.org/en/v1/manual/strings/#Regular-Expressions-1).
# They are very useful to parse text files.

line = ">sp|O43521|B2L11_HUMAN Bcl-2-like protein 11 OS=Homo sapiens OX=9606 GN=BCL2L11 PE=1 SV=1"

#-

regex = r"^>\w+\|(\w+)\|"

#-

m = match(regex, line)

#-

if m !== nothing
    println(m[1])
end

# ## String interpolation
#
# You can [interpolate values](https://docs.julialang.org/en/v1/manual/strings/#string-interpolation-1),
# single variables or the result of more complex expressions, into
# strings using `$`:

a = 1
b = 2
"$a + $b is $(a + b)"

max_stars_repo_name: diegozea/JuliaForBioinformatics

max_issues_repo_name: diegozea/JuliaForBioinformatics

max_forks_repo_name: diegozea/JuliaForBioinformatics


text:
@testset "Fluxes" begin
    At = CrossSectionalArea(ts,rand(2))
    cvt = AlongChannelVelocity(ts,randn(2))
    Q = At*cvt
    T = Temperature(ts,rand(2))
    F1 = T*Q
    F2 = Q*T
    @test typeof(F1) == Flux{Float64,Temperature{Float64}}
    @test quantity(F1) == quantity(F2)
end


text:
using Documenter
using InplaceArrays

pages = [
  "Home" => "index.md",
  "Gridap" => "Gridap.md",
  "Gridap.Helpers" => "Helpers.md",
  "Gridap.Inference" => "Inference.md",
  "Gridap.TensorValues" => "TensorValues.md",
  "Gridap.Arrays" => "Arrays.md",
  "Gridap.Fields" => "Fields.md",
  "Gridap.Polynomials" => "Polynomials.md",
  "Gridap.ReferenceFEs" => "ReferenceFEs.md",
 ]

makedocs(
    sitename = "InplaceArrays",
    format = Documenter.HTML(),
    modules = [InplaceArrays],
    pages = pages
)

deploydocs(
    repo = "github.com/fverdugo/InplaceArrays.jl.git",
)


text:
using OrdinaryDiffEq, DifferentialEquations

tspan = (0.0f0, 50.0f0)

"""
Van der Pol oscillator model from wikipedia
"""
function vdp_system(du, u, p, t;
                     U_func_time=U_func_time_trivial,
                     U_func_space=U_func_space_trivial)
    Ft = U_func_time(t)
    Fs = U_func_space(u)
    μ = p[1]
    x, y = u

    du[1] = μ*(x - (x^3)/3 - y)
    du[2] = x/μ

    # Add controllers in
    du .+= Ft .+ Fs
end

# Generate data, with a time component
p = [1.0]
u0 = [1.0, 1.0]
ts = range(tspan[1], tspan[2], length=5001)

function solve_vdp_system(;U_func_time=U_func_time_trivial,
                              U_func_space=U_func_space_trivial)
    prob = ODEProblem((du, u, p, t)->vdp_system(du, u, p, t,
                                U_func_time=U_func_time,
                                U_func_space=U_func_space),
                    u0, tspan, p)
    sol = solve(prob, Tsit5(), saveat=ts);
    return sol
end

#####
##### Turing function for Bayesian parameter estimation
#####
# Unforced version for use with Turing.jl
# function vdp_system(u, p, t)
#     u = convert.(eltype(p),u)
#     beta, gamma = p
#     S, I, R = u
#
#     N = sum([S,I,R])
#     tmp = beta*I*S/N
#     du = [-tmp; tmp - gamma*I; gamma*I]
#     return du
# end

#####
##### True model in SINDy syntax
#####
μ = p
#      x  y  c  xx xy yy xxx  xxy xyy yyy
 A = [[μ -μ  0  0  0  0  -μ/3 0   0   0];
      [1/μ 0 0  0  0  0  0    0   0   0]]
n = size(A, 1)
sindy_library = Dict(
    "cross_terms"=>[2,3],
    "constant"=>nothing
);
core_dyn_true = sindyModel(ts, A,
                            convert_string2function(sindy_library),
                            ["x", "y"])

# Actually export
export ts, solve_vdp_system, vdp_system, core_dyn_true


text:
using Base.Test
using LogClustering.Clustering

@testset "test soft & hard clustering on small sample data" begin
    data = [
        1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0;
        1.1 2.1 3.1 4.1 5.1 6.1 7.1 8.1 9.1;
        1.2 2.2 3.2 4.2 5.2 6.2 7.2 8.2 9.2;
    ]
    L, N = size(data)

    Cs, used = knn_clustering(data, k = 2, mode=:soft)
    @show length(Cs), Cs
    @test length(Cs) == N
    @show used
    # @test length(used) in 4:N
    @test length(used) == N

    Cs, used = knn_clustering(data, k = 2)
    @show length(Cs), Cs
    @test length(Cs) in 2:N
    @show used
    # @test length(used) in 4:N
    @test length(used) == N

    Vs = map(c->map(i->data[:,i],c),Cs)
    @show Vs
end

@testset "test soft & hard clustering on small diff sample data" begin
    data = [
        1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0;
        1.0 1.0 1.0 1.0 1.0000003 1.0000004 1.0000005 1.0000006 1.0000007 1.0000008 1.0000009;
    ]
    L, N = size(data)

    Cs, used = knn_clustering(data, k = 4, mode=:soft)
    @show length(Cs), Cs
    @test length(Cs) == N
    @show used
    # @test length(used) in 4:N
    @test length(used) == N

    Cs, used = knn_clustering(data, k = 4)
    @show length(Cs), Cs
    @test length(Cs) in 2:N
    @show used
    # @test length(used) in 4:N
    @test length(used) == N

    Vs = map(c->map(i->data[:,i],c),Cs)
    @show Vs
end

@testset "Validate Clustering with BetaCV" begin
    using LogClustering.Validation

    L = 3
    N = 1000
    data = randn(L,N)
    min = minimum(data)
    max = maximum(data)
    @show min, max

    C, used = knn_clustering(data)

    @show length(C)
    @show length(used) / N
    @show mean(map(length, values(C)))
    @show median(map(length, values(C)))

    Vs = map(c->data[:,c],C)
    bcv = @time betacv(Vs)
    @show bcv
    @test bcv > 0.7
    @test bcv < 0.8
end

max_stars_repo_name: swiesend/FluxMO

max_issues_repo_name: swiesend/FluxMO

max_forks_repo_name: swiesend/FluxMO


text:
include("canvas.jl")
include("crashes.jl")
include("compare.jl")


text:
module Geom

using Color
using Compose
using DataArrays
using Gadfly

import Compose.combine # Prevent DataFrame.combine from taking over.
import Gadfly: render, element_aesthetics, inherit, escape_id,
               default_statistic, setfield!, set
import Iterators
import Iterators: cycle, product, distinct, take, chain

# Geometry that renders nothing.
immutable Nil <: Gadfly.GeometryElement
end

const nil = Nil

function render(geom::Nil, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
end

# Catchall
function default_statistic(::Gadfly.GeometryElement)
    Gadfly.Stat.identity()
end

include("geom/bar.jl")
include("geom/boxplot.jl")
include("geom/errorbar.jl")
include("geom/hexbin.jl")
include("geom/hline.jl")
include("geom/label.jl")
include("geom/line.jl")
include("geom/point.jl")
include("geom/rectbin.jl")
include("geom/subplot.jl")
include("geom/vline.jl")
include("geom/ribbon.jl")

end # module Geom


text:
#=
my_tree:
- Julia version: 1.3.1
- Author: erathorn
- Date: 2019-05-07
=#

#TODO: Automate export of automatically genereated funtions

"""
    add_child!(mother_node::Node, child::Node, child_position::Union{Int64, Missing}=missing)

This function adds a child to the mother node.
The arity of the mother node is increased by `1` and the root
status of the child is set to `False`.

* `mother_node` : Node to add a child to.

* `child` : Node to add to mother_node.children.

* `child_position` : index at which to add the new child node; optional.
"""
function add_child!(mother_node::Node, child::Node, child_position::Union{Int64, Missing}=missing)
    if ismissing(child_position)
        push!(mother_node.children, child)
    else
        insert!(mother_node.children, child_position, child)
    end # if/else
    child.mother = mother_node
    mother_node.nchild += 1
    child.root = false
    mother_node.initialized = true
end # function add_child

"""
    remove_child!(mother_node::Node, left::Bool)::Node

This function removes a child from the list of nodes which are daughters of this
node. An input of "True" removes the left child,
while "False" removes the right child.

Returns the removed node.

* `mother_node` : Node from which to remove a child.

* `left` : boolean value determining which child of mother_node to remove.
"""
function remove_child!(mother_node::N, left::Bool)::N where N<:GeneralNode
    if left
        rv = popfirst!(mother_node.children)
        rv.mother = missing
    else
        rv = pop!(mother_node.children)
        rv.mother= missing
    end # end if
    mother_node.nchild -= 1
    return rv
end # function

"""
    remove_child!(mother_node::Node, child::Node)::Node

This function removes a child from the list of nodes which are daughters of this
node.

The removed node is returned.

* `mother_node` : Node from which to remove a child.

* `child` : specific Node to remove.
"""
function remove_child!(mother_node::N, child::N)::N where N<:GeneralNode
    ind = findfirst(x->x==child, mother_node.children)
    deleteat!(mother_node.children, ind)
    child.mother = missing

    mother_node.nchild -= 1
    return child
end # function

"""
    delete_node!(node::Node)::Nothing

This functions deletes node from a tree and assigns all its children to its
mother node.

* `node` : Node to be deleted.
"""
function delete_node!(node::T)::Nothing where T<:AbstractNode
    if node.root == true
        throw(ArgumentError("Cannot remove root node"))
    end
    mother = node.mother
    for child in node.children
        add_child!(mother, child, findfirst(x -> x == node, mother.children))
    end
    remove_child!(mother, node)
    return nothing
end

"""
    insert_node!(mother::Node, children::Vector{T})::T where T<:AbstractNode

This function inserts a node into a tree after a mother node and gains
a subset of the mother's children as its children.

Returns the inserted node.

* `mother` : Node under which to add the newly-inserted node.

* `children` : Children of node referenced by "mother" to reassign as children of the newly-inserted node.
"""
function insert_node!(mother::T, children::Vector{T})::T where T<:AbstractNode
    @assert length(children) >= 1
    inserted_node = Node()
    for child in children
        @assert child in mother.children
    end # for
    index = findfirst(x -> x in children, mother.children)
    for child in children
        remove_child!(mother, child)
        add_child!(inserted_node, child)
    end # for
    add_child!(mother, inserted_node, index)
    return inserted_node
end

function tree_from_leaves(leaf_nodes::Vector{String},node_size::Int, final_length::Int64)::Tuple{Vector{Node}, Int}
    my_node_list::Array{Node,1} = []

    # first create a list of leaf nodes
    for node_name in leaf_nodes
        nn = Node(node_name, data=zeros(Float64, (2, node_size)))

        push!(my_node_list,nn)
    end # for

    # Internal nodes are created using integers as names.
    temp_name::Int = length(my_node_list)+1

    # shuffle the node list to get a random tree
    Random.shuffle!(my_node_list)

    while length(my_node_list) > final_length
        # get two nodes
        # create a new mother node to which the two first nodes are added as children
        # add the new mother node to the list and reshuffle
        first_child::Node = pop!(my_node_list)
        first_child.inc_length = rand()#*0.1
        second_child::Node = pop!(my_node_list)
        second_child.inc_length = rand()
        curr_node::Node = Node(string(temp_name), data=zeros(Float64, (2, node_size)))

        add_child!(curr_node, first_child)
        add_child!(curr_node, second_child)
        push!(my_node_list, curr_node)
        temp_name += 1
        Random.shuffle!(my_node_list)
    end # while

    return my_node_list, temp_name
end

"""
    create_tree_from_leaves(leaf_nodes::Vector{T})::Node

This function creates a  random binary tree from a list of leaf nodes.

The root node as access point for the tree is returned.

* `leaf_nodes` : vector of leaf nodes.
"""
function create_tree_from_leaves_bin(leaf_nodes::Vector{String}, node_size::Int)::Node

    my_node_list, temp_name = tree_from_leaves(leaf_nodes, node_size ,2)

    root::Node = Node(string(temp_name), data=zeros(Float64, (2, node_size)))

    lchild = pop!(my_node_list)
    lchild.inc_length = rand()

    rchild = pop!(my_node_list)
    rchild.inc_length = rand()
    add_child!(root, lchild)
    add_child!(root, rchild)

    set_binary!(root)
    number_nodes!(root)

    return root
end # function create_tree_from_leaves

"""
    create_tree_from_leaves(leaf_nodes::Vector{T})::Node

This function creates a  random binary tree from a list of leaf nodes.

The root node as access point for the tree is returned.

* `leaf_nodes` : Names of nodes to create.

* `node_size` : used to initialize node.data

* `cu` : ??
"""
function create_tree_from_leaves(leaf_nodes::Vector{String}, node_size::Int64 = 1; cu::Bool=false)#::N where N <: GeneralNode

    my_node_list, temp_name = tree_from_leaves(leaf_nodes, node_size, 3)

    root::Node = Node(string(temp_name), data=zeros(Float64, (2, node_size)))
    lchild = pop!(my_node_list)
    lchild.inc_length = rand()
    mchild = pop!(my_node_list)
    mchild.inc_length = rand()
    rchild = pop!(my_node_list)
    rchild.inc_length = rand()
    add_child!(root, lchild)
    add_child!(root, rchild)
    add_child!(root, mchild)

    set_binary!(root)
    number_nodes!(root)

    return root
end # function create_tree_from_leaves

function create_tree_from_leaves_cu(leaf_nodes::Vector{String}, node_size::Int64 = 1)::Node_cu
    my_node_list::Array{Node_cu,1} = []

    # first create a list of leaf nodes
    for node_name in leaf_nodes
        nn =  Node_cu(node_name, zeros(Float64, (2, node_size)),missing, missing, missing, 0, true, 0.0, "0", 0, 0.0)
        push!(my_node_list,nn)
    end # for

    # Internal nodes are created using integers as names.
    temp_name::Int = length(my_node_list)+1

    # shuffle the node list to get a random tree
    Random.shuffle!(my_node_list)

    while length(my_node_list) != 1
        # get two nodes
        # create a new mother node to which the two first nodes are added as children
        # add the new mother node to the list and reshuffle
        first_child::Node = pop!(my_node_list)
        first_child.inc_length = rand(Uniform(0,1))
        second_child::Node = pop!(my_node_list)
        second_child.inc_length = rand(Uniform(0,1))
        curr_node::Node = Node_cu(string(temp_name), zeros(Float64, (2, node_size)), missing, missing, missing, 0, true, 0.0, "0", 0,0.0)
        add_child!(curr_node, first_child, true)
        add_child!(curr_node, second_child, false)
        push!(my_node_list, curr_node)
        temp_name += 1
        Random.shuffle!(my_node_list)
    end # while
    root = pop!(my_node_list)
    set_binary!(root)
    number_nodes!(root)

    return root
end # function create_tree_from_leaves

# legacy wrapper
function rescale_length(root::T) where T<:GeneralNode
    force_ultrametric(root)
end

"""
    force_ultrametric!(root::T) where T<:GeneralNode

Force an ultrametric version of the tree.
"""
function force_ultrametric!(root::T) where T<:GeneralNode
    po::Vector{T} = post_order(root)
    node2max_depth = zeros(UInt32, length(po))
    for node in po
        if node.nchild != 0
            mv = -1
            for child in node.children
                if node2max_depth[child.num] > mv
                    mv = node2max_depth[child.num]
                end # if
            end # for
            node2max_depth[node.num] = mv+1
        else
            node2max_depth[node.num] = 1
        end # if
    end # for

    node2dist = zeros(Float64, size(node2max_depth))
    tl = tree_height(root)
    nblv = zeros(Float64,length(po))
    for node in level_order(root)
        if node.root != true
            nv = (tl - node2dist[node.mother.num])/node2max_depth[node.num]

            nblv[node.num] = nv
            node2dist[node.num] = nv + node2dist[node.mother.num]
        end # if
    end # for

    set_branchlength_vector!(root, nblv)
end # function force_ultrametric!

#################### Tree length & height ####################

"""
    tree_length(root::T)::Float64  where T<:GeneralNode

This function calculates the tree length.
"""
function tree_length(root::T)::Float64  where T<:GeneralNode
    return tree_length(root, 0.0)
end # function tree_length

"""
    tree_length(root::T, tl::Float64)::Float64 where T<:GeneralNode

This function does the internal tree length recursion
"""
function tree_length(root::T, tl::Float64)::Float64 where T<:GeneralNode

    #if length(root.children) != 0
    for child in root.children
        tl = tree_length(child, tl)
    end
    #end # if
    if root.root !== true
        tl += root.inc_length
    end

    tl
end # function tree_length

"""
    tree_height(root::T)::Float64  where T<:GeneralNode

This function calculates the tree height.
"""
function tree_height(root::T)::Float64  where T<:GeneralNode
    node_height(root)
    return root.height
end

"""
    node_height(root::T, mv::Float64)::Float64  where T<:GeneralNode

Calculate the height of a node.
"""
function node_height(root::T)  where T<:GeneralNode

    if root.nchild != 0
        for node in root.children
            node_height(node)
        end
        root.height = maximum([child.inc_length+child.height for child in root.children])
    else
        root.height = 0.0
    end
end # function node_height

function node_height_vec(root::T, vec::Vector{N})  where {T<:GeneralNode, N<:Real}

    if root.nchild != 0
        for node in root.children
            node_height_vec(node, vec)
        end
        root.height = maximum([child.inc_length+child.height for child in root.children])
    else
        root.height = 0.0
    end
    vec[root.num] = root.height
end # function node_height

function node_height_vec(root::T)::Vector{Float64} where T<:GeneralNode
    t = zeros(length(post_order(root)))
    node_height_vec(root, t)
    t
end # function node_height

"""
    node_depth(node::T)::Int64 where T<:GeneralNode

Calculate the depth of a node.
"""
function node_depth(node::T)::Int64 where T<:GeneralNode
    return length(split(node.binary, ",")) - 1
end

function node_distance(tree::T, node1::T, node2::T)::Float64 where T<:GeneralNode
    lca = find_lca(tree, node1, node2)
    path_length(lca, node1)+path_length(lca,node2)
end

function get_path(ancestor::T, descendant::T)::Vector{Int64} where T<:GeneralNode
    path::Vector{Int64} = []
    while descendant.num != ancestor.num
        push!(path, descendant.num)
        descendant = descendant.mother
    end
    path
end

"""
    path_length(ancestor::T, descendant::T)::Float64  where T<:GeneralNode

Note: The function assumes there is an ancestral relationship between the two nodes.

This function calculates the length of the path separating the ancestor from the
offspring node. The function follows the path specified through the binary
description of the node.
"""
function path_length(ancestor::T, descendant::T)::Float64  where T<:GeneralNode
    l::Float64 = 0.0

    while descendant != ancestor
        l += descendant.inc_length
        descendant = descendant.mother
    end # while
    return l
end #function path_length

"""
    get_sister(node::T)::T  where T<:GeneralNode

This function gets the sister of `node`. It does so by looking for the respective
binary representation of the sister.
"""
@inline function get_sister(node::T)::T  where T<:GeneralNode
    node.mother.children[findfirst(y-> y!=node, node.mother.children)]
end # function

"""
    function get_mother(node::T)::T  where T<:GeneralNode

This function gets the mother of `node`. It does so by looking for the respective
binary representation of the mother node.
"""
@inline function get_mother(node::T)::T  where T<:GeneralNode
    return node.mother
end # function

"""
    set_binary!(root::Node)

Assign a binary representation to each node, which specifies the path from the
root to this node via the binary representation of the node.
A left turn is a 1 in binary and a right turn a 0.
"""
function set_binary!(root::T)  where T<:GeneralNode
    if root.root
        root.binary = "1"
    end # if
    if root.nchild != 0
        for (ind, node) in enumerate(root.children)
            ind -= 1
            node.binary = string(root.binary,",", ind)
            set_binary!(node)
        end
    end # if
end # function set_binary

"""
    number_nodes!(root::T)::Nothing  where T<:GeneralNode

This function assigns a unique, sequential number to each node. Leaves are numbered first
in alphabetical order.
"""
function number_nodes!(root::T)::Nothing  where T<:GeneralNode
    tips = [n.name for n in get_leaves(root)]
    sort!(tips)
    running = length(tips)
    for node in post_order(root)
        if node.nchild == 0
            node.num = findfirst(x-> x==node.name, tips)
        else
            node.num = running + 1
            running += 1
        end
    end # for
end # fuction number_nodes

"""
    function random_node(root::T)::T  where T<:GeneralNode

This function returns a random node from the tree.
"""
function random_node(root::T)::T  where T<:GeneralNode
    post_order_trav = post_order(root)
    return rand(post_order_trav)
end # function random_node

#################### Vector of branch lengths: get & set ####################

"""
    get_branchlength_vector(root::N)::Vector{T}  where {N <:GeneralNode, T<:Real}

Get the vector of branch lengths of the tree.
"""
function get_branchlength_vector(root::N)::Vector{Float64}  where {N <:GeneralNode}
    if length(root.blv) == 0
        root.blv = zeros(length(post_order(root))-1)
    end
    get_branchlength_vector(root, root.blv)
    return root.blv
end # function get_branchlength_vector

"""
    get_branchlength_vector(t::TreeStochastic)

Get the vector of branch lengths of the tree.
"""
function get_branchlength_vector(t::TreeStochastic)::Vector{Float64}
    get_branchlength_vector(t.value)
end # function

"""
    get_branchlength_vector(root::N, out_vec::Vector{T}) where {N<:GeneralNode, T<:Real}

Do post order traversal to retrieve a vector of branch lengths.
"""
function get_branchlength_vector(root::N, out_vec::Vector{T})::Nothing where {N<:GeneralNode, T<:Real}
    for child in root.children
        get_branchlength_vector(child, out_vec)
    end
    if !root.root
        out_vec[root.num] = root.inc_length
    end
    nothing
end

"""
    set_branchlength_vector!(t::TreeStochastic, blenvec::Array{T}) where T <: Real

Get the vector of branch lengths of the tree.
"""
function set_branchlength_vector!(t::TreeStochastic, blenvec::Array{T}) where T <: Real
    set_branchlength_vector!(t.value, blenvec)
end # function

"""
    set_branchlength_vector!(t::TreeStochastic, blenvec::ArrayStochastic)

Get the vector of branch lengths of the tree.
"""
function set_branchlength_vector!(t::TreeStochastic, blenvec::ArrayStochastic)
    set_branchlength_vector!(t.value, blenvec.value)
end # function

function set_branchlength_vector!(t::N, blenvec::ArrayStochastic) where N<:Node
    set_branchlength_vector!(t, blenvec.value)
end # function

"""
    set_branchlength_vector!(root::N, blenvec::Array{T}) where {N<:GeneralNode, T<:Real}

This function sets the branch lengths of a tree to the values specified in blenvec.
"""
function set_branchlength_vector!(root::N, blenvec::Array{T}) where {N<:GeneralNode, T<:Real}
    any(0 .> blenvec) && throw("this should never happen")
    for child in root.children
        set_branchlength_vector!(child, blenvec)
    end

    @views if root.root !== true
        root.inc_length = blenvec[root.num]
    end
    nothing
end # function set_branchlength_vector!

"""
    get_sum_seperate_length!(root::T)::Vector{Float64}  where T<:GeneralNode

This function gets the sum of the branch lengths of the internal branches and the
branches leading to the leave nodes.
"""
function get_sum_seperate_length!(root::T)::Vector{Float64}  where T<:GeneralNode
    return get_sum_seperate_length!(post_order(root))
end # function get_sum_seperate_length!

"""
    get_sum_seperate_length!(post_order::Vector{T})::Vector{Float64}  where T<:GeneralNode

This function gets the sum of the branch lengths of the internal branches and the
branches leading to the leave nodes.
"""
function get_sum_seperate_length!(post_order::Vector{T})::Vector{Float64}  where T<:GeneralNode
    res_int::Float64 = 0.0
    res_leave::Float64 = 0.0
    res_int_log::Float64 = 0.0
    res_leave_log::Float64 = 0.0
    @simd for node in post_order
        if node.nchild !== 0
            # internal branches
            if !node.root
                res_int += node.inc_length
                res_int_log += log(node.inc_length)
            end
        else
            # branches leading to leaves
            res_leave += node.inc_length
            res_leave_log += log(node.inc_length)
        end # if
    end # for
    return [res_int, res_leave, res_int_log, res_leave_log]
end # function get_sum_seperate_length!

function internal_external_map(t::TreeStochastic)::Vector{Int64}
    internal_external_map(t.value)
end

function internal_external_map(root::T)::Vector{Int64}  where T<:GeneralNode
    internal_external_map(post_order(root))
end

function internal_external_map(post_order::Vector{T})::Vector{Int64}  where T<:GeneralNode
    my_map::Vector{Int64} = zeros(Int64, length(post_order)-1)
    for node in post_order
        if !node.root
            if node.nchild != 0
                my_map[node.num] = 1
            end
        end
    end
    return my_map
end

function internal_external(root::T)::Vector{Int64}  where T<:GeneralNode
    v = root.IntExtMap
    length(v) != 0 && return v
    v = internal_external_map(root)
    root.IntExtMap = v

    return v
end

function find_lca(tree::T, node_l::Array{String, 1})::T  where T<:GeneralNode
    find_lca(tree, [find_by_name(tree, i) for i in node_l])
end

function find_lca(tree::T, node_l::Array{T})::T  where T<:GeneralNode
    @assert length(node_l) > 0
    if length(node_l) === 1
        return node_l[1]
    else
        n1 = popfirst!(node_l)
        n2 = popfirst!(node_l)
        lca = find_lca(tree, n1, n2)
        while length(node_l) !== 0
            n1 = popfirst!(node_l)
            lca = find_lca(tree, lca, n1)
        end
        return lca
    end
end

function find_lca(tree::T, node1::T, node2::T)::T  where T<:GeneralNode
    nb = lcp(node1.binary, node2.binary)
    if nb == node1.binary
        return node1
    end
    if nb == node2.binary
        return node2
    end
    find_binary(tree, nb)
end

"""
    check_binary(root::Node)::Bool

checks to see if given tree is binary; returns true if properly formatted and false otherwise
"""
function check_binary(root::Node)::Bool
    if root.root
        if root.nchild != 2 && root.nchild != 3
            return false
        end #if
    else
        if root.nchild != 0 && root.nchild != 2
            return false
        end #if
    end #else
    res::Bool = true
    for child in root.children
        res &= check_binary(child)
    end #for
    return res
end #function

"""
    check_leafsets(trees::Vector{T})::Nothing where T<:AbstractNode

Checks if an array of trees shares the same leafset (based on the leaf names)
"""
function check_leafsets(trees::Vector{T})::Nothing where T<:AbstractNode
    leaveset = Set([n.name for n in get_leaves(trees[1])])
    count = 0
    for (index, tree) in enumerate(trees[2:end])
        leaveset2 = Set([n.name for n in get_leaves(tree)])
        if leaveset != leaveset2
            println("Tree #$index has a different set of leaves than the first tree")
            count += 1
        end # if
    end # for
    if count != 0
        throw(ArgumentError("$count different trees have a different set of leaves than the first tree"))
    end # if
end # function


text:
using Base.Test
using PolynomialRings
using FGb: FGb_with, groebner

@testset "FGb" begin
    R = @ring! ℤ[x,y,z]

    FGb_with(R) do FGbPolynomial
        f = FGbPolynomial(x^2*y)
        g = FGbPolynomial(x)
        G = groebner([f,g])

        println(G)
        println(map(g->convert(R,g), G))
    end

    S = @ring! ℤ[c[]]
    c1,c2,c3 = c[]

    FGb_with(S, 3) do FGbPolynomial
        f = FGbPolynomial(c1*c2)
        g = FGbPolynomial(c1)
        G = groebner([f,g])

        println(G)
        println(map(g->convert(S,g), G))
    end
end

@testset "Groebner testset" begin
    include(Pkg.dir("PolynomialRings", "test", "Groebner.jl"))
end


text:
using BenchmarkTools
using PlainLinearAlgebra
import PlainLinearAlgebra: diagonal, altdiagonal

using LinearAlgebra: diag

BenchmarkTools.DEFAULT_PARAMETERS.overhead = BenchmarkTools.estimate_overhead()

const SUITE = BenchmarkGroup()

function compactstring(fun)
    io = IOBuffer()
    print(IOContext(io, :compact => true), fun)
    return String(take!(io))
end

###
### diagonal
###

g = addgroup!(SUITE, "diagonal", [])

m2 = rand(2, 2)
m10 = rand(10, 10)
m1000 = rand(1000, 1000)

r1 = addgroup!(g, "reduction", [])

dimstring(a) = string("n=", size(a, 1))
funop(fun, op) = string(compactstring(fun),"(", compactstring(op), "(a))")

for a in (m2, m10, m1000)
    for fun in (sum, prod)
        for diagfun in (diag, diagonal, altdiagonal)
            r1[funop(fun, diagfun), string(eltype(a)), dimstring(a)] = @benchmarkable $fun($diagfun($a))
        end
    end
end

# If a cache of tuned parameters already exists, use it, otherwise, tune and cache
# the benchmark parameters. Reusing cached parameters is faster and more reliable
# than re-tuning `SUITE` every time the file is included.
paramspath = joinpath(dirname(@__FILE__), "params.json")

if isfile(paramspath)
    loadparams!(SUITE, BenchmarkTools.load(paramspath)[1], :evals);
else
    tune!(SUITE)
    BenchmarkTools.save(paramspath, params(SUITE));
end

result = run(SUITE, verbose = true)


text:
using TunedSmoother,BenchmarkTools,Plots,DataStructures
using TunedSmoother: p₀

# Generate WaterLily simulation data
data = create_waterlily()
# Optimize the psuedo-inverse functions
begin # this block takes around 30 minutes. the next block has the results
    smooth! = pseudo!
    scaleloss = zeros(Float32,6,length(data))
    scaleloss[1,:] .= [loss(test;p=zeros(Float32,5),smooth!) for (name,test) ∈ data]
    scaleloss[2,:] .= [loss(test;p=p₀,smooth!) for (name,test) ∈ data]
    opt = OrderedDict(name=>p₀ for name ∈ keys(data))
    for scale ∈ 1:4
        @show scale
        scaledata = scale==4 ? data : create_waterlily(p=scale+2,cases=keys(data))
        opt = OrderedDict(name=>fit(train,opt[name]) for (name,train) ∈ scaledata)
        scaleloss[scale+2,:] .= [loss(test;p=opt[name],smooth!) for (name,test) ∈ data]
    end
end

begin # output from the block above
    opt = OrderedDict(
        circle => [-0.00727588, -0.102282, -0.0178857, 0.760556, -0.586334, 1.78245, 1.14314, 0.167734],
        TGV    => [-0.0617478, 0.00940791, 0.00780108, 0.391629, -0.260235, -1.16419, -0.5565, -0.0615989],
        donut  => [-0.134641, -0.0313802, 0.00234106, 0.286056, -0.159517, -0.021158, -0.213197, -0.0317727],
        wing   => [-0.229121, -0.253223, -0.042952, 0.782714, -0.604938, -0.00708945, 0.0287928, 0.00663962],
        shark  => [-0.279771, -0.281626, -0.0476559, 0.821296, -0.649919, -0.00109436, 0.055443, 0.00988565])

    scaleloss= [-0.46811  -0.63303  -0.64559  -0.52364  -0.58142
                -1.75374  -1.52343  -1.48713  -1.37156  -1.45274
                -1.7893   -1.69872  -1.63804  -1.42732  -1.46177
                -1.81289  -1.73545  -1.64428  -1.42765  -1.48817
                -1.80914  -1.73749  -1.65463  -1.43041  -1.49582
                -1.82367  -1.73825  -1.65887  -1.4328   -1.50249]
end

# plot loss across examples and scales
begin
    using StatsPlots,CategoricalArrays
    cats = ["Jacobi","transfer","⅛","¼","½","1"]
    ctg = CategoricalArray(repeat(cats,inner=length(data)))
    levels!(ctg,cats)
    groupedbar(scaleloss', size = (500,400),
                group=ctg, legend=:outerbottomright, palette=:Greens_6,
                yaxis=("log₁₀ residual reduction",:flip),
                xaxis=("cases",(1:length(data),keys(data))),)
end
savefig("scaleloss.pdf")

begin
    scaledata = create_waterlily(p=4)
    scaledata[wing] = filter(d->itcount(d,GS!)<32,scaledata[wing])
    opt2 = OrderedDict(name=>fit(train,p₀;it=2) for (name,train) ∈ scaledata)
end

begin
    opt2 = OrderedDict(
        circle => [0.143873, 0.114321, 0.0243589, -0.166979, -0.211476],
        TGV    => [-0.445295, -0.175371, -0.0100856, -0.147759, 0.0],
        donut  => [0.121331, 0.0259214, 0.00808, -0.149346, -0.0750463],
        wing   => [0.0259344, -0.123088, -0.0268142, -0.174569, -0.542905],
        shark  => [0.0673015, -0.104415, -0.024194, -0.182964, -0.298915]
    )
    # c = [avecount(test,pseudo!;p=opt2[name]) for (name,test) ∈ data]
    c = [2.0 2.11 2.99 5.99 6.9]
end

# time a single Vcycle
begin
    temp!(st;kw...) = mg!(st;inner=2,mxiter=1,kw...)
    jacobi_time = @belapsed temp!(st;smooth! = Jacobi!) setup=(st=state($data[shark][1]...)) #
    gauss_time = @belapsed temp!(st;smooth! = GS!) setup=(st=state($data[shark][1]...)) #
    sor_time = @belapsed temp!(st;smooth! = SOR!) setup=(st=state($data[shark][1]...)) #
    pseudo_time = @belapsed temp!(st) setup=(st=state($data[shark][1]...)) #
    gauss_time /= jacobi_time
    sor_time /= jacobi_time
    pseudo_time /= jacobi_time
end

# Count the number of cycles needed
crosscount = [[avecount(d,GS!) for (_,d) ∈ data]'.*gauss_time
              [avecount(d,SOR!) for (_,d) ∈ data]'.*sor_time
              [avecount(d,pseudo!) for (_,d) ∈ data]'.*pseudo_time
              [avecount(d,pseudo!;p=opt2[name]) for (name,d) ∈ data]'.*pseudo_time]
crosscount = [  7.0   7.56    10.5     18.585   20.825
                6.4   9.632    9.632   17.92    20.128
                2.74  4.11     4.1511   8.7954  10.1106
                2.74  2.8907   4.0963   8.2063   9.453
            ]

begin
    using StatsPlots,CategoricalArrays
    n = length(data)
    cats = ["Gauss-Sidel","SOR","Ã⁻¹ transfer","Ã⁻¹ tuned-¼"]
    colors = repeat([palette(:default)[2],palette(:default)[4],:lightgreen,palette(:default)[3]],inner=n)
    ctg = CategoricalArray(repeat(cats,inner=n))
    levels!(ctg,cats)
    groupedbar(crosscount',size=(400,400),
                group=ctg, legend=:topleft,c=colors,
                yaxis=("relative time"),
                xaxis=("cases",(1:length(data),keys(data))),)
end
savefig("crosscount.pdf")

# plot best pseudo-inverse functions across examples
begin
    using Plots
    p₀ = Float32[-0.104447, -0.00238399, 0.00841367, -0.158046, -0.115103]
    pmodels(p) = (D->1+p[1]+D*(p[2]+D*p[3]),L->L*(p[4]*(L-2)+p[5]*(L-1)))
    a,_ = pmodels(p₀)
    plot(-6:0.1:0,a,label="transfer",size=(400,400),
         xaxis=("aᵢᵢ"),yaxis=("fᵢᵢ"),legend=:bottomleft)
    for (name,p) in opt2
        x = name ∈ (TGV,donut) ? (-6:0.1:0) : (-4:0.1:0)
        a,_ = pmodels(p)
        plot!(x,a,label=name)
    end
    display(plot!(-6:0,i->1,label="jacobi",c=:grey))
    savefig("diag_fun.png")

    x = 0:0.02:1
    _,a = pmodels(p₀)
    plot(x,a,label="transfer",size=(400,400),
         xaxis=("aᵢⱼ"),yaxis=("fᵢⱼ"),legend=:bottomright)
    for (name,p) in opt2
        _,a = pmodels(p)
        plot!(x,a,label=name)
    end
    display(plot!(x,i->0,label="jacobi",c=:grey))
    savefig("lower_fun.png")
end

begin
    using GeometricMultigrid: Vcycle!
    for (name,h) in ((circle,160),(wing,300),(shark,160))
    # for (name,h) in ((TGV,300),)
    # for (name,h) in ((donut,160),)
            smooth! = pseudo!
        # st = state(data[name][25]...;smooth!)
        st = state(data[name][end]...;smooth!)
        x = copy(st.x.data)
        x .-= st.x[1]
        x[st.A.D .> -1.5] .= NaN
        r = copy(st.r)
        f(x) = @. log10(clamp(abs(x),1e-6,Inf64))
        Vcycle!(st;smooth!)
        plot(heatmap(x[st.x.R]',legend=false,clims=(-1,1)),
        heatmap(f(r.data[st.x.R]'),legend=false,c=:Reds,clims=(-6,-1)),
        heatmap(f(st.r.data[st.x.R]'),legend=false,c=:Reds,clims=(-6,-1)),
        # plot(heatmap(x[st.x.R[:,:,end÷4]]',legend=false,clims=(-1,1)),
        # plot(heatmap(x[st.x.R[:,:,end÷4]]',legend=false),
        # heatmap(f(r.data[st.x.R[:,:,end÷4]]'),legend=false,c=:Reds,clims=(-6,-1)),
        # heatmap(f(st.r.data[st.x.R[:,:,end÷4]]'),legend=false,c=:Reds,clims=(-6,-1)),
        layout = (1,3),size=(900,h),axis=nothing)
        savefig(string(name)*"triple.png")
    end
end
begin
    name = shark
    anim = @animate for d in data[name]
        st = state(d...)
        x = copy(st.x.data)
        x .-= st.x[1]
        x[st.A.D .> -1.5] .= NaN
        r = copy(st.r)
        f(x) = @. log10(clamp(abs(x),1e-6,Inf64))
        plot(heatmap(x[st.x.R]', legend=false, c=:bluesreds, clims=(-1,1)),
             heatmap(f(r.data[st.x.R]'), legend=false, c=:Reds, clims=(-6,-1)),
             layout=(1,2), size=(1000,250), axis=nothing)
    end
    gif(anim,string(name)*".gif")
end


text:
using StatsBase

StatsBase.describe(X::DataValueVector) = StatsBase.describe(STDOUT, X)

function StatsBase.describe(io::IO, X::DataValueVector{T}) where {T <: Real}
    nullcount = sum(X.isna)
    pnull = 100nullcount / length(X)
    if pnull != 100 # describe will fail if dropna returns an empty vector
        describe(io, dropna(X))
    else
        println(io, "Summary Stats:")
        println(io, "Type:           $(eltype(X))")
    end
    println(io, "Number Missing: $(nullcount)")
    @printf(io, "%% Missing:      %.6f\n", pnull)
    return
end

function StatsBase.describe(io::IO, X::DataValueVector)
    nullcount = sum(X.isna)
    pnull = 100nullcount / length(X)
    println(io, "Summary Stats:")
    println(io, "Length:         $(length(X))")
    println(io, "Type:           $(eltype(X))")
    println(io, "Number Unique:  $(length(unique(X)))")
    println(io, "Number Missing: $(nullcount)")
    @printf(io, "%% Missing:      %.6f\n", pnull)
    return
end

function StatsBase.describe(io::IO, X::AbstractVector{DataValue{T}}) where {T <: Real}
    nullcount = sum(isna, X)
    pnull = 100nullcount / length(X)
    if pnull != 100 # describe will fail if dropna returns an empty vector
        describe(io, dropna(X))
    else
        println(io, "Summary Stats:")
        println(io, "Type:           $(eltype(X))")
    end
    println(io, "Number Missing: $(nullcount)")
    @printf(io, "%% Missing:      %.6f\n", pnull)
    return
end

function StatsBase.describe(io::IO, X::AbstractVector{T}) where {T <: DataValue}
    nullcount = sum(isna, X)
    pnull = 100nullcount / length(X)
    println(io, "Summary Stats:")
    println(io, "Length:         $(length(X))")
    println(io, "Type:           $(eltype(X))")
    println(io, "Number Unique:  $(length(unique(X)))")
    println(io, "Number Missing: $(nullcount)")
    @printf(io, "%% Missing:      %.6f\n", pnull)
    return
end


text:
module CompactBasisFunctions

    using ContinuumArrays
    using StaticArrays

    import Base: *
    import ContinuumArrays: Mul, QMul2, QuasiAdjoint, (..), @simplify

    export Derivative, ℵ₁

    include("vandermonde_matrix.jl")

    export Basis, basis, nbasis, grid, degree, order

    include("basis.jl")

    export Bernstein, BernsteinDerivative,
           Chebyshev, ChebyshevDerivative,
           ChebyshevT, ChebyshevTDerivative,
           ChebyshevU, ChebyshevUDerivative,
           Lagrange, LagrangeDerivative,
           LagrangeGauß, LagrangeLobatto,
           Legendre, LegendreDerivative

    include("bernstein.jl")
    include("chebyshev.jl")
    include("lagrange.jl")
    include("legendre.jl")

end


text:
module SunbergTypes

using Random

using CPUTime
using POMDPs
using POMDPPolicies
using DiscreteValueIteration
using ParticleFilters
using POMDPModelTools
using QMDP
using MCTS
using BasicPOMCP
using LaserTag
using LinearAlgebra

using Parameters

import POMDPs: initialstate, generate_s, generate_o, generate_sor, support, discount, isterminal, generate_sr
import POMDPs: actions, n_actions, actionindex, action, dimensions, isterminal_act
import POMDPs: states, n_states, stateindex, transition
import POMDPs: observations, observation, n_observations, obsindex, isterminal_obs
import POMDPs: initialstate, initialstate_distribution
import POMDPs: updater, update
import POMDPs: reward
import POMDPs: convert_s, convert_a, convert_o
import POMDPs: solve

import POMDPs: max_possible_weight, new_particle, reset_distribution

export
    ObsAdaptiveParticleFilter,
    GBMDPSolver,
    MeanRewardBeliefMDP,

    PDPWSolver,

    reset_distribution,
    max_possible_weight,
    new_particle

include("heuristics.jl")
include("mr_belief_mdp.jl")
include("obs_adaptive_pf.jl")
include("gbmdp_solver.jl")
include("rollout.jl")

end # module


text:
using ImageQuilting
using Meshes
using GeoStatsBase
using GeoStatsImages
using Statistics
using Plots; gr(size=(600,400))
using ReferenceTests, ImageIO
using Test, Random

# workaround GR warnings
ENV["GKSwstype"] = "100"

# environment settings
isCI = "CI" ∈ keys(ENV)
islinux = Sys.islinux()
visualtests = !isCI || (isCI && islinux)
datadir = joinpath(@__DIR__,"data")

@testset "ImageQuilting.jl" begin
  @testset "Basic checks" begin
    # the output of a homogeneous image is also homogeneous
    TI = ones(20,20,20)
    reals = iqsim(TI, (10,10,10), size(TI))
    @test reals[1] == TI

    # categories are obtained from training image only
    ncateg = 3; TI = rand(0:ncateg, 20, 20, 20)
    reals = iqsim(TI, (10,10,10), size(TI))
    @test Set(reals[1]) ⊆ Set(TI)
  end

  @testset "Soft data" begin
    # trends with soft data
    TI = [zeros(10,20,1); ones(10,20,1)]
    trend = [zeros(20,10,1) ones(20,10,1)]
    reals = iqsim(TI, (10,10,1), size(TI), soft=[(trend,TI)], tol=1)
    @test mean(reals[1][:,1:10,:]) ≤ mean(reals[1][:,11:20,:])

    # no side effects with soft data
    TI = ones(20,20,20)
    TI[:,5,:] .= NaN
    aux = fill(1.0, size(TI))
    iqsim(TI, (10,10,10), size(TI), soft=[(aux,aux)])
    @test aux == fill(1.0, size(TI))

    # auxiliary variable with integer type
    TI = ones(20,20,20)
    aux = [i for i in 1:20, j in 1:20, k in 1:20]
    iqsim(TI, (10,10,10), size(TI), soft=[(aux,aux)])
    @test aux == [i for i in 1:20, j in 1:20, k in 1:20]
  end

  @testset "Hard data" begin
    # hard data is honored everywhere
    TI = ones(20,20,20)
    obs = rand(size(TI)...)
    data = Dict(CartesianIndex(i,j,k)=>obs[i,j,k] for i=1:20, j=1:20, k=1:20)
    reals = iqsim(TI, (10,10,10), size(TI), hard=data)
    @test reals[1] == obs

    # multiple realizations with hard data
    TI = ones(20,20,20)
    data = Dict(CartesianIndex(20,20,20)=>10)
    reals = iqsim(TI, (10,10,10), size(TI), hard=data, nreal=3)
    for real in reals
      @test real[20,20,20] == 10
    end
  end

  @testset "Masked grids" begin
    # masked simulation domain
    TI = ones(20,20,20)
    shape = Dict{CartesianIndex{3},Real}()
    active = trues(size(TI))
    for i=1:20, j=1:20, k=1:20
      if (i-10)^2 + (j-10)^2 + (k-10)^2 < 25
        push!(shape, CartesianIndex(i,j,k)=>NaN)
        active[i,j,k] = false
      end
    end
    reals = iqsim(TI, (10,10,10), size(TI), hard=shape)
    @test all(isnan.(reals[1][.!active]))
    @test all(.!isnan.(reals[1][active]))

    # masked training image
    TI = ones(20,20,20)
    TI[:,5,:] .= NaN
    reals = iqsim(TI, (10,10,10), size(TI))
    @test reals[1] == fill(1.0, size(TI))
    TI[1,5,:] .= 0
    reals = iqsim(TI, (10,10,10), size(TI))
    @test reals[1] == fill(1.0, size(TI))

    # masked domain and masked training image
    TI = ones(20,20,20)
    TI[:,5,:] .= NaN
    aux = fill(1.0, size(TI))
    shape = Dict(CartesianIndex(i,j,k)=>NaN for i=1:20, j=5, k=1:20)
    reals = iqsim(TI, (10,10,10), size(TI), hard=shape)
    @test all(isnan.(reals[1][:,5,:]))
    @test all(reals[1][:,1:4,:] .== 1)
    @test all(reals[1][:,6:20,:] .== 1)
    reals = iqsim(TI, (10,10,10), size(TI), hard=shape, soft=[(aux,aux)])
    @test all(isnan.(reals[1][:,5,:]))
    @test all(reals[1][:,1:4,:] .== 1)
    @test all(reals[1][:,6:20,:] .== 1)
  end

  @testset "Minimum error cut" begin
    # 3D cut
    TI = ones(20,20,20)
    _, _, voxs = iqsim(TI, (10,10,10), overlap=(1/3,1/3,1/3), debug=true)
    @test 0 ≤ voxs[1] ≤ 1

    A = ones(20,20)
    B = ones(20,20)
    C = ImageQuilting.boykov_kolmogorov_cut(A, B, 1)
    @test all(C[1:end-1,:] .== 1)
    @test all(C[end,:] .== 0)
    C = ImageQuilting.boykov_kolmogorov_cut(A, B, 2)
    @test all(C[:,1:end-1] .== 1)
    @test all(C[:,end] .== 0)
  end

  @testset "Simulation paths" begin
    # different simulation paths
    for kind in [:raster,:dilation,:random]
      rng = MersenneTwister(123)
      path = ImageQuilting.genpath(rng, (10,10,10), kind, Int[])
      @test length(path) == 1000
    end

    # data is visited first if present
    rng = MersenneTwister(123)
    path = ImageQuilting.genpath(rng, (10,10,10), :data, [1,1000])
    @test path[1:2] == [1,1000] || path[1:2] == [1000,1]
  end

  @testset "Voxel reuse" begin
    # mean voxel reuse is in range [0,1]
    TI = rand(20,20,20)
    μ, σ = voxelreuse(TI, (10,10,10), nreal=1)
    @test 0 ≤ μ ≤ 1
  end

  if visualtests
    for (TIname,var) in [("Strebelle",:facies),("StoneWall",:Z)]
      rng = MersenneTwister(2017)
      img = geostatsimage(TIname)
      TI = asarray(img, var)[1:50,1:50,:]
      reals = iqsim(TI, (30,30,1), size(TI), nreal=4, rng=rng)
      ps = [heatmap(real[:,:,1]) for real in reals]
      @test_reference "data/Reals$(TIname).png" plot(ps...)
    end

    for TIname in ["StoneWall","WalkerLake"]
      rng = MersenneTwister(2017)
      img = geostatsimage(TIname)
      TI = asarray(img, :Z)[1:20,1:20,:]
      @test_reference "data/Voxel$(TIname).png" voxelreuseplot(TI, rng=rng)
    end
  end

  @testset "GeoStats.jl API" begin
    sdata   = georef((facies=[1.,0.,1.],), [25. 50. 75.; 25. 75. 50.])
    sdomain = CartesianGrid(100,100)
    problem = SimulationProblem(sdata, sdomain, :facies, 3)

    rng = MersenneTwister(2017)
    trainimg = geostatsimage("Strebelle")
    inactive = [CartesianIndex(i,j) for i in 1:30 for j in 1:30]
    solver = IQ(:facies => (trainimg=trainimg, tilesize=(30,30), inactive=inactive), rng=rng)

    solution = solve(problem, solver)
    @test length(solution) == 3
    @test size(domain(solution[1])) == (100,100)

    incomplete_solver = IQ()
    @test_throws ErrorException solve(problem, incomplete_solver)

    if visualtests
      @test_reference "data/GeoStatsAPI.png" plot(solution,size=(900,300))
    end
  end
end


text:
using Dates

const wday = Dates.Fri
const lo = 1
const hi = 12

print("\nThis script will print the last ", Dates.dayname(wday))
println("s of each month of the year given.")
println("(Leave input empty to quit.)")

while true
    print("\nYear> ")
    y = chomp(readline())
    0 < length(y) || break
    y = try
        parseint(y)
    catch
        println("Sorry, but \"", y, "\" does not compute as a year.")
        continue
    end
    println()
    for m in Date(y, lo):Month(1):Date(y, hi)
        println("    ", tolast(m, wday))
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module SymbolicWedderburn

using LinearAlgebra
using SparseArrays
using Primes

using Cyclotomics
using GroupsCore
using PermutationGroups
using StarAlgebras

export symmetry_adapted_basis, WedderburnDecomposition
export basis,
    # degree, # too common name to export
    direct_summands,
    invariant_vectors,
    issimple,
    multiplicity

macro spawn_compat(expr)
    @static if VERSION < v"1.3.0"
        return :(@async $(esc(expr)))
    else
        return :(Threads.@spawn $(esc(expr)))
    end
end

include("Characters/Characters.jl")
using .Characters
import .Characters: row_echelon_form!
import .Characters.FiniteFields

include("util.jl")
include("actions.jl")
include("action_characters.jl")
include("matrix_projections.jl")
include("minimal_projections.jl")
include("direct_summands.jl")
include("sa_basis.jl")
include("wedderburn_decomposition.jl")

end # module


text:
using StochasticDiffEq, Test, Random, DiffEqDevTools
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_2Dlinear
Random.seed!(100)
prob = prob_sde_2Dlinear

## Solve and plot
println("Solve and Plot")
sol = solve(prob,EM(),dt=1/2^(3))
sol = solve(prob,LambaEM(),dt=1/2^(3))
sol = solve(prob,LambaEulerHeun(),dt=1/2^(3))
sol = solve(prob,RKMil(),dt=1/2^(3))
sol = solve(prob,SRI(),dt=1/2^(3))
sol = solve(prob,SRIW1(),dt=1/2^(3))
sol = solve(prob,SRIW2(),dt=1/2^(3))
sol = solve(prob,SOSRI(),dt=1/2^(3))
sol = solve(prob,SOSRI2(),dt=1/2^(3))
sol = solve(prob,ImplicitEM(),dt=1/2^(3))
sol = solve(prob,ImplicitEM(nlsolve=StochasticDiffEq.NLFunctional()),dt=1/2^(3),adaptive=false)
sol = solve(prob,ImplicitEM(autodiff=false),dt=1/2^(3))
sol = solve(prob,ImplicitRKMil(),dt=1/2^(3))

sol = solve(prob,SRIW1(),dt=1/2^(3),save_everystep=false)
#sol = solve(prob,SRIW1(),dt=1/2^(3),progress=true,progress_steps=1)

sol = solve(prob,SRIW1(),dt=1/2^(3),seed=1)
sol2 = solve(prob,SRI(error_terms=2),dt=1/2^(3),seed=1)
sol.t ≈ sol2.t

sol2 = solve(prob,SRI(tableau=StochasticDiffEq.constructSRIOpt1()),dt=1/2^(3),seed=1)
sol = solve(prob,SOSRI(),dt=1/2^(3),seed=1)
sol.t ≈ sol2.t

## Convergence Testing
println("Convergence Test on 2D Linear")
dts = (1/2) .^ (7:-1:4) #14->7 good plot

sim = test_convergence(dts,prob,EM(),numMonte=1000)
@test abs(sim.𝒪est[:l2]-.5) < 0.1

sim = test_convergence(dts,prob,LambaEM(),numMonte=1000)
@test abs(sim.𝒪est[:l2]-.5) < 0.1

sim = test_convergence(dts,prob,ImplicitEM(),numMonte=100)
@test abs(sim.𝒪est[:l2]-.5) < 0.1

sim = test_convergence(dts,prob,ImplicitEM(theta=1),numMonte=100)
@test abs(sim.𝒪est[:l2]-.5) < 0.1

sim = test_convergence(dts,prob,ImplicitEM(symplectic=true),numMonte=100)
@test abs(sim.𝒪est[:l2]-.5) < 0.1

sim = test_convergence(dts,prob,ImplicitEM(symplectic=true,autodiff=false),numMonte=100)
@test abs(sim.𝒪est[:l2]-.5) < 0.1

sim = test_convergence(dts,prob,ISSEM(),numMonte=100)
@test abs(sim.𝒪est[:l2]-.5) < 0.1

sim = test_convergence(dts,prob,ImplicitRKMil(),numMonte=100)
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim = test_convergence(dts,prob,ImplicitRKMil(theta=1),numMonte=100)
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim = test_convergence(dts,prob,ImplicitRKMil(theta=1,autodiff=false),numMonte=200)
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim = test_convergence(dts,prob,ImplicitRKMil(symplectic=true),numMonte=100)
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim = test_convergence(dts,prob,ImplicitRKMil(symplectic=true,autodiff=false),numMonte=100)
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim = test_convergence(dts,prob,ImplicitRKMil(nlsolve=StochasticDiffEq.NLFunctional()),numMonte=100)
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim2 = test_convergence(dts,prob,RKMil(),numMonte=100)
@test abs(sim2.𝒪est[:l∞]-1) < 0.2

sim3 = test_convergence(dts,prob,SRI(),numMonte=10)
@test abs(sim3.𝒪est[:final]-1.5) < 0.3

sim4 = test_convergence(dts,prob,SRIW1(),numMonte=100)
@test abs(sim4.𝒪est[:final]-1.5) < 0.3

sim4 = test_convergence(dts,prob,SRIW2(),numMonte=100)
@test abs(sim4.𝒪est[:final]-1.5) < 0.3

sim4 = test_convergence(dts,prob,SOSRI(),numMonte=100)
@test abs(sim4.𝒪est[:final]-1.5) < 0.3

sim4 = test_convergence(dts,prob,SOSRI2(),numMonte=100)
@test abs(sim4.𝒪est[:final]-1.5) < 0.3

# 2D oop
f_oop(u,p,t) = u
g_oop(u,p,t) = u
prob = SDEProblem(f_oop, g_oop, ones(2, 2), (0., 1.))
@test_nowarn solve(prob, ImplicitEM())


text:
@testset "group" begin
    @test Dictionaries.group(identity, 11:20) == HashDictionary(11:20, (x->[x]).(11:20))
    @test Dictionaries.group(iseven, 1:10) == dictionary(true => [2,4,6,8,10], false => [1,3,5,7,9])

    @test Dictionaries.group(iseven, x -> x*2, 1:10) == dictionary(true => [4,8,12,16,20], false => [2,6,10,14,18])

    @test Dictionaries.group((x,y) -> iseven(x+y), (x,y) -> x, 1:10, [1,3,4,2,5,6,4,2,3,9]) == dictionary(true => [1,4,5,6,8,9], false => [2,3,7,10])
end

#@testset "groupunique" begin
    #@test Dictionaries.groupunique(identity, 11:20)::HashIndices{Int} == HashIndices(11:20)
    #@test Dictionaries.groupunique(identity, 11:20)::HashIndices{Int} == HashIndices(11:20)
#end
#=
@testset "groupinds" begin
    @test Dictionaries.groupinds(identity, 11:20) == dictionary(Pair.(11:20, (x->[x]).(1:10)))
    @test Dictionaries.groupinds(iseven, 11:20) == dictionary(true => [2,4,6,8,10], false => [1,3,5,7,9])
end

@testset "groupview" begin
    @test Dictionaries.groupview(identity, 11:20)::Groups == group(identity, 11:20)::dictionary
    @test Dictionaries.groupview(iseven, 11:20)::Groups == group(iseven, 11:20)::dictionary
end
=#
@testset "groupreduce" begin
    @test Dictionaries.groupreduce(identity, +, 1:10) == dictionary(Pair.(1:10, 1:10))
    @test Dictionaries.groupreduce(iseven, +, 1:10) == dictionary(true => 30, false => 25)

    @test Dictionaries.groupreduce(iseven, x -> x*2, +, 1:10) == dictionary(true => 60, false => 50)

    @test Dictionaries.groupreduce(iseven, x -> x*2, +, 1:10; init=10) == dictionary(true => 70, false => 60)

    @test Dictionaries.groupreduce((x,y) -> iseven(x+y), (x,y) -> x+y, +, 1:10, 1:10; init=10) == dictionary(true => 120)
    @test Dictionaries.groupreduce((x,y) -> iseven(x+y), (x,y) -> x+y, +, 1:10, [1,3,4,2,5,6,4,2,3,9]; init=10) == dictionary(true => 62, false => 52)

    @test Dictionaries.groupcount(iseven, 1:10) == dictionary(true => 5, false => 5)
    @test Dictionaries.groupsum(iseven, 1:10) == dictionary(true => 30, false => 25)
    @test Dictionaries.groupprod(iseven, 1:10) == dictionary(true => 2*4*6*8*10, false => 1*3*5*7*9)
end


text:
include("electrostatics.jl")
include("Partridge_Schwenke.jl")
include("utilities.jl")

mutable struct TTM21F <: TTM_Potential
    num_waters::Int

    M_site_coords::Vector{SVector{3, Float64}}
    constants::TTM21_Constants
    elec_data::Electrostatics
end

TTM21F(num_waters::Int) = TTM21F(num_waters, [@SVector zeros(3) for _ in 1:(num_waters > 0 ? 4*num_waters : 1)], TTM21_Constants(), Electrostatics(zeros(Int(num_waters * 4)), TTM21_Constants()))
TTM21F() = TTM21F(0)

function evaluate!(ttm21f::TTM21F, coords::AbstractMatrix{Float64}, grads::Union{Matrix{Float64}, Nothing}=nothing)
    """
    Evaluates the ttm21f energy of the system containing water called coords.
    """
    @assert isinteger(size(coords, 2) / 3) "Number of coordinates not divisble by three. Is this water?"
    if (ttm21f.num_waters != Int(size(coords, 2) / 3))
        temp_ttm = TTM21F(Int(size(coords, 2) / 3))
        ttm21f.num_waters = temp_ttm.num_waters
        ttm21f.M_site_coords = temp_ttm.M_site_coords
        ttm21f.elec_data = temp_ttm.elec_data
    end

    compute_M_site_crd!(ttm21f, coords)
    #Rij_vectors = distance_vectors(M_site_coords)
    # @SPEED factor all of the gradient intermediates into a different function or block
    # which is only executed if the gradients are actually requested
    grads_E = [@MVector zeros(3) for _ in 1:3 * 4 * ttm21f.num_waters] # O H H M
    grads_q = [@MArray zeros(3,3,3) for _ in 1:ttm21f.num_waters]
    # calculate the INTRA-molecular energy and Dipole Moment
    # Surface (Partridge-Schwenke)
    E_int = PS_energies_and_gradients!(ttm21f, coords, grads_q, grads, false)

    ttm21f.elec_data.q *= chargecon()
    if grads !== nothing
        grads_q *= chargecon()
    end

    #---------------------------------------------------------------!
    # Calculate the vdW interactions for all atoms                  !
    #---------------------------------------------------------------!

    E_vdw::Float64 = 0.0
    # Loop through all of the waters for van der Waal's interactions
    for i in 1:ttm21f.num_waters-1
        for j in i+1:ttm21f.num_waters
            # Calculate Oxygen-Oxygen interactions
            @inbounds @views Rij = coords[:, 3*(i-1)+1] - coords[:, 3*(j-1)+1]
            expon = ttm21f.constants.vdwD * exp(-ttm21f.constants.vdwE * norm(Rij))

            # Add Evdw van der Waal energy for O-O interaction
            x::Float64 = 1 / (Rij ⋅ Rij)
            E_vdw += expon + x * x * x * (ttm21f.constants.vdwC + x * x * (ttm21f.constants.vdwB + x * ttm21f.constants.vdwA))

            if (grads !== nothing)
                tmp::Float64 = (ttm21f.constants.vdwE * expon / norm(Rij)
                + x*x*x*x*(6 * ttm21f.constants.vdwC
                + x*x*(10 * ttm21f.constants.vdwB + 12 * x * ttm21f.constants.vdwA)))

                # add forces to indices for the M site
                @inbounds @views grads[:, 3*(i-1)+1] -= tmp * Rij
                @inbounds @views grads[:, 3*(j-1)+1] += tmp * Rij
            end
        end
    end

    # if no gradients then just return energy here.
    if (grads === nothing)
        return E_int + E_vdw + electrostatics(ttm21f.elec_data, ttm21f.M_site_coords, nothing)
    end

    #-------------------------------------------------------------------------!
    # Calculate the remaining part of the derivatives                         !
    #-------------------------------------------------------------------------!

    E_elec::Float64 = electrostatics(ttm21f.elec_data, ttm21f.M_site_coords, grads_E)

    # this function call might be slowing things down.
    #
    add_electrostatic_and_dms_derivatives!(ttm21f, grads, grads_E, grads_q)

    return E_int + E_vdw + E_elec
end

max_stars_repo_name: heindelj/MolecularPotentials

max_issues_repo_name: heindelj/MolecularPotentials

max_forks_repo_name: heindelj/MolecularPotentials


text:
# tests for the global system jacobian

"""
  Compare finite difference and complex step jacobians
"""
function test_jac_res()
# check that finite differencing and complex stepping the residual agree

  @testset "----- Testing Jacobian -----" begin
    fname = "input_vals_8el.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    jac_fd = zeros(Float64, 3, 3, mesh.numEl)
    for el = 1:mesh.numEl
      println("----- Doing Finite Differences -----")
  #    jac_fd = zeros(Float64, 3,3)
      eps_fd = 1e-7
      # calculate jacobian of the first element

      fill!(eqn.res, 0.0)
      AdvectionEquationMod.evalResidual(mesh, sbp, eqn, opts)
      res_0 = copy(reshape(eqn.res[1, :, el], 3))
      for i=1:3
        eqn.q[1, i, el] += eps_fd
        fill!(eqn.res, 0.0)
        AdvectionEquationMod.evalResidual(mesh, sbp, eqn, opts)
        res_i = reshape(eqn.res[1, :, el], 3)
        for j=1:3
          jac_fd[j, i, el] = (res_i[j] - res_0[j])/eps_fd
        end

        #undo perturbation
        eqn.q[1, i, el] -= eps_fd
      end
    end

    # now do complex step
    println("----- Doing Complex step -----")
    include(fname)
    arg_dict["jac_method"] = 2
    f = open("input_vals_8elc.jl", "w")
    println(f, arg_dict)
    close(f)
    fname = "input_vals_8elc.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    jac_c = zeros(Float64, 3,3, mesh.numEl)
    eps_c = complex(0, 1e-20)
    for el=1:mesh.numEl
      for i=1:3
        eqn.q[1, i, el] += eps_c
        fill!(eqn.res, 0.0)
        AdvectionEquationMod.evalResidual(mesh, sbp, eqn, opts)
        res_i = reshape(eqn.res[1, :, el], 3)
        for j=1:3
          jac_c[j, i, el] = imag(res_i[j])/abs(eps_c)
        end

        #undo perturbation
        eqn.q[1, i, el] -= eps_c
      end

      @test isapprox( jac_c[:, :, el], jac_fd[:, :, el]) atol=1e-6
    end

    @test isapprox( jac_c, jac_fd) atol=1e-6

    # back to finite differences
    println("----- Testing Finite Difference Jacobian -----")
    fname = "input_vals_8el.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

    # needed for calls to NewtonData below
    Tsol = eltype(eqn.q)
    Tres = eltype(eqn.res)

    # now test full jacobian
    fill!(eqn.res, 0.0)
    AdvectionEquationMod.evalResidual(mesh, sbp, eqn, opts)
    array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
    res_3d0 = copy(eqn.res)
    res_0 = copy(eqn.res_vec)
    jac = zeros(Float64, mesh.numDof, mesh.numDof)
    eps_fd = 1e-7
    fill!(eqn.res, 0.0)
    Jacobian.calcJacFD(mesh, sbp, eqn, opts, AdvectionEquationMod.evalResidual, res_0, eps_fd, jac)

  #  jac_sparse = SparseMatrixCSC(mesh.sparsity_bounds, Float64)
    jac_sparse = SparseMatrixCSC(mesh.sparsity_bnds, Float64)
    println("create jac_sparse")
    fill!(eqn.res, 0.0)
    array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
    Jacobian.calcJacobianSparse(mesh, sbp, eqn, opts, AdvectionEquationMod.evalResidual, res_3d0, eps_fd, jac_sparse)

    jac_sparsefull = full(jac_sparse)
    jac_diff = jac - jac_sparsefull
    for i=1:mesh.numDof
      for j=1:mesh.numDof
        @test isapprox( abs(jac_diff[j, i]), 0.0) atol=1e-6
      end
    end

    # back to complex step
    println("----- Testing Complex Step Jacobian -----")
    fname = "input_vals_8elc.jl"
    arg_dict["jac_method"] = 2  # something screwy is going on because this is necessary
    mesh, sbp, eqn, opts = solvePDE(fname)
    array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

    # now test full jacobian
    jac_c = zeros(Complex128, mesh.numDof, mesh.numDof)
    eps_c = complex(0, 1e-20)
    fill!(eqn.res, 0.0)
    Jacobian.calcJacobianComplex(mesh, sbp, eqn, opts, AdvectionEquationMod.evalResidual, eps_c, jac_c)

  #  jac_csparse = SparseMatrixCSC(mesh.sparsity_bounds, Float64)
    jac_csparse = SparseMatrixCSC(mesh.sparsity_bnds, Float64)
    fill!(eqn.res, 0.0)
    res_3d0 = Array{Float64}(0, 0, 0)
    array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
    Jacobian.calcJacobianSparse(mesh, sbp, eqn, opts, AdvectionEquationMod.evalResidual, res_3d0, eps_c, jac_csparse)

    jac_csparsefull = full(jac_csparse)
    jac_diff = jac_c - jac_csparsefull
    for i=1:mesh.numDof
      for j=1:mesh.numDof
        @test isapprox( abs(jac_diff[j, i]), 0.0) atol=1e-12
      end
    end

  end  # end facts block

  return nothing
end

#test_jac_res()
add_func1!(AdvectionTests, test_jac_res, [TAG_SHORTTEST])

"""
  Test the various methods of calculating the jacobian
"""
function test_jac_calc()
  @testset "----- Testing Jacobian calculation -----" begin
    # back to finite differences
    println("----- Testing Finite Difference Jacobian -----")
    fname = "input_vals_8el.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

    # needed for calls to NewtonData below
    Tsol = eltype(eqn.q)
    Tres = eltype(eqn.res)

    # now test full jacobian
    fill!(eqn.res, 0.0)
    AdvectionEquationMod.evalResidual(mesh, sbp, eqn, opts)
    array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
    res_3d0 = copy(eqn.res)
    res_0 = copy(eqn.res_vec)
    jac = zeros(Float64, mesh.numDof, mesh.numDof)
    eps_fd = 1e-7
    fill!(eqn.res, 0.0)
    Jacobian.calcJacFD(mesh, sbp, eqn, opts, AdvectionEquationMod.evalResidual, res_0, eps_fd, jac)

  #  jac_sparse = SparseMatrixCSC(mesh.sparsity_bounds, Float64)
    jac_sparse = SparseMatrixCSC(mesh.sparsity_bnds, Float64)
    fill!(eqn.res, 0.0)
    array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
    Jacobian.calcJacobianSparse(mesh, sbp, eqn, opts, AdvectionEquationMod.evalResidual, res_3d0, eps_fd, jac_sparse)

    jac_sparsefull = full(jac_sparse)
    jac_diff = jac - jac_sparsefull
    for i=1:mesh.numDof
      for j=1:mesh.numDof
        @test isapprox( abs(jac_diff[j, i]), 0.0) atol=1e-6
      end
    end

    # back to complex step
    println("----- Testing Complex Step Jacobian -----")
    fname = "input_vals_8elc.jl"
    arg_dict["run_type"] = 5  # something screwy is going on because this is necessary
    arg_dict["jac_method"] = 2
    mesh, sbp, eqn, opts = solvePDE(fname)
    array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

    # now test full jacobian
    jac_c = zeros(Complex128, mesh.numDof, mesh.numDof)
    eps_c = complex(0, 1e-20)
    fill!(eqn.res, 0.0)
    Jacobian.calcJacobianComplex(mesh, sbp, eqn, opts, AdvectionEquationMod.evalResidual, eps_c, jac_c)

  #  jac_csparse = SparseMatrixCSC(mesh.sparsity_bounds, Float64)
    jac_csparse = SparseMatrixCSC(mesh.sparsity_bnds, Float64)
    fill!(eqn.res, 0.0)
    res_3d0 = Array{Float64}(0, 0, 0)
    array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
    Jacobian.calcJacobianSparse(mesh, sbp, eqn, opts, AdvectionEquationMod.evalResidual, res_3d0, eps_c, jac_csparse)

    jac_csparsefull = full(jac_csparse)
    jac_diff = jac_c - jac_csparsefull
    for i=1:mesh.numDof
      for j=1:mesh.numDof
        @test isapprox( abs(jac_diff[j, i]), 0.0) atol=1e-12
      end
    end

    # now check FD vs Complex step
    for i=1:mesh.numDof
      for j=1:mesh.numDof
        @test isapprox( abs(jac_c[i, j] - jac[i,j]), 0.0) atol=1e-6
      end
    end

  end  # end facts block

  return nothing
end

#test_jac_calc()
add_func1!(AdvectionTests, test_jac_calc, [TAG_SHORTTEST])


text:
using Catlab
using Catlab.Theories

@signature ZXCategory{Ob,Hom} <: DaggerCompactCategory{Ob,Hom} begin
    # Argument α is the phase, usually <: Real
    zphase(A::Ob, α)::(A → A)
    zcopy(A::Ob, α)::(A → (A⊗A))
    zdelete(A::Ob, α)::(A → munit())
    zmerge(A::Ob, α)::((A⊗A) → A)
    zcreate(A::Ob, α)::(munit() → A)

    xphase(A::Ob, α)::(A → A)
    xcopy(A::Ob, α)::(A → (A⊗A))
    xdelete(A::Ob, α)::(A → munit())
    xmerge(A::Ob, α)::((A⊗A) → A)
    xcreate(A::Ob, α)::(munit() → A)

    hadamard(A::Ob)::(A → A)
end

# Convenience methods for phaseless spiders.
zcopy(A) = zcopy(A,0)
zdelete(A) = zdelete(A,0)
zmerge(A) = zmerge(A,0)
zcreate(A) = zcreate(A,0)

xcopy(A) = xcopy(A,0)
xdelete(A) = xdelete(A,0)
xmerge(A) = xmerge(A,0)
xcreate(A) = xcreate(A,0);

import Catlab.Theories.Ob
@syntax ZXCalculus{ObExpr,HomExpr} ZXCategory begin
    # otimes(A::Ob, B::Ob) = associate_unit(new(A,B), munit)
    # otimes(f::Hom, g::Hom) = associate(new(f,g))
    # compose(f::Hom, g::Hom) = associate(new(f,g; strict=true))
end

using Metatheory, Metatheory.EGraphs
@metatheory_init

# Custom type APIs for the GATExpr
using TermInterface
TermInterface.operation(t::ObExpr) = :call
TermInterface.arguments(t::ObExpr) = [head(t), t.args...]
TermInterface.operation(t::HomExpr) = :call
TermInterface.arguments(t::HomExpr) = [head(t), t.args...]

abstract type CatType end
struct ObType <: CatType
    ob
    mod
end
struct HomType <: CatType
    dom
    codom
    mod
end

# Type information will be stored in the metadata
function TermInterface.metadata(t::HomExpr)
    return HomType(t.type_args[1], t.type_args[2], typeof(t).name.module)
end
TermInterface.metadata(t::ObExpr) = ObType(t, typeof(t).name.module)
TermInterface.istree(t::GATExpr) = true
TermInterface.arity(t::GATExpr) = length(arguments(t))

struct CatlabAnalysis <: AbstractAnalysis end
function EGraphs.make(an::Type{CatlabAnalysis}, g::EGraph, n::ENode{T}) where T
    !(T <: GATExpr) && return t
    return metadata(n)
end
EGraphs.join(an::Type{CatlabAnalysis}, from, to) = from
EGraphs.islazy(x::Type{CatlabAnalysis}) = false

function infer(t::GATExpr)
    g = EGraph(t)
    analyze!(g, CatlabAnalysis)
    getdata(geteclass(g, g.root), CatlabAnalysis)
end

function EGraphs.extractnode(g::EGraph, n::ENode{T}, extractor::Function) where {T <: ObExpr}
    @assert n.head == :call
    return metadata(n).ob
end

function EGraphs.extractnode(g::EGraph, n::ENode{T}, extractor::Function) where {T <: HomExpr}
    @assert n.head == :call
    nargs = extractor.(n.args)
    nmeta = metadata(n)
    return nmeta.mod.Hom{nargs[1]}(nargs[2:end], GATExpr[nmeta.dom, nmeta.codom])
end

# function EGraphs.instantiateterm(g::EGraph, pat::PatTerm,  T::Type{H{K}}, sub::Sub, rule::Rule) where {H <: GATExpr, K}
# # TODO
# end

t = Metatheory.@theory begin
    compose(hadamard(A), hadamard(A)) |>
    begin
        d = getdata(A, Main.CatlabAnalysis)
        return d.mod.id(d.ob)
    end
    compose(f, id(B)) |>
    begin
        bd = getdata(B, CatlabAnalysis)
        fd = getdata(f, CatlabAnalysis)
        if bd.ob == fd.codom
            return f
        else
            error("TYPE ERROR!")
            return _lhs_expr
        end
    end
    compose(id(A), f) |>
    begin
        ad = getdata(A, CatlabAnalysis)
        fd = getdata(f, CatlabAnalysis)
        if ad.ob == fd.dom
            return f
        else
            error("TYPE ERROR!")
            return _lhs_expr
        end
    end
end

t[2]

A = Ob(ZXCalculus.Ob, :A)
B = Ob(ZXCalculus.Ob, :B)
f = Hom(:f, A, B)
h = hadamard(A)
c = h ⋅ h
G = EGraph(c)
infer(zdelete(A)).codom == A

analyze!(G, CatlabAnalysis)
saturate!(G, t)
ex = extract!(G, astsize)
ex == id(A)

G = EGraph(f ⋅ id(B))
analyze!(G, CatlabAnalysis)
saturate!(G, t)
ex = extract!(G, astsize)
ex == f

x = id(A) ⋅ f ⋅ id(B)
G = EGraph(x)
analyze!(G, CatlabAnalysis)
saturate!(G, t)
ex = extract!(G, astsize)
ex == f

using Catlab, Catlab.Theories
using Catlab.WiringDiagrams, Catlab.Graphics
using Catlab.Syntax

A, B, C, D, E = Ob(FreeBiproductCategory, :A, :B, :C, :D, :E)
f = Hom(:f, A, B)
g = Hom(:g, B, C)
h = Hom(:h, B, A)
k = Hom(:k, C, B)
x = id(A) ⋅ f ⋅ id(B)

z = x ⊗ f ⊗ ((f ⊗ g) ⋅ braid(B,C) ⋅ (k ⊗ h) ⋅ (delete(B) ⊗ f))
to_composejl(z; orientation=LeftToRight)

drop = munit(FreeCompactClosedCategory.Ob)
delete()


text:
#
# WIP: work in progress on deepvalidate. This code is not expected to work at this time
#

#
# This file copies much of the same functionality from varextract.jl,
# with important modifications to write out test traces to a file -
# "traces.dat" - for subsequent classification modeling.
#

using Pkg

# packages = ['Distributions','DelimitedFiles','Cassette','Test']

# for package in packages
#     haskey(Pkg.installed(),package) || Pkg.add(package)

using Distributions
using DelimitedFiles
using Cassette
using Test
Cassette.@context TraceCtx

"""    varname(ir::Core.CodeInfo, sym::Symbol)
look up the name of a slot from the codeinfo slotnames.
see also: ir.slotnames
"""
function varname(ir::Core.CodeInfo, sym::Symbol)
    s = string(sym)[2:end]
    i = parse(Int,s)
    varname = ir.slotnames[i]
    return varname
end

"""    Extraction
- ir: a CodeInfo object that we are extracting from
- varnames: Variable names used as left hand sides
- funccalls: Tuples of (returnvar, funcname)
- literals: Literal values used as right hand sides in assignment
- SSAassigns: Locations used for the storage of subexpression results ie (2+2*4), the value 8 is stored in an SSAassign
"""
struct Extraction
    ir::Core.CodeInfo
    varnames::Vector{Any}
    funccalls::Vector{Any}
    literals::Vector{Any}
    SSAassigns::Vector{Any}
end

"""     Extraction(ir)
construct an Extraction object from a piece of code info
"""
function Extraction(ir)
    return Extraction(ir, Symbol[], Any[], Any[], Any[])
end

function findvars(ext, ir, expr)
    @info "Finding Variables"
    # @show ir
    # dump(expr)
    # if typeof(expr) <: SSAValue
    #     return
    # end

    vars = Any[]
    try
        args = expr.args
        for arg in args
            # @show arg
            if typeof(arg) <: Core.SlotNumber
                push!(vars, varname(ir, Symbol(arg)))
            elseif typeof(arg) <: GlobalRef
                # @show arg
                continue
            # elseif typeof(arg) <: SSAValue
            #     continue
            else
                push!(vars, findvars(ext,ir,arg))
            end
        end
    catch
        dump(expr)
        return varname(ir, Symbol(expr))
    end

    return vars
end

# add an expression to the Extraction struct by parsing out the relevant info.
function Base.push!(ext::Extraction, expr::Expr)
    ir = ext.ir
    if expr.head == :(=)
        # @show expr
        try
            # @show expr.args
            sym = Symbol(expr.args[1])
            vn = varname(ir, sym)
            # vntree = findvars(ext, ir, expr)
            # push!(ext.varnames, vntree)
            push!(ext.varnames, vn)
            if isa(expr.args[2], Expr)
                if expr.args[2].head == :(call)
                    fname = expr.args[2].args[1]
                    push!(ext.funccalls, (vn, fname))
                else
                    @warn "No method to handle Non-call Expr as RHS: $(expr.args[2])"
                end
            elseif isa(expr.args[2], Core.SSAValue)
                @info "hit an SSAValue $vn = $(expr.args[2])"
                push!(ext.SSAassigns, expr)
            else
                @info "Residual Clause: $(expr): $(expr.args[2])"
                push!(ext.literals, expr)
            end
        catch ex
            @show ex
            # @warn "could not find slotname for $(expr.args[1])"
            # @show slotnames
            # @show expr.args[1], s, i
        end
    end
end

function Base.show(ext::Extraction)
    if length(ext.varnames) > 0
        @show ext.varnames
    end
    if length(ext.funccalls) > 0
        @show ext.funccalls
    end
    if length(ext.literals) > 0
        @show ext.literals
    end
    if length(ext.SSAassigns) > 0
        @show ext.SSAassigns
    end
end

"""    extractpass(::Type{<:TraceCtx}, reflection::Cassette.Reflection)
is a Cassette pass to log the varnames and function calls to build the dynamic code graph
part of the SemanticModels knowledge graph.
"""
function extractpass(::Type{<:TraceCtx}, reflection::Cassette.Reflection)
    ir = reflection.code_info
    slotnames = ir.slotnames
    vn = slotnames[end]
    ext = Extraction(ir)
    s = ""
    i = 1
    modname = ir.linetable[1].mod
    methname = ir.linetable[1].method
    for expr in ir.code
        if expr == nothing
            continue
        end

        if !isa(expr, Expr) && expr != nothing
            @debug "Found non expression. $expr"
            continue
        end
        push!(ext, expr)
    end
    if length(ext.varnames) > 0
        @info "Working with method: $(modname).$(methname)"
        show(ext)
        # @show ext.ir
        @show ext.ir.slotnames
        println((modul=modname, func=methname, varnames=ext.ir.slotnames, callees=ext.funccalls))
    end
    # TODO: do something with functions without explicit assignment.
    # if length(varnames) == 0
    #     @show ir
    # end
    return ir

end

const ExtractPass = Cassette.@pass extractpass

# defines the set of methods we do NOT want to descend into.
function Cassette.canrecurse(ctx::TraceCtx,
                             f::Union{typeof(+), typeof(*), typeof(/), typeof(-),typeof(Base.iterate),
                                      typeof(Base.sum),
                                      typeof(Base.mapreduce),
                                      typeof(Base.Broadcast.copy),
                                      typeof(Base.Broadcast.instantiate),
                                      typeof(Base.Math.throw_complex_domainerror),
                                      typeof(Base.Broadcast.broadcasted)},
                             args...)
    return false
end

#
# This function has been modified from the version in varextract.jl
# in the first line, where instead of "@show f, args" we write these
# objects to our output file "traces.dat"
#
# Not elegant, but various other approaches did not write out the
# correct trace information. This should be updated with a more elegant
# solution when possible.
#

function Cassette.overdub(ctx::TraceCtx,
                          f,
                          args...)
    open("traces.dat", "a") do file
        write(file, string(f))
        write(file, string(args))
    end

    # if we are supposed to descend, we call Cassette.recurse
    if Cassette.canrecurse(ctx, f, args...)
        subtrace = (Any[],Any[])
        push!(ctx.metadata[1], (f, args) => subtrace)
        newctx = Cassette.similarcontext(ctx, metadata = subtrace)
        retval = Cassette.recurse(newctx, f, args...)
        # push!(ctx.metadata[2], subtrace[2])
    else
        retval = Cassette.fallback(ctx, f, args...)
        push!(ctx.metadata[1], :t)
        push!(ctx.metadata[2], retval)
    end
    # @info "returning"
    # @show retval
    return retval
end

function add(a, b)
    c = a + b
    return c
end

#
# This testset produces 30,000 traces on these basic arithmetic
# functions, and writes them to a file called "traces.dat"
#

@testset "TraceExtract" begin
    g(x) = begin
        y = add(x.*x, -x)
        z = 1
        v = y .- z
        s = sum(v)
        return s
    end
    h(x) = begin
        z = g(x)
        zed = sqrt(z)
        return zed
    end

    open("traces.dat", "w") do f
        write(f, "")
    end

    # Error conditions happen when our inputs are sufficiently small, so
    # Normal(0,2) gives us a good range of values to generate a reasonable
    # percentage of "bad" traces on which to train. Empirically the share
    # of "bad" traces is about 15-17%.

    seeds = rand(Normal(0,2),3000,3)

    for i=1:size(seeds,1)
        ctx = TraceCtx(pass=ExtractPass, metadata = (Any[], Any[]))
        try
            result = Cassette.overdub(ctx, h, seeds[i,:])
        catch DomainError
            dump(ctx.metadata)
        finally
            open("traces.dat", "a") do f
                write(f, "\n")
            end
        end
        if i%1000 == 0
            @info string(i)
        end
    end
end


text:
using TopologyPreprocessing
using Test
using Random
Random.seed!(1234)

# using MatrixOrganization

@testset "MatrixOrganization.jl -> matrix pooling" begin
	in_vector = [1 2 3 4 3 2 1]
	in_matrix = [1 2 3; 5 6 7; 8 9 0]
	in_matrix2 = [1 2 3; 5 6 7; 8 9 0]

	sqr_matrix0 = [1 2 3; 2 6 7; 3 7 0]
	resulting_matrix0_2 = [1 2 3; 2 6 7; 3 7 0]

	sqr_matrix1 = [	0	1	13	4	5	9;
					1	0	2	14	6	10;
					13	2	0	3	7	11;
					4	14	3	0	8	12;
					5	6	7	8	0	15;
					9	10	11	12	15	0]
	resulting_matrix1_2 = [0   1  14  14  10  10;
						 1   0  14  14  10  10;
						 14  14   0   3  12  12;
						 14  14   3   0  12  12;
						 10  10  12  12   0  15;
						 10  10  12  12  15   0]

	 sqr_matrix2 = [ 0	1	113	4	5	9	13	17	81	82	83	84;
						1	0	2	114	6	10	14	18	85	86	87	88;
						113	2	0	3	7	11	15	19	89	90	91	92;
						4	114	3	0	8	12	16	20	93	94	95	96;
						5	6	7	8	0	21	115	24	29	30	31	32;
						9	10	11	12	21	0	22	116	33	34	35	36;
						13	14	15	16	115	22	0	23	37	38	39	40;
						17	18	19	20	24	116	23	0	41	42	43	44;
						81	85	89	93	29	33	37	41	0	25	117	28;
						82	86	90	94	30	34	38	42	25	0	26	118;
						83	87	91	95	31	35	39	43	117	26	0	27;
						84	88	92	96	32	36	40	44	28	118	27	0]
		result_matrix2_2 = [   0    1  114  114  10  10   18   18   86   86   88   88;
								  1    0  114  114  10  10   18   18   86   86   88   88;
								114  114    0    3  12  12   20   20   94   94   96   96;
								114  114    3    0  12  12   20   20   94   94   96   96;
								 10   10   12   12   0  21  116  116   34   34   36   36;
								 10   10   12   12  21   0  116  116   34   34   36   36;
							  	 18   18   20   20 116 116    0   23   42   42   44   44;
								 18   18   20   20 116 116   23    0   42   42   44   44;
								 86   86   94   94  34  34   42   42    0   25  118  118;
								 86   86   94   94  34  34   42   42   25    0  118  118;
								 88   88   96   96  36  36   44   44  118  118    0   27;
								 88   88   96   96  36  36   44   44  118  118   27    0]
		result_matrix2_3 = [ 0	1	113	114	114	114	89	89	89	92	92	92;
								1	0	2	114	114	114	89	89	89	92	92	92;
								113	2	0	114	114	114	89	89	89	92	92	92;
								114	114	114	0	8	12	116	116	116	96	96	96;
								114	114	114	8	0	21	116	116	116	96	96	96;
								114	114	114	12	21	0	116	116	116	96	96	96;
								89	89	89	116	116	116	0	23	37	117	117	117;
								89	89	89	116	116	116	23	0	41	117	117	117;
								89	89	89	116	116	116	37	41	0	117	117	117;
								92	92	92	96	96	96	117	117	117	0	26	118;
								92	92	92	96	96	96	117	117	117	26	0	27;
								92	92	92	96	96	96	117	117	117	118	27	0]
		result_matrix2_4 = [ 0	1	114	114	20	20	20	20	96	96	96	96;
								1	0	114	114	20	20	20	20	96	96	96	96;
								114	114	0	3	20	20	20	20	96	96	96	96;
								114	114	3	0	20	20	20	20	96	96	96	96;
								20	20	20	20	0	21	116	116	44	44	44	44;
								20	20	20	20	21	0	116	116	44	44	44	44;
								20	20	20	20	116	116	0	23	44	44	44	44;
								20	20	20	20	116	116	23	0	44	44	44	44;
								96	96	96	96	44	44	44	44	0	25	118	118;
								96	96	96	96	44	44	44	44	25	0	118	118;
								96	96	96	96	44	44	44	44	118	118	0	27;
								96	96	96	96	44	44	44	44	118	118	27	0]

	@test matrix_poling(in_vector) !== in_vector
	@test matrix_poling(in_vector, method="max_pooling") == [4 4 4 4 4 4 4]
	# @test matrix_poling!(in_vector) === in_vector
	# @test matrix_poling!(in_vector) == [4 4 4 4 4 4 4]

	@test matrix_poling(in_matrix) !== in_matrix
	@test matrix_poling(in_matrix, method="max_pooling") == 9 .* ones(size(in_matrix))
	# @test matrix_poling!(in_matrix) === in_matrix
	# @test matrix_poling!(in_matrix) == 9 .* ones(size(in_matrix))

	@test matrix_poling(in_matrix2[1:2,1:2], method="max_pooling") == 6 .* ones(size(in_matrix2[1:2,1:2]))
	@test matrix_poling(in_matrix2[1:2,1:2], method="max_pooling") != in_matrix2[1:2,1:2]

	# ====
	# Subsampling matrix

	# Function is supposed to work only on upper half, and here the upper half is too small, so there are no operations
	@test subsample_matrix(sqr_matrix0, subsamp_size=2, method="max_pooling") == resulting_matrix0_2

	@test subsample_matrix(sqr_matrix1, subsamp_size=2, method="max_pooling") == resulting_matrix1_2

	@test subsample_matrix(sqr_matrix2, subsamp_size=2, method="max_pooling") == result_matrix2_2
	@test subsample_matrix(sqr_matrix2, subsamp_size=3, method="max_pooling") == result_matrix2_3
	@test subsample_matrix(sqr_matrix2, subsamp_size=4, method="max_pooling") == result_matrix2_4
end

@testset "MatrixOrganization.jl -> add_random_patch" begin
	# TODO set seed for add_random_path
	# TODO Seed has to be set for this test
	in_vector = [1, 2, 3, 4, 3, 2, 1]

	sqr_matrix0 = [ 1 2 3;
	 				2 6 7;
					3 7 0]
	sqr_matrix1 = [1 2 3  4  5;
					 2 1 6  7  8;
					 3 6 1  9  10;
					 4 7 9  1  11;
					 5 8 10 11 1]
	sqr_matrix2 = [	0	1	13	4	5	9;
					1	0	2	14	6	10;
					13	2	0	3	7	11;
					4	14	3	0	8	12;
					5	6	7	8	0	15;
					9	10	11	12	15	0]
	 sqr_matrix3 = [ 0	1	113	4	5	9	13	17	81	82	83	84;
 					1	0	2	114	6	10	14	18	85	86	87	88;
 					113	2	0	3	7	11	15	19	89	90	91	92;
 					4	114	3	0	8	12	16	20	93	94	95	96;
 					5	6	7	8	0	21	115	24	29	30	31	32;
 					9	10	11	12	21	0	22	116	33	34	35	36;
 					13	14	15	16	115	22	0	23	37	38	39	40;
 					17	18	19	20	24	116	23	0	41	42	43	44;
 					81	85	89	93	29	33	37	41	0	25	117	28;
 					82	86	90	94	30	34	38	42	25	0	26	118;
 					83	87	91	95	31	35	39	43	117	26	0	27;
 					84	88	92	96	32	36	40	44	28	118	27	0]
	function get_unchanged_indices(input_matrix,ind)
		indices = CartesianIndices(size(input_matrix))
		indices = findall(x->x!=ind,indices)
		for i = ind
			indices = indices[findall(x->x!=i,indices)]
		end
		return indices
	end

    out_m, ind = add_random_patch(sqr_matrix0)
	indices = get_unchanged_indices(sqr_matrix0,ind)
	@test size(ind) == (1,2)
	@test sqr_matrix0[indices] == out_m[indices]

	big_sqr_matrix0 = sqr_matrix0 .*100
	out_m, ind = add_random_patch(big_sqr_matrix0, patch_size=1,total_patches=2)
	indices = get_unchanged_indices(big_sqr_matrix0,ind)
	@test size(ind) == (2,2)
	@test big_sqr_matrix0[indices] == out_m[indices]
	@test sum(big_sqr_matrix0[ind] .!= out_m[ind]) == 4
	@test sum(big_sqr_matrix0[ind] .== out_m[ind]) == 0

	out_m, ind = add_random_patch(sqr_matrix1, patch_size=1,total_patches=2)
	indices = get_unchanged_indices(sqr_matrix1,ind)
	@test size(ind) == (2,2)
	@test sqr_matrix1[indices] == out_m[indices]

	# TODO those 2 tests fails when random value is the equal to one that is replaced
	# @test sum(sqr_matrix1[ind] .!= out_m[ind]) == 4
	# @test sum(sqr_matrix1[ind] .== out_m[ind]) == 0

	# ===
	input_matrix = sqr_matrix1
	function test_adding_rand_patch(input_matrix, t_patches,p_size)
		out_m, ind = add_random_patch(input_matrix, patch_size=p_size, total_patches=t_patches)
		indices = get_unchanged_indices(input_matrix,ind)
		@test size(ind) == (t_patches*p_size^2,2)
		@test input_matrix[indices] == out_m[indices]
		# For values from range, tests below does not make sense:
		# @test sum(input_matrix[ind] .!= out_m[ind]) == length(ind)
		# @test sum(input_matrix[ind] .== out_m[ind]) == 0
	end
	t_patches = 1
	p_size = 2
		test_adding_rand_patch(sqr_matrix0, t_patches,p_size)
		test_adding_rand_patch(sqr_matrix1, t_patches,p_size)
		test_adding_rand_patch(sqr_matrix2, t_patches,p_size)
		test_adding_rand_patch(sqr_matrix3, t_patches,p_size)

	t_patches = 1
	p_size = 3
		test_adding_rand_patch(sqr_matrix0, t_patches,p_size)
		test_adding_rand_patch(sqr_matrix1, t_patches,p_size)
		test_adding_rand_patch(sqr_matrix2, t_patches,p_size)
		test_adding_rand_patch(sqr_matrix3, t_patches,p_size)

	t_patches = 1
	p_size = 4
	correct_error = 3
    # TODO change this into test_throws
	try
		add_random_patch(sqr_matrix0, patch_size=p_size, total_patches=t_patches)
	catch err
		my_error = 0
		if isa(err, DomainError)
			println("DomainError")
			my_error = 1
		elseif isa(err, DimensionMismatch)
			println("DimensionMismatch")
			my_error = 2
		else
			println("Unknow error")
			my_error = 3
		end
		@test my_error == correct_error
	end

	test_adding_rand_patch(sqr_matrix1, t_patches, p_size)
	test_adding_rand_patch(sqr_matrix2, t_patches, p_size)
	test_adding_rand_patch(sqr_matrix3, t_patches, p_size)

	t_patches = 3
	p_size = 5
	test_adding_rand_patch(sqr_matrix2, t_patches,p_size)
	test_adding_rand_patch(sqr_matrix3, t_patches,p_size)

	# ===
	# locations
	locations1 = [CartesianIndex(1,2), CartesianIndex(2,3)]
	locations2 = [CartesianIndex(1,2), CartesianIndex(99,3)]

	t_patches = 1
	p_size = 1
	out_m, ind = add_random_patch(sqr_matrix1, patch_size=p_size, total_patches=t_patches,locations=locations1)
		indices = get_unchanged_indices(sqr_matrix1,ind)
		@test ind[:,1] == locations1
		@test size(ind) == (size(locations1)[1]*p_size^2,2)
		@test sqr_matrix1[indices] == out_m[indices]
		# The number of
		@test sum(sqr_matrix1[locations1] .!= out_m[locations1]) == length(locations1)
		@test sum(sqr_matrix1[locations1] .== out_m[locations1]) == 0

		correct_error = 0
		try
			out_m, ind = add_random_patch(sqr_matrix1, patch_size=p_size, total_patches=t_patches,locations=locations2)
		catch err
			# global correct_error
			if isa(err, DomainError)
				correct_error = 1
			else
				correct_error = 2
			end
		finally
			# global correct_error
			@test correct_error == 2
		end
	# TODO test for index below diagonal
	# TODO too many indices
end


text:
const StructTypes = JSON3.StructTypes

# StructTypes.StructType(::Type{T}) where T <: HasFields = StructTypes.DictType()
# StructTypes.construct(::Type{T}, x::AbstractDict; kw...) where T <: HasFields = T(_symbol_dict(x); kw...)
# StructTypes.StructType(::Type{T}) where T <: Union{Template} = StructTypes.Struct()
# StructTypes.StructType(::Type{T}) where T <: Union{Plot,PlotConfig} = StructTypes.Mutable()
# StructTypes.omitempties(::Type{PlotConfig}) = fieldnames(PlotConfig)

StructTypes.StructType(::Type{<:PlotlyBase.Plot}) = JSON3.RawType()
JSON3.rawbytes(plot::PlotlyBase.Plot) = codeunits(JSON.json(plot))


text:
struct AccelDesang
    a::Float64
    ϕdes::Float64
end
Base.show(io::IO, a::AccelDesang) = @printf(io, "AccelDesang(%6.3f,%6.3f)", a.a, a.ϕdes)
Base.length(::Type{AccelDesang}) = 2
Base.convert(::Type{AccelDesang}, v::Vector{Float64}) = AccelDesang(v[1], v[2])
function Base.copyto!(v::Vector{Float64}, a::AccelDesang)
    v[1] = a.a
    v[2] = a.ϕdes
    v
end
function propagate(veh::Vehicle, action::AccelDesang, roadway::Roadway, Δt::Float64; n_integration_steps::Int=4)

    a = action.a # accel
    ϕdes = action.ϕdes # desired heading angle

    x = veh.state.posG.x
    y = veh.state.posG.y
    θ = veh.state.posG.θ
    v = veh.state.v

    δt = Δt/n_integration_steps

    for i in 1 : n_integration_steps

        posF = Frenet(VecSE2(x, y, θ), roadway)
        ω = ϕdes - posF.ϕ

        x += v*cos(θ)*δt
        y += v*sin(θ)*δt
        θ += ω*δt
        v += a*δt
    end

    posG = VecSE2(x, y, θ)
    VehicleState(posG, roadway, v)
end

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


text:
module TestSolid

import Modia3D
using  Modia3D.Test
using  Modia3D.LinearAlgebra
using StaticArrays
EYE3() = @SMatrix[1.0  0.0  0.0;
                0.0   1.0    0.0 ;
                0.0   0.0    1.0]

# Test Sphere
r = 2.0
m = 3.0
shape   = Modia3D.Sphere(diameter=2r)
mass  = Modia3D.MassProperties(shape,m)
solid = Modia3D.Solid(shape=shape, massProperties=mass)
obj   = Modia3D.Object3D(feature = solid)

rref  = SVector(1.0,2.0,3.0)
Tref  = EYE3()
obj.r_abs = rref
obj.R_abs = Tref
eref  = rref/norm(rref)
sp    = Modia3D.supportPoint(obj,eref)
AABB1 = Modia3D.BoundingBox()
AABB2 = Modia3D.BoundingBox(rref[1]-r, rref[1]+r, rref[2]-r, rref[2]+r, rref[3]-r, rref[3]+r)
Modia3D.boundingBox!(obj, AABB1)

@testset "Modia3D.Solid: Test Sphere" begin
   @test mass.rCM == [0.0, 0.0, 0.0]
   @test isapprox(Modia3D.volume(shape), 4/3*pi*r^3)
   @test isapprox(mass.I             , 2/5*m*r^2*EYE3())
   @test isapprox(sp                 , rref + r*eref)
   @test isapprox(AABB1.x_min        , AABB2.x_min)
   @test isapprox(AABB1.x_max        , AABB2.x_max)
   @test isapprox(AABB1.y_min        , AABB2.y_min)
   @test isapprox(AABB1.y_max        , AABB2.y_max)
   @test isapprox(AABB1.z_min        , AABB2.z_min)
   @test isapprox(AABB1.z_max        , AABB2.z_max)
end

# Test Ellipsoid
a = 1.0
b = 2.0
c = 3.0
m = 4.0
shape  = Modia3D.Ellipsoid(lengthX=2a, lengthY=2b, lengthZ=2c)
mass = Modia3D.MassProperties(shape,m)
solid = Modia3D.Solid(shape=shape, massProperties=mass)
obj = Modia3D.Object3D(feature = solid)

rref  = SVector(1.0,2.0,3.0)
Tref  = EYE3()
eref  = SVector(0.0, 0.0, 1.0)
obj.r_abs = rref
obj.R_abs = Tref
sp    = Modia3D.supportPoint(obj,eref)
AABB1 = Modia3D.BoundingBox()
AABB2 = Modia3D.BoundingBox(rref[1]-a, rref[1]+a, rref[2]-b, rref[2]+b, rref[3]-c, rref[3]+c)
Modia3D.boundingBox!(obj, AABB1)

@testset "Modia3D.Solid: Test Ellipsoid" begin
   @test mass.rCM == [0.0, 0.0, 0.0]
   @test isapprox(Modia3D.volume(shape), 4/3*pi*a*b*c)
   @test isapprox(mass.I             , m/5*Diagonal([b^2+c^2, c^2+a^2, a^2+b^2]))
   @test isapprox(sp                 , rref + c*eref)
   @test isapprox(AABB1.x_min        , AABB2.x_min)
   @test isapprox(AABB1.x_max        , AABB2.x_max)
   @test isapprox(AABB1.y_min        , AABB2.y_min)
   @test isapprox(AABB1.y_max        , AABB2.y_max)
   @test isapprox(AABB1.z_min        , AABB2.z_min)
   @test isapprox(AABB1.z_max        , AABB2.z_max)
end

# Test Box and MassProperties
a = 1.0
b = 2.0
c = 3.0
m = 4.0
d = 2700.0
solidMaterial = Modia3D.SolidMaterial(density=d)
shape   = Modia3D.Box(lengthX=a, lengthY=b, lengthZ=c)
mass1 = Modia3D.MassProperties(shape,m)
mass2 = Modia3D.MassProperties(shape,solidMaterial)
mass3 = Modia3D.MassProperties(shape,"Aluminium")
solid = Modia3D.Solid(shape=shape, massProperties=mass1)
obj = Modia3D.Object3D(feature = solid)

V     = Modia3D.volume(shape)
m2    = mass2.m

rref  = SVector(1.0,2.0,3.0)
Tref  = EYE3()
eref  = SVector(0.0, 0.0, 1.0)
obj.r_abs = rref
obj.R_abs = Tref
sp    = Modia3D.supportPoint(obj,eref)
AABB1 = Modia3D.BoundingBox()
AABB2 = Modia3D.BoundingBox(rref[1]-a/2-solid.collisionSmoothingRadius, rref[1]+a/2+solid.collisionSmoothingRadius, rref[2]-b/2-solid.collisionSmoothingRadius, rref[2]+b/2+solid.collisionSmoothingRadius, rref[3]-c/2-solid.collisionSmoothingRadius, rref[3]+c/2+solid.collisionSmoothingRadius)
Modia3D.boundingBox!(obj, AABB1)

@testset "Modia3D.Solid: Test Box and MassProperties" begin
   @test mass1.rCM == [0.0, 0.0, 0.0]
   @test isapprox(V          , a*b*c)
   @test isapprox(mass1.I    , m/12*Diagonal([b^2+c^2, c^2+a^2, a^2+b^2]))
   @test isapprox(mass2.m    , d*V)
   @test isapprox(mass2.I    , m2/12*Diagonal([b^2+c^2, c^2+a^2, a^2+b^2]))
   @test isapprox(mass3.m    , mass2.m)
   @test isapprox(mass3.I    , mass2.I)
  # @test isapprox(AABB1.x_min, AABB2.x_min)
  # @test isapprox(AABB1.x_max, AABB2.x_max)
  # @test isapprox(AABB1.y_min, AABB2.y_min)
  # @test isapprox(AABB1.y_max, AABB2.y_max)
  # @test isapprox(AABB1.z_min, AABB2.z_min)
  # @test isapprox(AABB1.z_max, AABB2.z_max)
end

# Test MassProperties
m   = 2.0
rCM = [1,2,3]
II  = [5 1 2;
       1 6 3;
       2 3 7]

mass1 = Modia3D.MassProperties(m,rCM,II)
mass2 = Modia3D.MassProperties(mass=m,centerOfMass=rCM,Ixx=5,Iyy=6,Izz=7,Ixy=1,Ixz=2,Iyz=3)

@testset "Modia3D.Solid: Test MassProperties" begin
   @test isapprox(mass1.m, m)
   @test isapprox(mass2.m, m)

   @test mass1.rCM == rCM
   @test mass2.rCM == rCM

   @test mass1.I == II
   @test mass2.I == II
end

end