text:
abstract type Widget{T} end

gtype(::Type{Widget}) = GtkWidget
gtype(::Type{<: Widget{T}}) where {T} = T

# --------------------------------------------------------------------------------
# ----------------------------------| Containers |--------------------------------
# --------------------------------------------------------------------------------

abstract type Container{T} <: Widget{T} end

mutable struct NullContainer <: Container{Gtk.GtkNullContainer}
    widget::Gtk.GtkWidget
    NullContainer() = new(Gtk.GtkNullContainer())
end

mutable struct Window <: Container{GtkWindow}
    widget::GtkWidget
end

function Window(title::String, width::Int, height::Int, children::Widget...; props...)
    widget = @container GtkWindow(title, width, height)
    return @gcpreserve Window(widget)
end

"""
A container used to pack widgets together.
"""
mutable struct Box <: Container{GtkBox}
    widget::GtkWidget
end
"""
        Box(direction, children...; props...) -> Box

Creates a [`Box`](@ref) widget.

The `direction` parameter determines the direction the widgets
are packed. `:v` for vertical or `:h` for horizontal.
"""
function Box(direction::Symbol, children::Widget...; props...)
    widget = @container GtkBox(direction)
    return @gcpreserve Box(widget)
end
"""
A widget with two adjustable panes.
"""
mutable struct Paned <: Container{GtkPaned}
    widget::GtkWidget
end
"""
        Paned(direction, children...; props...) -> Paned

Creates a [`Paned`](@ref) widget.

The `direction` parameter determines how the two panes are
arranged, either horizontally `:h`, or vertically `:v`.
"""
function Paned(direction::Symbol, children::Widget...; props...)
    widget = @container GtkPaned(direction)
    return @gcpreserve Paned(widget)
end
"""
Adds scrollbars to its child widget.
"""
mutable struct Scrolled <: Container{GtkScrolledWindow}
    widget::GtkWidget
end
"""
        Scrolled(children::Widget...; props...) -> Scrolled

Creates a [`Scrolled`](@ref) widget.
"""
function Scrolled(children::Widget...; props...)
    widget = @container GtkScrolledWindow()
    return @gcpreserve Scrolled(widget)
end
"""
A bin with a decorative frame and optional label
"""
mutable struct Frame <: Container{GtkFrame}
    widget::GtkWidget
end

function Frame(children::Widget...; props...)
    widget = @container GtkFrame()
    return @gcpreserve Frame(widget)
end

function Frame(label::String, children::Widget...; props...)
    widget = @container GtkFrame(label)
    return @gcpreserve Frame(widget)
end

function Frame(label::Widget, children::Widget...; props...)
    widget = @container GtkFrame()
    this = @gcpreserve Frame(widget)
    this["label-widget"] = gwidget(label)
    return this
end

mutable struct Notebook <: Container{GtkNotebook}
    widget::GtkWidget
end

function Notebook(children::Widget...; props...)
    widget = @container GtkNotebook()
    return @gcpreserve Notebook(widget)
end

mutable struct TreeView <: Container{GtkTreeView}
    widget::GtkWidget
end

function TreeView(children...; props...)
    widget = @container GtkTreeView(types...)
    return @gcpreserve TreeView(widget)
end

mutable struct Expander <: Container{GtkExpander}
    widget::GtkWidget
end

function Expander(label::AbstractString, children::Widget...; props...)
    widget = @container GtkExpander(label)
    return @gcpreserve Expander(widget)
end

# ----------------------------------------------------------------------------
# ----------------------------------| Inputs |--------------------------------
# ----------------------------------------------------------------------------

abstract type Input{T, S} <: Widget{S} end
abstract type RangeWidget{T, S} <: Input{T, S} end

"""
A slider widget for selecting a value from a range
"""
mutable struct Slider{T} <: RangeWidget{T, GtkScale}
    widget::GtkWidget
    value::Observable{T}
    adjustment::GtkAdjustment
end

function Slider(range::AbstractRange{T}, vertical::Bool = false; start::Real = middle(range), props...) where {T}
    widget = @widget GtkScale(vertical, range)
    adj = G.adjustment(widget)
    this = @gcpreserve Slider{T}(widget, Observable{T}(start), adj)
    G.value(adj, start)
    signal_connect(adj, "value-changed") do args...
        this.value[] = G.value(adj)
    end
    return this
end
"""
Retrieve an integer or floating-point number from the user
"""
mutable struct SpinButton{T} <: RangeWidget{T, GtkSpinButton}
    widget::GtkWidget
    value::Observable{T}
    adjustment::GtkAdjustment
end

function SpinButton(range::AbstractRange{T}; start::Real = middle(range), props...) where {T}
    widget = @widget GtkSpinButton(range)
    adj = G.adjustment(widget)
    this = @gcpreserve SpinButton{T}(widget, Observable{T}(start), adj)
    G.value(adj, start)
    signal_connect(adj, "value-changed") do args...
        this.value[] = G.value(adj)
    end
    return this
end
"""
A single line text entry field.
"""
mutable struct TextField <: Input{String, GtkEntry}
    widget::GtkWidget
    value::Observable{String}
end

function TextField(; props...)
    widget = @widget GtkEntry()
    text = get(props, :text, "")
    this = @gcpreserve TextField(widget, Observable(text))
    onevent("changed", this) do args...
        this.value[] = Gtk.bytestring(G.text(widget))
    end
    return this
end

"""
A button to launch a color selection dialog
"""
mutable struct ColorButton{T <: Colorant} <: Input{T, GtkColorButton}
    widget::GtkWidget
    value::Observable{T}
end

function ColorButton(color::T; props...) where {T <: Colorant}
    widget = @widget GtkColorButton(convert(GdkRGBA, color))
    this = @gcpreserve ColorButton{T}(widget, Observable(color))
    onevent("notify::color", widget) do args...
        this.value[] = this["rgba", GdkRGBA]
    end
    return this
end

ColorButton(color::String; props...) = ColorButton(parse(Colorant, color); props...)

abstract type Activable{W} <: Input{Bool, W} end

"""
Create buttons which retain their state
"""
mutable struct ToggleButton <: Activable{GtkToggleButton}
    widget::GtkWidget
    value::Observable{Bool}
end

function ToggleButton(text::String; props...)
    widget = @widget GtkToggleButton(text)
    this = @gcpreserve ToggleButton(widget, Observable(false))
    onevent("toggled", widget) do args...
        this.value[] = G.active(widget)
    end
    return this
end

"""
Create widgets with a discrete toggle button
"""
mutable struct CheckBox <: Activable{GtkCheckButton}
    widget::GtkWidget
    value::Observable{Bool}
end

function CheckBox(checked::Bool; props...)
    widget = @widget GtkCheckButton()
    this = @gcpreserve CheckBox(widget, Observable(checked))
    G.active(widget, checked)
    onevent("toggled", widget) do args...
        this.value[] = G.active(widget)
    end
    return this
end

"""
A “light switch” style toggle
"""
mutable struct Switch <: Activable{GtkSwitch}
    widget::GtkWidget
    value::Observable{Bool}
end

function Switch(checked::Bool; props...)
    widget = @widget GtkSwitch()
    this = @gcpreserve Switch(widget, Observable(checked))
    this["active"] = true
    onevent("notify::active", this) do args...
        setindex!(this.value, this["active", Bool])
    end
    return this
end

mutable struct Dropdown <: Input{String, GtkComboBoxText}
    widget::GtkWidget
    choices::Tuple{Vararg{String}}
    index::Observable{Int}
    value::Observable{String}
end

"""
A widget used to choose from a list of items
"""
function Dropdown(choices::String...; active::Int = 1, props...)
    widget = @widget GtkComboBoxText()
    index = Observable(active)
    for choice in choices
        push!(widget, choice)
    end
    this = @gcpreserve Dropdown(widget, choices, index, Observable(choices[active]))
    this["active"] = active - 1
    onevent("changed", this) do args...
        this.index[] = getprop(widget, :active, Int)
        this.value[] = Gtk.bytestring( G.active_text(widget) )
    end
    return this
end

# ------------------------------------------------------------------------------------
# ----------------------------------| Common Widgets |--------------------------------
# ------------------------------------------------------------------------------------
"""
A widget that emits a signal when clicked on
"""
mutable struct Button <: Widget{GtkButton}
    widget::GtkWidget
end

Button(text::String; props...) = Button(@widget GtkButton(text))

function Button(icon::AbstractString, size::Symbol; props...)
    gsize::Cuint = getfield(Gtk.GtkIconSize, Symbol(uppercase(string(size))))
    widget = @widget icon_button(icon, gsize)
    return @gcpreserve Button(widget)
end

function Button(clicked::Function, text::String; props...)
    button = Button(text; props...)
    onevent(clicked, "clicked", button)
    return button
end

function Button(clicked::Function, icon::AbstractString, size::Symbol; props...)
    button = Button(icon, size; props...)
    onevent(clicked, "clicked", button)
    return button
end

"""
A widget that displays a small to medium amount of text
"""
mutable struct Label <: Widget{GtkLabel}
    widget::GtkWidget
    function Label(text::String; props...)
        widget = @widget GtkLabel(text)
        return @gcpreserve new(widget)
    end
end

"""
A widget which indicates progress visually

Use the [`fill!`](@ref) and the [`pulse!`](@ref) functions
to control the progress bar.
"""
mutable struct ProgressBar <: Widget{GtkProgressBar}
    widget::GtkWidget
end

function ProgressBar(; props...)
    widget = @widget GtkProgressBar()
    return @gcpreserve ProgressBar(widget)
end

function ProgressBar(fill::Float64; props...)
    this = ProgressBar(; props...)
    fill!(this, fill)
    return this
end

"""
A widget for displaying an image.
"""
mutable struct ImageView <: Widget{GtkImage}
    widget::GtkWidget
    function ImageView(; props...)
        widget = @widget GtkImage()
        return @gcpreserve new(widget)
    end
    function ImageView(file::AbstractString; props...)
        widget = @widget GtkImage(file)
        return @gcpreserve new(widget)
    end
end

"""
Custom drawing with OpenGL.
"""
mutable struct GLArea <: Widget{GtkGLArea}
    widget::GtkWidget
    function GLArea(; props...)
        widget = @widget GtkGLArea()
        return @gcpreserve new(widget)
    end
end


text:
# ===========================================================================
# Redirect info, warnings, and errors to the logger
out = open("runtests.log", "a")
logger = SimpleLogger(out)
global_logger(logger)
@info("stdout redirect and logging")

# Set some keyword defaults
SeisIO.KW.comp = 0x00
has_restricted = safe_isdir(path * "/SampleFiles/Restricted/")
keep_log = false
keep_samples = true

ext: jl


text:
module StdLib

using ..CxxWrap

abstract type CppBasicString <: AbstractString end

# These are defined in C++, but the functions need to exist to add methods
function append end
function cppsize end
function cxxgetindex end
function cxxsetindex! end
function push_back end
function resize end

@wrapmodule(CxxWrap.libcxxwrap_julia_stl)

function __init__()
  @initcxx
end

# Pass-through for fundamental types
_append_dispatch(v::StdVector,a::Vector,::Type{CxxWrap.IsNormalType}) = append(v,a)
# For C++ types, convert the array to an array of references, so the pointers can be read directly from a contiguous array on the C++ side
_append_dispatch(v::StdVector{T}, a::Vector{<:T},::Type{CxxWrap.IsCxxType}) where {T} = append(v,CxxWrap.CxxRef.(a))
# Choose the correct append method depending on the type trait
append(v::StdVector{T}, a::Vector{<:T}) where {T} = _append_dispatch(v,a,CxxWrap.cpp_trait_type(T))

Base.ncodeunits(s::CppBasicString)::Int = cppsize(s)
Base.codeunit(s::StdString) = UInt8
Base.codeunit(s::StdWString) = Cwchar_t == Int32 ? UInt32 : UInt16
Base.codeunit(s::CppBasicString, i::Integer) = s[i]
Base.isvalid(s::CppBasicString, i::Integer) = (0 < i <= ncodeunits(s))
function Base.iterate(s::CppBasicString, i::Integer=1)
  if !isvalid(s,i)
    return nothing
  end
  return(convert(Char,codeunit(s,i)),i+1)
end

function StdWString(s::String)
  char_arr = transcode(Cwchar_t, s)
  StdWString(char_arr, length(char_arr))
end

function StdVector(v::Vector{T}) where {T}
  if (CxxWrap.cpp_trait_type(T) == CxxWrap.IsCxxType)
    return StdVector(CxxRef.(v))
  end
  result = StdVector{T}()
  append(result, v)
  return result
end

function StdVector(v::Vector{CxxRef{T}}) where {T}
  result = isconcretetype(T) ? StdVector{supertype(T)}() : StdVector{T}()
  append(result, v)
  return result
end

function StdVector(v::Vector{Bool})
  result = StdVector{CxxBool}()
  append(result, convert(Vector{CxxBool}, v))
  return result
end

Base.IndexStyle(::Type{<:StdVector}) = IndexLinear()
Base.size(v::StdVector) = (Int(cppsize(v)),)
Base.getindex(v::StdVector, i::Int) = cxxgetindex(v,i)[]
Base.setindex!(v::StdVector{T}, val, i::Int) where {T} = cxxsetindex!(v, convert(T,val), i)

function Base.push!(v::StdVector, x)
  push_back(v, x)
  return v
end

function Base.resize!(v::StdVector, n::Integer)
  resize(v, n)
  return v
end

Base.empty!(v::StdVector) = Base.resize!(v, 0)

function Base.append!(v::StdVector, a::Vector)
  append(v, a)
  return v
end

function Base.append!(v::StdVector{CxxBool}, a::Vector{Bool})
  append(v, convert(Vector{CxxBool}, a))
  return v
end

# Make sure functions taking a C++ string as argument can also take a Julia string
CxxWrap.map_julia_arg_type(x::Type{<:StdString}) = AbstractString
Base.convert(::Type{T}, x::String) where {T<:StdString} = StdString(x)
Base.cconvert(::Type{CxxWrap.ConstCxxRef{StdString}}, x::String) = StdString(x)
Base.unsafe_convert(::Type{CxxWrap.ConstCxxRef{StdString}}, x::StdString) = ConstCxxRef(x)

function StdValArray(v::Vector{T}) where {T}
  return StdValArray{T}(v, length(v))
end

Base.IndexStyle(::Type{<:StdValArray}) = IndexLinear()
Base.size(v::StdValArray) = (Int(cppsize(v)),)
Base.getindex(v::StdValArray, i::Int) = cxxgetindex(v,i)[]
Base.setindex!(v::StdValArray{T}, val, i::Int) where {T} = cxxsetindex!(v, convert(T,val), i)

end


text:
"""
    equispaced_line_mesh(
    x1min::FT,
    x1max::FT,
    n1,
    per = false,
) where {FT<:AbstractFloat}

This function generates an equispaced line mesh with user provided coordinates for left and right
domain vertices, number of mesh elements and periodicity information.
"""
function equispaced_line_mesh(
    x1min::FT,
    x1max::FT,
    n1,
    per = false,
) where {FT<:AbstractFloat}
    x1c = range(x1min, x1max; length = n1 + 1)
    return Mesh1D(x1c, per)
end

function Mesh1D(x1c, per = false)
    I = Int
    FT = eltype(x1c)
    nverts = length(x1c)
    nelems = nverts - 1
    nbndry = per ? 1 : 3
    coordinates = Array{FT}(undef, nverts, 1)
    coordinates[:, 1] .= x1c
    vertex_neighbors = zeros(I, nverts, 4)
    tag_names = (:interior, :left, :right)
    if per
        unique_verts = Array{I}(1:nverts-1)
        boundary_tag_names = tuple(tag_names[1])
        boundary_tags = [0]
        vertex_neighbors[1, :] .= [nelems, 2, 1, 1]
        vertex_neighbors[nverts, :] .= [nelems, 2, 1, 1]
        vertex_boundary = Vector(1:nverts-1)
        vertex_boundary_offset = [1, nverts]
    else
        unique_verts = Array{I}(1:nverts)
        vertex_neighbors[1, :] .= [0, 0, 1, 1]
        vertex_neighbors[nverts, :] .= [nelems, 2, 0, 0]
        if nelems > 1
            boundary_tag_names = tuple(tag_names...)
            boundary_tags = [0, 1, 2]
            vertex_boundary = vcat(Vector(2:nverts-1), 1, nverts)
            vertex_boundary_offset = [1, nverts - 1, nverts, nverts + 1]
        else
            boundary_tag_names = tuple(tag_names[2:end]...)
            boundary_tags = [1, 2]
            vertex_boundary = vcat(1, nverts)
            vertex_boundary_offset = [1, nverts, nverts + 1]
        end
    end
    for vt = 2:nverts-1
        vertex_neighbors[vt, :] .= [vt - 1, 2, vt, 1]
    end
    elem_verts = Array{I}(undef, nelems, 2)
    for el = 1:nelems
        elem_verts[el, :] .= [el, el + 1]
    end
    return Mesh1D(
        nverts,
        nelems,
        nbndry,
        coordinates,
        unique_verts,
        vertex_neighbors,
        vertex_boundary,
        boundary_tags,
        boundary_tag_names,
        vertex_boundary_offset,
        elem_verts,
    )
end

max_stars_repo_name: CliMA/GCMMesh

max_issues_repo_name: CliMA/GCMMesh

max_forks_repo_name: CliMA/GCMMesh


text:
# 2D points in some minority class.
p1_x = 0.0
p1_y = 2.0
p2_x = 1.0
p2_y = 3.0

data = [p1_x p2_x;
        p1_y p2_y]

@test Resample._detect_k(data) == 1

n = 2
new_points = smote(data, n)
@test Resample._npoints(new_points) == n

for p in eachcol(new_points)
    @test Resample._is_in_between(p, [p1_x, p1_y], [p2_x, p2_y]) == true
end

n = 100
ndims = 30
data = rand(ndims, 10_000)
new_points = smote(data, n)
@test Resample._npoints(new_points) == n
@test Resample._ndims(new_points) == ndims

n = 4
ndims = 2
data = (; x=[p1_x, p2_x], y=[p1_y, p2_y])
new_points = smote(data, n)
@test Resample._npoints(new_points) == n
@test Resample._ndims(new_points) == ndims

# 10D
data = rand(10, 2)
n = 6
new_points = smote(data, n)
for p in eachcol(new_points)
    p1, p2 = Resample._point.(Ref(data), [1, 2])
    @test Resample._is_in_between(p, p1, p2)
end

@testset "integers" begin
    data = (; X=1:2, Y=2:3)
    n = 3
    @test length(smote(data, n).X) == n
end

@testset "auto balance" begin
    data = (; X=[1, 2, 1, 2, 1], class=[1, 2, 1, 2, 1])
    if v"1.7" ≤ VERSION
        @inferred Resample._col2int(data, :class)
    end

    out = smote(data, :class)
    @test out.X == [2.0]
    @test out.class == [2.0]
end

@testset "dataframes" begin
    data = DataFrame(; X=[1, 2, 1, 2, 1], class=[1, 2, 1, 2, 1])
    out = DataFrame(smote(data, :class))
    combined = vcat(data, out)
    @test combined.class == [1, 2, 1, 2, 1, 2]
end

@testset "rng" begin
    rng() = StableRNG(3)
    data = (; X=rand(5), Y=[1, 1, 1, 2, 2])
    @test smote(rng(), data, :Y).X == smote(rng(), data, :Y).X
end


text:
##### NewtonsMethod

export NewtonsMethod

"""
    NewtonsMethod(
            f!::F!,
            j!::J!,
            x_init::A,
        ) where {F! <: Function, J! <: Function, A <: AbstractArray}

A non-linear system of equations type.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct NewtonsMethod{FT, F!, J!, A, JA} <: AbstractNonlinearSolverMethod{FT}
    "Function to find the root of"
    f!::F!
    "Jacobian of `f!`"
    j!::J!
    "Initial guess"
    x_init::A
    "Storage"
    x1::A
    "Storage"
    J::JA
    "Storage"
    J⁻¹::JA
    "Storage"
    F::A
    function NewtonsMethod(
        f!::F!,
        j!::J!,
        x_init::A,
    ) where {F! <: Function, J! <: Function, A <: AbstractArray}
        x1 = similar(x_init)
        J = similar(x_init, (length(x_init), length(x_init)))
        J⁻¹ = similar(J)
        F = similar(x_init)
        JA = typeof(J)
        FT = eltype(x_init)
        return new{FT, F!, J!, A, JA}(f!, j!, x_init, x1, J, J⁻¹, F)
    end
end

method_args(m::NewtonsMethod) = (m.x_init, m.x1, m.f!, m.F, m.j!, m.J, m.J⁻¹)
function solve!(
    ::NewtonsMethod,
    x0::AT,
    x1::AT,
    f!::F!,
    F::FA,
    j!::J!,
    J::JA,
    J⁻¹::J⁻¹A,
    soltype::SolutionType,
    tol::AbstractTolerance{FT},
    maxiters::Int,
) where {FA, J⁻¹A, JA, F! <: Function, J! <: Function, AT, FT}

    x_history = init_history(soltype, AT)
    F_history = init_history(soltype, AT)
    if soltype isa VerboseSolution
        f!(F, x0)
        j!(J, x0)
        push_history!(x_history, x0, soltype)
        push_history!(F_history, F, soltype)
    end
    for i in 1:maxiters
        f!(F, x0)
        j!(J, x0)
        x1 .= x0 .- J \ F
        push_history!(x_history, x1, soltype)
        push_history!(F_history, F, soltype)
        if tol(x0, x1, F)
            return SolutionResults(
                soltype,
                x1,
                true,
                F,
                i,
                x_history,
                F_history,
            )
        end
        x0 = x1
    end
    return SolutionResults(
        soltype,
        x0,
        false,
        F,
        maxiters,
        x_history,
        F_history,
    )
end


text:
# ] add JSON, HTTP;

include("src/Pluto.jl")
@async Pluto.run(1234);


text:
import AlphaZero.GI
using Random

# Initialize the dimension constants
const NUM_SUITS = 4
const NUM_VALUES= 13
const NUM_ARRAYS = 18
const NUM_CARDS = 52

# Initialize the player constants
const Player = UInt8
const WHITE = 0x01
const BLACK = 0x02
switch_players(p::Player) = 0x03 - p

# Initialize the board type
const Cell = UInt8
const EMPTY = 0x00
const Board = Array{Cell, 3}

# Initialize the deck
deck = zeros(UInt8, NUM_VALUES, NUM_SUITS, NUM_CARDS)
for i in 1:NUM_CARDS
  deck[:,:,i] = setindex!(deck[:,:,i],1,i)
end

#####
##### Define the Game
#####

mutable struct Game <: GI.AbstractGame
  board :: Board
  curplayer :: Player
  finished :: Bool
  trick_winner :: Player
  amask :: Vector{Bool} # actions mask
  # Actions history, which uniquely identifies the current board position
  # Used by external solvers
  history :: Union{Nothing, Vector{Int}}
end

function Game()
  board = random_starting_state(deck)
  curplayer = WHITE
  finished = false
  trick_winner = 0x00
  amask = get_legal_actions(board)
  history = Int[]
  Game(board, curplayer, finished, trick_winner, amask, history)
end

function Game(state)
  board = state.board
  curplayer = state.curplayer
  trick_winner = state.trick_winner
  amask = get_legal_actions(board)
  finished = !any(amask)
  history = Int[]
  Game(board, curplayer, finished, trick_winner, amask, history)
end

function GI.play!(g::Game, action)
  g.board = copy(g.board)
  isnothing(g.history) || push!(g.history, action)
  g.board[:,:,1] = setindex!(g.board[:,:,1], 0x00, action)
  g.board[:,:,6] = setindex!(g.board[:,:,6], maximum(g.board[:,:,6]) + 1, action)
  g.board[:,:,1:4] = circshift(g.board[:,:,1:4], (0,0,-1))
  update_status!(g)
  g.curplayer = switch_players(g.curplayer)
end

function GI.game_terminated(g::Game)
  return g.finished
end

GI.State(::Type{Game}) = typeof((board=zeros(UInt8, NUM_VALUES, NUM_SUITS, NUM_ARRAYS), curplayer=WHITE, trick_winner=WHITE))
GI.Action(::Type{Game}) = Int
GI.two_players(::Type{Game}) = true
const ACTIONS = collect(1:NUM_CARDS)
GI.actions(::Type{Game}) = ACTIONS
history(g::Game) = g.history
GI.actions_mask(g::Game) = g.amask
GI.current_state(g::Game) = (board=g.board, curplayer=g.curplayer, trick_winner=g.trick_winner)
GI.white_playing(::Type{Game}, state) = state.curplayer == WHITE

function GI.white_reward(g::Game)
  g.trick_winner == WHITE && (return  1.)
  g.trick_winner == BLACK && (return -1.)
  return 0.
end

#####
##### Helper Functions
#####

function random_starting_state(deck::Array{UInt8,3})
  # Initialize the hand
  board = zeros(UInt8, NUM_VALUES, NUM_SUITS, NUM_ARRAYS)

  # Shuffle the Deck
  deck = deck[:,:,shuffle(1:end)]

  # Create the hands
  for i in 1:4
    board[:,:,i] = sum(deck[:,:,(13*(i-1))+1:13*i], dims=3)
  end

  # Pick one of the five possible trump suits
  trump_suit = rand(1:5)
  if trump_suit == 5
    board[:,:,5] = ones(UInt8, NUM_VALUES, NUM_SUITS)
  else
    board[:,trump_suit,5] = ones(UInt8, NUM_VALUES)
  end

  return board
end

function get_legal_actions(hand::Board)
  return vec(Array{Bool}(hand[:,:,1]))
end

function update_status!(g::Game)
  g.amask = get_legal_actions(g.board)
  g.finished = !any(g.amask)
  if maximum(g.board[:,:,6]) == 4
    g.trick_winner = calculate_winner(g.board)
    if !g.finished
      g.board[:,:,19-sum(g.board[:,:,1])] = g.board[:,:,6]
      g.board[:,:,6] = zeros(UInt8, NUM_VALUES, NUM_SUITS)
    end
  else
    g.trick_winner = 0x00
  end
end

function calculate_winner(board::Board)
  # Analyze the leading card
  leading_card = findfirst(x -> x == 1, board[:,:,6])
  winning_value = leading_card[1]
  winning_suit = leading_card[2]
  winning_player = 0x01
  if isnothing(findfirst(iszero, board[:,:,5]))
    trump_suit = 0x05
  else
    trump_suit = findfirst(x -> x == 1, board[:,:,5])[2]
  end

  # Loop through the other cards to calculate the winner
  for i in 2:4
    next_card = findfirst(x -> x == i, board[:,:,6])
    next_value = next_card[1]
    next_suit = next_card[2]
    if next_suit == trump_suit && winning_suit != trump_suit
      winning_suit = next_suit
      winning_value = next_value
      winning_player = i
    elseif next_suit == winning_suit && next_value > winning_value
      winning_value = next_value
      winning_player = i
    end
  end

  # Return the reward corresponding to the correct winner
  if winning_player == 1 || winning_player == 3
    return 0x01
  else
    return 0x02
  end
end

function Base.copy(g::Game)
  history = isnothing(g.history) ? nothing : copy(g.history)
  Game(g.board, g.curplayer, g.finished, g.trick_winner, copy(g.amask), history)
end

function Base.copy(state)
  return (board=copy(state.board), curplayer=state.curplayer, trick_winner=state.trick_winner)
end

#####
##### Simple heuristic for minmax
#####

function GI.heuristic_value(g::Game)
  # Find the trump suit
  if isnothing(findfirst(iszero, g.board[:,:,5]))
    trump_suit = 0x05
  else
    trump_suit = findfirst(x -> x == 1, g.board[:,:,5])[2]
  end

  # Add all of the card values (as higher values are better)
  mine = sum([card[1] for card in findall(x -> x == 0x01, g.board[:,:,1])]) + sum([card[1] for card in findall(x -> x == 0x01, g.board[:,:,3])])
  yours = sum([card[1] for card in findall(x -> x == 0x01, g.board[:,:,2])]) + sum([card[1] for card in findall(x -> x == 0x01, g.board[:,:,4])])

  # If the trump suit is a single suit, count those values again
  if trump_suit != 0x05
    mine += sum([card[1] for card in findall(x -> x == 0x01, g.board[:,trump_suit,1])]) + sum([card[1] for card in findall(x -> x == 0x01, g.board[:,trump_suit,3])])
    yours += sum([card[1] for card in findall(x -> x == 0x01, g.board[:,trump_suit,2])]) + sum([card[1] for card in findall(x -> x == 0x01, g.board[:,trump_suit,4])])
  end

  hvalue = mine - yours
  return convert(Float64, hvalue)
end

#####
##### ML interface
#####

function GI.vectorize_state(::Type{Game}, state)
  return convert(Array{Float32, 3}, state.board)
end

#####
##### Symmetries
#####

function flipped_board(board)
  return board[:,reverse(1:NUM_SUITS),:]
end

function GI.symmetries(::Type{Game}, state)
  symb = flipped_board(state.board)
  σ = vec(reshape(collect(1:NUM_CARDS), (NUM_VALUES, NUM_SUITS))[:,reverse(1:NUM_SUITS)])
  syms = (board=Board(symb), curplayer=state.curplayer, trick_winner=state.trick_winner)
  return [(syms, σ)]
end

#####
##### User interface
#####

GI.action_string(::Type{Game}, a) = string(a)

function GI.parse_action(g::Game, str)
  try
    p = parse(Int, str)
    1 <= p <= NUM_CARDS ? p : nothing
  catch
    nothing
  end
end

function GI.render(g::Game)
  if isnothing(findfirst(iszero, g.board[:,:,5]))
    trump_suit = 5
  else
    trump_suit = findfirst(x -> x == 1, g.board[:,:,5])[2]
  end
  println("Trump:" * string(trump_suit))
  println("Previous Reward: " * string(GI.white_reward(g)))
  println("Current Trick: " * string(findfirst(x -> x == 1, g.board[:,:,6])) * ", " * string(findfirst(x -> x == 2, g.board[:,:,6])) * ", " * string(findfirst(x -> x == 3, g.board[:,:,6])) * ", " * string(findfirst(x -> x == 4, g.board[:,:,6])))
  println("Hand 1:" * string(findall(x -> x == 0x01, g.board[:,:,1])))
  println("Hand 2:" * string(findall(x -> x == 0x01, g.board[:,:,2])))
  println("Hand 3:" * string(findall(x -> x == 0x01, g.board[:,:,3])))
  println("Hand 4:" * string(findall(x -> x == 0x01, g.board[:,:,4])))
end

function GI.read_state(::Type{Game})
  return nothing
end

max_stars_repo_name: btheodorou/Bridge-NN-Research

max_issues_repo_name: btheodorou/Bridge-NN-Research

max_forks_repo_name: btheodorou/Bridge-NN-Research


text:
using MultivariateStats
using Base.Test

srand(34568)

## PCA with zero mean

X = randn(5, 10)
Y = randn(3, 10)

P = qr(randn(5, 5))[1][:, 1:3]
pvars = [5., 4., 3.]
M = PCA(Float64[], P, pvars, 15.0)

@test indim(M) == 5
@test outdim(M) == 3
@test mean(M) == zeros(5)
@test projection(M) == P
@test principalvars(M) == pvars
@test principalvar(M, 2) == pvars[2]
@test tvar(M) == 15.0
@test tprincipalvar(M) == 12.0
@test tresidualvar(M) == 3.0
@test principalratio(M) == 0.8

@test transform(M, X[:,1]) ≈ P'X[:,1]
@test transform(M, X) ≈ P'X

@test reconstruct(M, Y[:,1]) ≈ P * Y[:,1]
@test reconstruct(M, Y) ≈ P * Y

## PCA with non-zero mean

mv = rand(5)
M = PCA(mv, P, pvars, 15.0)

@test indim(M) == 5
@test outdim(M) == 3
@test mean(M) == mv
@test projection(M) == P
@test principalvars(M) == pvars
@test principalvar(M, 2) == pvars[2]
@test tvar(M) == 15.0
@test tprincipalvar(M) == 12.0
@test tresidualvar(M) == 3.0
@test principalratio(M) == 0.8

@test transform(M, X[:,1]) ≈ P' * (X[:,1] .- mv)
@test transform(M, X) ≈ P' * (X .- mv)

@test reconstruct(M, Y[:,1]) ≈ P * Y[:,1] .+ mv
@test reconstruct(M, Y) ≈ P * Y .+ mv

## prepare training data

d = 5
n = 1000

R = qr(randn(d, d))[1]
@test R'R ≈ eye(5)
scale!(R, sqrt.([0.5, 0.3, 0.1, 0.05, 0.05]))

X = R'randn(5, n) .+ randn(5)
mv = vec(mean(X, 2))
Z = X .- mv

C = cov(X, 2)
pvs0 = sort(eigvals(C); rev=true)
tv = sum(diag(C))

## pcacov (default using cov when d < n)

M = fit(PCA, X)
P = projection(M)
pvs = principalvars(M)

@test indim(M) == 5
@test outdim(M) == 5
@test mean(M) == mv
@test P'P ≈ eye(5)
@test C*P ≈ P*Diagonal(pvs)
@test issorted(pvs; rev=true)
@test pvs ≈ pvs0
@test tvar(M) ≈ tv
@test sum(pvs) ≈ tvar(M)
@test reconstruct(M, transform(M, X)) ≈ X

M = fit(PCA, X; mean=mv)
@test projection(M) ≈ P

M = fit(PCA, Z; mean=0)
@test projection(M) ≈ P

M = fit(PCA, X; maxoutdim=3)
P = projection(M)

@test indim(M) == 5
@test outdim(M) == 3
@test P'P ≈ eye(3)
@test issorted(pvs; rev=true)

M = fit(PCA, X; pratio=0.85)

@test indim(M) == 5
@test outdim(M) == 3
@test P'P ≈ eye(3)
@test issorted(pvs; rev=true)

## pcastd

M = fit(PCA, X; method=:svd)
P = projection(M)
pvs = principalvars(M)

@test indim(M) == 5
@test outdim(M) == 5
@test mean(M) == mv
@test P'P ≈ eye(5)
@test isapprox(C*P, P*Diagonal(pvs), atol=1.0e-3)
@test issorted(pvs; rev=true)
@test isapprox(pvs, pvs0, atol=1.0e-3)
@test isapprox(tvar(M), tv, atol=1.0e-3)
@test sum(pvs) ≈ tvar(M)
@test reconstruct(M, transform(M, X)) ≈ X

M = fit(PCA, X; method=:svd, mean=mv)
@test projection(M) ≈ P

M = fit(PCA, Z; method=:svd, mean=0)
@test projection(M) ≈ P

M = fit(PCA, X; method=:svd, maxoutdim=3)
P = projection(M)

@test indim(M) == 5
@test outdim(M) == 3
@test P'P ≈ eye(3)
@test issorted(pvs; rev=true)

M = fit(PCA, X; method=:svd, pratio=0.85)

@test indim(M) == 5
@test outdim(M) == 3
@test P'P ≈ eye(3)
@test issorted(pvs; rev=true)

# test that fit works with Float32 values
X2 = convert(Array{Float32,2}, X)
# Float32 input, default pratio
M = fit(PCA, X2; maxoutdim=3)
# Float32 input, specified Float64 pratio
M = fit(PCA, X2, pratio=0.85)
# Float32 input, specified Float32 pratio
M = fit(PCA, X2, pratio=0.85f0)
# Float64 input, specified Float32 pratio
M = fit(PCA, X, pratio=0.85f0)


text:
include("default_params.jl")
include("model_def.jl")

sim_data = CSV.read((@__DIR__)*"/data/sims.csv",DataFrame)
x_vals=range(0,1.6;length=n_err_evals)

function cdf(x::AbstractVector,data::DataFrame)
    sort!(data,:X)

    x2 = data.X[data.φ.==2]
    x4 = data.X[data.φ.==4]

    c2 = 1
    c4 = 1

    F2 = zeros(length(x))
    F4 = zeros(length(x))

    for (c,xi) in enumerate(x)
        tmp2 = findfirst(y->y>xi,x2[c2:end])
        if !(tmp2===nothing)
            c2 += tmp2-1
        else
            c2 = length(x2)+1
        end
        F2[c] = c2-1

        tmp4 = findfirst(y->y>xi,x4[c4:end])
        if !(tmp4===nothing)
            c4 += tmp4-1
        else
            c4 = length(x4)+1
        end
        F4[c] = c4-1
    end
    return F2./length(data.X), F4./length(data.X)
end

F2,F4 = cdf(x_vals,sim_data)

CSV.write((@__DIR__)*"/data/sims_cdf_evaluated.csv",DataFrame(phase_2=F2,phase_4=F4))

function bootstrap(sim_data,x_vals,n_boot)
    rng = StableRNGs.StableRNG(1412)
    n_sims = length(sim_data.X)
    for n in 1:n_boot
        ind = rand(rng,1:n_sims,n_sims)
        boot_data = sim_data[ind,:]
        F2,F4 = cdf(x_vals,boot_data)
        CSV.write(
            (@__DIR__)*"/data/bootstrap_first_return/sample_"*string(n)*".csv",
            DataFrame(phase_2=F2,phase_4=F4),
        )
    end
    return nothing
end

bootstrap(sim_data,x_vals,n_boot)


text:
"""
    qrdelete!(Q, R, k)

Delete the left-most column of F = Q[:, 1:k] * R[1:k, 1:k] by updating Q and R.
Only Q[:, 1:(k-1)] and R[1:(k-1), 1:(k-1)] are valid on exit.
"""
function qrdelete!(Q::AbstractMatrix, R::AbstractMatrix, k::Int)
  n, m = size(Q)
  m == LinearAlgebra.checksquare(R) || throw(DimensionMismatch())
  1 ≤ k ≤ m || throw(ArgumentError())

  # apply Givens rotations
  for i in 2:k
      g = first(givens(R, i - 1, i, i))
      lmul!(g, R)
      rmul!(Q, g')
  end

  # move columns of R
  @inbounds for j in 1:(k-1)
    for i in 1:(k-1)
      R[i, j] = R[i, j + 1]
    end
  end

  Q, R
end

"""
    qradd!(Q, R, v, k)

Replace the right-most column of F = Q[:, 1:k] * R[1:k, 1:k] with v by updating Q and R.
This implementation modifies vector v as well. Only Q[:, 1:k] and R[1:k, 1:k] are valid on
exit.
"""
function qradd!(Q::AbstractMatrix, R::AbstractMatrix, v::AbstractVector, k::Int)
  n, m = size(Q)
  n == length(v) || throw(DimensionMismatch())
  m == LinearAlgebra.checksquare(R) || throw(DimensionMismatch())
  1 ≤ k ≤ m || throw(ArgumentError())

  @inbounds for i in 1:(k-1)
    q = view(Q, :, i)
    r = dot(q, v)

    R[i, k] = r
    axpy!(-r, q, v)
  end

  @inbounds begin
    d = norm(v)
    R[k, k] = d
    @. Q[:, k] = v / d
  end

  Q, R
end

function qradd!(Q::AbstractMatrix, R::AbstractMatrix, v::Number, k::Int)
  1 == LinearAlgebra.checksquare(Q) == LinearAlgebra.checksquare(R) ||
    throw(DimensionMismatch())
  k == 1 || throw(ArgumentError())

  R[1, 1] = abs(v)
  Q[1, 1] = one(v)

  Q, R
end

isnewton(nlsolver::NLSolver) = isnewton(nlsolver.cache)
isnewton(nlcache::Union{NLNewtonCache,NLNewtonConstantCache}) = true
isnewton(nlcache::AbstractNLSolverCache) = false

set_new_W!(nlsolver::NLSolver, val::Bool) = set_new_W!(nlsolver.cache, val)
set_new_W!(nlcache::NLNewtonCache, val::Bool) = (nlcache.new_W = val; nothing)
set_new_W!(nlcache::AbstractNLSolverCache, val::Bool) = nothing

get_W(nlsolver::NLSolver) = get_W(nlsolver.cache)
get_W(nlcache::Union{NLNewtonCache,NLNewtonConstantCache}) = nlcache.W

set_W!(nlsolver::NLSolver, W) = set_W!(nlsolver.cache, W)
set_W!(nlcache::Union{NLNewtonCache,NLNewtonConstantCache}, W) = (nlcache.W = W; W)

function get_κtol(nlalg::Union{NLAnderson,NLFunctional,NLNewton}, uTolType, reltol)
  κ = nlalg.κ === nothing ? uTolType(1//100) : uTolType(nlalg.κ)
  tol = nlalg.tol === nothing ? uTolType(min(0.03, first(reltol)^(0.5))) : uTolType(nlalg.tol)
  κ * tol
end

DiffEqBase.@def iipnlsolve begin
  # define additional fields of cache of non-linear solver
  z = similar(u); dz = similar(u); tmp = similar(u); b = similar(u)
  k = zero(rate_prototype)

  # adapt options of non-linear solver to current integration problem
  uTolType = real(uBottomEltypeNoUnits)
  κtol = get_κtol(alg.nlsolve, uTolType, reltol)

  # create cache of non-linear solver
  if alg.nlsolve isa NLNewton
    nf = nlsolve_f(f, alg)

    # check if `nf` is linear
    islin = f isa Union{ODEFunction,SplitFunction} && islinear(nf.f)

    if islin
      # get the operator
      J = nf.f
      W = WOperator(f.mass_matrix, dt, J, true)
    else
      if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype !== nothing
        W = WOperator(f, dt, true)
        J = nothing # is J = W.J better?
      else
        J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
        W = similar(J)
      end
    end

    nlcache = NLNewtonCache(true,W)
  elseif alg.nlsolve isa NLFunctional
    z₊ = similar(z)

    nlcache = NLFunctionalCache(z₊)
  elseif alg.nlsolve isa NLAnderson
    z₊ = similar(z)

    max_history = min(alg.nlsolve.max_history, alg.nlsolve.max_iter, length(z))
    Δz₊s = [zero(z) for i in 1:max_history]
    Q = Matrix{uEltypeNoUnits}(undef, length(z), max_history)
    R = Matrix{uEltypeNoUnits}(undef, max_history, max_history)
    γs = Vector{uEltypeNoUnits}(undef, max_history)
    dzold = zero(z)
    z₊old = zero(z)

    nlcache = NLAndersonCache(z₊,dzold,z₊old,Δz₊s,Q,R,γs,alg.nlsolve.aa_start,alg.nlsolve.droptol)
  end

  # create non-linear solver
  nlsolver = NLSolver{true,typeof(z),typeof(k),uTolType,typeof(γ),typeof(c),typeof(nlcache)}(z,dz,tmp,b,k,one(uTolType),κtol,γ,c,alg.nlsolve.max_iter,10000,nlcache)

  # define additional fields of cache
  fsalfirst = zero(rate_prototype)
  if alg.nlsolve isa NLNewton
    if islin
      du1 = rate_prototype
      uf = nothing
      jac_config = nothing
      linsolve = alg.linsolve(Val{:init},nf,u)
    else
      du1 = zero(rate_prototype)
      # if the algorithm specializes on split problems the use `nf`
      uf = DiffEqDiffTools.UJacobianWrapper(nf,t,p)
      jac_config = build_jac_config(alg,nf,uf,du1,uprev,u,tmp,dz)
      linsolve = alg.linsolve(Val{:init},uf,u)
    end
  else
    J = nothing
    W = nothing
    du1 = rate_prototype
    uf = nothing
    jac_config = nothing
    linsolve = nothing
  end
end

DiffEqBase.@def oopnlsolve begin
  # define additional fields of cache of non-linear solver (all aliased)
  z = uprev; dz = z; tmp = z; b = z; k = rate_prototype

  # define tolerances
  uTolType = real(uBottomEltypeNoUnits)
  κtol = get_κtol(alg.nlsolve, uTolType, reltol)

  # create cache of non-linear solver
  if alg.nlsolve isa NLNewton
    nf = nlsolve_f(f, alg)

    # only use `nf` if the algorithm specializes on split eqs
    uf = DiffEqDiffTools.UDerivativeWrapper(nf,t,p)

    islin = f isa Union{ODEFunction,SplitFunction} && islinear(nf.f)
    if islin || DiffEqBase.has_jac(f)
      # get the operator
      J = islin ? nf.f : f.jac(uprev, p, t)
      if !isa(J, DiffEqBase.AbstractDiffEqLinearOperator)
        J = DiffEqArrayOperator(J)
      end
      W = WOperator(f.mass_matrix, dt, J, false)
    else
      # https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/pull/672
      if u isa StaticArray
        # get a "fake" `J`
        J = if u isa AbstractMatrix && size(u, 1) > 1 # `u` is already a matrix
          u
        elseif size(u, 1) == 1 # `u` is a row vector
          vcat(u, u)
        else # `u` is a column vector
          hcat(u, u)
        end
        W = lu(J)
      else
        W = u isa Number ? u : LU{LinearAlgebra.lutype(uEltypeNoUnits)}(Matrix{uEltypeNoUnits}(undef, 0, 0),
                                                                        Vector{LinearAlgebra.BlasInt}(undef, 0),
                                                                        zero(LinearAlgebra.BlasInt))
      end
    end

    nlcache = NLNewtonConstantCache(W)
  elseif alg.nlsolve isa NLFunctional
    uf = nothing

    nlcache = NLFunctionalConstantCache()
  elseif alg.nlsolve isa NLAnderson
    uf = nothing

    max_history = min(alg.nlsolve.max_history, alg.nlsolve.max_iter, length(z))
    Δz₊s = Vector{typeof(z)}(undef, max_history)
    Q = Matrix{uEltypeNoUnits}(undef, length(z), max_history)
    R = Matrix{uEltypeNoUnits}(undef, max_history, max_history)
    γs = Vector{uEltypeNoUnits}(undef, max_history)

    nlcache = NLAndersonConstantCache(Δz₊s,Q,R,γs,alg.nlsolve.aa_start,alg.nlsolve.droptol)
  end

  # create non-linear solver
  nlsolver = NLSolver{false,typeof(z),typeof(k),uTolType,typeof(γ),typeof(c),typeof(nlcache)}(z,dz,tmp,b,k,one(uTolType),κtol,γ,c,alg.nlsolve.max_iter,10000,nlcache)
end


text:
@testset "Hash equality" begin
    text = "this is a text"
    path = "/a/b/c"
    @test hash(DocumentMetadata()) == hash(DocumentMetadata())
    @test hash(StringDocument(text)) == hash(StringDocument(text))
    @test hash(TokenDocument(text)) == hash(TokenDocument(text))
    @test hash(NGramDocument(text)) == hash(NGramDocument(text))
    @test hash(FileDocument(path)) == hash(FileDocument(path))
    @test hash(Corpus([StringDocument(text)])) ==
        hash(Corpus([StringDocument(text)]))
end


text:
struct DateTimeVariate <: AbstractDateTimeVariate
    name::String
    data::AbstractVector{Int64}
end

Base.length(dt::DateTimeVariate) = length(dt.data)

function DateTimeVariate(data::AbstractVector{Int64})
    DateTimeVariate("", data)
end

function DateTimeVariate(name::String, length::Integer, datpath::String)
    data = Vector{Int64}(undef, length)
    open(datpath) do f
        read!(f, data)
    end
    DateTimeVariate(name, data)
end

function slice(dtvariate::DateTimeVariate, fromobs::Integer, toobs::Integer, slicelength::Integer)
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slice(dtvariate.data, fromobs, toobs, slicelength)
end

function Base.map(dtvariate::DateTimeVariate, dataframe::AbstractDataFrame)
    dataframe[getname(dtvariate)]
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
"""
    CosmoParticles.matrix_rotate(vals::AbstractMatrix{<:Number}, rotmat::AbstractMatrix{<:Real})

Returns `vals` rotated by the rotation matrix `rotmat` by regular matrix multiplication.
Works for any dimensions and is optimized for unitful arrays.

This is not exported.
"""
matrix_rotate(vals::AbstractMatrix{<:Number}, rotmat::AbstractMatrix{<:Real}) = rotmat *ᵘ vals

"""
    CosmoParticles.matrix_rotate!(vals::AbstractMatrix{<:Number}, rotmat::AbstractMatrix{<:Real})

In-place version of [`CosmoParticles.matrix_rotate`](@ref), overwriting `vals`.

This is not exported.
"""
function matrix_rotate!(vals::AbstractMatrix{<:Number}, rotmat::AbstractMatrix{<:Real})
    vals .= matrix_rotate(vals, rotmat)
    return vals
end

# Without this special method for lazy arrays, matrix multiplication with unitful arrays leads to bizarre errors
# where only one row or certain diagonal parts of the array or properly rotated.
function matrix_rotate!(vals::ApplyMatrix{<:Number}, rotmat::AbstractMatrix{<:Real})
    matrix_rotate!.(vals.args, (rotmat,))
    return vals
end

"""
    CosmoParticles.matrix_rotate!(dst::AbstractMatrix{<:Number}, vals::AbstractMatrix{<:Number}, rotmat::AbstractMatrix{<:Real})

In-place version of [`CosmoParticles.matrix_rotate`](@ref), overwriting `dst`.

This is not exported.
"""
function matrix_rotate!(dst::AbstractMatrix{<:Number}, vals::AbstractMatrix{<:Number}, rotmat::AbstractMatrix{<:Real})
    uproduct!(dst, rotmat, vals)
end

"""
    rotate(p::AbstractParticles, rotmat::AbstractMatrix{<:Real}, prop)
    rotate(p::AbstractParticles, rotmat::AbstractMatrix{<:Real}, props=[:pos, :vel])

Rotates the specified properties `props` of the particles `p` by the rotation matrix `rotmat`.

Creates a copy of the particles with only new pointers to the rotated properties
(by default position and velocity). This means that mutations to other properties will
affect both the original particles as well as the rotated copied particles.
The properties should be given as a single `Symbol`, or a vector of `Symbol`s. Only
the existing quantities will be rotated (e.g., this method does not throw an error if
velocities are not given for the particles).
"""
function rotate(p::AbstractParticles, rotmat::AbstractMatrix{<:Real}, props=[:pos, :vel])
    p = copy(p)

    # only rotate existing quantities
    Threads.@threads for prop in intersect(keys(p), props) |> collect # collect used for compatibility with threads
        p[prop] = matrix_rotate(p[prop], rotmat)
    end

    return p
end

rotate(p::AbstractParticles, rotmat::AbstractMatrix{<:Real}, prop::Symbol) = rotate(p, rotmat, (prop,))

"""
    rotate(pc::AbstractParticleCollection, rotmat::AbstractMatrix{<:Real}, prop)
    rotate(pc::AbstractParticleCollection, rotmat::AbstractMatrix{<:Real}, props=[:pos, :vel])

Rotates the specified properties `props` of the particles in the collection by the rotation matrix `rotmat`.

Creates a copy of the particle collection with only new pointers to the rotated particles
(by default position and velocity). The properties should be given as a single `Symbol`, or a vector of
`Symbol`s.
"""
function rotate(pc::AbstractParticleCollection, rotmat::AbstractMatrix{<:Real}, props=[:pos, :vel])
    pc = copy(pc)

    Threads.@threads for ptype in keys(pc) |> collect # collect used for compatibility with threads
        pc[ptype] = rotate(pc[ptype], rotmat, props)
    end

    return pc
end

function rotate(pc::AbstractParticleCollection, rotmat::AbstractMatrix{<:Real}, prop::Symbol)
    rotate(pc, rotmat, (prop,))
end

"""
    LinearAlgebra.rotate!(p::AbstractParticles, rotmat::AbstractMatrix{<:Real}, props=[:pos, :vel])
    LinearAlgebra.rotate!(p::AbstractParticles, rotmat::AbstractMatrix{<:Real}, prop)
    LinearAlgebra.rotate!(pc::AbstractParticleCollection, rotmat::AbstractMatrix{<:Real}, props=[:pos, :vel])
    LinearAlgebra.rotate!(pc::AbstractParticleCollection, rotmat::AbstractMatrix{<:Real}, prop)

In-place version of [`rotate`](@ref).

This function is reexported.
"""
function LinearAlgebra.rotate!(p::AbstractParticles, rotmat::AbstractMatrix{<:Real}, props=[:pos, :vel])
    # only rotate existing quantities
    Threads.@threads for prop in intersect(keys(p), props) |> collect # collect used for compatibility with threads
        matrix_rotate!(p[prop], rotmat)
    end

    return p
end

LinearAlgebra.rotate!(p::AbstractParticles, rotmat::AbstractMatrix{<:Real}, prop::Symbol) =
    rotate!(p, rotmat, (prop,))

function LinearAlgebra.rotate!(
    pc::AbstractParticleCollection,
    rotmat::AbstractMatrix{<:Real},
    props=[:pos, :vel],
)
    Threads.@threads for ptype in keys(pc) |> collect # collect used for compatibility with threads
        rotate!(pc[ptype], rotmat, props)
    end

    return pc
end

LinearAlgebra.rotate!(pc::AbstractParticleCollection, rotmat::AbstractMatrix{<:Real}, prop::Symbol) =
    rotate!(pc, rotmat, (prop,))

"""
    translate(p::AbstractParticles, Δx, prop::Symbol=:pos)

Translates the specified property of the particles `p` by `Δx`.

Creates a copy of the particles with only new pointers to the translated property.
Typically, `Δx` will be an `AbstractVector`, e.g., for positions or velocities.
"""
function translate(p::AbstractParticles, Δx, prop::Symbol=:pos)
    p = copy(p)

    if haskey(p, prop)
        p[prop] = p[prop] .+ Δx
    end

    return p
end

"""
    translate(pc::AbstractParticleCollection, Δx::AbstractVector{<:Number}, prop::Symbol=:pos)

Translates the specified property of the particles in the collection `pc` by `Δx`.

Creates a copy of the particle collection with only new pointers to the translated properties for the particles.
Typically, `Δx` will be an `AbstractVector`, e.g., for positions or velocities.
"""
function translate(pc::AbstractParticleCollection, Δx, prop::Symbol=:pos)
    pc = copy(pc)

    Threads.@threads for ptype in keys(pc) |> collect # collect used for compatibility with threads
        pc[ptype] = translate(pc[ptype], Δx, prop)
    end

    return pc
end

"""
    translate!(p::AbstractParticles, Δx, prop::Symbol=:pos)
    translate!(pc::AbstractParticleCollection, Δx, prop::Symbol=:pos)

In-place version of [`translate`](@ref).
"""
function translate!(p::AbstractParticles, Δx, prop::Symbol=:pos)
    if haskey(p, prop)
        p[prop] .+= Δx
    end

    return p
end

function translate!(pc::AbstractParticleCollection, Δx, prop::Symbol=:pos)
    Threads.@threads for ptype in keys(pc) |> collect # collect used for compatibility with threads
        translate!(pc[ptype], Δx, prop)
    end

    return pc
end

"""
    _translate_periodic(x::Number, x₀::Number, period::Number, period_half::Number=1//2*period)

Returns the coordinate `x` shifted by the `period` if reference coordinate `x₀` is on the other side.

For periodic simulation boxes, the period corresponds to the boxlength.

This is not exported.
"""
@inline function _translate_periodic(
    x::Number,
    x₀::Number,
    period::Number,
    period_half::Number=1 // 2 * period,
)
    if x - x₀ > period_half
        return x - period
    elseif x₀ - x > period_half
        return x + period
    end
    return x
end

"""
    translate_periodic(p::AbstractParticles, centerpos, period, prop::Symbol=:pos)
    translate_periodic(pc::AbstractParticleCollection, centerpos, period, prop::Symbol=:pos)

Translates the specified property of the particles or collection by the `period` such that all particles are the closest to `centerpos` as possible.

Creates a copy of the particles or the collection with only new pointers to the translated property.
Typically, `centerpos` will be an `AbstractVector`, e.g., for positions, and `period` can be a `Number` or
`AbstractVector`, depending on if the period is the same for all dimensions or if they are different.
"""
function translate_periodic(p::AbstractParticles, centerpos, period, prop::Symbol=:pos)
    p = copy(p)

    if haskey(p, prop)
        p[prop] = _translate_periodic.(p[prop], centerpos, period, 0.5period)
    end

    return p
end

function translate_periodic(pc::AbstractParticleCollection, centerpos, period, prop::Symbol=:pos)
    pc = copy(pc)

    Threads.@threads for ptype in keys(pc) |> collect # collect used for compatibility with threads
        pc[ptype] = translate_periodic(pc[ptype], centerpos, period, prop)
    end

    return pc
end

"""
    translate_periodic!(p::AbstractParticles, centerpos, period, prop::Symbol=:pos)
    translate_periodic!(pc::AbstractParticleCollection, centerpos, period, prop::Symbol=:pos)

In-place version of [`translate_periodic`](@ref).
"""
function translate_periodic!(p::AbstractParticles, centerpos, period, prop::Symbol=:pos)
    if haskey(p, prop)
        p[prop] .= _translate_periodic.(p[prop], centerpos, period, 0.5period)
    end

    return p
end

function translate_periodic!(pc::AbstractParticleCollection, centerpos, period, prop::Symbol=:pos)
    Threads.@threads for ptype in keys(pc) |> collect # collect used for compatibility with threads
        translate_periodic!(pc[ptype], centerpos, period, prop)
    end

    return pc
end

"""
    translate_periodic_to_center(p::AbstractParticles, centerpos, period, prop::Symbol=:pos)
    translate_periodic_to_center(pc::AbstractParticleCollection, centerpos, period, prop::Symbol=:pos)

Translates the specified property of the particles or collection to place `centerpos` at the origin.

The particles are translated periodically to place them closest to `centerpos`.

Creates a copy of the particles or the collection with only new pointers to the translated property.
Typically, `centerpos` will be an `AbstractVector`, e.g., for positions, and `period` can be a `Number` or
`AbstractVector`, depending on if the period is the same for all dimensions or if they are different.

Also see [`translate_periodic`](@ref).
"""
function translate_periodic_to_center(p::AbstractParticles, centerpos, period, prop::Symbol=:pos)
    p = copy(p)

    if haskey(p, prop)
        p[prop] = _translate_periodic.(p[prop], centerpos, period, 0.5period) .- centerpos
    end

    return p
end

function translate_periodic_to_center(pc::AbstractParticleCollection, centerpos, period, prop::Symbol=:pos)
    pc = copy(pc)

    Threads.@threads for ptype in keys(pc) |> collect # collect used for compatibility with threads
        pc[ptype] = translate_periodic_to_center(pc[ptype], centerpos, period, prop)
    end

    return pc
end

"""
    translate_periodic_to_center!(p::AbstractParticles, centerpos, period, prop::Symbol=:pos)
    translate_periodic_to_center!(pc::AbstractParticleCollection, centerpos, period, prop::Symbol=:pos)

In-place version of [`translate_periodic_to_center`](@ref).
"""
function translate_periodic_to_center!(p::AbstractParticles, centerpos, period, prop::Symbol=:pos)
    if haskey(p, prop)
        p[prop] .= _translate_periodic.(p[prop], centerpos, period, 0.5period) .- centerpos
    end

    return p
end

function translate_periodic_to_center!(pc::AbstractParticleCollection, centerpos, period, prop::Symbol=:pos)
    Threads.@threads for ptype in keys(pc) |> collect # collect used for compatibility with threads
        translate_periodic_to_center!(pc[ptype], centerpos, period, prop)
    end

    return pc
end

"""
    factor_to_comoving(z::Real, n::Real)

Returns ``(1 + z)^n``.

This is not exported.
"""
function factor_to_comoving(z::Real, n::Real)
    (1 + z)^n
end

"""
    factor_to_physical(z::Real, n::Real)

Returns ``(1 + z)^{-n}``.

This is not exported.
"""
function factor_to_physical(z::Real, n::Real)
    1 / (1 + z)^n
end

"""
    to_comoving(p::AbstractParticles, z::Real; propexp=[(:pos, 1), (:vel, 1)])
    to_comoving(pc::AbstractParticleCollection; propexp=[(:pos, 1), (:vel, 1)])

Create new particles or collection with particle properties converted from physical to comoving.

The properties and the positional exponent of their units (e.g., 1 for positions, 3 for volumes, and -3 for
densities) are passed as the keyword argument `propexp` as a vector of tuples.
The redshift `z` is obtained from [`redshift`](@ref) for the particle collection.

The properties are multiplied by ``(1 + z)^n`` where ``n`` is the positional exponent.
"""
function to_comoving end

"""
    to_comoving!(p::AbstractParticles, z::Real; propexp=[(:pos, 1), (:vel, 1)])
    to_comoving!(pc::AbstractParticleCollection; propexp=[(:pos, 1), (:vel, 1)])

Convert particle properties from physical to comoving in-place.

In-place version of [`to_comoving`](@ref).
"""
function to_comoving! end

"""
    to_physical(p::AbstractParticles, z::Real; propexp=[(:pos, 1), (:vel, 1)])
    to_physical(pc::AbstractParticleCollection; propexp=[(:pos, 1), (:vel, 1)])

Create new particles or collection with particle properties converted from comoving to physical.

The properties and the positional exponent of their units (e.g., 1 for positions, 3 for volumes, and -3 for
densities) are passed as the keyword argument `propexp` as a vector of tuples.
The redshift `z` is obtained from [`redshift`](@ref) for the particle collection.

The properties are multiplied by ``(1 + z)^{-n}`` where ``n`` is the positional exponent.
"""
function to_physical end

"""
    to_physical!(p::AbstractParticles, z::Real; propexp=[(:pos, 1), (:vel, 1)])
    to_physical!(pc::AbstractParticleCollection; propexp=[(:pos, 1), (:vel, 1)])

Convert particle properties from comoving to physical in-place.

In-place version of [`to_physical`](@ref).
"""
function to_physical! end

for (name, factorexpr) in zip(["to_comoving", "to_physical"], [:(1 + z), :(1 / (1 + z))])
    quote
        function $(Symbol(name))(p::AbstractParticles, z::Real; propexp=[(:pos, 1), (:vel, 1)])
            p = copy(p)

            factor = $factorexpr
            Threads.@threads for (prop, n) in propexp
                if haskey(p, prop)
                    p[prop] = product_preserve_type(p[prop], factor^n)
                end
            end

            return p
        end

        function $(Symbol(name))(pc::AbstractParticleCollection; propexp=[(:pos, 1), (:vel, 1)])
            pc = copy(pc)

            Threads.@threads for ptype in keys(pc) |> collect # collect used for compatibility with threads
                pc[ptype] = $(Symbol(name))(pc[ptype], redshift(pc); propexp=propexp)
            end

            return pc
        end

        function $(Symbol(name, "!"))(p::AbstractParticles, z::Real; propexp=[(:pos, 1), (:vel, 1)])
            factor = $factorexpr
            Threads.@threads for (prop, n) in propexp
                if haskey(p, prop)
                    product_preserve_type!(p[prop], factor^n)
                end
            end

            return p
        end

        function $(Symbol(name, "!"))(pc::AbstractParticleCollection; propexp=[(:pos, 1), (:vel, 1)])
            Threads.@threads for ptype in keys(pc) |> collect # collect used for compatibility with threads
                $(Symbol(name, "!"))(pc[ptype], redshift(pc); propexp=propexp)
            end

            return pc
        end
    end |> eval
end


text:
module Visualizations
export
    plot_tfbdry,
    wiggle,
    wiggle!

using
    Plots,
    Statistics
using
    ..Utils

"""
    treenodes_matrix(x)

Given a `BitVector` of nodes in a binary tree, output the matrix representation of the
nodes.

# Arguments
- `x::BitVector`: BitVector representing a binary tree, where an input is 1 if the
  corresponding node exists and has children, and 0 if the corresponding node does not
  exist/does not have children.

# Returns
`::BitMatrix`: BitMatrix representation of the tree, where each column corresponds to a
level of the binary tree. The inputs in the matrix are 1 if the corresponding node exists
and has children, and 0 if the corresponding node does not exist/does not have children.

# Visualization
```
# Binary tree
      x
     / \\
    x   o
   / \\
  x   o
 / \\
o   o

# BitVector representation
[1,1,0,1,0,0,0]

# BitMatrix representation
[1 1 1 1;
 1 1 0 0;
 1 0 0 0]
```

# Examples
```@repl
using Wavelets, WaveletsExt

tree = maketree(8, 3, :dwt)
WaveletsExt.Visualizations.treenodes_matrix(tree)
```
"""
function treenodes_matrix(x::BitVector)
    n = (length(x) + 1) >> 1
    L = floor(Integer, log2(n))

    result = falses(L+1, n)
    i = 1
    for lvl in 0:L
        c = 1 << lvl                  # number of nodes at current level
        n₀ = nodelength(n, lvl)       # length of nodes at current level
        for node in 0:(c - 1)
            rng = (node * n₀ + 1):((node + 1) * n₀)
            result[lvl+1, rng] = repeat([x[i]], n₀)
            i += 1
        end
    end

    return result
end

"""
    plot_tfbdry(tree[; start, nd_col, ln_col, bg_col])

Given a tree, output a visual representation of the leaf nodes, user will have the option to
start the node count of each level with 0 or 1.

# Arguments
- `tree::BitVector`: Tree for plotting the leaf nodes. Comes in the form of a `BitVector`.
- `depth::Integer`: (Default: `log2(length(tree)+1)-1 |> Int`) Maximum depth to be
  displayed.

# Keyword Arguments
- `start::Integer`: (Default: `0`) Whether to zero-index or one-index the root of the tree.
- `nd_col::Symbol`: (Default: `:white`) Color of the leaf nodes.
- `ln_col::Symbol`: (Default: `:white`) Color of lines in plot.
- `bg_col::Symbol`: (Default: `:black`) Color of background.

# Returns
`::Plots.Plot`: Plot object with the visual representation of the leaf nodes.

# Examples
```julia
using Wavelets, WaveletsExt

# Build a tree using Wavelets `maketree`
tree = maketree(128, 7, :dwt)

# Plot the leaf nodes
plot_tfbdry(tree)
```
"""
function plot_tfbdry(tree::BitVector,
                     depth::Integer = log2(length(tree)+1)-1 |> Int;
                     start::Integer = 0,
                     nd_col::Symbol = :white,
                     ln_col::Symbol = :white,
                     bg_col::Symbol = :black)

    @assert 0 <= start <= 1
    leaf = getleaf(tree,:binary)

    ncol = 1 << depth
    nrow = depth+1
    mat = treenodes_matrix(leaf[1:(1<<(depth+1)-1)])

    p = heatmap(0:(ncol-1), start:(nrow+start-1), mat, color = [bg_col, nd_col],
                legend = false, background_color = bg_col)

    plot!(p, xlims = (start-0.5, ncol+start-0.5), ylims = (-0.5, nrow-0.5), yticks = 0:nrow)
    # plot horizontal lines
    for j in 0:(nrow-1)
        x_rng = (start-0.5):(ncol+start-0.5)
        y_val = (j-0.5)*ones(ncol+1)
        plot!(p, x_rng, y_val, color = ln_col, legend = false)
    end

    # plot vertical lines
    for j in 1:(nrow-1)
        for jj in 1:2^(j-1)
            x_val = (ncol/2^j)*(2*jj-1)+start-0.5;
            y_rng = j-0.5:nrow-0.5
            plot!(p, x_val*ones(nrow-j+1), y_rng, color = ln_col)
        end
    end
    plot!(p, (ncol+start-0.5)*ones(nrow+1), -0.5:nrow-0.5, color = ln_col)
    plot!(p, (start-0.5)*ones(nrow+1), -0.5:nrow-0.5, color = ln_col)
    plot!(p, yaxis = :flip)

    return p
end

"""
    wiggle(wav[; taxis, zaxis, sc, EdgeColor, FaceColor, Orient, Overlap, ZDir])
    wiggle(plt, wav[; taxis, zaxis, sc, EdgeColor, FaceColor, Orient, Overlap, ZDir])

Plots a set of shaded wiggles.

# Arguments
- `plt::Plots.Plot`: Input plot to plot shaded wiggles.
- `wav::AbstractArray{<:Number,2}`: Matrix of waveform columns.

# Keyword Arguments
- `taxis::AbstractVector`: (Default: `1:size(wav,1)`) Time axis vector
- `zaxis::AbstractVector`: (Default: `1:size(wav,2)`) Space axis vector
- `sc::Real`: (Default: `1`) Scale factor/magnification.
- `EdgeColor::Symbol`: (Default: `:black`) Sets edge of wiggles color.
- `FaceColor::Symbol`: (Default: `:black`) Sets shading color of wiggles.
- `Overlap::Bool`: (Default: `true`) How signals are scaled.
    - `true`  - Signals overlap (default);
    - `false` - Signals are scaled so they do not overlap.
- `Orient::Symbol`: (Default: `:across`) Controls orientation of wiggles.
    - `:across` - from left to right
    - `:down`   - from top to down
- `ZDir::Symbol`: (Default: `:normal`) Direction of space axis.
    - `:normal`  - First signal at bottom (default)
    - `:reverse` - First signal at top.

# Returns
`::Plots.Plot`: Shaded wiggles on top of current plot object.

# Examples
```julia
using Plots, WaveletsExt

# Generate random signals
x = randn(16, 5)

# ----- Build wiggles -----
# Method 1
wiggle(x)

# Method 2
p = Plot()
wiggle(p, x)
```

Translated by Nicholas Hausch -- MATLAB file provided by Naoki Saito. The previous MATLAB
version contributors are Anthony K. Booer (SLB) and Bradley Marchand (NSWC-PC).

Revised by Naoki Saito, Feb. 05, 2018. Maintained by Zeng Fung Liew for newest Julia version
compatibility.

**See also:** [`wiggle!`](@ref)
"""
wiggle(args...; kwargs...) = wiggle!(Plots.Plot(), args...; kwargs...)
wiggle(plt::Plots.Plot, args...; kwargs...) = wiggle!(deepcopy(plt), args...; kwargs...)

"""
    wiggle!(wav[; taxis, zaxis, sc, EdgeColor, FaceColor, Orient, Overlap, ZDir])
    wiggle!(plt, wav[; taxis, zaxis, sc, EdgeColor, FaceColor, Orient, Overlap, ZDir])

Plot a set of shaded wiggles on the current displayed graphics or on top of `plt`. If there
are no displayed graphics currently available, a new `Plots.Plot` object is generated to
plot the shaded wiggles.

# Arguments
- `plt::Plots.Plot`: Input plot to plot shaded wiggles.
- `wav::AbstractArray{<:Number,2}`: Matrix of waveform columns.

# Keyword Arguments
- `taxis::AbstractVector`: (Default: `1:size(wav,1)`) Time axis vector
- `zaxis::AbstractVector`: (Default: `1:size(wav,2)`) Space axis vector
- `sc::Real`: (Default: `1`) Scale factor/magnification.
- `EdgeColor::Symbol`: (Default: `:black`) Sets edge of wiggles color.
- `FaceColor::Symbol`: (Default: `:black`) Sets shading color of wiggles.
- `Overlap::Bool`: (Default: `true`) How signals are scaled.
    - `true`  - Signals overlap (default);
    - `false` - Signals are scaled so they do not overlap.
- `Orient::Symbol`: (Default: `:across`) Controls orientation of wiggles.
    - `:across` - from left to right
    - `:down`   - from top to down
- `ZDir::Symbol`: (Default: `:normal`) Direction of space axis.
    - `:normal`  - First signal at bottom (default)
    - `:reverse` - First signal at top.

# Returns
`::Plots.Plot`: Shaded wiggles on top of current plot object.

# Examples
```julia
using Plots, WaveletsExt

# Generate random signals
x = randn(16, 5)

# ----- Build wiggles -----
# Build onto existing plot
plt = plot()
wiggle!(x)

# Build onto a specified plot
wiggle!(plt, x)
```

Translated by Nicholas Hausch -- MATLAB file provided by Naoki Saito. The previous MATLAB
version contributors are Anthony K. Booer (SLB) and Bradley Marchand (NSWC-PC).

Revised by Naoki Saito, Feb. 05, 2018. Maintained by Zeng Fung Liew for newest Julia version
compatibility.

**See also:** [`wiggle`](@ref)
"""
function wiggle!(args...; kwargs...)
    @nospecialize
    local plt
    try
        plt = current()
    catch
        return wiggle(args..., kwargs...)
    end
    wiggle!(current(), args...; kwargs...)
end

function wiggle!(plt::Plots.Plot,
                 wav::AbstractArray{T,2};
                 taxis::AbstractVector=1:size(wav,1),
                 zaxis::AbstractVector=1:size(wav,2),
                 sc::Real=1,
                 EdgeColor::Symbol=:black,
                 FaceColor::Symbol=:black,
                 Overlap::Bool=true,
                 Orient::Symbol=:across,
                 ZDir::Symbol=:normal) where T<:Number
    # Set axes
    (n, m) = size(wav)

    # Sanity check
    @assert Orient ∈ [:across, :down]
    @assert ZDir ∈ [:normal, :reverse]
    if length(taxis) != n
        error("Inconsistent taxis dimension!")
    end
    if length(zaxis) != m
        error("Inconsistent zaxis dimension!")
    end

    # For calculation purposes
    maxrow = zeros(m)
    minrow = zeros(m)
    for (i, wavᵢ) in enumerate(eachcol(wav))
        maxrow[i] = maximum(wavᵢ)
        minrow[i] = minimum(wavᵢ)
    end

    # Scale the data for plotting
    dz = mean(diff(zaxis))
    if Overlap
        wamp = 2 * dz * (sc/maximum(maxrow-minrow)) * wav
    else
        wmax = maximum(maxrow) <= 0 ? 0 : maximum(maxrow)
        wmin = minimum(minrow) >= 0 ? 0 : minimum(minrow)
        wamp = sc*wav/(wmax-wmin)
    end

    # Set initial plot
    t0 = minimum(taxis)
    t1 = maximum(taxis)
    z0 = minimum(zaxis)
    z1 = maximum(zaxis)
    if Orient == :down
        plot!(plt, xlims=(z0-dz,z1+dz), ylims=(t0,t1), yflip=true, legend=:none)
    else
        plot!(plt, xlims=(t0,t1), ylims=(z0-dz,z1+dz), legend=:none)
    end
    if ZDir == :reverse
        wamp = reverse(wamp, dims=2)
    end

    # Plot each wavelet
    for (i, wampᵢ) in enumerate(eachcol(wamp))
        t = deepcopy(taxis)
        w_sign = sign.(wampᵢ)
        for j in 1:(n-1)
            if (w_sign[j]!=w_sign[j+1] && w_sign[j]!=0 && w_sign[j+1]!=0)
                wampᵢ = [wampᵢ; 0]
                t = [t; t[j]-wampᵢ[j]*(t[j+1]-t[j])/(wampᵢ[j+1]-wampᵢ[j])]
            end
        end
        ix = sortperm(t)
        t = t[ix]
        wampᵢ = wampᵢ[ix]
        len = length(t)
        indperm = [1:len; len:-1:1]
        inputx = t[indperm]
        inputy = zaxis[i] .+ [wampᵢ; min.(wampᵢ[len:-1:1],0)]
        # In the plot! functions below, theoretically speaking, either fillrange = zaxis[k]
        # or fillrange=[zaxis[k], zaxis[k]+dz] should be used. However, those do not
        # generate the desired plots as of O3/19/2018. Somehow, the relative value of 0,
        # i.e., fillrange=0, works well, which is used temporarily.
        if Orient == :down
            plot!(plt, inputy, inputx, fillrange=0, fillalpha=0.75,
                  fillcolor=FaceColor, linecolor=EdgeColor, orientation=:v)
        else
            plot!(plt, inputx, inputy, fillrange=0, fillalpha=0.75,
                  fillcolor=FaceColor, linecolor=EdgeColor)
        end
    end
    return plt
end

end # end module


text:
# Schamber's Vector quant algorithm - It foregoes the a weighted least squares fit for a simple
# homoskedastic fit.  In doing this it can precompute a fitting matrix which requires nothing
# more than a single matrix multiplication to perform the fit.  This makes this mechanism
# extremely quick. This makes it ideal for processing in real-time or HyperSpectrum objects.

struct VectorQuant
    # Vector(label[1], roi[2], charonly[3], sum(charonly)[4], scale[5])
    references::Vector{Tuple{ReferenceLabel,UnitRange,Vector{Float64},Float64,Float64}}
    vectors::Matrix{Float64}

    """
        VectorQuant(frefs::Vector{FilteredReference}, filt::TopHatFilter)

    Constructs a structure used to perform accelerated filtered spectrum fits based on the specified
    collection of `FilteredReference`(s), and a `TopHatFilter`.
    """
    function VectorQuant(frefs::Vector{FilteredReference}, filt::TopHatFilter)
        refs = [
            (fref.label, fref.roi, fref.charonly, sum(fref.charonly), fref.scale) for
            fref in frefs
        ]
        x = zeros(Float64, (length(filt.filters), length(frefs)))
        for (c, fref) in enumerate(frefs)
            x[fref.ffroi, c] = fref.filtered
        end
        # ((ch × ne)T * (ch × ne))^(-1) * (ch × ne) * (ch × ch) => (ne × ch)
        xTxIxf =
            pinv(transpose(x) * x) *
            transpose(x) *
            NeXLSpectrum.filterdata(filt, 1:length(filt.filters))
        return new(refs, xTxIxf)
    end

    VectorQuant(ffrs::FilterFitPacket) = VectorQuant(ffrs.references, ffrs.filter)
end

NeXLCore.minproperties(::VectorQuant) = (:BeamEnergy, :TakeOffAngle, :)

Base.show(io::IO, vq::VectorQuant) = print(
    io,
    "VectorQuant[\n" * join(map(r -> "\t" * repr(r[1]), vq.references), ",\n") * "\n]",
)

function fit_spectrum(
    vq::VectorQuant,
    spec::Spectrum,
    zero = x -> max(0.0, x),
)::FilterFitResult
    raw = counts(spec, Float64)
    krs = zero.(vq.vectors * raw)
    spsc = dose(spec)
    residual = copy(raw)
    for (i, (_, roi, co, _, _)) in enumerate(vq.references)
        residual[roi] -= krs[i] * co
    end
    peakback = Dict{ReferenceLabel,NTuple{2,Float64}}()
    dkrs = zeros(Float64, length(vq.references))
    for (i, (lbl, roi, _, ico, _)) in enumerate(vq.references)
        ii, bb = krs[i] * ico, sum(residual[roi])
        peakback[lbl] = (ii, bb)
        dkrs[i] = sqrt(max(0.0, ii + bb)) / ico
    end
    kratios = uvs(
        map(ref -> ref[1], vq.references), #
        map(i -> krs[i] / (vq.references[i][5] * spsc), eachindex(krs)), #
        map(i -> (dkrs[i] / (vq.references[i][5] * spsc))^2, eachindex(krs)),
    )
    return FilterFitResult(
        UnknownLabel(spec),
        kratios,
        1:length(raw),
        raw,
        residual,
        peakback,
    )
end

function fit_spectrum(
    vq::VectorQuant,
    hs::HyperSpectrum,
    zero = x -> max(0.0, x),
)::Array{KRatios}
    krs = zeros(Float32, length(vq.references), size(hs)...)
    vecs = vq.vectors[:, 1:depth(hs)]
    scales = [ vq.references[i][5] for i in eachindex(vq.references) ]
    data = counts(hs)
    # @threads seems to slow this (maybe cache misses??)
    for ci in CartesianIndices(hs)
        @inbounds @avx krs[:, ci] = (vecs * data[:, ci]) ./ (dose(hs,ci) * scales)
    end
    # ensure positive...
    map!(zero, krs, krs)
    res = KRatios[]
    for i in filter(ii -> vq.references[ii][1] isa CharXRayLabel, eachindex(vq.references))
        k, lbl = krs[i], vq.references[i][1]
        rprops = properties(spectrum(lbl))
        push!(
            res,
            KRatios(xrays(lbl), properties(hs), rprops, rprops[:Composition], krs[i, :, :]),
        )
    end
    return res
end


text:
module JLBoostMLJ

using JLBoost

export AUC, gini, feature_importance

include("mlj.jl")

end # module


text:
using DateLib
using Base.Test

for i=1:100
	FOY=DateClass(1,1,1890+i);
	@test 1==FOY.day;
	@test 1==FOY.month;
	@test 1890+i==FOY.year;
end
println("Simple Test Passed")

println("Starting Domain Test")
for i=0:100
	@test_throws(ErrorException, DateClass(1,1,1900+i*1.0000001))
end
for i=0:11
	@test_throws(ErrorException, DateClass(1,1+i*1.0000001,1900))
end
for i=0:20
	@test_throws(ErrorException, DateClass(1+i*1.0000001,1,1900+i))
end
println("Domain Test Passed")

println("Starting Last of February Test")

FileLeapYears=open("LeapYears.txt");
LeapYearStrings = readlines(FileLeapYears);
LeapYears=Array{Int64}(length(LeapYearStrings));
for i=1:length(LeapYearStrings)
	LeapYears[i]=parse(Int64,LeapYearStrings[i]);
end

ii=1;
for i=1904:2076
	if i==LeapYears[ii]
		@test DateClass(29,2,i)==DateClass(29,2,i);
		@test true==isLeapYear(LeapYears[ii])
		ii+=1
	else
		@test_throws(ErrorException, DateClass(29,2,i))
		@test false==isLeapYear(i)
	end
end
println("Test Last of February Passed")


text:
"""

abstract type SpeciesDynamics{F<:Species} end

Each SpeciesDynamics struct needs a functor:
    (s::SomeDynamics)(n::Neuron, variable, currents)

"""

get_parameters(::SpeciesDynamics) = nothing
get_states(::SpeciesDynamics) = nothing
default_params(::SpeciesDynamics, a, b, c) = Dict{Num,Float64}()
default_states(::SpeciesDynamics, a, b, c) = Dict{Num,Float64}()

struct BasicVoltageDynamics <: SpeciesDynamics{Voltage} end

function (b::BasicVoltageDynamics)(n::Neuron, vars, varnames, flux)
    Cₘ, = get_parameters(b)
    V = vars[findfirst(x -> x == Voltage, varnames)]
    return D(V) ~ (1 / Cₘ) * (flux) #non-standard convention, sum of fluxes already has negative sign because of the (E-V) in currents
end

get_parameters(::BasicVoltageDynamics) = @parameters Cₘ
default_params(v::BasicVoltageDynamics, n::Neuron, vars, varnames) = Dict(get_parameters(v)... => capacitance(n.geometry))


text:
function grayscale( filename::String; path::String=pwd(), ftyp=Float32, ityp=Int16 )

	occursin( '/', filename ) && ( path = "" )

	img = FileIO.load( path*filename )

	# Temporary array, used to store the grayscale values of "img" as Floats

	imgf = zeros( ftyp, size( img ) )

	for idx in 1:length(img)
		imgf[idx] = ftyp( ColorTypes.Gray( img[idx] ) )
	end

	# The grayscale image is stored in "imgG" as Integer values

	imgG = Array{ ityp, 2 }(undef, size(img) )
	for li in 1:length( imgG )
		imgG[li] = 1 + round( ityp, imgf[li]* 254 )
	end

	return imgG
end

function grayscale!( filename::String, img::Array{<:Integer,2}; path::String=pwd(), ftyp=Float32 )

	occursin( '/', filename ) && ( path = "" )
	( isfile( path*filename) ) || error( "file does not exists" )

	img_ = FileIO.load( path*filename )

	( size( img_ ) == size(img) ) || error( "input image should have dimensions: ", size(img_) )

	T = eltype( img );

	for li in 1:length( img )
		img[li] = 1 + round( T, ftyp( ( img_[li].val ) )*254 )
	end
end

function grayscale255( filename::String; path::String=pwd(), typ=Float32 )

	occursin( '/', filename ) && ( path = "" )

	img = FileIO.load( path*filename )

	imgf = zeros( typ, size( img ) )

	for idx in 1:length(img)
		imgf[idx] = typ( ColorTypes.Gray( img[idx] ) )
	end

	mn, mx = extrema( imgf )

	im255 = Array{ Int64, 2 }(undef, size(img) )
	for e in 1:length( img )
		im255[e] = 1 + round( Int64, ( imgf[e] - mn )/( mx - mn ) * 254 )
	end

	return im255
end

function grayscale255rev( filename::String; path::String=pwd() )

	occursin( '/', filename ) && ( path = "" )

	img = load( path*filename )

	imgf = zeros( Float64, size( img ) )

	for idx in 1:length(img)
		imgf[idx] = Float64( ColorTypes.Gray( img[idx] ) )
	end

	mn, mx = extrema( imgf )

	im255 = Array{ Int64, 2 }(undef, size(imgf) )
	for e in 1:length( img )
		im255[e] = 256 - 1 - round( Int64, ( imgf[e] - mn )/( mx - mn ) * 254 )
	end

	return im255
end

function RGBfloat( filename::String; path::String=pwd(), typ=Float32 )

	occursin( '/', filename ) && ( path = "" )

	data = load( path*filename )

	RGBimg = Array{ typ, 3 }( undef, size( data, 1 ), size( data, 2 ), 3 )

	Goff = length( data );
	Boff = Goff*2;

	@inbounds for li in 1:length(data)
	    pixel = data[li]
	    RGBimg[   li  ] = typ( data[li].r )
	    RGBimg[li+Goff] = typ( data[li].g )
	    RGBimg[li+Boff] = typ( data[li].b )
	end

	return RGBimg
end


text:
function detectCollision(
    dlnParams::DislocationParameters,
    network::DislocationNetwork,
    skipSegs = Vector{Tuple{Int, Int}}(),
)
    collisionDistSq = dlnParams.collisionDistSq
    label = network.label
    numNode = network.numNode[1]
    numSeg = network.numSeg[1]
    coord = network.coord
    nodeVel = network.nodeVel
    links = network.links
    bVec = network.bVec
    connectivity = network.connectivity

    elemT = eltype(coord)
    smallestMinDist = zero(elemT)
    collision = false
    collisionType = :null
    n1s1, n2s1, n1s2, n2s2, s1, s2 = 0, 0, 0, 0, 0, 0
    L1, L2 = zero(elemT), zero(elemT)
    mobileNodes = (intMobDln, srfMobDln)

    # Hinge condition.
    for i in 1:numNode
        # Only hinge collisions with mobile nodes.
        label[i] ∉ mobileNodes && continue

        for j in 1:connectivity[1, i]
            # Find connected node.
            link1, rowColLink1, connectedNode1 = findConnectedNode(network, i, j)
            # Make t vector.
            tVec = SVector{3, elemT}(
                coord[1, i] - coord[1, connectedNode1],
                coord[2, i] - coord[2, connectedNode1],
                coord[3, i] - coord[3, connectedNode1],
            )
            # Square of the norm.
            tVecN1 = tVec ⋅ tVec

            # Other connections.
            for k in (j + 1):connectivity[1, i]
                link2, rowColLink2, connectedNode2 = findConnectedNode(network, i, k)
                tVec = SVector{3, elemT}(
                    coord[1, i] - coord[1, connectedNode2],
                    coord[2, i] - coord[2, connectedNode2],
                    coord[3, i] - coord[3, connectedNode2],
                )
                tVecN2 = tVec ⋅ tVec

                # We collapse one of the lines into a point, we want to make sure we collapse the smallest distance so the distance calculation is the actual minimum.
                tVecN1 < tVecN2 ?
                begin
                    link2,
                    rowColLink2,
                    connectedNode2,
                    tVecN2,
                    link1,
                    rowColLink1,
                    connectedNode1,
                    tVecN1 = link1,
                    rowColLink1,
                    connectedNode1,
                    tVecN1,
                    link2,
                    rowColLink2,
                    connectedNode2,
                    tVecN2
                end : nothing

                # tVecN1 should now be bigger than tVecN2. i is the hinge node.
                x0 = SVector{3, elemT}(coord[1, i], coord[2, i], coord[3, i])
                vx0 = SVector{3, elemT}(nodeVel[1, i], nodeVel[2, i], nodeVel[3, i])

                x1 = SVector{3, elemT}(
                    coord[1, connectedNode1],
                    coord[2, connectedNode1],
                    coord[3, connectedNode1],
                )
                vx1 = SVector{3, elemT}(
                    nodeVel[1, connectedNode1],
                    nodeVel[2, connectedNode1],
                    nodeVel[3, connectedNode1],
                )

                # We collapse the other segment into a point.
                y0 = SVector{3, elemT}(
                    coord[1, connectedNode2],
                    coord[2, connectedNode2],
                    coord[3, connectedNode2],
                )
                vy0 = SVector{3, elemT}(
                    nodeVel[1, connectedNode2],
                    nodeVel[2, connectedNode2],
                    nodeVel[3, connectedNode2],
                )

                skip = false
                for k in 1:length(skipSegs)
                    s1Tmp, s2Tmp = skipSegs[k]
                    if link1 == s1Tmp && link2 == s2Tmp || link1 == s2Tmp && link2 == s1Tmp
                        skip = true
                        break
                    end
                end

                skip && continue

                distSq, dDistSqDt, L1, L2 =
                    minimumDistance(x0, x1, y0, y0, vx0, vx1, vy0, vy0)

                if distSq < collisionDistSq && dDistSqDt < -eps(elemT)
                    smallestMinDistTmp = 1 / distSq
                    if smallestMinDistTmp > smallestMinDist
                        smallestMinDist = smallestMinDistTmp
                        collision = true
                        collisionType = :hinge
                        n1s1 = i # Hinge node
                        n2s1 = connectedNode1
                        n1s2 = connectedNode2
                        n2s2 = n1s2
                        s1 = link1
                        s2 = link2
                    end
                end
            end
        end
    end

    collision && return collision, collisionType, n1s1, n2s1, n1s2, n2s2, s1, s2, L1, L2

    # Collision of unconnected segments, twoline collision.
    for i in 1:numSeg
        n1s1_i = links[1, i]
        n2s1_i = links[2, i]

        label1 = label[n1s1_i]
        label2 = label[n2s1_i]

        (label1 ∉ mobileNodes || label2 ∉ mobileNodes) && continue

        x0 = SVector{3, elemT}(coord[1, n1s1_i], coord[2, n1s1_i], coord[3, n1s1_i])
        vx0 = SVector{3, elemT}(nodeVel[1, n1s1_i], nodeVel[2, n1s1_i], nodeVel[3, n1s1_i])

        x1 = SVector{3, elemT}(coord[1, n2s1_i], coord[2, n2s1_i], coord[3, n2s1_i])
        vx1 = SVector{3, elemT}(nodeVel[1, n2s1_i], nodeVel[2, n2s1_i], nodeVel[3, n2s1_i])

        bi = SVector{3, elemT}(bVec[1, i], bVec[2, i], bVec[3, i])

        for j in (i + 1):numSeg
            n1s2_j = links[1, j]
            n2s2_j = links[2, j]

            # Exclude hinges.
            if n1s1_i != n1s2_j && n1s1_i != n2s2_j && n2s1_i != n1s2_j && n2s1_i != n2s2_j
                label1 = label[n1s2_j]
                label2 = label[n2s2_j]

                (label1 ∉ mobileNodes || label2 ∉ mobileNodes) && continue

                y0 = SVector{3, elemT}(coord[1, n1s2_j], coord[2, n1s2_j], coord[3, n1s2_j])
                vy0 = SVector{3, elemT}(
                    nodeVel[1, n1s2_j],
                    nodeVel[2, n1s2_j],
                    nodeVel[3, n1s2_j],
                )

                y1 = SVector{3, elemT}(coord[1, n2s2_j], coord[2, n2s2_j], coord[3, n2s2_j])
                vy1 = SVector{3, elemT}(
                    nodeVel[1, n2s2_j],
                    nodeVel[2, n2s2_j],
                    nodeVel[3, n2s2_j],
                )

                bj = SVector{3, elemT}(bVec[1, j], bVec[2, j], bVec[3, j])

                # Stop superdislocations from forming.
                (
                    (
                        (bi - bj) ⋅ (bi - bj) < eps(elemT) &&
                        (x1 - x0) ⋅ (y1 - y0) > eps(elemT)
                    ) || (
                        (bi + bj) ⋅ (bi + bj) < eps(elemT) &&
                        (x1 - x0) ⋅ (y1 - y0) < eps(elemT)
                    )
                ) && continue

                skip = false
                for k in 1:length(skipSegs)
                    s1Tmp, s2Tmp = skipSegs[k]
                    if i == s1Tmp && j == s2Tmp || i == s2Tmp && j == s1Tmp
                        skip = true
                        break
                    end
                end

                skip && continue

                distSq, dDistSqDt, L1, L2 =
                    minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)

                if distSq < collisionDistSq && dDistSqDt < -eps(elemT)
                    smallestMinDistTmp = 1 / distSq
                    if smallestMinDistTmp > smallestMinDist
                        smallestMinDist = smallestMinDistTmp
                        collision = true
                        collisionType = :twoLine
                        n1s1 = n1s1_i
                        n2s1 = n2s1_i
                        n1s2 = n1s2_j
                        n2s2 = n2s2_j
                        s1 = i
                        s2 = j
                    end
                end
            end
        end
    end

    return collision, collisionType, n1s1, n2s1, n1s2, n2s2, s1, s2, L1, L2
end

function resolveCollision!(
    dlnParams::DislocationParameters,
    matParams::MaterialParameters,
    network::DislocationNetwork,
    collisionType,
    n1s1,
    n2s1,
    n1s2,
    n2s2,
    s1,
    s2,
    L1,
    L2,
)
    collisionDist = dlnParams.collisionDist
    collisionDistSq = dlnParams.collisionDistSq
    minSegLenSq = dlnParams.minSegLenSq
    coord = network.coord
    nodeVel = network.nodeVel
    segForce = network.segForce
    connectivity = network.connectivity
    linksConnect = network.linksConnect
    elemT = eltype(coord)

    if collisionType == :twoLine
        # Identify closeness to nodes of segments 1 and 2.
        mergeNode1, network = findSplitNode!(dlnParams, network, n1s1, n2s1, s1, L1)
        mergeNode2, network = findSplitNode!(dlnParams, network, n1s2, n2s2, s2, L2)
    elseif collisionType == :hinge
        # The hinge node is n1s1. We merge node n2s1, into node n1s2.
        mergeNode1 = n1s2
        mergeNode2 = n2s1
    end

    powerPreCollision = calcPowerDissipation(network, mergeNode1)
    powerPreCollision += calcPowerDissipation(network, mergeNode2)
    powerPreCollision *= 1.05

    collisionPoint = findCollisionPoint(network, mergeNode1, mergeNode2)
    coord[:, mergeNode1] = collisionPoint
    nodeVel[:, mergeNode1] .= 0

    nodeKept, network = mergeNode!(network, mergeNode1, mergeNode2)

    # If the node still exists.
    if nodeKept > 0
        connectivity = network.connectivity
        # Calculate mobility of nodes connected to it.
        links, missing, connectedNodes =
            findConnectedNode(network, nodeKept, 1:connectivity[1, nodeKept])
        calcSegForce!(dlnParams, matParams, network, links)
        dlnMobility!(dlnParams, matParams, network, connectedNodes)
        # Calculate mobility of the node.
        dlnMobility!(dlnParams, matParams, network, nodeKept)
    end

    return powerPreCollision, network
end

function findSplitNode!(dlnParams, network, n1, n2, s, L)
    minSegLenSq = dlnParams.minSegLenSq
    collisionDist = dlnParams.collisionDist

    coord = network.coord
    nodeVel = network.nodeVel
    linksConnect = network.linksConnect
    elemT = eltype(network.coord)

    tVec = SVector{3, eltype(coord)}(
        coord[1, n2] - coord[1, n1],
        coord[2, n2] - coord[2, n1],
        coord[3, n2] - coord[3, n1],
    )
    tVecSq = tVec ⋅ tVec

    close_n1 = ((L^2 * tVecSq) < minSegLenSq)
    close_n2 = (((1 - L)^2 * tVecSq) < minSegLenSq)

    if close_n1 && L <= 0.5
        mergeNode1 = n1
    elseif close_n2
        mergeNode1 = n2
    else
        splitNode = n1
        splitConnect = linksConnect[1, s]

        if close_n1
            L = collisionDist / sqrt(tVec)
        else
            L = 1 - collisionDist / sqrt(tVec)
        end

        midCoord =
            SVector{3, elemT}(coord[1, n1], coord[2, n1], coord[3, n1]) * (1 - L) +
            SVector{3, elemT}(coord[1, n2], coord[2, n2], coord[3, n2]) * L

        midVel =
            SVector{3, elemT}(nodeVel[1, n1], nodeVel[2, n1], nodeVel[3, n1]) * (1 - L) +
            SVector{3, elemT}(nodeVel[1, n2], nodeVel[2, n2], nodeVel[3, n2]) * L

        network = splitNode!(network, splitNode, splitConnect, midCoord, midVel)
        mergeNode = network.numNode[1]
    end

    return mergeNode, network
end

function calcPowerDissipation(network, node)
    nodeVel = network.nodeVel
    segForce = network.segForce
    connectivity = network.connectivity
    elemT = eltype(segForce)

    force = SVector(zeros(elemT, 3))
    c = connectivity[1, node]
    for i in 1:c
        link = connectivity[2 * i, node]
        pos = connectivity[2 * i + 1, node]
        force += SVector{3, elemT}(
            segForce[1, 3 - pos, link],
            segForce[2, 3 - pos, link],
            segForce[3, 3 - pos, link],
        )
    end
    nodeVelTmp = SVector{3, elemT}(nodeVel[1, node], nodeVel[2, node], nodeVel[3, node])
    powerDissipation = nodeVelTmp ⋅ force

    return powerDissipation
end

function findCollisionPoint(network::DislocationNetwork, node1, node2)
    bVec = network.bVec
    label = network.label
    coord = network.coord
    nodeVel = network.nodeVel
    connectivity = network.connectivity
    elemT = eltype(coord)

    p1 = SVector{3, elemT}(coord[1, node1], coord[2, node1], coord[3, node1])
    p2 = SVector{3, elemT}(coord[1, node2], coord[2, node2], coord[3, node2])
    label1 = label[node1]
    label2 = label[node2]

    if label1 == intFixDln
        collisionPoint = p1
        return collisionPoint
    elseif label2 == intFixDln
        collisionPoint = p2
        return collisionPoint
    end

    N = 0
    matN = zeros(elemT, MMatrix{3, 3})
    vecN = zeros(elemT, MVector{3})

    N = findPlaneCondition!(matN, vecN, N, network, node1, p1)
    N < 3 ? findPlaneCondition!(matN, vecN, N, network, node2, p2) : nothing

    A = SMatrix{6, 6, elemT}([I(3) matN[:, 1:3]'; matN[:, 1:3] zeros(3, 3)])
    b = [(p1 + p2) / 2; vecN]
    x = A \ b
    collisionPoint = x[1:3]

    return collisionPoint
end

function findPlaneCondition!(matN, vecN, N, network, node, point)
    bVec = network.bVec
    coord = network.coord
    nodeVel = network.nodeVel
    connectivity = network.connectivity
    elemT = eltype(coord)

    newPlaneCondition = 0.875

    for i in 1:connectivity[1, node]
        if N < 3
            link, missing, connectedNode = findConnectedNode(network, node, i)
            tVec =
                point - SVector{3, elemT}(
                    coord[1, connectedNode],
                    coord[2, connectedNode],
                    coord[3, connectedNode],
                )
            tVec = normalize(tVec)

            b = SVector{3, elemT}(bVec[1, link], bVec[2, link], bVec[3, link])
            vel = SVector{3, elemT}(nodeVel[1, node], nodeVel[2, node], nodeVel[3, node])

            n1 = tVec × b
            n2 = tVec × vel

            n1Sq = n1 ⋅ n1
            n2Sq = n2 ⋅ n2

            if n1Sq > eps(elemT)
                plane = normalize(n1)
            elseif n2Sq > eps(elemT)
                plane = normalize(n2)
            end

            if n1Sq > eps(elemT) || n2Sq > eps(elemT)
                if N == 0
                    condition = true
                elseif N == 1
                    condition = (matN[:, 1] ⋅ plane)^2 < newPlaneCondition^2
                else
                    detN = det(
                        SMatrix{3, 3, elemT}(
                            matN[1, 1],
                            matN[2, 1],
                            matN[3, 1],
                            matN[1, 2],
                            matN[2, 2],
                            matN[3, 2],
                            plane[1],
                            plane[2],
                            plane[3],
                        ),
                    )
                    condition = detN^2 < (1 - newPlaneCondition)^2
                end

                if condition
                    N += 1
                    matN[:, N] = plane
                    vecN[N] = plane ⋅ point
                end
            end
        end
    end
    return N
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD