text:
# ===================================================================================
# AbstractTypes

abstract type AbstractTree{N} end
abstract type AbstractNode{T} end
abstract type AbstractBinaryNode{T} <: AbstractNode{T} end

# ------------------------------------------------------------------------------------
# Nodes

"""
    NullNode()

Empty node for leaves and roots.
"""
struct NullNode <: AbstractNode{Nothing} end

Base.getproperty(::NullNode, ::Symbol) = nothing
Base.getproperty(::NullNode, ::Nothing) = nothing

"""
    SimpleBinaryNode{T} <: AbstractBinaryNode{T}
    const SBN = SimpleBinaryNode
    SBN(data)
    SBN(::Type{T}, data)

Binary nodes without height information.

* `data`: stored data.
* `parent`: parent node.
* `left`: left child node.
* `right`: right child node.
"""
mutable struct SimpleBinaryNode{T} <: AbstractBinaryNode{T}
    data::T
    parent::Union{SimpleBinaryNode{T}, NullNode}
    left::Union{SimpleBinaryNode{T}, NullNode}
    right::Union{SimpleBinaryNode{T}, NullNode}

    # Initial constructor
    SimpleBinaryNode{T}(data) where T = new{T}(convert(T, data), NullNode(), NullNode(), NullNode())
end

const SBN = SimpleBinaryNode

"""
    HeightBinaryNode{T} <: AbstractBinaryNode{T}
    const HBN = HeightBinaryNode
    HBN(data)
    HBN(::Type{T}, data)

Binary nodes with height information.

* `height`: if this node is a leaf of tree, it is 0; otherwise, it is the number of nodes between this node and the deepest node (not including itself).

For other fields, please see `SimpleBinaryNode`.
"""
mutable struct HeightBinaryNode{T} <: AbstractBinaryNode{T}
    data::T
    parent::Union{HeightBinaryNode{T}, NullNode}
    left::Union{HeightBinaryNode{T}, NullNode}
    right::Union{HeightBinaryNode{T}, NullNode}
    height::Int

    # Initial constructor
    HeightBinaryNode{T}(data) where T = new{T}(convert(T, data), NullNode(), NullNode(), NullNode(), 0)
end

const HBN = HeightBinaryNode

"""
    RedBlackBinaryNode{T} <: AbstractBinaryNode{T}
    const RBN = RedBlackBinaryNode
    RBN(data)
    RBN(::Type{T}, data)

Binary nodes for red black tree.

* `red`: true for red, false for black.

For other fields, please see `SimpleBinaryNode`.
"""
mutable struct RedBlackBinaryNode{T} <: AbstractBinaryNode{T}
    data::T
    parent::Union{RedBlackBinaryNode{T}, NullNode}
    left::Union{RedBlackBinaryNode{T}, NullNode}
    right::Union{RedBlackBinaryNode{T}, NullNode}
    red::Bool # true for red, false for black

    # Initial constructor
    RedBlackBinaryNode{T}(data) where T = new{T}(convert(T, data), NullNode(), NullNode(), NullNode(), false)
end

const RBN = RedBlackBinaryNode

const NSBN = Union{NullNode, SimpleBinaryNode}

for Node in [:SBN, :HBN, :RBN]
    @eval begin
        $Node(data) = $Node{typeof(data)}(data)
        $Node(::Type{T}, data) where T = $Node{T}(data)
    end
end

# ------------------------------------------------------------------------------------
# Trees

"""
    BinarySearchTree <: AbstractTree

Binary search tree. Nodes can be `SimpleBinaryNode` or `HeightBinaryNode`.

# Constructors
    const BST = BinarySearchTree
    BST(; height::Bool = false)
    BST(datas...; height::Bool = false)
    BST(::Type{T}; height::Bool = false)
    BST(::Type{T}, datas...; height::Bool = false)

# Fields
* `root`: the root node.
* `size`: number of nodes in this tree.
* `height`: the height of root. Only valid when using `HeightBinaryNode`.
"""
mutable struct BinarySearchTree{N <: AbstractBinaryNode} <: AbstractTree{N}
    root::Union{N, NullNode}
    size::Int
    height::Int # Only valid when using HBN

    # Initial constructor
    BinarySearchTree{N}() where N = new{N}(NullNode(), 0, -1)
    # W/ data constructor
    BinarySearchTree{N}(data) where N = new{N}(N(data), 1, 0)
end

const BST = BinarySearchTree

#BST(data; height::Bool = false) = BST(typeof(data), data; height)

BST(datas...; height::Bool = false) = BST(promote(datas...)...; height)
BST(; height::Bool = false) = BST(Any; height)

function BST(datas::Vararg{T}; height::Bool = false) where T
    data = first(datas)
    tree = BST(T, data; height)
    insert!(tree, datas[2:end]...)
    tree
end

BST(::Type{T}; height::Bool = false) where T = height ? BST{HBN{T}}() : BST{SBN{T}}()
BST(::Type{T}, data; height::Bool = false) where T = height ? BST{HBN{T}}(data) : BST{SBN{T}}(data)

function BST(::Type{T}, datas...; height::Bool = false) where T
    data = first(datas)
    tree = BST(T, data; height)
    insert!(tree, datas[2:end]...)
    tree
end

"""
    AVLTree <: AbstractTree

AVL tree. Nodes must be `HeightBinaryNode` to be able to do rotations.

# Constructors
    const AVL = AVLTree
    AVL()
    AVL(datas...)
    AVL(::Type{T})
    AVL(::Type{T}), datas...)

# Fields
* `root`: the root node.
* `size`: number of nodes in this tree.
* `height`: the height of root.
"""
mutable struct AVLTree{N <: HBN} <: AbstractTree{N}
    root::Union{N, NullNode}
    size::Int
    height::Int

    # Initial constructor
    AVLTree{N}() where N = new{N}(NullNode(), 0, -1)
    # W/ data constructor
    AVLTree{N}(data) where N = new{N}(N(data), 1, 0)
end

const AVL = AVLTree

"""
    SplayTree <: AbstractTree

Splay tree. Nodes should be `SimpleBinaryNode`, though `HeightBinaryNode` is ok in theory.

# Constructors
    const Splay = SplayTree
    Splay()
    Splay(datas...)
    Splay(::Type{T})
    Splay(::Type{T}, datas...)

# Fields
* `root`: the root node.
* `size`: number of nodes in this tree.
"""
mutable struct SplayTree{N <: SBN} <: AbstractTree{N}
    root::Union{N, NullNode}
    size::Int

    # Initial constructor
    SplayTree{N}() where N = new{N}(NullNode(), 0)
    # W/ data constructor
    SplayTree{N}(data) where N  = new{N}(N(data), 1)
end

const Splay = SplayTree

"""
    RedBlackTree <: AbstractTree

Red black tree. Nodes should be `RedBlackBinaryNode`.

# Constructors
    const RBT = RedBlackTree
    RBT()
    RBT(datas...)
    RBT(::Type{T})
    RBT(::Type{T}, datas...)

# Fields
* `root`: the root node.
* `size`: number of nodes in this tree.
"""
mutable struct RedBlackTree{N <: RBN} <: AbstractTree{N}
    root::Union{N, NullNode}
    size::Int

    # Initial constructor
    RedBlackTree{N}() where N = new{N}(NullNode(), 0)
    # W/ data constructor
    RedBlackTree{N}(data) where N = new{N}(N(data), 1)
end

const RBT = RedBlackTree

# ------------------------------------------------------------------------------------------

canonicalnode(::Type{<: AbstractTree}) = SBN
canonicalnode(::Type{<: AVL}) = HBN
canonicalnode(::Type{<: RBT}) = RBN

for Tree in [:AVL, :Splay, :RBT]
    @eval begin
        #$Tree(data) = $Tree{canonicalnode($Tree){typeof(data)}}(data)

        $Tree(datas...) = $Tree(promote(datas...)...)
        $Tree() = $Tree{canonicalnode($Tree){Any}}()

        function $Tree(datas::Vararg{T}) where T
            data = first(datas)
            tree = $Tree(T, data)
            insert!(tree, datas[2:end]...)
            tree
        end

        $Tree(::Type{T}) where T = $Tree{canonicalnode($Tree){T}}()
        $Tree(::Type{T}, data) where T = $Tree{canonicalnode($Tree){T}}(data)

        function $Tree(::Type{T}, datas...) where T
            data = first(datas)
            tree = $Tree(T, data)
            insert!(tree, datas[2:end]...)
            tree
        end
    end
end


text:
export SymbolicAmplitudeAwarePermutation

"""
    SymbolicAmplitudeAwarePermutation(; τ = 1, m = 3, A = 0.5, lt = Entropies.isless_rand) <: PermutationProbabilityEstimator

See docstring for [`SymbolicPermutation`](@ref).
"""
struct SymbolicAmplitudeAwarePermutation{F} <: PermutationProbabilityEstimator
    τ::Int
    m::Int
    A::Float64
    lt::F
end
function SymbolicAmplitudeAwarePermutation(; τ::Int = 1, m::Int = 2, A::Real = 0.5,
        lt::F = isless_rand) where {F <: Function}
    2 ≤ m || error("Need m ≥ 2, otherwise no dynamical information is encoded in the symbols.")
    0 ≤ A ≤ 1 || error("Weighting factor A must be on interval [0, 1]. Got A=$A.")
    SymbolicAmplitudeAwarePermutation{F}(τ, m, A, lt)
end

"""
    AAPE(x, A::Real = 0.5, m::Int = length(a))

Encode relative amplitude information of the elements of `a`.
- `A = 1` emphasizes only average values.
- `A = 0` emphasizes changes in amplitude values.
- `A = 0.5` equally emphasizes average values and changes in the amplitude values.
"""
function AAPE(x; A::Real = 0.5, m::Int = length(x))
    (A/m)*sum(abs.(x)) + (1-A)/(m-1)*sum(abs.(diff(x)))
end

function probabilities(x::AbstractDataset{m, T}, est::SymbolicAmplitudeAwarePermutation) where {m, T}
    πs = symbolize(x, SymbolicPermutation(m = m, lt = est.lt)) # motif length controlled by dimension of input data
    wts = AAPE.(x.data, A = est.A, m = est.m)

    Probabilities(symprobs(πs, wts, normalize = true))
end

function probabilities(x::AbstractVector{T}, est::SymbolicAmplitudeAwarePermutation) where {T<:Real}
    τs = tuple([est.τ*i for i = 0:est.m-1]...)
    emb = genembed(x, τs)
    πs = symbolize(emb, SymbolicPermutation(m = est.m, lt = est.lt))  # motif length controlled by estimator m
    wts = AAPE.(emb.data, A = est.A, m = est.m)
    p = symprobs(πs, wts, normalize = true)
    Probabilities(p)
end

ext: jl


text:
module Errors

import Base: show_method_candidates, ExceptionStack, InterpreterIP

import Term: highlight, truncate, reshape_text, load_code_and_highlight

import ..Layout:
    hLine, rvstack, cvstack, rvstack, vstack, vLine, Spacer, hstack, lvstack, pad
import ..Renderables: RenderableText, AbstractRenderable
import ..Panels: Panel

export install_term_stacktrace

include("_errors.jl")

const ErrorsExplanations = Dict(
    ArgumentError => "The parameters to a function call do not match a valid signature.",
    AssertionError => "comes up when an assertion's check fails (e.g., `@assert 1==2`)",
    BoundsError => "comes up when trying to acces a container at invalid position (e.g., a string a='abcd' with 4 characters cannot be accessed as a[5]).",
    DimensionMismatch => "comes up when trying to perform an operation on objects which don't have matching dimensionality (e.g., summing matrixes of different size).",
    DivideError => "comes up when attempting integer division with 0 as denominator. {blue}2/0=Inf{/blue} is okay, but {orange1}div(2, ){/orange1} will give an error",
    DomainError => "comes up when the argument to a function is outside its domain (e.g., √(-1))",
    ErrorException => "is a generic error type",
    KeyError => "comes up when attempting to access a non-existing {blue}Dict{/blue} key.",
    InexactError => "comes up when a type cannot exactly be converted to another (e.g. Int(2.5) cannot convert Float64 to Int64, but Int(round(2.5)) will work)",
    LoadError => "occurs when another comes up while evaluating 'include', 'require' or 'using' statements",
    MethodError => "comes up when to method can be found with a given name and for a given set of argument types.",
    StackOverflowError => "usually comes up when functions call each other recursively.",
    TypeError => "is a type assertion failure, or calling an intrinsic function with an incorrect argument type.",
    UndefKeywordError => "comes up when a function has a keyword argument with no default value and no value is passed to a function call",
    UndefVarError => "comes up when a variable is used which is either not defined, or, which is not visible in the current variables scope (e.g.: variable defined in function A and used in function B)",
)

_width() = min(console_width(stderr), 100)

# ----------------------- error type specific messages ----------------------- #

# ! ARGUMENT ERROR
error_message(er::ArgumentError) = er.msg, ""

# ! ASSERTION ERROR
error_message(er::AssertionError) = return er.msg, ""

# ! BOUNDS ERROR
function error_message(er::BoundsError)
    # @info "bounds error" er fieldnames(typeof(er))
    main_msg = "Attempted to access `$(er.a)` at index $(er.i)"

    additional_msg = ""

    if isdefined(er, :a)
        if er.a isa AbstractString
            nunits = ncodeunits(er.a)
            additional_msg = "S\ntring has $nunits codeunits, $(length(er.a)) characters."
        end
    else
        additional_msg = "\n{red}Variable is not defined!.{/red}"
    end
    return main_msg, additional_msg
end

# ! Domain ERROR
function error_message(er::DomainError)
    # @info "err exceprion" er fieldnames(DomainError) er.val
    # msg = split(er.msg, " around ")[1]
    return er.msg, "\nThe invalid value is: $(er.val)."
end

# ! DimensionMismatch
error_message(er::DimensionMismatch) = er.msg, ""

# ! DivideError
error_message(er::DivideError) = "Attempted integer division by {bold}0{/bold}", ""

# ! EXCEPTION ERROR
function error_message(er::ErrorException)
    # @info "err exceprion" er fieldnames(ErrorException) er.msg
    msg = split(er.msg, " around ")[1]
    return msg, ""
end

# !  KeyError
function error_message(er::KeyError)
    # @info "err KeyError" er fieldnames(KeyError)
    msg = "Key `$(er.key)` not found!"
    return msg, ""
end

# ! InexactError
function error_message(er::InexactError)
    # @info "load error message"  fieldnames(InexactError)
    msg = "Cannot convert $(er.val) to type ::$(er.T)"
    subm = "\nConversion error in function: $(er.func)"
    return msg, subm
end

# ! LoadError
function error_message(er::LoadError)
    # @info "load error message"  fieldnames(LoadError)
    msg = "At {grey62 underline}$(er.file){/grey62 underline} line {bold}$(er.line){/bold}"
    subm = "The cause is an error of type: {bright_red}$(string(typeof(er.error)))"
    return msg, subm
end

# ! METHOD ERROR
method_error_regex = r"(?<group>\!Matched\:\:(\w|\.)+)"
function method_error_candidate(fun, candidate)
    # highlight non-matched types
    candidate = replace(
        candidate,
        method_error_regex => SubstitutionString("{red}" * s"\g<0>" * "{/red}"),
    )
    # remove
    candidate = replace(candidate, "!Matched" => "")

    # highlight fun
    candidate = replace(candidate, string(fun) => "{bold yellow}$(fun){/bold yellow}")
    return candidate
end

function error_message(er::MethodError; kwargs...)
    f = er.f
    ft = typeof(f)
    name = ft.name.mt.name
    kwargs = ()
    if endswith(string(ft.name.name), "##kw")
        f = er.args[2]
        ft = typeof(f)
        name = ft.name.mt.name
        arg_types_param = arg_types_param[3:end]
        kwargs = pairs(er.args[1])
        # er = MethodError(f, er.args[3:end::Int])
    end

    # get main error message
    _args = join(
        map(
            a ->
                "   {dim bold}($(a[1])){/dim bold} $(truncate(highlight("::"*string(typeof(a[2]))), 30))",
            enumerate(er.args),
        ),
        "\n",
    )
    main_line = "No method matching `$name` with arguments types:" / _args

    # get recomended candidates
    _candidates = split(sprint(show_method_candidates, er), "\n")[3:(end - 1)]

    if length(_candidates) > 0
        _candidates = map(c -> split(c, " at ")[1], _candidates)
        candidates = map(c -> method_error_candidate(name, c), _candidates)
        main_line =
            main_line /
            lvstack("" / "Alternative candidates:", RenderableText.(candidates)...;)
    else
        main_line = main_line / " " / "{dim}No alternative candidates found"
    end
    return string(main_line), ""
end

# ! StackOverflowError
error_message(er::StackOverflowError) = "Stack overflow error: too many function calls.", ""

# ! TYPE ERROR
function error_message(er::TypeError)
    # @info "type err" er fieldnames(typeof(er)) er.func er.context er.expected er.got
    # var = string(er.var)
    msg = "In `$(er.func)` > `$(er.context)` got"
    msg *= " {orange1 bold}$(er.got){/orange1 bold}(::$(typeof(er.got))) but expected argument of type ::$(er.expected)"
    return msg, ""
end

# ! UndefKeywordError
function error_message(er::UndefKeywordError)
    # @info "UndefKeywordError" er er.var typeof(er.var) fieldnames(typeof(er.var))
    var = string(er.var)
    return "Undefined function keyword argument: `$(er.var)`.", ""
end

# ! UNDEFVAR ERROR
function error_message(er::UndefVarError)
    # @info "undef var error" er er.var typeof(er.var)
    var = string(er.var)
    return "Undefined variable `$(er.var)`.", ""
end

# ! STRING INDEX ERROR
function error_message(er::StringIndexError)
    # @info er typeof(er) fieldnames(typeof(er))
    m1 = "attempted to access a String at index $(er.index)\n"
    return m1, ""
end

# ! catch all other errors
function error_message(er)
    # @debug "Error message type doesnt have a specialized method!" er typeof(er) fieldnames(
    #     typeof(er)
    # )
    if hasfield(typeof(er), :error)
        # @info "nested error" typeof(er.error)
        m1, m2 = error_message(er.error)
        msg = "\n{bold red}LoadError:{/bold red}\n" * m1
    else
        msg = if hasfield(typeof(er), :msg)
            er.msg
        else
            "no message for error of type $(typeof(er)), sorry."
        end
        m2 = ""
    end
    return msg, m2
end

# ---------------------------------------------------------------------------- #
#                              INSTALL STACKTRACE                              #
# ---------------------------------------------------------------------------- #
function install_term_stacktrace(; reverse_backtrace::Bool = true, max_n_frames::Int = 30)
    @eval begin
        function Base.showerror(io::IO, er, bt; backtrace = true)
            (length(bt) == 0 && !isa(er, StackOverflowError)) && return nothing
            try
                println("\n")
                ename = string(typeof(er))
                error =
                    hLine("{default bold red}$ename{/default bold red}"; style = "dim red")
                if length(bt) > 0
                    rendered_bt = render_backtrace(
                        bt;
                        reverse_backtrace = $(reverse_backtrace),
                        max_n_frames = $(max_n_frames),
                    )
                    error /= rendered_bt
                    W = rendered_bt.measure.w
                else
                    W = 88
                end
                err, _ = error_message(er)
                msg =
                    "" / Panel(
                        "{#aec2e8}$(err){/#aec2e8}";
                        width = W,
                        title = "{bold red default underline}$(typeof(er)){/bold red default underline}",
                        padding = (2, 2, 1, 1),
                        style = "dim red",
                        title_justify = :center,
                        # fit = true,
                        # width=W
                    )
                error /= msg
                print(error)
            catch err
                @error "ERROR: " exception = err
                @warn "Term.jl: failed to render error message" err
                Base.showerror(io, er)
            end
        end
    end
end

end


text:
addprocs(3)
@everywhere using NFP, DataFrames
@everywhere const sStart_s = "01/01/15" # start out of sample
@everywhere const iSymbol = :NFP # dependent variable
@everywhere const vSymbol = [:Date, :NFP, :NFP_bb_median] # remove from dataset (non-numerical and dep. var.)
@everywhere const H = [1,2] # horizons
@everywhere const iBest = 3 # Check the best x variables among the two criteria
@everywhere const ncomb_load = 3 # max comb to use for the forecast

@everywhere mae(vX::Vector,vY::Vector) = mean(abs.(vX-vY))      ## MAE loss function
@everywhere rmse(vX::Vector,vY::Array) = sqrt(mean((vX-vY).^2)) ## RMSE loss function

@everywhere const fLoss = [mae rmse] ## Vector of loss functions

@everywhere dfData = readtable(joinpath(Pkg.dir("NFP"),"test","data.csv"), header = true)
@everywhere const iStart = find(dfData[:Date] .== sStart_s)[1]

# For NFP level forecast
l_plot,r = sforecast(dfData,vSymbol,iSymbol,H,iStart,iBest,ncomb_load,fLoss)
l_plot,r= fforecast(dfData,vSymbol,iSymbol,H,iStart,iBest,ncomb_load,fLoss)

# See the forecasts plot
l_plot

rmprocs(2:4)


text:
using GLM

@testset "L0_reg normal" begin
	# Since my code seems to work, putting in some output as they can be verified by comparing with simulation

	#simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = Normal
	l = canonicallink(d())

	#set random seed
	Random.seed!(1111)

	#construct SnpArraym, snpmatrix, and non genetic covariate (intercept)
	x = simulate_random_snparray(n, p, "test1.bed")
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	z = ones(n, 1)

	# simulate response, true model b, and the correct non-0 positions of b
	y, true_b, correct_position = simulate_random_response(x, xbm, k, d, l)

	#run result
	result = L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false, init=false, use_maf=false)

	@test length(result.beta) == 10000
	@test findall(!iszero, result.beta) == [2384;3352;3353;4093;5413;5609;7403;8753;9089;9132]
	@test all(result.beta[findall(!iszero, result.beta)] .≈ [-1.2601406011046452;
	 				-0.2674202492177914; 0.14120810664715883; 0.289955803600036;
	  				0.3666894767520663; -0.1371805027382694; -0.3082545756160329;
	  				0.3328814701200445; 0.9645980728400257; -0.5094607091364866])
	@test result.c[1] == 0.0
	@test result.k == 10
	@test result.logl ≈ -1407.2533232402275
end

@testset "L0_reg Bernoulli" begin
	# Since my code seems to work, putting in some output as they can be verified by comparing with simulation

	#simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = Bernoulli
	l = canonicallink(d())

	#set random seed
	Random.seed!(1111)

	#construct SnpArraym, snpmatrix, and non genetic covariate (intercept)
	x = simulate_random_snparray(n, p, "test1.bed")
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	z = ones(n, 1)

	# simulate response, true model b, and the correct non-0 positions of b
	y, true_b, correct_position = simulate_random_response(x, xbm, k, d, l)

	#run result
	result = L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false, init=false, use_maf=false)

	@test length(result.beta) == 10000
	@test findall(!iszero, result.beta) == [1733;1816;2384;5413;7067;8753;8908;9089;9132;9765]
	@test all(result.beta[findall(!iszero, result.beta)] .≈ [-0.2787326116508012;
					  0.3113511410050774;-1.1292096054341005;0.5001816459301949;
					 -0.32694130827328116;0.4134742776599116;-0.3275424847038566;
					  0.8619785898062307;-0.5068258295825918;-0.32972421733995294])
	@test result.c[1] == 0.0
	@test result.k == 10
	@test result.logl ≈ -489.8770526620568
end

@testset "L0_reg Poisson" begin
	# Since my code seems to work, putting in some output as they can be verified by comparing with simulation

	#simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = Poisson
	l = canonicallink(d())

	#set random seed
	Random.seed!(1111)

	#construct SnpArraym, snpmatrix, and non genetic covariate (intercept)
	x = simulate_random_snparray(n, p, "test1.bed")
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	z = ones(n, 1)

	# simulate response, true model b, and the correct non-0 positions of b
	y, true_b, correct_position = simulate_random_response(x, xbm, k, d, l)

	#run result
	result = L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false, init=false, use_maf=false)
	@test length(result.beta) == 10000
	@test findall(!iszero, result.beta) == [298; 606; 2384; 5891; 7067; 8753; 8755; 8931; 9089; 9132]
	@test all(result.beta[findall(!iszero, result.beta)] .≈ [0.10999211487301704;
		-0.09628969787009399; -0.3660582504298778;  0.11767397809862554;
		0.09686501699067837;  0.11419451741236888;  0.12373749347128933;
		0.11916107757737655;  0.2904980599350941; -0.12920302008477738])
	@test result.c[1] == 0.0
	@test result.k == 10
	@test result.logl ≈ -1293.4456256102478
	@test result.iter == 14
end

@testset "L0_reg NegativeBinomial" begin
	# Since my code seems to work, putting in some output as they can be verified by comparing with simulation

	#simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = NegativeBinomial
	l = LogLink()

	#set random seed
	Random.seed!(1111)

	#construct SnpArraym, snpmatrix, and non genetic covariate (intercept)
	x = simulate_random_snparray(n, p, "test1.bed")
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	z = ones(n, 1)

	# simulate response, true model b, and the correct non-0 positions of b
	y, true_b, correct_position = simulate_random_response(x, xbm, k, d, l)

	#run result
	result = L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false, init=false, use_maf=false)

	@test length(result.beta) == 10000
	@test findall(!iszero, result.beta) == [1245; 1774; 1982; 2384; 5413; 5440; 5614; 7166; 9089; 9132;]
	@test all(result.beta[findall(!iszero, result.beta)] .≈ [-0.13251205076442696;
		-0.16906821875893546;  0.12865324984770152; -0.30791709019019947;
		0.1202449253579259; 0.12545318690748591;  0.10252799982767402;
		0.12937785947321034;  0.27113364351607033;-0.19797373419860373])
	@test result.c[1] == 0.0
	@test result.k == 10
	@test result.logl ≈ -1387.341396480908
	@test result.iter == 9
end

@testset "L0_reg with non-genetic covariates" begin
	# Since my code seems to work, putting in some output as they can be verified by comparing with simulation

	#simulat data with k true predictors, from distribution d and with link l.
	n = 1000
	p = 10000
	k = 10
	d = Normal
	l = canonicallink(d())

	#set random seed
	Random.seed!(1111)

	#construct SnpArraym, snpmatrix, and non genetic covariate (intercept)
	x = simulate_random_snparray(n, p, "test1.bed")
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	z = ones(n, 2) # the intercept
	z[:, 2] .= randn(n)

	#define true_b and true_c
	true_b = zeros(p)
	true_b[1:k-2] = randn(k-2)
	shuffle!(true_b)
	correct_position = findall(!iszero, true_b)
	true_c = [3.0; 3.5]

	#simulate phenotype
	prob = GLM.linkinv.(l, xbm * true_b .+ z * true_c)
	y = [rand(d(i)) for i in prob]

	#run result
	result = L0_reg(x, xbm, z, y, 1, k, d(), l, debias=false, init=false, use_maf=false)

	@test length(result.beta) == 10000
	@test length(result.c) == 2
	@test findall(!iszero, result.beta) == [2984;4147;4604;6105;6636;7575;8271;9300]
	@test findall(!iszero, result.c) == [1;2]
	@test all(result.beta[findall(!iszero, result.beta)] .≈ [  1.0813690725698488;
 		-0.21167725381808253; -0.3320058055771102; 0.40296120706642935;
 		-0.1320095894527104;  1.6694054572246806;  0.3151337065844457;
 		-1.5967591054279828])
	@test all(result.c .≈ [2.931093515105184; 3.4674609432737893])
	@test result.k == 10
	@test result.logl ≈ -1372.8963328607294

	#clean up
	rm("test1.bed", force=true)
end

@testset "L0_reg with correlated predictors and double sparsity" begin
	# Since my code seems to work, putting in some output as they can be verified by comparing with simulation

	#simulat data with k true predictors, from distribution d and with link l.
    n = 1000
    p = 10000
    d = Normal
	l = canonicallink(d())
    block_size = 20
    num_blocks = Int(p / block_size)

	#set random seed
	Random.seed!(1111)

    # assign group membership
    membership = collect(1:num_blocks)
    g = zeros(Int64, p + 1)
    for i in 1:length(membership)
        for j in 1:block_size
            cur_row = block_size * (i - 1) + j
            g[block_size*(i - 1) + j] = membership[i]
        end
    end
    g[end] = membership[end]

    #simulate correlated snparray
    x = simulate_correlated_snparray(n, p, "tmp.bed")
    z = ones(n, 1) # the intercept
    x_float = convert(Matrix{Float64}, x, model=ADDITIVE_MODEL, center=true, scale=true)

    #simulate true model, where 5 groups each with 3 snps contribute
    true_b = zeros(p)
    true_groups = randperm(num_blocks)[1:5]
    within_group = [randperm(block_size)[1:3] randperm(block_size)[1:3] randperm(block_size)[1:3] randperm(block_size)[1:3] randperm(block_size)[1:3]]
    correct_position = zeros(Int64, 15)
    for i in 1:5
        cur_group = block_size * (true_groups[i] - 1)
        cur_group_snps = cur_group .+ within_group[:, i]
        correct_position[3*(i-1)+1:3i] .= cur_group_snps
    end
    for i in 1:15
        true_b[correct_position[i]] = rand(-1:2:1) * 0.2
    end
    sort!(correct_position)

    # simulate phenotype
    if d == Normal || d == Bernoulli || d == Poisson
        prob = GLM.linkinv.(l, x_float * true_b)
        clamp!(prob, -20, 20)
        y = [rand(d(i)) for i in prob]
    elseif d == NegativeBinomial
        nn = 10
        μ = GLM.linkinv.(l, x_float * true_b)
        clamp!(μ, -20, 20)
        prob = 1 ./ (1 .+ μ ./ nn)
        y = [rand(d(nn, i)) for i in prob] #number of failtures before nn success occurs
    end
    y = Float64.(y)

    #run IHT without groups
    k = 15
    ungrouped = L0_reg(x_float, z, y, 1, k, d(), l, debias=false)

    #run IHT with groups
    J = 5
    k = 3
    grouped = L0_reg(x_float, z, y, J, k, d(), l, debias=false, group=g)

    #clean up
    rm("tmp.bed", force=true)

    @test ungrouped.iter == 36
    @test ungrouped.logl ≈ -1339.457450752207
    @test all(findall(!iszero, ungrouped.beta) .== [ 3674; 5969; 5975; 7006; 7015;
    					7020; 7764; 7766; 7769; 8423; 8424; 8429; 9321; 9323; 9337])
    @test all(ungrouped.beta[findall(!iszero, ungrouped.beta)] .≈ [
		 -0.11410368504207032;
		 -0.12255816354649987;
		  0.20533252358870932;
		  0.18063333820398803;
		 -0.17263255338849545;
		 -0.28933626047953254;
		  0.20523417123190285;
		 -0.2138908635675793;
		 -0.18086336715675141;
		  0.17944252853814854;
		  0.17589384672160724;
		  0.2303348657532707;
		  0.27590559980983664;
		 -0.17119356771598235;
		 -0.12463707920692581])

    @test grouped.iter == 31
    @test grouped.logl ≈ -1339.4622452299625
    @test all(findall(!iszero, grouped.beta) .== [5964;5969;5975;7006;7015;7020;
    						7764;7766; 7769; 8423; 8424; 8428; 9321; 9326; 9338])
    @test all(grouped.beta[findall(!iszero, grouped.beta)] .≈ [
    	 -0.1211609170236838
		 -0.14680300260861762
		  0.1662368258482021
		  0.18915447812538866
		 -0.17408389637789895
		 -0.2829678623871821
		  0.19735076031573381
		 -0.21902842957653598
		 -0.16680170315905823
		  0.12147615146809417
		  0.1842115499470037
		  0.23792933900511254
		  0.25490356768515327
		 -0.1815202361606688
		 -0.25255175830728194])
end


text:
module PictureSegmentation

using Images
using Makie
using ColorVectorSpace

abstract type SegmentationAlgorithm end
struct GrowCut <: SegmentationAlgorithm end

include("common.jl")
include("seedpixel.jl")
include("growcut.jl")
include("growcut3d.jl")
include("display.jl")

export
	# main functions
	segment_image,
	set_seed_pixels!,
	GrowCut,
	display_3D,
	help_rule_template,
	create_rule,
	convert_labels,
	make_contour
end # module


text:
module CentroidBased
    using Pluto

    run_notebook_k_means_pp() = Pluto.run(notebook=joinpath(@__DIR__, "../notebook/k_means_pp.jl"))
    run_notebook_k_medoids() = Pluto.run(notebook=joinpath(@__DIR__, "../notebook/k_medoids.jl"))
    run_notebook_medoid_median_mode() = Pluto.run(notebook=joinpath(@__DIR__, "../notebook/medoid_median_mode.jl"))
end


text:
int_rules_7_3_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 7.3.3*Exponentials*of*inverse*hyperbolic*tangent =#
    @apply_utils Antiderivative(ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => Antiderivative((1 + ~a * ~x) ^ ((~n + 1) / 2) / ((1 - ~a * ~x) ^ ((~n - 1) / 2) * sqrt(1 - (~a) ^ 2 * (~x) ^ 2)), ~x) <-- FreeQ(~a, ~x) && IntegerQ((~n - 1) / 2)
    @apply_utils Antiderivative((~x) ^ ~(m') * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => Antiderivative((~x) ^ ~m * ((1 + ~a * ~x) ^ ((~n + 1) / 2) / ((1 - ~a * ~x) ^ ((~n - 1) / 2) * sqrt(1 - (~a) ^ 2 * (~x) ^ 2))), ~x) <-- FreeQ([~a, ~m], ~x) && IntegerQ((~n - 1) / 2)
    @apply_utils Antiderivative(ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => Antiderivative((1 + ~a * ~x) ^ (~n / 2) / (1 - ~a * ~x) ^ (~n / 2), ~x) <-- FreeQ([~a, ~n], ~x) && Not(IntegerQ((~n - 1) / 2))
    @apply_utils Antiderivative((~x) ^ ~(m') * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => Antiderivative(((~x) ^ ~m * (1 + ~a * ~x) ^ (~n / 2)) / (1 - ~a * ~x) ^ (~n / 2), ~x) <-- FreeQ([~a, ~m, ~n], ~x) && Not(IntegerQ((~n - 1) / 2))
    @apply_utils Antiderivative((~c + ~(d') * ~x) ^ ~(p') * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (~c) ^ ~n * Antiderivative((~c + ~d * ~x) ^ (~p - ~n) * (1 - (~a) ^ 2 * (~x) ^ 2) ^ (~n / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ(~a * ~c + ~d, 0) && (IntegerQ((~n - 1) / 2) && IntegerQ(2 * ~p)))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~c + ~(d') * ~x) ^ ~(p') * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (~c) ^ ~n * Antiderivative((~e + ~f * ~x) ^ ~m * (~c + ~d * ~x) ^ (~p - ~n) * (1 - (~a) ^ 2 * (~x) ^ 2) ^ (~n / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~m, ~p], ~x) && (EqQ(~a * ~c + ~d, 0) && (IntegerQ((~n - 1) / 2) && ((IntegerQ(~p) || (EqQ(~p, ~n / 2) || EqQ((~p - ~n / 2) - 1, 0))) && IntegerQ(2 * ~p))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x) ^ ~(p') * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative((~u * (1 + (~d * ~x) / ~c) ^ ~p * (1 + ~a * ~x) ^ (~n / 2)) / (1 - ~a * ~x) ^ (~n / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~a) ^ 2 * (~c) ^ 2 - (~d) ^ 2, 0) && (IntegerQ(~p) || GtQ(~c, 0)))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x) ^ ~(p') * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => Antiderivative((~u * (~c + ~d * ~x) ^ ~p * (1 + ~a * ~x) ^ (~n / 2)) / (1 - ~a * ~x) ^ (~n / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~a) ^ 2 * (~c) ^ 2 - (~d) ^ 2, 0) && Not(IntegerQ(~p) || GtQ(~c, 0)))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / ~x) ^ ~(p') * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (~d) ^ ~p * Antiderivative((~u * (1 + (~c * ~x) / ~d) ^ ~p * ℯ ^ (~n * atanh(~a * ~x))) / (~x) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~c) ^ 2 - (~a) ^ 2 * (~d) ^ 2, 0) && IntegerQ(~p))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / ~x) ^ ~p * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => (-1) ^ (~n / 2) * (~c) ^ ~p * Antiderivative((~u * (1 + ~d / (~c * ~x)) ^ ~p * (1 + 1 / (~a * ~x)) ^ (~n / 2)) / (1 - 1 / (~a * ~x)) ^ (~n / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ((~c) ^ 2 - (~a) ^ 2 * (~d) ^ 2, 0) && (Not(IntegerQ(~p)) && (IntegerQ(~n / 2) && GtQ(~c, 0))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / ~x) ^ ~p * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => Antiderivative((~u * (~c + ~d / ~x) ^ ~p * (1 + ~a * ~x) ^ (~n / 2)) / (1 - ~a * ~x) ^ (~n / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ((~c) ^ 2 - (~a) ^ 2 * (~d) ^ 2, 0) && (Not(IntegerQ(~p)) && (IntegerQ(~n / 2) && Not(GtQ(~c, 0)))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / ~x) ^ ~p * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (((~x) ^ ~p * (~c + ~d / ~x) ^ ~p) / (1 + (~c * ~x) / ~d) ^ ~p) * Antiderivative((~u * (1 + (~c * ~x) / ~d) ^ ~p * ℯ ^ (~n * atanh(~a * ~x))) / (~x) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~c) ^ 2 - (~a) ^ 2 * (~d) ^ 2, 0) && Not(IntegerQ(~p)))
    @apply_utils Antiderivative(ℯ ^ (~n * atanh(~(a') * ~x)) / (~c + ~(d') * (~x) ^ 2) ^ (3 / 2), ~x) => ((~n - ~a * ~x) * ℯ ^ (~n * atanh(~a * ~x))) / (~a * ~c * ((~n) ^ 2 - 1) * sqrt(~c + ~d * (~x) ^ 2)) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => ((~n + 2 * ~a * (~p + 1) * ~x) * (~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * atanh(~a * ~x))) / (~a * ~c * ((~n) ^ 2 - 4 * (~p + 1) ^ 2)) - ((2 * (~p + 1) * (2 * ~p + 3)) / (~c * ((~n) ^ 2 - 4 * (~p + 1) ^ 2))) * Antiderivative((~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * atanh(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (LtQ(~p, -1) && (Not(IntegerQ(~n)) && (NeQ((~n) ^ 2 - 4 * (~p + 1) ^ 2, 0) && IntegerQ(2 * ~p)))))
    @apply_utils Antiderivative(ℯ ^ (~(n') * atanh(~(a') * ~x)) / (~c + ~(d') * (~x) ^ 2), ~x) => ℯ ^ (~n * atanh(~a * ~x)) / (~a * ~c * ~n) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && Not(IntegerQ(~n / 2)))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative((1 - (~a) ^ 2 * (~x) ^ 2) ^ (~p - ~n / 2) * (1 + ~a * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (IntegerQ(~p) && (IGtQ((~n + 1) / 2, 0) && Not(IntegerQ(~p - ~n / 2)))))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative((1 - (~a) ^ 2 * (~x) ^ 2) ^ (~p + ~n / 2) / (1 - ~a * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (IntegerQ(~p) && (ILtQ((~n - 1) / 2, 0) && Not(IntegerQ(~p - ~n / 2)))))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative((1 - ~a * ~x) ^ (~p - ~n / 2) * (1 + ~a * ~x) ^ (~p + ~n / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (IntegerQ(~p) || GtQ(~c, 0)))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => (~c) ^ (~n / 2) * Antiderivative((~c + ~d * (~x) ^ 2) ^ (~p - ~n / 2) * (1 + ~a * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (Not(IntegerQ(~p) || GtQ(~c, 0)) && IGtQ(~n / 2, 0)))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => (1 / (~c) ^ (~n / 2)) * Antiderivative((~c + ~d * (~x) ^ 2) ^ (~p + ~n / 2) / (1 - ~a * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (Not(IntegerQ(~p) || GtQ(~c, 0)) && ILtQ(~n / 2, 0)))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (((~c) ^ IntPart(~p) * (~c + ~d * (~x) ^ 2) ^ FracPart(~p)) / (1 - (~a) ^ 2 * (~x) ^ 2) ^ FracPart(~p)) * Antiderivative((1 - (~a) ^ 2 * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atanh(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && Not(IntegerQ(~p) || GtQ(~c, 0)))
    @apply_utils Antiderivative((~x * ℯ ^ (~n * atanh(~(a') * ~x))) / (~c + ~(d') * (~x) ^ 2) ^ (3 / 2), ~x) => ((1 - ~a * ~n * ~x) * ℯ ^ (~n * atanh(~a * ~x))) / (~d * ((~n) ^ 2 - 1) * sqrt(~c + ~d * (~x) ^ 2)) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative(~x * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => ((~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * atanh(~a * ~x))) / (2 * ~d * (~p + 1)) - ((~a * ~c * ~n) / (2 * ~d * (~p + 1))) * Antiderivative((~c + ~d * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atanh(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (LtQ(~p, -1) && (Not(IntegerQ(~n)) && IntegerQ(2 * ~p))))
    #= Antiderivative((~x)*((~c)+(~d')*(~x)^2)^(~p)*ℯ^((~n)*atanh((~a')*(~x))),~x) :=  -(2*(p+1)+a*n*x)*(c+d*x^2)^(p+1)*ℯ^(n*atanh(a*x))/(d*(n^2-4*(p+1)^2) ) - n*(2*p+3)/(a*c*(n^2-4*(p+1)^2))*Antiderivative((c+d*x^2)^(p+1)*ℯ^(n*atanh(a* x)),x) <-- FreeQ([a,c,d,n],x) && EqQ(a^2*c+d,0) && LeQ(p,-1) &&  NeQ(n^2-4*(p+1)^2,0) && Not(IntegerQ(n)) =#
    @apply_utils Antiderivative((~x) ^ 2 * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => ((1 - ~a * ~n * ~x) * (~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * atanh(~a * ~x))) / (~a * ~d * ~n * ((~n) ^ 2 - 1)) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (EqQ((~n) ^ 2 + 2 * (~p + 1), 0) && Not(IntegerQ(~n))))
    @apply_utils Antiderivative((~x) ^ 2 * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => (-((~n + 2 * (~p + 1) * ~a * ~x)) * (~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * atanh(~a * ~x))) / (~a * ~d * ((~n) ^ 2 - 4 * (~p + 1) ^ 2)) + (((~n) ^ 2 + 2 * (~p + 1)) / (~d * ((~n) ^ 2 - 4 * (~p + 1) ^ 2))) * Antiderivative((~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * atanh(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (LtQ(~p, -1) && (Not(IntegerQ(~n)) && (NeQ((~n) ^ 2 - 4 * (~p + 1) ^ 2, 0) && IntegerQ(2 * ~p)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative((~x) ^ ~m * (1 - (~a) ^ 2 * (~x) ^ 2) ^ (~p - ~n / 2) * (1 + ~a * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && ((IntegerQ(~p) || GtQ(~c, 0)) && (IGtQ((~n + 1) / 2, 0) && Not(IntegerQ(~p - ~n / 2)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative(((~x) ^ ~m * (1 - (~a) ^ 2 * (~x) ^ 2) ^ (~p + ~n / 2)) / (1 - ~a * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && ((IntegerQ(~p) || GtQ(~c, 0)) && (ILtQ((~n - 1) / 2, 0) && Not(IntegerQ(~p - ~n / 2)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative((~x) ^ ~m * (1 - ~a * ~x) ^ (~p - ~n / 2) * (1 + ~a * ~x) ^ (~p + ~n / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (IntegerQ(~p) || GtQ(~c, 0)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => (~c) ^ (~n / 2) * Antiderivative((~x) ^ ~m * (~c + ~d * (~x) ^ 2) ^ (~p - ~n / 2) * (1 + ~a * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (Not(IntegerQ(~p) || GtQ(~c, 0)) && IGtQ(~n / 2, 0)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => (1 / (~c) ^ (~n / 2)) * Antiderivative(((~x) ^ ~m * (~c + ~d * (~x) ^ 2) ^ (~p + ~n / 2)) / (1 - ~a * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (Not(IntegerQ(~p) || GtQ(~c, 0)) && ILtQ(~n / 2, 0)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (((~c) ^ IntPart(~p) * (~c + ~d * (~x) ^ 2) ^ FracPart(~p)) / (1 - (~a) ^ 2 * (~x) ^ 2) ^ FracPart(~p)) * Antiderivative((~x) ^ ~m * (1 - (~a) ^ 2 * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atanh(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (Not(IntegerQ(~p) || GtQ(~c, 0)) && Not(IntegerQ(~n / 2))))
    @apply_utils Antiderivative(~u * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative(~u * (1 - ~a * ~x) ^ (~p - ~n / 2) * (1 + ~a * ~x) ^ (~p + ~n / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (IntegerQ(~p) || GtQ(~c, 0)))
    @apply_utils Antiderivative(~u * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => (((~c) ^ IntPart(~p) * (~c + ~d * (~x) ^ 2) ^ FracPart(~p)) / ((1 - ~a * ~x) ^ FracPart(~p) * (1 + ~a * ~x) ^ FracPart(~p))) * Antiderivative(~u * (1 - ~a * ~x) ^ (~p - ~n / 2) * (1 + ~a * ~x) ^ (~p + ~n / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (Not(IntegerQ(~p) || GtQ(~c, 0)) && IntegerQ(~n / 2)))
    @apply_utils Antiderivative(~u * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (((~c) ^ IntPart(~p) * (~c + ~d * (~x) ^ 2) ^ FracPart(~p)) / (1 - (~a) ^ 2 * (~x) ^ 2) ^ FracPart(~p)) * Antiderivative(~u * (1 - (~a) ^ 2 * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atanh(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (Not(IntegerQ(~p) || GtQ(~c, 0)) && Not(IntegerQ(~n / 2))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (~d) ^ ~p * Antiderivative((~u / (~x) ^ (2 * ~p)) * (1 - (~a) ^ 2 * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atanh(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~c + (~a) ^ 2 * ~d, 0) && IntegerQ(~p))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / (~x) ^ 2) ^ ~p * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative(~u * (1 - 1 / (~a * ~x)) ^ ~p * (1 + 1 / (~a * ~x)) ^ ~p * ℯ ^ (~n * atanh(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ(~c + (~a) ^ 2 * ~d, 0) && (Not(IntegerQ(~p)) && (IntegerQ(~n / 2) && GtQ(~c, 0))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / (~x) ^ 2) ^ ~p * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => (((~x) ^ (2 * ~p) * (~c + ~d / (~x) ^ 2) ^ ~p) / ((1 - ~a * ~x) ^ ~p * (1 + ~a * ~x) ^ ~p)) * Antiderivative((~u / (~x) ^ (2 * ~p)) * (1 - ~a * ~x) ^ ~p * (1 + ~a * ~x) ^ ~p * ℯ ^ (~n * atanh(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~c + (~a) ^ 2 * ~d, 0) && (Not(IntegerQ(~p)) && (IntegerQ(~n / 2) && Not(GtQ(~c, 0)))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * atanh(~(a') * ~x)), ~x) => (((~x) ^ (2 * ~p) * (~c + ~d / (~x) ^ 2) ^ ~p) / (1 + (~c * (~x) ^ 2) / ~d) ^ ~p) * Antiderivative((~u / (~x) ^ (2 * ~p)) * (1 + (~c * (~x) ^ 2) / ~d) ^ ~p * ℯ ^ (~n * atanh(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~c + (~a) ^ 2 * ~d, 0) && (Not(IntegerQ(~p)) && Not(IntegerQ(~n / 2))))
    @apply_utils Antiderivative(ℯ ^ (~(n') * atanh(~(c') * (~a + ~(b') * ~x))), ~x) => Antiderivative((1 + ~a * ~c + ~b * ~c * ~x) ^ (~n / 2) / ((1 - ~a * ~c) - ~b * ~c * ~x) ^ (~n / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~n * atanh(~(c') * (~a + ~(b') * ~x))), ~x) => (4 / (~n * (~b) ^ (~m + 1) * (~c) ^ (~m + 1))) * Subst(Antiderivative(((~x) ^ (2 / ~n) * ((-1 - ~a * ~c) + (1 - ~a * ~c) * (~x) ^ (2 / ~n)) ^ ~m) / (1 + (~x) ^ (2 / ~n)) ^ (~m + 2), ~x), ~x, (1 + ~c * (~a + ~b * ~x)) ^ (~n / 2) / (1 - ~c * (~a + ~b * ~x)) ^ (~n / 2)) <-- FreeQ([~a, ~b, ~c], ~x) && (ILtQ(~m, 0) && LtQ(-1, ~n, 1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ℯ ^ (~(n') * atanh(~(c') * (~a + ~(b') * ~x))), ~x) => Antiderivative(((~d + ~e * ~x) ^ ~m * (1 + ~a * ~c + ~b * ~c * ~x) ^ (~n / 2)) / ((1 - ~a * ~c) - ~b * ~c * ~x) ^ (~n / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x + ~(e') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atanh(~a + ~(b') * ~x)), ~x) => (~c / (1 - (~a) ^ 2)) ^ ~p * Antiderivative(~u * ((1 - ~a) - ~b * ~x) ^ (~p - ~n / 2) * (1 + ~a + ~b * ~x) ^ (~p + ~n / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (EqQ(~b * ~d - 2 * ~a * ~e, 0) && (EqQ((~b) ^ 2 * ~c + ~e * (1 - (~a) ^ 2), 0) && (IntegerQ(~p) || GtQ(~c / (1 - (~a) ^ 2), 0))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x + ~(e') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atanh(~a + ~(b') * ~x)), ~x) => ((~c + ~d * ~x + ~e * (~x) ^ 2) ^ ~p / (((1 - (~a) ^ 2) - 2 * ~a * ~b * ~x) - (~b) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative(~u * (((1 - (~a) ^ 2) - 2 * ~a * ~b * ~x) - (~b) ^ 2 * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atanh(~a * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (EqQ(~b * ~d - 2 * ~a * ~e, 0) && (EqQ((~b) ^ 2 * ~c + ~e * (1 - (~a) ^ 2), 0) && Not(IntegerQ(~p) || GtQ(~c / (1 - (~a) ^ 2), 0))))
    @apply_utils Antiderivative(~(u') * ℯ ^ (~(n') * atanh(~(c') / (~(a') + ~(b') * ~x))), ~x) => Antiderivative(~u * ℯ ^ (~n * acoth(~a / ~c + (~b * ~x) / ~c)), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
    @apply_utils Antiderivative(~(u') * ℯ ^ (~n * acoth(~(a') * ~x)), ~x) => (-1) ^ (~n / 2) * Antiderivative(~u * ℯ ^ (~n * atanh(~a * ~x)), ~x) <-- FreeQ(~a, ~x) && IntegerQ(~n / 2)
    @apply_utils Antiderivative(ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => -(Subst(Antiderivative((1 + ~x / ~a) ^ ((~n + 1) / 2) / ((~x) ^ 2 * (1 - ~x / ~a) ^ ((~n - 1) / 2) * sqrt(1 - (~x) ^ 2 / (~a) ^ 2)), ~x), ~x, 1 / ~x)) <-- FreeQ(~a, ~x) && IntegerQ((~n - 1) / 2)
    @apply_utils Antiderivative((~x) ^ ~(m') * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => -(Subst(Antiderivative((1 + ~x / ~a) ^ ((~n + 1) / 2) / ((~x) ^ (~m + 2) * (1 - ~x / ~a) ^ ((~n - 1) / 2) * sqrt(1 - (~x) ^ 2 / (~a) ^ 2)), ~x), ~x, 1 / ~x)) <-- FreeQ(~a, ~x) && (IntegerQ((~n - 1) / 2) && IntegerQ(~m))
    @apply_utils Antiderivative(ℯ ^ (~n * acoth(~(a') * ~x)), ~x) => -(Subst(Antiderivative((1 + ~x / ~a) ^ (~n / 2) / ((~x) ^ 2 * (1 - ~x / ~a) ^ (~n / 2)), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~n], ~x) && Not(IntegerQ(~n))
    @apply_utils Antiderivative((~x) ^ ~(m') * ℯ ^ (~n * acoth(~(a') * ~x)), ~x) => -(Subst(Antiderivative((1 + ~x / ~a) ^ (~n / 2) / ((~x) ^ (~m + 2) * (1 - ~x / ~a) ^ (~n / 2)), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~n], ~x) && (Not(IntegerQ(~n)) && IntegerQ(~m))
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => -((~x) ^ ~m) * (1 / ~x) ^ ~m * Subst(Antiderivative((1 + ~x / ~a) ^ ((~n + 1) / 2) / ((~x) ^ (~m + 2) * (1 - ~x / ~a) ^ ((~n - 1) / 2) * sqrt(1 - (~x) ^ 2 / (~a) ^ 2)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~m], ~x) && (IntegerQ((~n - 1) / 2) && Not(IntegerQ(~m)))
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~n * acoth(~(a') * ~x)), ~x) => -((~x) ^ ~m) * (1 / ~x) ^ ~m * Subst(Antiderivative((1 + ~x / ~a) ^ (~n / 2) / ((~x) ^ (~m + 2) * (1 - ~x / ~a) ^ (~n / 2)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~m, ~n], ~x) && (Not(IntegerQ(~n)) && Not(IntegerQ(~m)))
    @apply_utils Antiderivative((~c + ~(d') * ~x) ^ ~(p') * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => ((1 + ~a * ~x) * (~c + ~d * ~x) ^ ~p * ℯ ^ (~n * acoth(~a * ~x))) / (~a * (~p + 1)) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~a * ~c + ~d, 0) && (EqQ(~p, ~n / 2) && Not(IntegerQ(~n / 2))))
    #= Antiderivative(((~c)+(~d')*(~x))^(~p')*ℯ^((~n')*acoth((~a')*(~x))),~x) :=  -(-a)^n*c^n*Subst(Antiderivative((d+c*x)^(p-n)*(1-x^2/a^2)^(n/2)/x^(p+2),x),x,1/ x) <-- FreeQ([a,c,d],x) && EqQ(a*c+d,0) && IntegerQ((n-1)/2) &&  IntegerQ(p) =#
    #= Antiderivative((~x)^(~m')*((~c)+(~d')*(~x))^(~p')*ℯ^((~n')*acoth((~a')*(~x))),~x) :=  -(-a)^n*c^n*Subst(Antiderivative((d+c*x)^(p-n)*(1-x^2/a^2)^(n/2)/x^(m+p+2),x),x, 1/x) <-- FreeQ([a,c,d],x) && EqQ(a*c+d,0) && IntegerQ((n-1)/2) &&  IntegerQ(m) && IntegerQ(p) =#
    #= Antiderivative(((~c)+(~d')*(~x))^(~p')*ℯ^((~n')*acoth((~a')*(~x))),~x) :=  -(-a)^n*c^n*sqrt(c+d*x)/(sqrt(x)*sqrt(d+c/x))*Subst(Antiderivative((d+c*x)^(p-n)* (1-x^2/a^2)^(n/2)/x^(p+2),x),x,1/x) <-- FreeQ([a,c,d],x) &&  EqQ(a*c+d,0) && IntegerQ((n-1)/2) && IntegerQ(p-1/2) =#
    #= Antiderivative((~x)^(~m')*((~c)+(~d')*(~x))^(~p')*ℯ^((~n')*acoth((~a')*(~x))),~x) :=  -(-a)^n*c^n*sqrt(c+d*x)/(sqrt(x)*sqrt(d+c/x))*Subst(Antiderivative((d+c*x)^(p-n)* (1-x^2/a^2)^(n/2)/x^(m+p+2),x),x,1/x) <-- FreeQ([a,c,d],x) &&  EqQ(a*c+d,0) && IntegerQ((n-1)/2) && IntegerQ(m) && IntegerQ(p-1/2) =#
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x) ^ ~(p') * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => (~d) ^ ~p * Antiderivative(~u * (~x) ^ ~p * (1 + ~c / (~d * ~x)) ^ ~p * ℯ ^ (~n * acoth(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * (~c) ^ 2 - (~d) ^ 2, 0) && (Not(IntegerQ(~n / 2)) && IntegerQ(~p)))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x) ^ ~p * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => ((~c + ~d * ~x) ^ ~p / ((~x) ^ ~p * (1 + ~c / (~d * ~x)) ^ ~p)) * Antiderivative(~u * (~x) ^ ~p * (1 + ~c / (~d * ~x)) ^ ~p * ℯ ^ (~n * acoth(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~a) ^ 2 * (~c) ^ 2 - (~d) ^ 2, 0) && (Not(IntegerQ(~n / 2)) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~c + ~(d') / ~x) ^ ~(p') * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => -((~c) ^ ~n) * Subst(Antiderivative(((~c + ~d * ~x) ^ (~p - ~n) * (1 - (~x) ^ 2 / (~a) ^ 2) ^ (~n / 2)) / (~x) ^ 2, ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ(~c + ~a * ~d, 0) && (IntegerQ((~n - 1) / 2) && ((IntegerQ(~p) || (EqQ(~p, ~n / 2) || EqQ(~p, ~n / 2 + 1))) && IntegerQ(2 * ~p))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') / ~x) ^ ~(p') * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => -((~c) ^ ~n) * Subst(Antiderivative(((~c + ~d * ~x) ^ (~p - ~n) * (1 - (~x) ^ 2 / (~a) ^ 2) ^ (~n / 2)) / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ(~c + ~a * ~d, 0) && (IntegerQ((~n - 1) / 2) && (IntegerQ(~m) && ((IntegerQ(~p) || (EqQ(~p, ~n / 2) || (EqQ(~p, ~n / 2 + 1) || LtQ(-5, ~m, -1)))) && IntegerQ(2 * ~p)))))
    @apply_utils Antiderivative((~c + ~(d') / ~x) ^ ~(p') * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => -((~c) ^ ~p) * Subst(Antiderivative(((1 + (~d * ~x) / ~c) ^ ~p * (1 + ~x / ~a) ^ (~n / 2)) / ((~x) ^ 2 * (1 - ~x / ~a) ^ (~n / 2)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~c) ^ 2 - (~a) ^ 2 * (~d) ^ 2, 0) && (Not(IntegerQ(~n / 2)) && (IntegerQ(~p) || GtQ(~c, 0))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') / ~x) ^ ~(p') * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => -((~c) ^ ~p) * Subst(Antiderivative(((1 + (~d * ~x) / ~c) ^ ~p * (1 + ~x / ~a) ^ (~n / 2)) / ((~x) ^ (~m + 2) * (1 - ~x / ~a) ^ (~n / 2)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~c) ^ 2 - (~a) ^ 2 * (~d) ^ 2, 0) && (Not(IntegerQ(~n / 2)) && ((IntegerQ(~p) || GtQ(~c, 0)) && IntegerQ(~m))))
    @apply_utils Antiderivative((~x) ^ ~m * (~c + ~(d') / ~x) ^ ~(p') * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => -((~c) ^ ~p) * (~x) ^ ~m * (1 / ~x) ^ ~m * Subst(Antiderivative(((1 + (~d * ~x) / ~c) ^ ~p * (1 + ~x / ~a) ^ (~n / 2)) / ((~x) ^ (~m + 2) * (1 - ~x / ~a) ^ (~n / 2)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ((~c) ^ 2 - (~a) ^ 2 * (~d) ^ 2, 0) && (Not(IntegerQ(~n / 2)) && ((IntegerQ(~p) || GtQ(~c, 0)) && Not(IntegerQ(~m)))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / ~x) ^ ~p * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => ((~c + ~d / ~x) ^ ~p / (1 + ~d / (~c * ~x)) ^ ~p) * Antiderivative(~u * (1 + ~d / (~c * ~x)) ^ ~p * ℯ ^ (~n * acoth(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~c) ^ 2 - (~a) ^ 2 * (~d) ^ 2, 0) && (Not(IntegerQ(~n / 2)) && Not(IntegerQ(~p) || GtQ(~c, 0))))
    @apply_utils Antiderivative(ℯ ^ (~(n') * acoth(~(a') * ~x)) / (~c + ~(d') * (~x) ^ 2), ~x) => ℯ ^ (~n * acoth(~a * ~x)) / (~a * ~c * ~n) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && Not(IntegerQ(~n / 2)))
    @apply_utils Antiderivative(ℯ ^ (~n * acoth(~(a') * ~x)) / (~c + ~(d') * (~x) ^ 2) ^ (3 / 2), ~x) => ((~n - ~a * ~x) * ℯ ^ (~n * acoth(~a * ~x))) / (~a * ~c * ((~n) ^ 2 - 1) * sqrt(~c + ~d * (~x) ^ 2)) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => ((~n + 2 * ~a * (~p + 1) * ~x) * (~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * acoth(~a * ~x))) / (~a * ~c * ((~n) ^ 2 - 4 * (~p + 1) ^ 2)) - ((2 * (~p + 1) * (2 * ~p + 3)) / (~c * ((~n) ^ 2 - 4 * (~p + 1) ^ 2))) * Antiderivative((~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * acoth(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (Not(IntegerQ(~n / 2)) && (LtQ(~p, -1) && (NeQ(~p, -3 / 2) && (NeQ((~n) ^ 2 - 4 * (~p + 1) ^ 2, 0) && (IntegerQ(~p) || Not(IntegerQ(~n))))))))
    @apply_utils Antiderivative((~x * ℯ ^ (~n * acoth(~(a') * ~x))) / (~c + ~(d') * (~x) ^ 2) ^ (3 / 2), ~x) => (-((1 - ~a * ~n * ~x)) * ℯ ^ (~n * acoth(~a * ~x))) / ((~a) ^ 2 * ~c * ((~n) ^ 2 - 1) * sqrt(~c + ~d * (~x) ^ 2)) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && Not(IntegerQ(~n)))
    @apply_utils Antiderivative(~x * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => ((2 * (~p + 1) + ~a * ~n * ~x) * (~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * acoth(~a * ~x))) / ((~a) ^ 2 * ~c * ((~n) ^ 2 - 4 * (~p + 1) ^ 2)) - ((~n * (2 * ~p + 3)) / (~a * ~c * ((~n) ^ 2 - 4 * (~p + 1) ^ 2))) * Antiderivative((~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * acoth(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (Not(IntegerQ(~n / 2)) && (LeQ(~p, -1) && (NeQ(~p, -3 / 2) && (NeQ((~n) ^ 2 - 4 * (~p + 1) ^ 2, 0) && (IntegerQ(~p) || Not(IntegerQ(~n))))))))
    @apply_utils Antiderivative((~x) ^ 2 * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => (-((~n + 2 * (~p + 1) * ~a * ~x)) * (~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * acoth(~a * ~x))) / ((~a) ^ 3 * ~c * (~n) ^ 2 * ((~n) ^ 2 - 1)) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (Not(IntegerQ(~n / 2)) && (EqQ((~n) ^ 2 + 2 * (~p + 1), 0) && NeQ((~n) ^ 2, 1))))
    @apply_utils Antiderivative((~x) ^ 2 * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => ((~n + 2 * (~p + 1) * ~a * ~x) * (~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * acoth(~a * ~x))) / ((~a) ^ 3 * ~c * ((~n) ^ 2 - 4 * (~p + 1) ^ 2)) - (((~n) ^ 2 + 2 * (~p + 1)) / ((~a) ^ 2 * ~c * ((~n) ^ 2 - 4 * (~p + 1) ^ 2))) * Antiderivative((~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * acoth(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (Not(IntegerQ(~n / 2)) && (LeQ(~p, -1) && (NeQ((~n) ^ 2 + 2 * (~p + 1), 0) && (NeQ((~n) ^ 2 - 4 * (~p + 1) ^ 2, 0) && (IntegerQ(~p) || Not(IntegerQ(~n))))))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => (-((-(~c)) ^ ~p) / (~a) ^ (~m + 1)) * Subst(Antiderivative((ℯ ^ (~n * ~x) * coth(~x) ^ (~m + 2 * (~p + 1))) / cosh(~x) ^ (2 * (~p + 1)), ~x), ~x, acoth(~a * ~x)) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (Not(IntegerQ(~n / 2)) && (IntegerQ(~m) && (LeQ(3, ~m, -2 * (~p + 1)) && IntegerQ(~p)))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => (~d) ^ ~p * Antiderivative(~u * (~x) ^ (2 * ~p) * (1 - 1 / ((~a) ^ 2 * (~x) ^ 2)) ^ ~p * ℯ ^ (~n * acoth(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (Not(IntegerQ(~n / 2)) && IntegerQ(~p)))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => ((~c + ~d * (~x) ^ 2) ^ ~p / ((~x) ^ (2 * ~p) * (1 - 1 / ((~a) ^ 2 * (~x) ^ 2)) ^ ~p)) * Antiderivative(~u * (~x) ^ (2 * ~p) * (1 - 1 / ((~a) ^ 2 * (~x) ^ 2)) ^ ~p * ℯ ^ (~n * acoth(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~a) ^ 2 * ~c + ~d, 0) && (Not(IntegerQ(~n / 2)) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => ((~c) ^ ~p / (~a) ^ (2 * ~p)) * Antiderivative((~u / (~x) ^ (2 * ~p)) * (-1 + ~a * ~x) ^ (~p - ~n / 2) * (1 + ~a * ~x) ^ (~p + ~n / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~c + (~a) ^ 2 * ~d, 0) && (Not(IntegerQ(~n / 2)) && ((IntegerQ(~p) || GtQ(~c, 0)) && IntegersQ(2 * ~p, ~p + ~n / 2))))
    @apply_utils Antiderivative((~c + ~(d') / (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => -((~c) ^ ~p) * Subst(Antiderivative(((1 - ~x / ~a) ^ (~p - ~n / 2) * (1 + ~x / ~a) ^ (~p + ~n / 2)) / (~x) ^ 2, ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~c + (~a) ^ 2 * ~d, 0) && (Not(IntegerQ(~n / 2)) && ((IntegerQ(~p) || GtQ(~c, 0)) && Not(IntegersQ(2 * ~p, ~p + ~n / 2)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') / (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => -((~c) ^ ~p) * Subst(Antiderivative(((1 - ~x / ~a) ^ (~p - ~n / 2) * (1 + ~x / ~a) ^ (~p + ~n / 2)) / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~c + (~a) ^ 2 * ~d, 0) && (Not(IntegerQ(~n / 2)) && ((IntegerQ(~p) || GtQ(~c, 0)) && (Not(IntegersQ(2 * ~p, ~p + ~n / 2)) && IntegerQ(~m)))))
    @apply_utils Antiderivative((~x) ^ ~m * (~c + ~(d') / (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => -((~c) ^ ~p) * (~x) ^ ~m * (1 / ~x) ^ ~m * Subst(Antiderivative(((1 - ~x / ~a) ^ (~p - ~n / 2) * (1 + ~x / ~a) ^ (~p + ~n / 2)) / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ(~c + (~a) ^ 2 * ~d, 0) && (Not(IntegerQ(~n / 2)) && ((IntegerQ(~p) || GtQ(~c, 0)) && (Not(IntegersQ(2 * ~p, ~p + ~n / 2)) && Not(IntegerQ(~m))))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * acoth(~(a') * ~x)), ~x) => (((~c) ^ IntPart(~p) * (~c + ~d / (~x) ^ 2) ^ FracPart(~p)) / (1 - 1 / ((~a) ^ 2 * (~x) ^ 2)) ^ FracPart(~p)) * Antiderivative(~u * (1 - 1 / ((~a) ^ 2 * (~x) ^ 2)) ^ ~p * ℯ ^ (~n * acoth(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~c + (~a) ^ 2 * ~d, 0) && (Not(IntegerQ(~n / 2)) && Not(IntegerQ(~p) || GtQ(~c, 0))))
    @apply_utils Antiderivative(~(u') * ℯ ^ (~n * acoth(~(c') * (~a + ~(b') * ~x))), ~x) => (-1) ^ (~n / 2) * Antiderivative(~u * ℯ ^ (~n * atanh(~c * (~a + ~b * ~x))), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && IntegerQ(~n / 2)
    @apply_utils Antiderivative(ℯ ^ (~(n') * acoth(~(c') * (~a + ~(b') * ~x))), ~x) => (((~c * (~a + ~b * ~x)) ^ (~n / 2) * (1 + 1 / (~c * (~a + ~b * ~x))) ^ (~n / 2)) / (1 + ~a * ~c + ~b * ~c * ~x) ^ (~n / 2)) * Antiderivative((1 + ~a * ~c + ~b * ~c * ~x) ^ (~n / 2) / (-1 + ~a * ~c + ~b * ~c * ~x) ^ (~n / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && Not(IntegerQ(~n / 2))
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~n * acoth(~(c') * (~a + ~(b') * ~x))), ~x) => (-4 / (~n * (~b) ^ (~m + 1) * (~c) ^ (~m + 1))) * Subst(Antiderivative(((~x) ^ (2 / ~n) * (1 + ~a * ~c + (1 - ~a * ~c) * (~x) ^ (2 / ~n)) ^ ~m) / (-1 + (~x) ^ (2 / ~n)) ^ (~m + 2), ~x), ~x, (1 + 1 / (~c * (~a + ~b * ~x))) ^ (~n / 2) / (1 - 1 / (~c * (~a + ~b * ~x))) ^ (~n / 2)) <-- FreeQ([~a, ~b, ~c], ~x) && (ILtQ(~m, 0) && LtQ(-1, ~n, 1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ℯ ^ (~(n') * acoth(~(c') * (~a + ~(b') * ~x))), ~x) => (((~c * (~a + ~b * ~x)) ^ (~n / 2) * (1 + 1 / (~c * (~a + ~b * ~x))) ^ (~n / 2)) / (1 + ~a * ~c + ~b * ~c * ~x) ^ (~n / 2)) * Antiderivative(((~d + ~e * ~x) ^ ~m * (1 + ~a * ~c + ~b * ~c * ~x) ^ (~n / 2)) / (-1 + ~a * ~c + ~b * ~c * ~x) ^ (~n / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && Not(IntegerQ(~n / 2))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x + ~(e') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acoth(~a + ~(b') * ~x)), ~x) => (~c / (1 - (~a) ^ 2)) ^ ~p * ((~a + ~b * ~x) / (1 + ~a + ~b * ~x)) ^ (~n / 2) * ((1 + ~a + ~b * ~x) / (~a + ~b * ~x)) ^ (~n / 2) * (((1 - ~a) - ~b * ~x) ^ (~n / 2) / (-1 + ~a + ~b * ~x) ^ (~n / 2)) * Antiderivative(~u * ((1 - ~a) - ~b * ~x) ^ (~p - ~n / 2) * (1 + ~a + ~b * ~x) ^ (~p + ~n / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (Not(IntegerQ(~n / 2)) && (EqQ(~b * ~d - 2 * ~a * ~e, 0) && (EqQ((~b) ^ 2 * ~c + ~e * (1 - (~a) ^ 2), 0) && (IntegerQ(~p) || GtQ(~c / (1 - (~a) ^ 2), 0)))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x + ~(e') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acoth(~a + ~(b') * ~x)), ~x) => ((~c + ~d * ~x + ~e * (~x) ^ 2) ^ ~p / (((1 - (~a) ^ 2) - 2 * ~a * ~b * ~x) - (~b) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative(~u * (((1 - (~a) ^ 2) - 2 * ~a * ~b * ~x) - (~b) ^ 2 * (~x) ^ 2) ^ ~p * ℯ ^ (~n * acoth(~a * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (Not(IntegerQ(~n / 2)) && (EqQ(~b * ~d - 2 * ~a * ~e, 0) && (EqQ((~b) ^ 2 * ~c + ~e * (1 - (~a) ^ 2), 0) && Not(IntegerQ(~p) || GtQ(~c / (1 - (~a) ^ 2), 0)))))
    @apply_utils Antiderivative(~(u') * ℯ ^ (~(n') * acoth(~(c') / (~(a') + ~(b') * ~x))), ~x) => Antiderivative(~u * ℯ ^ (~n * atanh(~a / ~c + (~b * ~x) / ~c)), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
end

max_stars_repo_path: src/rules/7 Inverse hyperbolic functions/7.3 Inverse hyperbolic tangent/7.3.3 Exponentials of inverse hyperbolic tangent.jl


text:
module Commutator

export commutator!

"""
    commutator!(alpha, A, B, beta, C)

Calculate `C = beta * C + alpha * (A * B - B * A)` in-place.
"""
function commutator!(alpha, A, B, beta, C)
    C .*= beta
    C .+= alpha * (A * B - B * A)
end

function commutator!(
        alpha::Number, A::AbstractMatrix, B::AbstractMatrix, beta::Number,
        C::AbstractMatrix)
    C[:,:] = beta * C[:,:] + alpha * (A * B - B * A)
end

for elty in (:Float64, :Float32, :Complex128, :Complex64)
    @eval begin
        function commutator!(
                alpha::($elty), A::StridedVecOrMat{$elty},
                B::StridedVecOrMat{$elty}, beta::($elty),
                C::StridedVecOrMat{$elty})
            BLAS.gemm!('N', 'N', alpha, A, B, beta, C)
            BLAS.gemm!('N', 'N', -alpha, B, A, one(beta), C)
        end
        function commutator!(
                alpha::Number, A::StridedVecOrMat{$elty},
                B::StridedVecOrMat{$elty}, beta::Number,
                C::StridedVecOrMat{$elty})
            alpha = convert($elty, alpha)
            beta = convert($elty, beta)
            BLAS.gemm!('N', 'N', alpha, A, B, beta, C)
            BLAS.gemm!('N', 'N', -alpha, B, A, one(beta), C)
        end
    end
end

function commutator!(
        alpha::Number, A::SparseMatrixCSC, B::AbstractMatrix,
        beta::Number, C::AbstractMatrix)
    if beta != 1
        C .*= beta
    end
    for j = 1: A.n
        for k = 1:A.n, i_val = A.colptr[k] : (A.colptr[k+1]-1)
            i = A.rowval[i_val]
            A_ik = A.nzval[i_val]
            C[i,j] = C[i,j] + alpha * A_ik * B[k,j]
        end
    end
    for j = 1:A.n, i_val = A.colptr[j] : (A.colptr[j+1]-1)
        k = A.rowval[i_val]
        A_kj = A.nzval[i_val]
        for i = 1: A.n
            C[i,j] = C[i,j] - alpha * B[i, k] * A_kj
        end
    end
end

end # module


text:
module func

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
"""
	function append_to_curve!
"""
function append_to_curve!(target::Curve, newstuff::Curve)
    s_end = target[end].s
    for c in newstuff
        push!(target, CurvePt(c.pos, c.s+s_end, c.k, c.kd))
    end
    return target
end

"""
    function get_new_angle
- Does the actual angle calculation based on the x y coordinates
"""
function get_new_angle(tangent_vector::Array{Float64})
    # it might be a problem when we switch quadrants
    # use signs of tangent vector to get the quadrant of the heading
    x = tangent_vector[1]
    y = tangent_vector[2]
    if x == 0. && y == 0.
        heading = 0.0
    elseif x == 0.
        heading = π/2 * sign(y)
    elseif y == 0.
        heading = convert(Float64, π) # this could be either pi or -pi, but just go with pi
    elseif sign(x) == 1 && sign(y) == 1 # first quadrant
        heading = atan(y, x)
    elseif sign(x) == -1 && sign(y) == 1 # second quadrant
        heading = atan(y, x)
    elseif sign(x) == 1 && sign(y) == -1 # fourth quadrant
        heading = atan(y, x)
    elseif sign(x) == -1 && sign(y) == -1 # third quadrant
        heading = atan(y, x)
    end
    # bound_heading doesn't end up getting called cause Julia takes care of it apparently
    bound_heading(heading)

    return heading
end

"""
    function bound_heading
- Make the angle range from 0 to pi instead of going beyond
"""
function bound_heading(heading::Float64)
    if heading > π # send this to be just near -pi
        heading = -π + (heading - π)    # if heading is 3.15, then the new angle will be (-pi + (3.15-pi)) = -3.13
    elseif heading < -π # send this to be just near pi
        heading = π + (heading + π)     # if heading is -3.15, then the new angle will be (pi + (-3.15+pi)) = 3.13
    end
    return heading
end

"""
    function append_headings
- Used create the angles and append them into the coordinates

# Examples
```julia
x_coods = [1089.07510, 1093.82626, 1101.19325, 1112.59899, 1123.96733, 1133.24150, 1146.47964]
y_coods = [936.31213, 936.92692,938.52419, 940.93865, 943.27882, 945.21039, 947.88488]
coods = hcat(x_coods,y_coods)
append_headings(coods)
```
"""
function append_headings(coordinates::Matrix{Float64})
    headings = ones(size(coordinates)[1])
    for i = 1:size(coordinates)[1]-1
        # need to make sure that math is right, and that bounds are kept
        tangent_vector = [coordinates[i+1,1]-coordinates[i,1], coordinates[i+1,2]-coordinates[i,2]]
        # @show tangent_vector
        current_heading = get_new_angle(tangent_vector)
        # @show current_heading
        headings[i] = current_heading
    end
    headings[end] = headings[end-1] # assume this is fine
    coordinates = hcat(coordinates, headings)
    return coordinates
end

# function: centerlines_txt2tracks_new
"""
    function centerlines_txt2tracks(filename)
- Reads a .txt file which contains x coords in col 1 and y coords in col2
- Returns a track, i.e. a `curve` from AutomotiveDrivingModels

# Examples
```julia
# See make_roadway_interaction()
```
"""
function centerlines_txt2tracks(filename)
    coods = readdlm(filename,',')
    coods_app = append_headings(coods) # Append with the angle

    mid_coods = coods_app'

    first_cood = VecSE2(mid_coods[1,1], mid_coods[2,1], mid_coods[3,1])
    second_cood = VecSE2(mid_coods[1,2], mid_coods[2,2], mid_coods[3,2])
    radius = 0.01
    nsamples = 20

    track = gen_bezier_curve(first_cood, second_cood, radius, radius, nsamples)

    nsamples = 20
    for i = 3:size(coods,1)
        turn1 = VecSE2(mid_coods[1, i-1], mid_coods[2, i-1], mid_coods[3, i-1])
        turn2 = VecSE2(mid_coods[1, i], mid_coods[2, i], mid_coods[3, i])
        curve = gen_bezier_curve(turn1, turn2, radius, radius, nsamples)
        append_to_curve!(track, curve)
    end

    return track
end

# function: read in the centerline text files and make `roadway_interaction`
"""
    function make_roadway_interaction()
# Overview
- Make the `DR_CHN_Merging` roadway by reading in the centerlines in `AutomotiveInteraction.jl/dataset`
- These centerlines have been modified from the original ones in the `centerlines_DR_CHN_Merging_ZS` folder
- The long lane being merged into by on ramp was split into before and after merge point in 2 separate txt files
- Eg: `output_centerline_1.txt` has been split into `centerlines_b1.txt` (before merge point) and `centerlines_b2.txt` (after).
- Similarly, `output_centerline_5.txt` has been split into `cetnerlines_g.txt` (before) and `centerlines_h.txt` (after)
- Finally, the files in `centerlines_DR_CHN_Merging_ZS` folder `output_centerline_<>.txt` had x coords in row 1 and y coords in row 2.
- The current `centerlines_<>.txt` have x in col 1 and y in col 2

# Details
- Reads in the centerline information from `AutomotiveInteraction.jl/dataset`
- Segment 1: Has 2 lanes: On ramp `a` is lane 1 and `b1` is lane 2
- Segment 2: Has 1 lane: `b1` is lane 1. Both `a` and `b1` connect into `b2`
- Segment 3: Has 2 lanes: `c` is lane 1 and `d` is lane 2
- Segment 4: Has 3 lanes: `g` is lane 1, `f1` is lane 2, `e1` is lane 3.
- Segment 5: Has 2 lanes: `f2` is lane 1, `e2` is lane 2. `e1` connects to `e2` and both `g` and `f1` to `e2`
    - Note that convention used here is that right most lane in direction of travel is given lane id 1

# Examples
```julia
roadway_interaction = make_roadway_interaction()
```
"""
function make_roadway_interaction()
    road = Roadway()
        # Make segment 1: the on ramp a and first part of lane: b1
    track_a = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_a.txt")); # Top most on ramp
    lane_a = Lane(LaneTag(1,1),track_a,boundary_left=LaneBoundary(:broken,:white))

    track_b1 = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_b1.txt")); # Top most on ramp
    lane_b1 = Lane(LaneTag(1,2),track_b1)

        # Make segment 2: second part of lane: b2. And connect both lanes of segment 1 into segment 2
    track_b2 = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_b2.txt"))
    lane_b2 = Lane(LaneTag(2,1),track_b2)
    connect!(lane_a,lane_b2)
    connect!(lane_b1,lane_b2)

    push!(road.segments,RoadSegment(1,[lane_a,lane_b1]))
    push!(road.segments,RoadSegment(2,[lane_b2]))

        # Make segment 3: c and d
    track_c = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_c.txt"))
    lane_c = Lane(LaneTag(3,1),track_c,boundary_left=LaneBoundary(:broken,:white))

    track_d = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_d.txt"))
    lane_d = Lane(LaneTag(3,2),track_d)

    push!(road.segments,RoadSegment(3,[lane_c,lane_d]))

            # Other side of the divider
        # Make segment 4: g,f1,e1
    track_g = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_g.txt"))
    lane_g = Lane(LaneTag(4,1),track_g,boundary_left=LaneBoundary(:broken,:white))

    track_f1 = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_f1.txt"))
    lane_f1 = Lane(LaneTag(4,2),track_f1,boundary_left=LaneBoundary(:broken,:white))

    track_e1 = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_e1.txt"))
    lane_e1 = Lane(LaneTag(4,3),track_e1)

        # Make segment 5: f2,e2. And connect 3 lanes of segment 4 into segment 5 two lanes
    track_f2 = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_f2.txt"))
    lane_f2 = Lane(LaneTag(5,1),track_f2,boundary_left=LaneBoundary(:broken,:white))

    track_e2 = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_e2.txt"))
    lane_e2 = Lane(LaneTag(5,2),track_e2)
    connect!(lane_e1,lane_e2)
    connect!(lane_f1,lane_f2)
    connect!(lane_g,lane_f2)

    push!(road.segments,RoadSegment(4,[lane_g,lane_f1,lane_e1]))
    push!(road.segments,RoadSegment(5,[lane_f2,lane_e2]))

    return road
end

"""
    function make_roadway_ngsim()

- Read in the ngsim 101 roadway from the provided roadway file `/dataset/ngsim_101.txt`

# Examples
```julia
roadway_ngsim = make_roadway_ngsim()
```
"""
function make_roadway_ngsim()
    roadway_ngsim = open(io->read(io, Roadway),joinpath(@__DIR__,"dataset/ngsim_101.txt"), "r")
    return roadway_ngsim
end

# function: To create road extensions
"""
- Takes a CurvePt `c` and a distance `l`
- Returns a curve of length `l` starting at `c` and preserving its angle

# Examples
```julia
c_new = straight_extension(c)
track_extension = gen_straight_segment(c,c_new)
append_to_curve(track_orig,track_extension)
```
"""
function straight_extension(c::CurvePt,l::Float64)
    x1 = c.pos.x
    y1 = c.pos.y
    t = c.pos.θ
    x2 = x1+l*cos(t)
    y2 = y1+l*sin(t)
    return CurvePt(VecSE2(x2,y2,t),c.s+l)
    #return gen_straight_curve(VecE2(x1,y1),VecE2(x2,y2),2)
end

# function: Take a track and straight extension it
"""
    function extend_track_straight!(track,l::Float64)
- Take a `track` and length `l` and extend straight after end

# Examples
```julia
extended_track_e2 = extend_track_straight(track_e2,50.)
```
"""
function extend_track_straight!(track,l::Float64)
    c_ext = straight_extension(track[end],l)
    push!(track,c_ext)
    return track
end

# function: roadway_interaction with straight extensions to avoid veh reversal
"""
- Make roadway interaction but with added extensions to avoid vehicle getting blocked at road end
"""
function make_roadway_interaction_with_extensions()
    road = Roadway()
        # Make segment 1: the on ramp a and first part of lane: b1
    track_a = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_a.txt")); # Top most on ramp
    lane_a = Lane(LaneTag(1,1),track_a,boundary_left=LaneBoundary(:broken,:white))

    track_b1 = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_b1.txt")); # Top most on ramp
    lane_b1 = Lane(LaneTag(1,2),track_b1)

        # Make segment 2: second part of lane: b2. And connect both lanes of segment 1 into segment 2
    track_b2_unext = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_b2.txt"))
    track_b2 = extend_track_straight!(track_b2_unext,100.)
    lane_b2 = Lane(LaneTag(2,1),track_b2)
    connect!(lane_a,lane_b2)
    connect!(lane_b1,lane_b2)

    push!(road.segments,RoadSegment(1,[lane_a,lane_b1]))
    push!(road.segments,RoadSegment(2,[lane_b2]))

        # Make segment 3: c and d
    track_c_unext = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_c.txt"))
    track_c = extend_track_straight!(track_c_unext,100.)
    lane_c = Lane(LaneTag(3,1),track_c,boundary_left=LaneBoundary(:broken,:white))

    track_d_unext = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_d.txt"))
    track_d = extend_track_straight!(track_d_unext,100.)
    lane_d = Lane(LaneTag(3,2),track_d)

    push!(road.segments,RoadSegment(3,[lane_c,lane_d]))

            # Other side of the divider
        # Make segment 4: g,f1,e1
    track_g = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_g.txt"))
    lane_g = Lane(LaneTag(4,1),track_g,boundary_left=LaneBoundary(:broken,:white))

    track_f1 = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_f1.txt"))
    lane_f1 = Lane(LaneTag(4,2),track_f1,boundary_left=LaneBoundary(:broken,:white))

    track_e1 = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_e1.txt"))
    lane_e1 = Lane(LaneTag(4,3),track_e1)

        # Make segment 5: f2,e2. And connect 3 lanes of segment 4 into segment 5 two lanes
    track_f2_unext = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_f2.txt"))
    track_f2 = extend_track_straight!(track_f2_unext,100.)
    lane_f2 = Lane(LaneTag(5,1),track_f2,boundary_left=LaneBoundary(:broken,:white))

    track_e2_unext = centerlines_txt2tracks(joinpath(@__DIR__,"dataset/centerlines_e2.txt"))
    track_e2 = extend_track_straight!(track_e2_unext,100.)
    lane_e2 = Lane(LaneTag(5,2),track_e2)
    connect!(lane_e1,lane_e2)
    connect!(lane_f1,lane_f2)
    connect!(lane_g,lane_f2)

    push!(road.segments,RoadSegment(4,[lane_g,lane_f1,lane_e1]))
    push!(road.segments,RoadSegment(5,[lane_f2,lane_e2]))

    return road
end

# function: Build roadways to test jumpy behavior in presence of discontinuity
"""
    function make_discont_roadway(;segment_length::Float64=20., separation::Float64=2.)
- Intended to test whether jumpy behavior is shown by vehicle

# Used by
`test_jumpy_vehicle`
"""
function make_discont_roadway_straight(;segment_length::Float64=20., separation::Float64=2.)
    track_a = gen_straight_curve(VecE2(0.,0.),VecE2(segment_length,0.),2)
    track_b = gen_straight_curve(VecE2(segment_length+separation,0.),VecE2(segment_length+separation+segment_length,0.),2)
    lane_a = Lane(LaneTag(1,1),track_a)
    lane_b = Lane(LaneTag(2,1),track_b)
    connect!(lane_a,lane_b)
    road = Roadway()
    push!(road.segments,RoadSegment(1,[lane_a]))
    push!(road.segments,RoadSegment(2,[lane_b]))
    return road
end

"""
    function make_discont_roadway_jagged(;segment_length::Float64=20., separation::Float64=2.)

- road with two horizontal segments but seperated in y to create a zig zag

# Examples
```julia
road_zigzag = make_discont_roadway_jagged()
```
"""
function make_discont_roadway_jagged(;segment_length::Float64=20., separation::Float64=2.)
    track_a = gen_straight_curve(VecE2(0.,0.),VecE2(segment_length,0.),2)
    track_b = gen_straight_curve(VecE2(segment_length+separation,10.),VecE2(segment_length+separation+segment_length,10.),2)
    lane_a = Lane(LaneTag(1,1),track_a)
    lane_b = Lane(LaneTag(2,1),track_b)
    connect!(lane_a,lane_b)
    road = Roadway()
    push!(road.segments,RoadSegment(1,[lane_a]))
    push!(road.segments,RoadSegment(2,[lane_b]))
    return road
end

# Function: highlight lanes by overlaying colors
"""
    function show_lane_overlays(road_ext::Roadway,traj_ext)
- Draw lane colored overlays on the road_ext

# Examples
```julia
show_lane_overlays(road_ext,traj_ext)
```
"""
function show_lane_overlays(road_ext::Roadway,traj_ext)
    lo_a = LaneOverlay(road_ext[LaneTag(1,1)],RGBA(0.,0.,1.,0.5))
    lo_b1 = LaneOverlay(road_ext[LaneTag(1,2)],RGBA(1.,0.,0.,0.5))
    lo_b2 = LaneOverlay(road_ext[LaneTag(2,1)],RGBA(0.,1.,0.,0.5))

    lo_c = LaneOverlay(road_ext[LaneTag(3,1)],RGBA(1.,1.,0.,0.5))
    lo_d = LaneOverlay(road_ext[LaneTag(3,2)],RGBA(0.,1.,1.,0.5))

    scene = traj_ext[1]
    return render(scene,road_ext,[lo_a,lo_b1,lo_b2,lo_c,lo_d])
end

# Special rendering to avoid lane boundary go till end
"""
struct MergingRoadway
- Wrapper around a roadway for special rendering avoiding lane boundary wierd in merge

# Example
```julia
road_ext = make_roadway_interaction_with_extensions()
roadway = MergingRoadway(road_ext)
render([roadway])
```
"""
struct MergingRoadway
    roadway::Roadway
end

"""
function AutomotiveVisualization.add_renderable!(rendermodel::RenderModel,mr::MergingRoadway)

- Dedicated render function for merging roadway
"""
function AutomotiveVisualization.add_renderable!(rendermodel::RenderModel,
    mr::MergingRoadway)
    #print("This is the (still code copied) rendering for merging lane boundary")
    lane_marking_width=0.15
    lane_dash_len=1.0
    lane_dash_spacing=2.0
    lane_dash_offset=0.00
    color_asphalt = colorant"0x708090"
    roadway = mr.roadway

    # render the asphalt between the leftmost and rightmost lane markers
    for seg in roadway.segments
        if !isempty(seg.lanes)
            laneR = seg.lanes[1]
            laneL = seg.lanes[end]

            pts = Array{Float64}(undef, 2, length(laneL.curve) + has_next(laneL) +
                                        length(laneR.curve) + has_next(laneR) +
                                        2*length(seg.lanes))
            pts_index = 0
            for pt in laneL.curve
                edgept = pt.pos + polar(laneL.width/2, pt.pos.θ + π/2)
                pts_index += 1
                pts[1, pts_index] = edgept.x
                pts[2, pts_index] = edgept.y
            end
            if has_next(laneL)
                pt = next_lane_point(laneL, roadway)
                edgept = pt.pos + polar(laneL.width/2, pt.pos.θ + π/2)
                pts_index += 1
                pts[1, pts_index] = edgept.x
                pts[2, pts_index] = edgept.y
            end
            for i in reverse(1:length(seg.lanes))
                lane = seg.lanes[i]
                if has_next(lane)
                    pt = next_lane_point(lane, roadway).pos
                else
                    pt = lane.curve[end].pos
                end
                pts_index += 1
                pts[1, pts_index] = pt.x
                pts[2, pts_index] = pt.y
            end

            if has_next(laneR)
                pt = next_lane_point(laneR, roadway)
                edgept = pt.pos + polar(laneR.width/2, pt.pos.θ - π/2)
                pts_index += 1
                pts[1, pts_index] = edgept.x
                pts[2, pts_index] = edgept.y
            end
            for j in length(laneR.curve) : -1 : 1
                pt = laneR.curve[j]
                edgept = pt.pos + polar(laneR.width/2, pt.pos.θ - π/2)
                pts_index += 1
                pts[1, pts_index] = edgept.x
                pts[2, pts_index] = edgept.y
            end
            for i in 1:length(seg.lanes)
                lane = seg.lanes[i]
                pt = lane.curve[1].pos
                pts_index += 1
                pts[1, pts_index] = pt.x
                pts[2, pts_index] = pt.y
            end

            add_instruction!(rendermodel, render_fill_region, (pts, color_asphalt))
        end
        # for lane in seg.lanes
        #     render!(rendermodel, lane, roadway)
        # end
    end

    # render the lane edges
    for seg in roadway.segments
        for lane in seg.lanes
            if lane.tag == LaneTag(1,1)
                # This is the top merge lane i.e. lane a. Make dashed line cut short
                curve_tweak = lane.curve[1:end-40]
                N = length(curve_tweak)
                halfwidth = lane.width/2

                pts_left = Array{Float64}(undef, 2, N)

                for (i,pt) in enumerate(curve_tweak)
                    p_left = pt.pos + polar(halfwidth, pt.pos.θ + π/2)
                    pts_left[1,i] = p_left.x
                    pts_left[2,i] = p_left.y
                end

                add_renderable!(rendermodel, lane.boundary_left, pts_left, lane_marking_width, lane_dash_len, lane_dash_spacing, lane_dash_offset)

            elseif lane.tag == LaneTag(4,1)
                # This is the bottom merge lane i.e. lane g maybe. Make dashed line cut short
                curve_tweak = lane.curve[1:end-80]
                N = length(curve_tweak)
                halfwidth = lane.width/2

                pts_left = Array{Float64}(undef, 2, N)

                for (i,pt) in enumerate(curve_tweak)
                    p_left = pt.pos + polar(halfwidth, pt.pos.θ + π/2)
                    pts_left[1,i] = p_left.x
                    pts_left[2,i] = p_left.y
                end

                add_renderable!(rendermodel, lane.boundary_left, pts_left, lane_marking_width, lane_dash_len, lane_dash_spacing, lane_dash_offset)

            else

                N = length(lane.curve)
                halfwidth = lane.width/2

                # always render the left lane marking
                pts_left = Array{Float64}(undef, 2, N)

                for (i,pt) in enumerate(lane.curve)

                    p_left = pt.pos + polar(halfwidth, pt.pos.θ + π/2)

                    pts_left[1,i] = p_left.x
                    pts_left[2,i] = p_left.y
                end
                if has_next(lane)
                    lane2 = next_lane(lane, roadway)
                    pt = lane2.curve[1]
                    p_left = pt.pos + polar(lane2.width/2, pt.pos.θ + π/2)
                    pts_left = hcat(pts_left, [p_left.x, p_left.y])
                end

                add_renderable!(rendermodel, lane.boundary_left, pts_left, lane_marking_width, lane_dash_len, lane_dash_spacing, lane_dash_offset)

                # only render the right lane marking if this is the first lane
                if lane.tag.lane == 1
                    pts_right = Array{Float64}(undef, 2, N)

                    for (i,pt) in enumerate(lane.curve)

                        p_right = pt.pos - polar(halfwidth, pt.pos.θ + π/2)

                        pts_right[1,i] = p_right.x
                        pts_right[2,i] = p_right.y
                    end

                    if has_next(lane)
                        lane2 = next_lane(lane, roadway)
                        pt = lane2.curve[1]
                        p_right = pt.pos - polar(lane2.width/2, pt.pos.θ + π/2)
                        pts_right = hcat(pts_right, [p_right.x, p_right.y])
                    end

                    add_renderable!(rendermodel, lane.boundary_right, pts_right, lane_marking_width, lane_dash_len, lane_dash_spacing, lane_dash_offset)
                end
            end # Ends the if else to check whether merging lane
        end
    end

    return rendermodel
end

max_stars_repo_name: sisl/AutomotiveInteraction

max_issues_repo_name: sisl/AutomotiveInteraction

max_forks_repo_name: sisl/AutomotiveInteraction


text:
"""
    abstract type AbstractConfig{T <: SSDFloat} end

Supertype of all detector/world/object configs.

User defined geometries must be subtype of `AbstractConfig{T}`.

There are a few functions which must be defined for a user config, e.g. `struct UserConfig{T} <: AbstractConfig{T}`:

For cylindrical grids:

* in(pt::cylindrical{T}, config::UserConfig{T})::Bool where {T <: SSDFloat}
* Grid(config::UserConfig{T})::Grid{T, 3, :cylindrical} where {T <: SSDFloat}
* get\\_ρ\\_and\\_ϵ(pt::cylindrical{T}, config::UserConfig{T})::Tuple{T, T} where {T <: SSDFloat}
* set\\_pointtypes\\_and\\_fixed\\_potentials!(pointtypes::Array{PointType, 3}, potential::Array{T, 3},
        grid::Grid{T, 3, :cylindrical}, config::UserConfig{T}; weighting\\_potential\\_channel\\_idx::Union{Missing, Int} = missing)::Nothing where {T <: SSDFloat}

For cartesian grids:

* in(pt::StaticVector{3, T}, config::UserConfig{T})::Bool
* Grid(config::UserConfig{T})::Grid{T, 3, :cartesian} where {T <: SSDFloat}
* get\\_ρ\\_and\\_ϵ(pt::StaticVector{3, T}, config::UserConfig{T})::Tuple{T, T} where {T <: SSDFloat}
* set\\_pointtypes\\_and\\_fixed\\_potentials!(pointtypes::Array{PointType, 3}, potential::Array{T, 3},
        grid::Grid{T, 3, :cartesian}, config::UserConfig{T}; weighting\\_potential\\_channel\\_idx::Union{Missing, Int} = missing)::Nothing where {T <: SSDFloat}

"""
abstract type AbstractConfig{T <: SSDFloat} end


text:
using SGSDist
using Documenter

DocMeta.setdocmeta!(SGSDist, :DocTestSetup, :(using SGSDist); recursive=true)

makedocs(
    sitename = "SGSDist.jl",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Robert West",
    modules = [SGSDist],
    pages = [
        "Home" => "index.md",
        "Examples" => "examples.md",
        "Library" => "library.md",
        "Index" => "api.md"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/brwst/SGSDist.jl.git",
)


text:
include("src/Blake3Hash.jl")

using Profile
using .Blake3Hash

DATA_SIZE = 2*8192

function run()
    h = Blake3Hash.Simd128.Blake3Ctx()
    d = rand(UInt8, 2*8192)
    o = Vector{UInt8}(undef, 32)

    Blake3Hash.Simd128.update!(h, d)
    Blake3Hash.Simd128.digest(h, o)
end

run()
Profile.clear_malloc_data()
run()


text:
using DataFrames, MLDataUtils
using Clustering, Distances
using CSV
using Random
using Logging

# Set up Logging - we recommend to use this command to avoid package warnings during the model training process.
logger = Logging.SimpleLogger(stderr,Logging.Warn);
global_logger(logger);

#### Set parameters for the learners
cr = :dunnindex
method = "ICOT_local"
warm_start = :none;
geom_search = false
threshold = 0.0
seed = 1
gridsearch = false
num_tree_restarts = 100
complexity_c = 0.0
min_bucket = 10
maxdepth = 5

###### Step 1: Prepare the data
# Read the data - recommend the use of the (deprecated) readtable() command to avoid potential version conflicts with the CSV package.
data = readtable("../data/ruspini.csv");

# Convert the dataset to a matrix
data_array = convert(Matrix{Float64}, data);
# Get the number of observations and features
n, p = size(data_array)
data_t = data_array';

##### Step 2: Fit K-means clustering on the dataset to generate a warm-start for ICOT
#Fix the seed
Random.seed!(seed);

# The ruspini dataset has pre-defined clusters, which we will use to select the cluster count (K) for the K-means algorithm.
# In an unsupervised setting (with no prior-known K), the number of clusters for K means can be selected using the elbow method.
K = length(unique(data_array[:,end]))

# Run k-means and save the assignments
kmeans_result = kmeans(data_t, K);
assignment = kmeans_result.assignments;

data_full = DataFrame(hcat(data, assignment, makeunique=true));
names!(data_full, [:x1, :x2, :true_labels, :kmean_assign]);

# Prepare data for ICOT: features are stored in the matrix X, and the warm-start labels are stored in y
X = data_full[:,1:2]; y = data_full[:,:true_labels];

##### Step 3a. Before running ICOT, start by testing the IAI license
lnr_oct = ICOT.IAI.OptimalTreeClassifier(localsearch = false, max_depth = maxdepth,
													 minbucket = min_bucket,
													 criterion = :misclassification
													 )
grid = ICOT.IAI.GridSearch(lnr_oct)
ICOT.IAI.fit!(grid, X, y)
ICOT.IAI.showinbrowser(grid.lnr)

##### Step 3b. Run ICOT

# Run ICOT with no warm-start:
warm_start= :none
lnr_ws_none = ICOT.InterpretableCluster(ls_num_tree_restarts = num_tree_restarts, ls_random_seed = seed, cp = complexity_c, max_depth = maxdepth,
	minbucket = min_bucket, criterion = cr, ls_warmstart_criterion = cr, kmeans_warmstart = warm_start,
	geom_search = geom_search, geom_threshold = threshold);
run_time_icot_ls_none = @elapsed ICOT.fit!(lnr_ws_none, X, y);

ICOT.showinbrowser(lnr_ws_none)

score_ws_none = ICOT.score(lnr_ws_none, X, y, criterion=:dunnindex);
score_al_ws_none = ICOT.score(lnr_ws_none, X, y, criterion=:silhouette);

# Run ICOT with an OCT warm-start: fit an OCT as a supervised learning problem with labels "y" and use this as the warm-start
warm_start= :oct
lnr_ws_oct = ICOT.InterpretableCluster(ls_num_tree_restarts = num_tree_restarts, ls_random_seed = seed, cp = complexity_c, max_depth = maxdepth,
	minbucket = min_bucket, criterion = cr, ls_warmstart_criterion = cr, kmeans_warmstart = warm_start,
	geom_search = geom_search, geom_threshold = threshold);
run_time_icot_ls_oct = @elapsed ICOT.fit!(lnr_ws_oct, X, y);

score_ws_oct = ICOT.score(lnr_ws_oct, X, y, criterion=:dunnindex);
score_al_ws_oct = ICOT.score(lnr_ws_oct, X, y, criterion=:silhouette);

##### Comments
# Note that in this example, the OCT tree and ICOT (OCT warm-start, or no warm-start) result in the same solution.
# This is not generally the case, but can occur when the data is easily separated (as in the ruspini dataset)

# The score printed in the browser view is negative to reflect the formulation as a minimization problem.
# The positive score returned by the ICOT.score() functions reflect the "correct" interpretation of the score,
# in which we seek to maximize the criterion (with a maximum value of 1).

# For larger datasets, we recommend setting warm_start = :oct and threshold = 0.99 to improve the solve time.

max_stars_repo_name: agniorf/ICOT

max_issues_repo_name: agniorf/ICOT

max_forks_repo_name: agniorf/ICOT


text:
Pkg.add("PyPlot")
using PyPlot

include("../src/LDPCsim.jl")
using LDPCsim

ldpcSim(1000,[-2.0:1:10],[50])
savefig("ldpcVsSNR_M50rHalf.png") # run from root directory

ldpcSim(1000,[3.5],[4,8,16,32,64,128])
savefig("ldpcVsM_snr3p5_rHalf.png") # run from root directory

max_stars_repo_name: christianpeel/LDPCsim

max_issues_repo_name: christianpeel/LDPCsim

max_forks_repo_name: christianpeel/LDPCsim


text:
using CoverageBase
using Coverage
using Base.Test

@test !isempty(CoverageBase.julia_top())
@test !isempty(testnames())
@test_throws ErrorException runtests(["junk"])
runtests([joinpath(CoverageBase.julia_top(), "test", "goto")])


text:
#### HOW TO CREATE MARKDOWN FILES ####

using Literate, PyPlot
cd(@__DIR__)
pygui(false) # no interactive plots

tmp_outputdir = "tmp"
mkpath(tmp_outputdir)
pdfdir = joinpath(@__DIR__, "../../../products/additional_figures/")
mkpath(pdfdir)

# Run Literate.jl over the files:
# (execute=true produces plots)
Literate.markdown("raw_data.jl", tmp_outputdir; execute=true, documenter=false)
Literate.markdown("calibration_plots.jl", tmp_outputdir; execute=true, documenter=false)
Literate.markdown("derived_quantities_extra.jl", tmp_outputdir; execute=true, documenter=false)
Literate.markdown("other_figures.jl", tmp_outputdir; execute=true, documenter=false)

# to convert .md files to pdfs, linux command, e.g.:
cd(tmp_outputdir)
for fl in readdir()
    if splitext(fl)[2]==".md"
        fll = splitext(fl)[1]
        try
            run(`pandoc -s -V geometry:margin=1in -o $fll.pdf $fll.md`)
            cp("$fll.pdf", joinpath(pdfdir, "$fll.pdf"), force=true)
        catch e
            println(e)
            println("Could not convert markdown-file $fll.md to pdf")
        end
    end
end

cd(@__DIR__)
#
rm(tmp_outputdir, recursive=true)


text:
@testset "$TEST $G" begin

##################################################################
#
#   Maximal cliques of undirected graph
#   Derived from Graphs.jl: https://github.com/julialang/Graphs.jl
#
##################################################################

if !@isdefined(setofsets)
    function setofsets(array_of_arrays)
        Set(map(Set, array_of_arrays))
    end

    function test_cliques(graph, expected)
        # Make test results insensitive to ordering
        setofsets(maximal_cliques(graph)) == setofsets(expected)
    end
end

g = G(3)
add_edge!(g, 1, 2)
@test test_cliques(g, Array[[1,2], [3]])
add_edge!(g, 2, 3)
@test test_cliques(g, Array[[1,2], [2,3]])

# Test for "pivotdonenbrs not defined" bug
h = G(6)
add_edge!(h, 1, 2)
add_edge!(h, 1, 3)
add_edge!(h, 1, 4)
add_edge!(h, 2, 5)
add_edge!(h, 2, 6)
add_edge!(h, 3, 4)
add_edge!(h, 3, 6)
add_edge!(h, 5, 6)

@test maximal_cliques(h) != []

# test for extra cliques bug

h = G(7)
add_edge!(h,1,3)
add_edge!(h,2,6)
add_edge!(h,3,5)
add_edge!(h,3,6)
add_edge!(h,4,5)
add_edge!(h,4,7)
add_edge!(h,5,7)
@test test_cliques(h, Array[[7,4,5], [2,6], [3,5], [3,6], [3,1]])

end # testset