text:
# StaticCPD
let
    df = DataFrame(a=randn(100))
    cpd = fit(StaticCPD{Normal}, df, :a)
    @test name(cpd) == :a
    @test parentless(cpd)
    @test parents(cpd) == NodeName[]
    @test disttype(cpd) <: Normal
    @test nparams(cpd) == 2
    @test isa(cpd(), Normal)
    @test pdf(cpd, :a=>0.5) > 0.2
    @test pdf(cpd, :a=>0.0) > pdf(cpd, :a=>0.5)
    @test logpdf(cpd, :a=>0.5) > log(0.2)
end

# CategoricalCPD
let

    # no parents
    let
        df = DataFrame(a=[1,2,1,2,3])
        cpd = fit(DiscreteCPD, df, :a)

        @test name(cpd) == :a
        @test parentless(cpd)
        @test nparams(cpd) == 3

        d = cpd(Assignment())
        @test isa(d, Categorical) && isa(d, disttype(cpd))
        @test pdf(d, 1) == 0.4
        @test pdf(d, 2) == 0.4
        @test pdf(d, 3) == 0.2

        df = DataFrame(a=randn(100))
        cpd = fit(CategoricalCPD{Normal}, df, :a)

        @test isa(cpd(), disttype(cpd))
        @test pdf(cpd, :a=>0.5) > 0.2
        @test pdf(cpd, :a=>0.0) > pdf(cpd, :a=>0.55)

        cpd = fit(DiscreteCPD, DataFrame(a=[1,2,1,2,2]), :a, ncategories=3)
        @test pdf(cpd, :a=>1) == 0.4
        @test pdf(cpd, :a=>2) == 0.6
        @test pdf(cpd, :a=>3) == 0.0
        @test nparams(cpd) == 3
    end

    # with parents
    let
        # Example with Categorical
        df = DataFrame(a=[1,2,1,2,3], b=[1,1,2,1,2])
        cpd = fit(DiscreteCPD, df, :b, [:a])

        @test name(cpd) == :b
        @test parents(cpd) == [:a]
        @test !parentless(cpd)
        @test nparams(cpd) == 6

        @test cpd(:a=>1).p == [0.5,0.5]
        @test cpd(:a=>2).p == [1.0,0.0]
        @test cpd(:a=>3).p == [0.0,1.0]

        cpd = fit(DiscreteCPD, df, :b, [:a], parental_ncategories=[3], target_ncategories=5)
        @test nparams(cpd) == 15

        # Example with Bernoulli and more than one parent
        df = DataFrame(a=[   1,    1,    1,    1,    2,    2,    2,    2],
                       b=[   1,    1,    2,    2,    1,    1,    2,    2],
                       c=[true, true,false,false, true,false,false, true])
        cpd = fit(CategoricalCPD{Bernoulli}, df, :c, [:a, :b])

        @test nparams(cpd) == 4

        @test isa(cpd(Assignment(:a=>1, :b=>1)), disttype(cpd))
        @test cpd(Assignment(:a=>1, :b=>1)).p == 1.0
        @test cpd(Assignment(:a=>1, :b=>2)).p == 0.0
        @test cpd(Assignment(:a=>2, :b=>1)).p == 0.5
        @test cpd(Assignment(:a=>2, :b=>2)).p == 0.5
    end
end

# Linear Gaussian
let

    # no parents
    let
        df = DataFrame(a=randn(100))
        cpd = fit(LinearGaussianCPD, df, :a, min_stdev=0.0)

        @test name(cpd) == :a
        @test parents(cpd) == NodeName[]
        @test parentless(cpd)
        @test disttype(cpd) <: Normal
        @test nparams(cpd) == 2
        @test pdf(cpd, :a=>0.5) > 0.2
        @test pdf(cpd, :a=>0.0) > pdf(cpd, :a=>0.55)
    end

    # with parents
    let
        a = randn(1000)
        b = randn(1000) .+ 2*a .+ 1

        df = DataFrame(a=a, b=b)
        cpd = fit(LinearGaussianCPD, df, :b, [:a])

        @test !parentless(cpd)
        @test parents(cpd) == [:a]
        @test nparams(cpd) == 3

        p = cpd(:a=>0.0)
        @test isapprox(p.μ, 1.0, atol=0.25)
        @test isapprox(p.σ, 2.0, atol=0.50)

        p = cpd(:a=>1.0)
        @test isapprox(p.μ, 3.0, atol=0.25)
        @test isapprox(p.σ, 2.0, atol=0.50)

        cpd = fit(LinearGaussianCPD, df, :b, [:a], min_stdev=10.0)
        @test cpd(:a=>1.0).σ == 10.0
    end
end

# ConditionalLinearGaussianCPD
let

    # no parents
    let
        df = DataFrame(a=randn(10))
        cpd = fit(ConditionalLinearGaussianCPD, df, :a)

        @test name(cpd) == :a
        @test parentless(cpd)
        @test nparams(cpd) == 2

        d = cpd()
        @test isa(d, Normal) && isa(d, disttype(cpd))
    end

    # with parents
    let
        df = DataFrame(a=[  1,   1,   1,   1,     2,   2,   2,   2],
                       b=[0.5, 1.0, 1.5, 1.0,   2.5, 3.0, 3.5, 3.0],
                       c=[0.55,1.05,1.53,1.02,  2.52,3.01,3.55,3.03])
        cpd = fit(ConditionalLinearGaussianCPD, df, :c, [:a, :b])

        @test name(cpd) == :c
        @test parents(cpd) == [:a, :b]
        @test !parentless(cpd)
        @test nparams(cpd) == 6

        d = cpd(Assignment(:a=>1, :b=>0.5))
        @test isapprox(d.μ, 0.5475, atol=0.001)
        @test isapprox(d.σ, 0.4003, atol=0.001)

        d = cpd(Assignment(:a=>2, :b=>3.0))
        @test isapprox(d.μ, 3.027, atol=0.001)
        @test isapprox(d.σ, 0.421, atol=0.001)
    end
end

# FunctionalCPD
let
    bn2 = BayesNet()
    push!(bn2, StaticCPD(:sighted, NamedCategorical([:bird, :plane, :superman], [0.40, 0.55, 0.05])))
    push!(bn2, FunctionalCPD{Bernoulli}(:happy, [:sighted], a->Bernoulli(a == :superman ? 0.95 : 0.2)))

    @test name(get(bn2, 2)) == :happy
    @test parents(get(bn2, 2)) == [:sighted]

    val = rand(bn2)
    @test in(val[:sighted], [:bird, :plane, :superman])
    @test in(val[:happy], [0,1])

    # named cat
    ncat = get(bn2, 1).d
    show(IOBuffer(), ncat)
    @test ncategories(ncat) == 3
    @test probs(ncat) == [0.40, 0.55, 0.05]
    @test params(ncat) == ([0.40, 0.55, 0.05],)
    @test pdf(ncat, :bird) == 0.4
    @test isapprox(logpdf(ncat, :bird), log(0.4))
    show(IOBuffer(), sampler(ncat))

    fcpd = FunctionalCPD{Normal}(:sad, a->Normal(0.0, 1.0))
    @test name(fcpd) == :sad
    @test isempty(parents(fcpd))
    @test fcpd(Assignment()) == Normal(0.0, 1.0)
end

# ParentFunctionalCPD
let
    a = StaticCPD(:a, Bernoulli(0.5))
    b = StaticCPD(:b, Bernoulli(0.6))
    p = [:a,:b]
    c = ParentFunctionalCPD{Bernoulli}(:c, p, (seq,par)->begin
                Bernoulli(mean(seq[k] for k in par))
            end
        )
    bn = BayesNet()
    push!(bn, a)
    push!(bn, b)
    push!(bn, c)
    @test mean(rand(bn, 20,:a=>0)[:c]) <= 0.6
    @test mean(rand(bn, 20, :a=>1, :b=>1)[:c]) == 1;
end


text:
using LMCLUS
using Test
using Statistics
using DelimitedFiles
using Random

@testset "Clustering" begin

	# Test histogram size calculations
	Random.seed!(24975754857)
	x = rand(1000)
	@test LMCLUS.gethistogrambins(x, 0.0, 10, 5) == 10
	@test LMCLUS.gethistogrambins(x, 0.0, 10, 15) == 15
	@test LMCLUS.gethistogrambins(x, 0.01, 10, 15) == 233

	# Test distance calculation
	basis = reshape([1.,0.,0.,0.,1.,0.],3,2)
	point = [1.0,1.0,1.0]
	@test distance_to_manifold(point, basis)[1] ≈ 1.0
	@test distance_to_manifold(point, basis)[2] ≈ sqrt(2.0)

	# Test sampling parameters
	p = LMCLUS.Parameters(3)
	p.number_of_clusters = 1
	@test LMCLUS.sample_quantity(1,3,1000,p,1) == 1
	p.number_of_clusters = 2
	p.sampling_heuristic = 1
	@test LMCLUS.sample_quantity(1,3,1000,p,1) == 13
	p.sampling_heuristic = 2
	@test LMCLUS.sample_quantity(1,3,1000,p,1) == 10
	p.sampling_heuristic = 3
	@test LMCLUS.sample_quantity(1,3,1000,p,1) == 10

	# Test clustering
	p = LMCLUS.Parameters(5)
	p.random_seed = 4572489057

	testDataFile = joinpath(dirname(@__FILE__),"testData")
	ds = readdlm(testDataFile, ',')
	data = ds[:,1:end-1]'

	# test separation calculations
	M = [297.654, -183.908, -164.718, -339.345, -53.5994, -142.06, -207.939, -180.871, 469.81, 190.212]
	B = [0.0114666 0.20954 -0.371882 0.0527596 0.429496 -0.0887774 0.0872738 0.734395 0.121809 -0.246463]' |> Matrix
	s = LMCLUS.find_separation(data, M, B, p, MersenneTwister())
	@test typeof(s) == LMCLUS.Separation
	@test criteria(s) ≈ 7.285358462818518
	@test threshold(s) ≈ 852 atol=1.0

	# run clustering
	res = lmclus(data, p)
	@test nclusters(res) >= 3
	@test sum(counts(res)) == size(ds, 1)

	cnts = counts(res)
	nc = nclusters(res)
	@testset "Label Match" for idxs in [(i,j) for i in 1:nc, j in 1:nc if i < j]
		i = idxs[1]
		j = idxs[2]
		@test length(symdiff(points(manifold(res,i)), points(manifold(res, j)))) == cnts[i] + cnts[j]
	end

	# Sampling
	p.histogram_sampling = true
	p.number_of_clusters = 3
	res = lmclus(data, p)
	@test nclusters(res) >= 3
	p.histogram_sampling = false

	# RNG seed
	p.random_seed = 0
	res = lmclus(data, p)
	@test nclusters(res) >= 3

	# MDL
	p.mdl_algo = LMCLUS.MDL.SizeIndependent
	p.mdl = true
	res = lmclus(data, p)
	@test nclusters(res) >= 3
	p.mdl_algo = LMCLUS.MDL.OptimalQuant
	p.mdl = false

	# adjust cluster bases
	p.basis_alignment = true
	res = lmclus(data, p)
	@test nclusters(res) >= 3
	@test threshold(manifold(res,1))[1] > 0.0
	@test threshold(manifold(res,1))[2] == Inf
	p.basis_alignment = false

	# calculate cluster bounsds
	p.basis_alignment = true
	p.bounded_cluster = true
	res = lmclus(data, p)
	@test nclusters(res) >= 3
	@test threshold(manifold(res,1))[1] > 0.0
	@test threshold(manifold(res,1))[2] > 0.0
	p.bounded_cluster = false
	p.basis_alignment = false

	# Otsu
	p.sep_algo = LMCLUS.Otsu
	p.min_cluster_size = 200
	p.basis_alignment = true
	res = lmclus(data, p)
	@test nclusters(res) >= 3
	p.sep_algo = LMCLUS.Kittler
	p.min_cluster_size = 20

	# Iterative refinement
	dfun = (X,m)  -> distance_to_manifold(X, mean(m), projection(m))
	efun = (X,ms) -> LMCLUS.MDL.calculate(LMCLUS.MDL.OptimalQuant, ms, X, 32, 16)
	res2 = refine(res, data, dfun, efun, bounds=true)
	@test nclusters(res) >= nclusters(res2)
	@test threshold(manifold(res2,1))[1] > 0.0
	@test threshold(manifold(res2,1))[2] > 0.0
	# no changes to original clustering
	@testset for (m1, m2) in zip(manifolds(res), manifolds(res))
		@test all(points(m1) .== points(m2))
	end

end

ext: jl


text:
function CategoricalPool{S, T <: Integer, R <: Integer}(index::Vector{S},
                                           invindex::Dict{S, T},
                                           order::Vector{R},
                                           ordered::Bool=false)
    invindex = convert(Dict{S, R}, invindex)
    CategoricalPool{S, R, CategoricalValue{S, R}}(index, invindex, order, ordered)
end

@compat (::Type{CategoricalPool{T, R}}){T, R}(ordered::Bool=false) =
    CategoricalPool(T[], Dict{T, R}(), R[], ordered)
@compat (::Type{CategoricalPool{T}}){T}(ordered::Bool=false) =
    CategoricalPool(T[], Dict{T, DefaultRefType}(), DefaultRefType[], ordered)

@compat function (::Type{CategoricalPool{T, R}}){T, R}(index::Vector,
                                                       ordered::Bool=false)
    invindex = buildinvindex(index, R)
    order = Vector{R}(1:length(index))
    CategoricalPool(index, invindex, order, ordered)
end

function CategoricalPool(index::Vector, ordered::Bool=false)
    invindex = buildinvindex(index)
    order = Vector{DefaultRefType}(1:length(index))
    return CategoricalPool(index, invindex, order, ordered)
end

function CategoricalPool{S, R <: Integer}(invindex::Dict{S, R},
                                          ordered::Bool=false)
    index = buildindex(invindex)
    order = Vector{DefaultRefType}(1:length(index))
    return CategoricalPool(index, invindex, order, ordered)
end

# TODO: Add tests for this
function CategoricalPool{S, R <: Integer}(index::Vector{S},
                                          invindex::Dict{S, R},
                                          ordered::Bool=false)
    order = Vector{DefaultRefType}(1:length(index))
    return CategoricalPool(index, invindex, order, ordered)
end

function CategoricalPool{T}(index::Vector{T},
                            levels::Vector{T},
                            ordered::Bool=false)
    invindex = buildinvindex(index)
    order = buildorder(invindex, levels)
    return CategoricalPool(index, invindex, order, ordered)
end

function CategoricalPool{S, R <: Integer}(invindex::Dict{S, R},
                                          levels::Vector{S},
                                          ordered::Bool=false)
    index = buildindex(invindex)
    order = buildorder(invindex, levels)
    return CategoricalPool(index, invindex, order, ordered)
end

Base.convert(::Type{CategoricalPool}, pool::CategoricalPool) = pool
Base.convert{T}(::Type{CategoricalPool{T}}, pool::CategoricalPool{T}) = pool
Base.convert{T, R}(::Type{CategoricalPool{T, R}}, pool::CategoricalPool{T, R}) = pool

Base.convert{S, T, R}(::Type{CategoricalPool{S}}, pool::CategoricalPool{T, R}) =
    convert(CategoricalPool{S, R}, pool)
Base.convert{T, R}(::Type{CategoricalPool}, pool::CategoricalPool{T, R}) =
    convert(CategoricalPool{T, R}, pool)

function Base.convert{S, R}(::Type{CategoricalPool{S, R}}, pool::CategoricalPool)
    indexS = convert(Vector{S}, pool.index)
    invindexS = convert(Dict{S, R}, pool.invindex)
    order = convert(Vector{R}, pool.order)
    return CategoricalPool(indexS, invindexS, order)
end

function Base.show{T, R}(io::IO, pool::CategoricalPool{T, R})
    @printf(io, "%s{%s,%s}([%s])", typeof(pool).name, T, R,
                join(map(repr, levels(pool)), ","))

    pool.ordered && print(io, " with ordered levels")
end

Base.length(pool::CategoricalPool) = length(pool.index)

Base.getindex(pool::CategoricalPool, i::Integer) = pool.valindex[i]
Base.get(pool::CategoricalPool, level::Any) = pool.invindex[level]

function Base.get!{T, R, V}(pool::CategoricalPool{T, R, V}, level)
    get!(pool.invindex, level) do
        i = length(pool) + 1
        push!(pool.index, level)
        push!(pool.order, i)
        push!(pool.levels, level)
        push!(pool.valindex, V(i, pool))
        i
    end
end

Base.push!(pool::CategoricalPool, level) = (get!(pool, level); pool)

# TODO: optimize for multiple additions
function Base.append!(pool::CategoricalPool, levels)
    for level in levels
        push!(pool, level)
    end
    return pool
end

function Base.delete!{S, R, V}(pool::CategoricalPool{S, R, V}, levels...)
    for level in levels
        levelS = convert(S, level)
        if haskey(pool.invindex, levelS)
            ind = pool.invindex[levelS]
            delete!(pool.invindex, levelS)
            splice!(pool.index, ind)
            ord = splice!(pool.order, ind)
            splice!(pool.levels, ord)
            splice!(pool.valindex, ind)
            for i in ind:length(pool)
                pool.invindex[pool.index[i]] -= 1
                pool.valindex[i] = V(i, pool)
            end
            for i in 1:length(pool)
                pool.order[i] > ord && (pool.order[i] -= 1)
            end
        end
    end
    return pool
end

function levels!{S, R, V}(pool::CategoricalPool{S, R, V}, newlevels::Vector)
    if !allunique(newlevels)
        throw(ArgumentError(string("duplicated levels found in newlevels: ",
                                   join(unique(filter(x->sum(newlevels.==x)>1, newlevels)), ", "))))
    end

    n = length(newlevels)

    # No deletions: can preserve position of existing levels
    if issubset(pool.index, newlevels)
        append!(pool, setdiff(newlevels, pool.index))
    else
        empty!(pool.invindex)
        resize!(pool.index, n)
        resize!(pool.valindex, n)
        resize!(pool.order, n)
        resize!(pool.levels, n)
        for i in 1:n
            v = newlevels[i]
            pool.index[i] = v
            pool.invindex[v] = i
            pool.valindex[i] = V(i, pool)
        end
    end

    buildorder!(pool.order, pool.invindex, newlevels)
    for (i, x) in enumerate(pool.order)
        pool.levels[x] = pool.index[i]
    end
    return newlevels
end

index(pool::CategoricalPool) = pool.index
levels(pool::CategoricalPool) = pool.levels
order(pool::CategoricalPool) = pool.order

ordered(pool::CategoricalPool) = pool.ordered
ordered!(pool::CategoricalPool, ordered) = pool.ordered = ordered


text:
"""
Cell linked list splits the vector space a regular grid with edge length `EPS`.

# Member variables
- `data::Array{T, 2} where {T <: Real}`: an array holding the original data
- `grid::Array{LinkedList{T}} where {T <: Int}`: the gridded space
- `EPS::Real`: the edge length of each grid cube
- `maxs::Vector{T} where {T <: Real}`: stores the maximal coordinate over all data
- `mins::Vector{T} where {T <: Real}`: stores the minimal coordinate over all data
- `prevquery::Real`: skip `cll_neighbor_cells` in `cll_inrange` if `Int(ceil(dist / edge)) == Int(ceil(prevquery / edge))`
- `dirs::Vector{T}`: cached vector of neighbor cells
"""
struct CellLinkedList
  data::Array{T, 2} where {T <: Real}
  grid::Array{LinkedList{T}} where {T <: Int}
  EPS::Real
  maxs::Vector{T} where {T <: Real}
  mins::Vector{T} where {T <: Real}
  prevquery::Real
  dirs::Vector{Int}
end

"""
Calculates the amount of cells in the grid.
A cell linked list with the resulting amount of cells will be able to hold any vector `v` where `v .<= maxs` and `v .>= mins`.
"""
function cll_cubecount(maxs::Vector{T}, mins::Vector{T}, EPS::Real) where {T <: Real}
  size(maxs, 1) != size(mins, 1) && throw(DimensionMismatch())
  return [Int(floor(maxs[i] / EPS) - floor(mins[i] / EPS)) for i in 1:size(maxs, 1)] + 1
end

cllCubeCount(a...) = warn("`cllCubeCount` is deprecated, use `cll_cubecount` instead.")

"""
Creates a grid of empty linked lists.
"""
function cll_initgrid(cells::Vector{T}) where {T <: Integer}
  return fill(nil(T), Tuple(cells))
end

cllCreateEmptyGrid(a...) = warn("`cllCreateEmptyGrid` is deprecated, use `cll_initgrid` instead.")

"""
Finds point's index in the grid of `cll`.
"""
function cll_index(cll::CellLinkedList, pt::Point)
  return cll_index(cll.mins, cll.EPS, pt)
end

cllIndex(a...) = warn("`cllIndex` is deprecated, use `cll_index` instead.")

"""
Finds point's index in the grid given by `mins` and `EPS`.
"""
function cll_index(mins::Vector{T}, EPS::Real, pt::Point) where {T <: Real}
  return Int.(floor.((pt - mins) / EPS)) + 1
end

"""
Some operations on the CLL might increase the size of the gridded space.
Instead of recreating a new CLL, we reuse the already created lists assuming that `EPS` hasn't changed.
"""
function cll_reusegrid!(cll::CellLinkedList, newgrid::Array{LinkedList{T}}) where {T <: Integer} # FIXME: Integer?
  for l in cll.grid
    l == nil() && continue
    pos = cll_index(cll.mins, cll.EPS, cll.data[:, head(l)])
    newgrid[pos...] = l
  end
  return nothing
end

cllCreateEmptyGrid!(a...) = warn("`cllCreateEmptyGrid` is deprecated, use `cll_reusegrid` instead.")

"""
Each point is in a single column.

# Arguments
- `data::Array{T} where {T <: Real}`: a 2d array with data, where each vector is in a single column,
- `EPS::Real`: the edge length of the regular grid.
"""
function CellLinkedList(data::Array{T, 2}, EPS::Real) where {T <: Real}
  dim = size(data, 1)
  cnt = size(data, 2)
  maxs = [maximum(data[i, :]) for i in 1:dim]
  mins = [minimum(data[i, :]) for i in 1:dim]
  grid::Array{LinkedList{Int}, dim} = cll_initgrid(cll_cubecount(maxs, mins, EPS))

  @inbounds for i in 1:cnt
    cur = data[:, i]
    pos = cll_index(mins, EPS, cur)
    grid[pos...] = cons(i, grid[pos...])
  end

  return CellLinkedList(data, grid, EPS, maxs, mins, 0., Vector{Int}())
end

"""
Adds a point to a CLL.
"""
function cll_add!(cll::CellLinkedList, pt::Point)
  size(pt, 1) != size(cll.data, 1) && throw(DimensionMismatch())
  newmax = max.(pt, cll.maxs)
  newmin = min.(pt, cll.mins)
  if newmax != cll.maxs || newmin != cll.mins
    newgrid::Array{LinkedList{Int}, size(pt, 1)} = cll_initgrid(cll_cubecount(newmax, newmin, cll.EPS))
    cll_reusegrid!(cll, newgrid)
    cll.grid = newgrid
    cll.maxs = newmax
    cll.mins = newmin
  end
  pos = cll_index(cll, pt)
  cll.data = hcat(cll.data, pt)
  cll.grid[pos...] = cons(size(cll.data, 2), cll.grid[pos...])
  return nothing
end

cllAdd!(a...) = warn("`cllAdd!` is deprecated, use `cll_add` instead.")

"""
Removes a point from the its cell.  It remains in `data`.
"""
function cll_removeindex!(cll::CellLinkedList, idx::Integer, pos)
  h = cll.grid[pos...]
  todelete = tail(h)
  if head(h) == idx
    cll.grid[pos...] = todelete
  else
    while todelete != nil()
      if head(todelete) == idx
        h.tail = tail(todelete)
        break
      else
        h = tail(h)
        todelete = tail(todelete)
      end
    end
  end
  return nothing
end

cllRemoveIdx!(a...) = warn("`cllRemoveIdx!` is deprecated, use `cll_removeindex!` instead.")

"""
Only removes it from the linked list, because array deletions are expensive.
Maybe add a deletion counter and perform a deletion when `counter > sqrt(length(data))`
"""
function cll_remove!(cll::CellLinkedList, idx::Integer)
  if idx < 1 || idx > size(data, 2)
    throw(BoundsError(cll.data, idx))
  end

  pos = cll.index(data[:, idx])
  cll_removeindex!(cll, idx, pos)
  return nothing
end

cllRemove!(a...) = warn("`cllRemove!` is deprecated, use `cll_remove!` instead")

"""
Changes `idx`th data to new coordinates.
"""
function cll_modify!(cll::CellLinkedList, idx::Integer, newpt::Point)
  oldPos = cll_index(cll, cll.data[:, idx])
  newPos = cll_index(cll, newpt)
  if oldPos == newPos
    return nothing
  end
  cll_removeindex!(cll, idx, oldPos)

  newmax = max.(newpt, cll.maxs)
  newmin = min.(newpt, cll.mins)
  if newmax != cll.maxs || newmin != cll.mins
    newgrid::Array{LinkedList{Int}, size(newpt, 1)} =
      cll_initgrid(cll_cubecount(newmax, newmin, cll.EPS))
    cll_reusegrid!(cll, newgrid)
    cll.grid = newgrid
    cll.maxs = newmax
    cll.mins = newmin
    newPos = cll_index(cll, newpt)
  end

  cll.grid[newPos...] = cons(idx, cll.grid[newPos...])
  return nothing
end

cllModify!(a...) = warn("`cllModify!` is deprecated, use `cll_modify!` instead.")

"""
Precalculates the cells, which need to be checked in cll_inrange.

### Example
`cll_neighbor_cells(e::Real, d::Real)`
"""
function cll_neighbor_cells(edge::Real, dist::Real, dim::Integer)
  dim < 1 && error("dimension has to be at least 1, but it was $dim")
  (dist < 0 || dist < 1e-9) && error("distance has to be positive, but it was $dist")
  (edge < 0 || edge < 1e-9) && error("edge length has to be positive, but it was $edge")
  r = Int(ceil(dist / edge)) + 1
  to = ntuple(i -> r, dim)
  from = .- to
  res::Vector{Vector{Int}} = []
  for c in CartesianRange(CartesianIndex(from), CartesianIndex(to))
    n = sqrt(sum(map(x -> x^2, c.I)))
    if r < n - n / 1e8
      continue
    end
    push!(res, collect(c.I))
  end
  return res
end

cllNeighborCells(a...) = warn("`cllNeighborCells` is deprecated, use `cll_neighbor_cells` instead.")

"""
Obtains the indices of `cll.data` of points, that are within `d` from `pt`.
"""
function cll_inrange(cll::CellLinkedList, pt::Point, d::Real = cll.EPS)
  if Int(ceil(d / cll.EPS)) == Int(ceil(cll.prevquery) / cll.EPS)
    return cll_inrange(cll, pt, cll.dirs, d)
  else
    neighbors = cll_neighbor_cells(cll.EPS, d, size(cll.mins, 1))
    #FIXME: cache prev query
    return cll_inrange(cll, pt, neighbors, d)
  end
end

function cll_inrange(cll::CellLinkedList, pt::Point, neighbor_dirs, d::Real = cll.EPS)
  ptCell = cll_index(cll, pt)
  res::Vector{Int} = []
  for c in neighbor_dirs
    cur = ptCell + c
    if any(cur .< 1)
      continue
    end
    if any(cur .> collect(size(cll.grid)))
      continue
    end
    for p in cll.grid[cur...]
      if norm(cll.data[:, p] - pt) < d + d / 1e8
        push!(res, p)
      end
    end
  end
  return res
end

cllInrange(a...) = warn("`cllInrange` is deprecated, use `cll_inrange` instead.")


text:
module FanBeam

export para2fan, fan2para

using ..Monads
using ..CTImages: CTSinogram
using ..Geometry
using ..Interpolation: interpolate, AbstractInterp2DOrNone
using ..Utils: linspace, ORI, _wrap_angle
using IntervalSets

"""
    para2fan(
        sinog_para::AbstractMatrix{T},
        fbg::FanBeamGeometry;
        <keyword arguments>
    ) where {T<:Real} -> fbg, sinog_fan

Convert given sinogram `sinog_para` from parallel beam geometry
to fan beam projections.

The function returns a tuple where `fbg` is the new geometry
with fan beam projections parameters and `sinog_fan` is the
converted sinogram.

# Arguments
- `sinog_para`: the input sinogram.
- `pbg`: the geometry parameters.
- `D`: focal spot to ISO distance.
- `D′=nothing`: focal spot to detectors distance; if not specified
    and also `γ` is not specified, is defined so that the total
    fan angle is 1 radian.
- `γ=nothing`: total fan angle in degrees.
- `δ=1`: detectors spacing (cell size).
- `background=nothing`: background value to be used, defaults to 0.
- `interpolation`: interpolation strategy; it should be a function
    taking a matrix as input and returning a function of the indices to
    get the interpolated value. By default it is a bilinear interpolation.

See Also: [`fan2para`](@ref)
"""
function para2fan(
    sinog_para::AbstractMatrix{T},
    fbg::FanBeamGeometry;
    background::Optional{Real} = nothing,
    interpolation::Optional{Interp} = nothing,
) where {T <: Real, Interp <: AbstractInterp2DOrNone}
    nd, nϕ = num_det(fbg), num_proj(fbg)

    @assert (nd, nϕ) == size(sinog_para) "Sinogram size $(size(sinog_para)) should match geometry ($nd,$nϕ)"

    D::T = f2iso(fbg)
    Δβ::T = deg2rad(scan_angle(fbg)) # This should be 2π?
    βs = linspace(T, ORI(0..Δβ), nϕ)
    β₀::T = deg2rad(start_angle(fbg)) + 1
    δϕi::T = nϕ / 2π # This should be α / nϕ?
    γ::T = fan_angle(fbg)
    δγ::T = γ / nd
    x′max::T = D * sin(γ/2)
    δx′i::T = nd / 2x′max
    center = center_channel(fbg)
    γ₀::T = center * δγ
    γs = linspace(T, ORI(-γ₀..γ₀), nd)
    x′₀::T = center + 1
    nthalf = nϕ + 1//2

    interpolation = maybe(interpolate, interpolation)
    interp = interpolation(sinog_para)
    z::T = maybe(zero(T), background)
    sinog_fan = similar(sinog_para, nd, nϕ)
    fill!(sinog_fan, z)

    @inline function compute_value(x′::T, ϕ::T)
        if ϕ < 0
            ϕ += π
            x′ = -x′
        end
        if ϕ > 2π
            ϕ -= π
            x′ = -x′
        end
        x′ = x′ * δx′i + x′₀
        ϕ = mod2pi(ϕ) * δϕi + 1 # need +1 in order to be in 1:nϕ
        return x′ ∈ 1..nd && ϕ ∈ 1..nϕ+1 ?
            interp(x′, _wrap_angle(ϕ, nϕ, nthalf)) : z
    end

    Threads.@threads for iγ ∈ 1:nd
        @inbounds γ′ = γs[iγ]
        x′ = D * sin(γ′)
        @inbounds @simd for iβ ∈ 1:nϕ
            β = βs[iβ]
            ϕ = β - γ′
            sinog_fan[iγ, iβ] = compute_value(x′, ϕ)
        end
    end

    fbg, CTSinogram(sinog_fan)
end

function para2fan(
    sinog_para::AbstractMatrix{T},
    pbg::ParallelBeamGeometry;
    D::Real = 500,
    D′::Optional{Real} = nothing,
    D1::Optional{Real} = nothing,
    γ::Optional{Real} = nothing,
    gamma::Optional{Real} = nothing,
    δ::Optional{Real} = one(T),
    dx::Optional{Real} = one(T),
    kwargs...,
) where {T<:Real}
    D′ = maybe(D′, D1)
    γ = maybe(γ, gamma)
    δ = maybe(δ, dx)
    fbg = FanBeamGeometry(
        T,
        pbg.ct;
        pbg.nϕ,
        pbg.nd,
        pbg.rows,
        pbg.cols,
        D,
        D′,
        γ,
        δ,
        pbg.α,
        pbg.α₀,
        pbg.center,
    )
    para2fan(sinog_para, fbg; kwargs...)
end

para2fan(g::AbstractGeometry; kwargs...) = x -> para2fan(x, g; kwargs...)

"""
    fan2para(
        sinog_fan::AbstractMatrix{T},
        fbg::FanBeamGeometry;
        <keyword arguments>
    ) where {T<:Real} -> pbg, sinog_para

Convert given sinogram `sinog_fan` from fan beam geometry
to parallel beam projections.

It is assumed that the detectors are arranged on a circle arc,
not flat.

The function returns a tuple where `fbg` is the new geometry
with fan beam projections parameters and `sinog_fan` is the
converted sinogram.

# Arguments
- `sinog_para`: the input sinogram.
- `pbg`: the parallel geometry parameters.
- `background=nothing`: background value to be used, defaults to 0.
- `interpolation`: interpolation strategy; it should be a function
    taking a matrix as input and returning a function of the indices to
    get the interpolated value. By default it is a bilinear interpolation.

See Also: [`para2fan`](@ref)
"""
function fan2para(
    sinog_fan::AbstractMatrix{T},
    fbg::FanBeamGeometry{T,DefaultTomograph};
    background::Optional{Real} = nothing,
    interpolation::Optional{Interp} = nothing,
) where {T <: Real, Interp <: AbstractInterp2DOrNone}
    nd, nϕ = num_det(fbg), num_proj(fbg)

    @assert (nd, nϕ) == size(sinog_fan) "Sinogram size $(size(sinog_fan)) should match geometry ($nd,$nϕ)"

    pbg = ParallelBeamGeometry(
        T,
        fbg.ct;
        nϕ,
        nd,
        fbg.rows,
        fbg.cols,
        fbg.α,
        fbg.center,
    )

    D::T = f2iso(fbg)
    γ::T = fan_angle(fbg)
    Δϕ::T = deg2rad(scan_angle(fbg))
    ϕs = linspace(T, ORI(0..Δϕ), nϕ)
    δβi::T = nϕ / 2π
    β₀::T = deg2rad(start_angle(fbg))
    δγi::T = nd / γ
    x′max::T = D * sin(γ/2)
    δx′::T = 2x′max / nd
    center::T = center_channel(pbg)
    x′₀::T = center * δx′
    x′₁::T = (nd - center) * δx′
    xs = linspace(T, -x′₀..x′₁, nd) # support for custom center channel
    @assert (x′₀ / D) <= 1 "Fan beam parameters incompatible |x′₀/D| = $(x′₀ / D) which should be less than 1"
    γ₀::T = center + 1 # support for custom center channel.

    interpolation = maybe(interpolate, interpolation)
    interp = interpolation(sinog_fan)
    z::T = maybe(zero(T), background)
    sinog_para = similar(sinog_fan, nd, nϕ)
    fill!(sinog_para, z)
    nthalf = nϕ + 1//2

    @inline function compute_value(γ′::T, β::T)
        if β < 0
            # β += 2γ′ + π
            # γ′ = -γ′
            β += π - 2γ′ # We use by default the symmetry at π instead of 2π
            γ′ = -γ′
        end
        if β > 2π
            # β += 2γ′ - π
            # γ′ = -γ′
            β -= π + 2γ′
            γ′ = -γ′
        end
        γ′′ = γ′ * δγi + γ₀
        β = mod2pi(β) * δβi + 1 # need +1 in order to be in 1:nβ
        return γ′′ ∈ 1..nd && β ∈ 1..nϕ+1 ?
            interp(γ′′, _wrap_angle(β, nϕ, nthalf)) : z
    end

    Threads.@threads for ix ∈ 1:nd
        @inbounds x′ = xs[ix]
        γ′ = asin(x′ / D)
        @inbounds @simd for iϕ ∈ 1:nϕ
            ϕ = ϕs[iϕ]
            #β::T = ϕ - γ′ + β₀ # Default tomograph geometry requires '-'
            β = ϕ + γ′ + β₀ # It does not!
            sinog_para[ix, iϕ] = compute_value(γ′, β)
        end
    end

    pbg, CTSinogram(sinog_para)
end

fan2para(g::AbstractGeometry; kwargs...) = x -> fan2para(x, g; kwargs...)

end # module


text:
"""
Internal data for the worker process of the parallel evaluator.
"""
struct ParallelEvaluatorWorker{P<:OptimizationProblem}
    problem::P

    ParallelEvaluatorWorker(problem::P) where {P<:OptimizationProblem} =
        new{P}(problem)
end

fitness(params::Individual, worker::ParallelEvaluatorWorker) =
    fitness(params, worker.problem)

const ChannelRef{T} = RemoteChannel{Channel{T}}
const ParallelEvaluatorWorkerRef{P} = ChannelRef{ParallelEvaluatorWorker{P}}

fitness(params::Individual, worker_ref::ParallelEvaluatorWorkerRef{P}) where P =
    fitness(params, fetch(fetch(worker_ref))::ParallelEvaluatorWorker{P})

"""
Current state of fitness function evaluation for the vector of candidates.
"""
mutable struct ParallelEvaluationState{FA, A}
    archive::A
    candidates::Vector{Candidate{FA}}  # candidates to calculate fitness for

    worker_busy::Vector{Bool}         # if the worker is busy calculating
    retry_queue::Vector{Int}          # queue to candidate indices to retry

    worker_finished::Condition        # gets notified each time worker is done
                                      # fitness calculation
    next_index::Int                   # index of the next candidate to calculate fitness

    function ParallelEvaluationState(archive::A, nworkers::Int) where {A<:Archive}
        FA = fitness_type(archive)
        new{FA,A}(archive, Vector{Candidate{FA}}(),
                  fill(false, nworkers), Vector{Int}(), Condition(), 0)
    end
end

is_stopped(estate::ParallelEvaluationState) = estate.next_index == 0

abort!(estate::ParallelEvaluationState) = (estate.next_index = 0)

"""
Reset the current `ParallelEvaluationState` and the vector
of candidates that need fitness evaluation.
"""
function reset!(estate::ParallelEvaluationState{FA}, candidates::Vector{Candidate{FA}}) where FA
    estate.candidates = candidates
    fill!(estate.worker_busy, false)
    empty!(estate.retry_queue)
    estate.next_index = 1
    return estate
end

"""
Get the index of the next candidate for evaluation
based on the Base.pmap() code.
"""
function next_candidate!(estate::ParallelEvaluationState, worker_ix::Int)
    @assert !estate.worker_busy[worker_ix]

    task_ix = 0
    if !is_stopped(estate) && estate.next_index <= length(estate.candidates)
        # find the next candidate with unevaluated fitness
        while !is_stopped(estate) && estate.next_index <= length(estate.candidates)
            if isnafitness(estate.candidates[estate.next_index].fitness, fitness_scheme(estate.archive))
                task_ix = estate.next_index
                estate.next_index += 1
                break;
            else
                estate.next_index += 1
            end
        end
    end
    if task_ix == 0
        if !isempty(estate.retry_queue)
            task_ix = popfirst!(estate.retry_queue)
        else
            # Handles the condition where we have finished processing the requested
            # lsts as well as any retryqueue entries, but there are still some jobs
            # active that may result in an error and have to be retried.
            while any(estate.worker_busy)
                wait(estate.worker_finished)
                if !isempty(estate.retry_queue)
                    task_ix = popfirst!(estate.retry_queue)
                    break
                end
            end
        end
    end
    estate.worker_busy[worker_ix] = task_ix != 0
    return task_ix
end

"""
Notify that the worker process is finished and reset its busy flag.
"""
function worker_finished!(estate::ParallelEvaluationState, worker_ix::Int)
    @assert estate.worker_busy[worker_ix]
    estate.worker_busy[worker_ix] = false
    notify(estate.worker_finished; all=true)
    return estate
end

"""
Fitness evaluator that distributes candidates fitness calculation
among several worker processes.
"""
mutable struct ParallelEvaluator{F, FA, FS, P<:OptimizationProblem, A<:Archive} <: Evaluator{P}
    problem::P
    archive::A
    num_evals::Int
    last_fitness::F
    arch_nafitness::FA  # NA fitness

    worker_refs::Vector{ParallelEvaluatorWorkerRef{P}}
    eval_state::ParallelEvaluationState{FA, A}
end

function ParallelEvaluator(
        problem::P, archive::A;
        pids::Vector{Int} = workers()) where {P<:OptimizationProblem, A<:Archive}
    fs = fitness_scheme(problem)
    F = fitness_type(fs)
    FA = fitness_type(archive)
    ParallelEvaluator{F, FA, typeof(fs), P, A}(problem,
        archive,
        0, nafitness(fs), nafitness(FA),
        [RemoteChannel(function ()
                     # create fake channel and put problem there
                     ch = Channel{ParallelEvaluatorWorker{P}}(1)
                     put!(ch, ParallelEvaluatorWorker(copy(problem)))
                     ch
                   end, pid) for pid in pids],
        ParallelEvaluationState(archive, length(pids))
    )
end

ParallelEvaluator(problem::OptimizationProblem;
        pids::Vector{Int} = workers(),
        archiveCapacity::Integer = 10) =
    ParallelEvaluator(problem, TopListArchive(fitness_scheme(problem), numdims(problem), archiveCapacity),
                      pids=pids)

num_evals(e::ParallelEvaluator) = e.num_evals

function update_fitness!(f::Any, e::ParallelEvaluator, candidates::Any;
                         force::Bool=false)
    # FIXME use force
    reset!(e.eval_state, candidates)

    # based on pmap() code
    @sync begin
        for (widx, wref) in enumerate(e.worker_refs)
            @async begin
                while (candi_ix = next_candidate!(e.eval_state, widx)) != 0
                    try
                        candi = candidates[candi_ix]
                        candi_fitness = remotecall_fetch(fitness, wref.where, candi.params, wref)
                        candi.fitness = archived_fitness(candi_fitness, e.archive)
                        e.last_fitness = candi_fitness
                        e.num_evals += 1
                        add_candidate!(e.archive, candi.fitness, candi.params, e.num_evals)
                        (f !== nothing) && f(candi)
                    catch ex
                        abort!(e.eval_state) # when one worker fails, the whole process is aborted
                        rethrow(ex)
                    finally
                        worker_finished!(e.eval_state, widx)
                    end
                end
            end
        end
    end

    return candidates
end

# FIXME it's not efficient to calculate fitness like that with `ParallelEvaluator`
update_fitness!(f::Any, e::ParallelEvaluator, candidate::Candidate;
                force::Bool=false) = update_fitness!(f, e, [candidate])

# FIXME it's not efficient to calculate fitness like that with `ParallelEvaluator`
function fitness(params::Individual, e::ParallelEvaluator{F, FA}, tag::Int=0) where {F, FA}
    candi = Candidate{FA}(params, -1, e.arch_nafitness, nothing, tag)
    update_fitness!(e, [candi])
    return candi.fitness
end


text:
mutable struct _ScatterProblem{O, S <: WavePacketState, WP <: _WavePacket{S}} <: WaveguideProblem
	H::O
	Ls::Array{O}
	σ::O
	system_state

	ψᵢ::WP
	ψₒ::Array{Mode}

	ts::StepRangeLen

	displace_output::Bool
end

function WaveguideProblem(sys, ψₒ, ts, disp_out::Bool = false)
	WaveguideProblem(sys, NoPacket(), ψₒ, ts, disp_out)
end

function WaveguideProblem(sys, ψᵢ::_WavePacket, ψₒ, ts, disp_out::Bool = false)
	ts isa StepRangeLen || (ts = _toRange(ts))

	ψₒ = ψₒ isa Mode ? [ψₒ] : ψₒ

	if length(sys) == 3
		return _ScatterProblem(sys[1], typeof(sys[1])[], sys[2], sys[3], ψᵢ, ψₒ, ts, disp_out)
	end
	return _ScatterProblem(sys[1], sys[2], sys[3], sys[4], ψᵢ, ψₒ, ts, disp_out)
end


text:
struct OLS
    X::Array{Float64,2}
    y::Array{Float64,1}
    betas::Array{Float64,1}
end

ols(X::Array{Float64,2}, y::Array{Float64,1})::OLS = OLS(X, y, qr(X, Val(true)) \ y)

function wls(X::Array{Float64,2}, y::Array{Float64,1}, wts::Array{Float64,1})
    W = Diagonal(sqrt.(wts))
    return ols(W*X, W*y)
end

residuals(ols::OLS) = ols.y .- ols.X * ols.betas

coef(ols::OLS) = ols.betas

predict(ols::OLS) = ols.X * ols.betas

predict(ols::OLS, X::Array{Float64,2}) = X * ols.betas

max_stars_repo_name: akadal/LinRegOutliers

max_issues_repo_name: akadal/LinRegOutliers

max_forks_repo_name: akadal/LinRegOutliers


text:
using HPAT

#HPAT.CaptureAPI.set_debug_level(3)
ParallelAccelerator.set_debug_level(3)
using CompilerTools
CompilerTools.OptFramework.set_debug_level(3)

@acc hpat function test1(file_name, file_name2)
    t1 = DataSource(DataTable{:userid=Int64, :val2=Float64}, HDF5, file_name)
    t2 = DataSource(DataTable{:userid=Int64, :val3=Float64}, HDF5, file_name2)
    t3 = join(t1,t2, :userid==:userid, :userid)
    t4 = aggregate(t3, :userid, :sumo2 = sum(:val2==1.1), :size_val3 = length(:val3))
    return t4[:sumo2], t4[:size_val3]
end

a,b = test1("test1_1.hdf5","test1_2.hdf5")

println(a)
println(b)


text:
function calculate_tgrad(sys::AbstractODESystem;
                         simplify=false)
  isempty(get_tgrad(sys)[]) || return get_tgrad(sys)[]  # use cached tgrad, if possible

  # We need to remove explicit time dependence on the state because when we
  # have `u(t) * t` we want to have the tgrad to be `u(t)` instead of `u'(t) *
  # t + u(t)`.
  rhs = [detime_dvs(eq.rhs) for eq ∈ equations(sys)]
  iv = get_iv(sys)
  xs = states(sys)
  rule = Dict(map((x, xt) -> xt=>x, detime_dvs.(xs), xs))
  rhs = substitute.(rhs, Ref(rule))
  tgrad = [expand_derivatives(ModelingToolkit.Differential(iv)(r), simplify) for r in rhs]
  reverse_rule = Dict(map((x, xt) -> x=>xt, detime_dvs.(xs), xs))
  tgrad = Num.(substitute.(tgrad, Ref(reverse_rule)))
  get_tgrad(sys)[] = tgrad
  return tgrad
end

function calculate_jacobian(sys::AbstractODESystem;
                            sparse=false, simplify=false)
    isempty(get_jac(sys)[]) || return get_jac(sys)[]  # use cached Jacobian, if possible
    rhs = [eq.rhs for eq ∈ equations(sys)]

    iv = get_iv(sys)
    dvs = states(sys)

    if sparse
        jac = sparsejacobian(rhs, dvs, simplify=simplify)
    else
        jac = jacobian(rhs, dvs, simplify=simplify)
    end

    get_jac(sys)[] = jac  # cache Jacobian
    return jac
end

function generate_tgrad(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys);
                        simplify=false, kwargs...)
    tgrad = calculate_tgrad(sys,simplify=simplify)
    return build_function(tgrad, dvs, ps, get_iv(sys); kwargs...)
end

function generate_jacobian(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys);
                           simplify=false, sparse = false, kwargs...)
    jac = calculate_jacobian(sys;simplify=simplify,sparse=sparse)
    return build_function(jac, dvs, ps, get_iv(sys); kwargs...)
end

@noinline function throw_invalid_derivative(dervar, eq)
    msg = "The derivative variable must be isolated to the left-hand " *
    "side of the equation like `$dervar ~ ...`.\n Got $eq."
    throw(InvalidSystemException(msg))
end

function check_derivative_variables(eq, expr=eq.rhs)
    istree(expr) || return nothing
    if operation(expr) isa Differential
        throw_invalid_derivative(expr, eq)
    end
    foreach(Base.Fix1(check_derivative_variables, eq), arguments(expr))
end

function generate_function(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys); kwargs...)
    # optimization
    #obsvars = map(eq->eq.lhs, observed(sys))
    #fulldvs = [dvs; obsvars]

    eqs = equations(sys)
    foreach(check_derivative_variables, eqs)
    # substitute x(t) by just x
    rhss = [deq.rhs for deq in eqs]
    #obss = [makesym(value(eq.lhs)) ~ substitute(eq.rhs, sub) for eq ∈ observed(sys)]
    #rhss = Let(obss, rhss)

    # TODO: add an optional check on the ordering of observed equations
    return build_function(rhss,
                          map(x->time_varying_as_func(value(x), sys), dvs),
                          map(x->time_varying_as_func(value(x), sys), ps),
                          get_iv(sys); kwargs...)
end

function time_varying_as_func(x, sys)
    # if something is not x(t) (the current state)
    # but is `x(t-1)` or something like that, pass in `x` as a callable function rather
    # than pass in a value in place of x(t).
    #
    # This is done by just making `x` the argument of the function.
    if istree(x) &&
        operation(x) isa Sym &&
        !(length(arguments(x)) == 1 && isequal(arguments(x)[1], independent_variable(sys)))
        return operation(x)
    end
    return x
end

function calculate_massmatrix(sys::AbstractODESystem; simplify=false)
    eqs = equations(sys)
    dvs = states(sys)
    M = zeros(length(eqs),length(eqs))
    for (i,eq) in enumerate(eqs)
        if eq.lhs isa Term && operation(eq.lhs) isa Differential
            j = findfirst(x->isequal(tosymbol(x),tosymbol(var_from_nested_derivative(eq.lhs)[1])),dvs)
            M[i,j] = 1
        else
            _iszero(eq.lhs) || error("Only semi-explicit constant mass matrices are currently supported. Faulty equation: $eq.")
        end
    end
    M = simplify ? ModelingToolkit.simplify.(M) : M
    # M should only contain concrete numbers
    M == I ? I : M
end

jacobian_sparsity(sys::AbstractODESystem) =
    jacobian_sparsity([eq.rhs for eq ∈ equations(sys)],
                      [dv for dv in states(sys)])

function isautonomous(sys::AbstractODESystem)
    tgrad = calculate_tgrad(sys;simplify=true)
    all(iszero,tgrad)
end

function DiffEqBase.ODEFunction(sys::AbstractODESystem, args...; kwargs...)
    ODEFunction{true}(sys, args...; kwargs...)
end

"""
```julia
function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys);
                                     version = nothing, tgrad=false,
                                     jac = false,
                                     sparse = false,
                                     kwargs...) where {iip}
```

Create an `ODEFunction` from the [`ODESystem`](@ref). The arguments `dvs` and `ps`
are used to set the order of the dependent variable and parameter vectors,
respectively.
"""
function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys), u0 = nothing;
                                     version = nothing, tgrad=false,
                                     jac = false,
                                     eval_expression = true,
                                     sparse = false, simplify=false,
                                     eval_module = @__MODULE__,
                                     kwargs...) where {iip}

    f_gen = generate_function(sys, dvs, ps; expression=Val{eval_expression}, expression_module=eval_module, kwargs...)
    f_oop,f_iip = eval_expression ? (@RuntimeGeneratedFunction(eval_module, ex) for ex in f_gen) : f_gen
    f(u,p,t) = f_oop(u,p,t)
    f(du,u,p,t) = f_iip(du,u,p,t)

    if tgrad
        tgrad_gen = generate_tgrad(sys, dvs, ps;
                                   simplify=simplify,
                                   expression=Val{eval_expression}, expression_module=eval_module, kwargs...)
        tgrad_oop,tgrad_iip = eval_expression ? (@RuntimeGeneratedFunction(eval_module, ex) for ex in tgrad_gen) : tgrad_gen
        _tgrad(u,p,t) = tgrad_oop(u,p,t)
        _tgrad(J,u,p,t) = tgrad_iip(J,u,p,t)
    else
        _tgrad = nothing
    end

    if jac
        jac_gen = generate_jacobian(sys, dvs, ps;
                                    simplify=simplify, sparse = sparse,
                                    expression=Val{eval_expression}, expression_module=eval_module, kwargs...)
        jac_oop,jac_iip = eval_expression ? (@RuntimeGeneratedFunction(eval_module, ex) for ex in jac_gen) : jac_gen
        _jac(u,p,t) = jac_oop(u,p,t)
        _jac(J,u,p,t) = jac_iip(J,u,p,t)
    else
        _jac = nothing
    end

    M = calculate_massmatrix(sys)

    _M = (u0 === nothing || M == I) ? M : ArrayInterface.restructure(u0 .* u0',M)

    observedfun = let sys = sys, dict = Dict()
        function generated_observed(obsvar, u, p, t)
            obs = get!(dict, value(obsvar)) do
                build_explicit_observed_function(sys, obsvar)
            end
            obs(u, p, t)
        end
    end

    ODEFunction{iip}(
                     f,
                     jac = _jac === nothing ? nothing : _jac,
                     tgrad = _tgrad === nothing ? nothing : _tgrad,
                     mass_matrix = _M,
                     jac_prototype = sparse ? similar(get_jac(sys)[],Float64) : nothing,
                     syms = Symbol.(states(sys)),
                     indepsym = Symbol(independent_variable(sys)),
                     observed = observedfun,
                    )
end

"""
```julia
function DiffEqBase.ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys);
                                     version = nothing, tgrad=false,
                                     jac = false,
                                     sparse = false,
                                     kwargs...) where {iip}
```

Create a Julia expression for an `ODEFunction` from the [`ODESystem`](@ref).
The arguments `dvs` and `ps` are used to set the order of the dependent
variable and parameter vectors, respectively.
"""
struct ODEFunctionExpr{iip} end

struct ODEFunctionClosure{O, I} <: Function
    f_oop::O
    f_iip::I
end
(f::ODEFunctionClosure)(u, p, t) = f.f_oop(u, p, t)
(f::ODEFunctionClosure)(du, u, p, t) = f.f_iip(du, u, p, t)

function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys), u0 = nothing;
                                     version = nothing, tgrad=false,
                                     jac = false,
                                     linenumbers = false,
                                     sparse = false, simplify=false,
                                     kwargs...) where {iip}

    f_oop, f_iip = generate_function(sys, dvs, ps; expression=Val{true}, kwargs...)
    fsym = gensym(:f)
    _f = :($fsym = ModelingToolkit.ODEFunctionClosure($f_oop, $f_iip))
    tgradsym = gensym(:tgrad)
    if tgrad
        tgrad_oop, tgrad_iip = generate_tgrad(sys, dvs, ps;
                                simplify=simplify,
                                expression=Val{true}, kwargs...)
        _tgrad = :($tgradsym = ModelingToolkit.ODEFunctionClosure($tgrad_oop, $tgrad_iip))
    else
        _tgrad = :($tgradsym = nothing)
    end

    jacsym = gensym(:jac)
    if jac
        jac_oop,jac_iip = generate_jacobian(sys, dvs, ps;
                                 sparse=sparse, simplify=simplify,
                                 expression=Val{true}, kwargs...)
        _jac = :($jacsym = ModelingToolkit.ODEFunctionClosure($jac_oop, $jac_iip))
    else
        _jac = :($jacsym = nothing)
    end

    M = calculate_massmatrix(sys)

    _M = (u0 === nothing || M == I) ? M : ArrayInterface.restructure(u0 .* u0',M)

    jp_expr = sparse ? :(similar($(get_jac(sys)[]),Float64)) : :nothing
    ex = quote
        $_f
        $_tgrad
        $_jac
        M = $_M
        ODEFunction{$iip}(
                          $fsym,
                          jac = $jacsym,
                          tgrad = $tgradsym,
                          mass_matrix = M,
                          jac_prototype = $jp_expr,
                          syms = $(Symbol.(states(sys))),
                          indepsym = $(QuoteNode(Symbol(independent_variable(sys)))),
                         )
    end
    !linenumbers ? striplines(ex) : ex
end

function process_DEProblem(constructor, sys::AbstractODESystem,u0map,parammap;
                           version = nothing, tgrad=false,
                           jac = false,
                           checkbounds = false, sparse = false,
                           simplify=false,
                           linenumbers = true, parallel=SerialForm(),
                           eval_expression = true,
                           kwargs...)
    dvs = states(sys)
    ps = parameters(sys)
    defs = defaults(sys)

    u0 = varmap_to_vars(u0map,dvs; defaults=defs)
    p = varmap_to_vars(parammap,ps; defaults=defs)

    if u0 !== nothing
        length(dvs) == length(u0) || throw(ArgumentError("States ($(length(dvs))) and initial conditions ($(length(u0))) are of different lengths."))
    end

    f = constructor(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,checkbounds=checkbounds,
                    linenumbers=linenumbers,parallel=parallel,simplify=simplify,
                    sparse=sparse,eval_expression=eval_expression,kwargs...)
    return f, u0, p
end

function ODEFunctionExpr(sys::AbstractODESystem, args...; kwargs...)
    ODEFunctionExpr{true}(sys, args...; kwargs...)
end

function DiffEqBase.ODEProblem(sys::AbstractODESystem, args...; kwargs...)
    ODEProblem{true}(sys, args...; kwargs...)
end

"""
```julia
function DiffEqBase.ODEProblem{iip}(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    simplify=false,
                                    linenumbers = true, parallel=SerialForm(),
                                    kwargs...) where iip
```

Generates an ODEProblem from an ODESystem and allows for automatically
symbolically calculating numerical enhancements.
"""
function DiffEqBase.ODEProblem{iip}(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();kwargs...) where iip
    f, u0, p = process_DEProblem(ODEFunction{iip}, sys, u0map, parammap; kwargs...)
    ODEProblem{iip}(f,u0,tspan,p;kwargs...)
end

"""
```julia
function DiffEqBase.ODEProblemExpr{iip}(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    linenumbers = true, parallel=SerialForm(),
                                    skipzeros=true, fillzeros=true,
                                    simplify=false,
                                    kwargs...) where iip
```

Generates a Julia expression for constructing an ODEProblem from an
ODESystem and allows for automatically symbolically calculating
numerical enhancements.
"""
struct ODEProblemExpr{iip} end

function ODEProblemExpr{iip}(sys::AbstractODESystem,u0map,tspan,
                             parammap=DiffEqBase.NullParameters();
                             kwargs...) where iip

    f, u0, p = process_DEProblem(ODEFunctionExpr{iip}, sys, u0map, parammap; kwargs...)
    linenumbers = get(kwargs, :linenumbers, true)

    ex = quote
        f = $f
        u0 = $u0
        tspan = $tspan
        p = $p
        ODEProblem(f,u0,tspan,p;$(kwargs...))
    end
    !linenumbers ? striplines(ex) : ex
end

function ODEProblemExpr(sys::AbstractODESystem, args...; kwargs...)
    ODEProblemExpr{true}(sys, args...; kwargs...)
end

### Enables Steady State Problems ###
function DiffEqBase.SteadyStateProblem(sys::AbstractODESystem, args...; kwargs...)
    SteadyStateProblem{true}(sys, args...; kwargs...)
end

"""
```julia
function DiffEqBase.SteadyStateProblem(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    linenumbers = true, parallel=SerialForm(),
                                    kwargs...) where iip
```
Generates an SteadyStateProblem from an ODESystem and allows for automatically
symbolically calculating numerical enhancements.
"""
function DiffEqBase.SteadyStateProblem{iip}(sys::AbstractODESystem,u0map,
                                            parammap=DiffEqBase.NullParameters();
                                            kwargs...) where iip
    f, u0, p = process_DEProblem(ODEFunction{iip}, sys, u0map, parammap; kwargs...)
    SteadyStateProblem(f,u0,p;kwargs...)
end

"""
```julia
function DiffEqBase.SteadyStateProblemExpr(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    skipzeros=true, fillzeros=true,
                                    linenumbers = true, parallel=SerialForm(),
                                    kwargs...) where iip
```
Generates a Julia expression for building a SteadyStateProblem from
an ODESystem and allows for automatically symbolically calculating
numerical enhancements.
"""
struct SteadyStateProblemExpr{iip} end

function SteadyStateProblemExpr{iip}(sys::AbstractODESystem,u0map,
                                    parammap=DiffEqBase.NullParameters();
                                    kwargs...) where iip
    f, u0, p = process_DEProblem(ODEFunctionExpr{iip}, sys, u0map, parammap; kwargs...)
    linenumbers = get(kwargs, :linenumbers, true)
    ex = quote
        f = $f
        u0 = $u0
        p = $p
        SteadyStateProblem(f,u0,p;$(kwargs...))
    end
    !linenumbers ? striplines(ex) : ex
end

function SteadyStateProblemExpr(sys::AbstractODESystem, args...; kwargs...)
    SteadyStateProblemExpr{true}(sys, args...; kwargs...)
end

isdifferential(expr) = istree(expr) && operation(expr) isa Differential
isdiffeq(eq) = isdifferential(eq.lhs)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# cheap LPS(Logical Program Synthesis)

function lps_prover(spec, steplimit, contralimit)
 cdx=readcore(spec)
 tdx=alltemplateof(cdx)
 gb=[lidsof(:C1, cdx)]
 conds = []
 nstep = 0;

 evalproc(cdx.proc)
 execs = execof(cdx)
 cdx, branch = srule(cdx, execs)

 while true
  ga=dostepgoals1(gb, cdx)
  nstep += 1
  conds = contradictionsof(cdx)
  if length(conds) >= contralimit;break end
  if nstep >= steplimit;break end
  gb = ga
 end
 return conds,cdx,branch,INVARS,OUTVARS
end

## S-rule
function fname(str)
  ix=findfirst("(", str)
  if isa(ix, Nothing)
    return Nothing
  else
    return Symbol(str[1:(ix[1]-1)])
  end
end

function execof(core)
  exdb = []
  for pd in core.proc
    fn = fname(pd)
    if !isa(fn, Nothing)
      push!(exdb, fn)
    end
  end
  return exdb
end

function isexec(sym, edb)
  return in(sym,edb)
end

function srule(core, execs)
 branch = Dict{Symbol,Array}()
 for cid in keys(core.cdb)
   lits = []
   exes = []
   for lid in lidsof(cid, core)
     psym = literalof(lid, core).body.args[2].args[1]
     if isexec(psym, execs)
       push!(exes, lid)
     else
       push!(lits, lid)
     end
   end
   core.clmap[cid]  = lits
   branch[cid] = exes
 end
 core, branch
end

### extraction
function traceof(rid, vars, core)
  if rid in keys(core.proof)
    step = core.proof[rid]
    vars = traceof(cidof(step.leftp,core), vars, core)
    vars = traceof(cidof(step.rightp,core), vars, core)
    vars = apply(ovarsof(step.leftp, step.rightp, core), vars, step.sigma)
    vars = apply(step.rename[1], vars, step.rename[2])
# else
#  printclause(rid, core)
  end
   return vars
end

function extract(core,ins,ons)
  einf = []
  for conid in contradictionsof(core)
    rins = traceof(conid, ins, core)
    rons = traceof(conid, ons, core)
    push!(einf, [rins,rons,conid])
  end
  return [ins,ons],einf
end

### coder for julia
function printExpr(io,e::Number)
  print(io,e)
end

function printExpr(io,e::Symbol)
  print(io,e)
end

function printExpr(io,ea::Array)
  for e in ea
    print(io,e);print(io," ")
  end
end

function printExpr(io,e::Expr)
  print(io,"$(e.args[1])(")
  if 1==length(e.args)
    print(io,")")
    return
  end
  for ix in 2:length(e.args)
    ag = e.args[ix]
    printExpr(io,ag)
    if ix != length(e.args); print(io,",") end
  end
  print(io,")")
end

function printExprVec(io,ev)
  for ix in 1:length(ev)
    printExpr(io,ev[ix])
    if ix != length(ev); print(io,",") end
  end
end

function printaline(io, core, fname, iov, einf, bra)
  ea=apply(einf[1], einf[2], iov[1])
  printExprVec(io,ea)
end

function econdsof(rid, core, brdb, pabr)
  bras = []
  crid=rid
  if iscid(crid)
    bras = union(brdb[crid], pabr)
    return bras
  else # has parents
    ps = stepof(rid, core)
    lbrs=econdsof(cidof(ps.leftp,core), core, brdb, pabr)
    rbrs=econdsof(cidof(ps.rightp,core), core, brdb, pabr)
    bras = union(lbrs, rbrs)
  end
  bras = intersect(bras,bras)
  return bras
end

function terms2julia(core, fname, econds, einf)
  jterm = []
  for elid in econds
    lf2 = literalof(elid, core)
    vars = varsof(cidof(lf2.lid, core),core)
    eterm = lf2.body
    if eterm.args[1] == :+
      eterm.args[1] = :!
    else
      eterm = eterm.args[2]
    end
    fterm = apply(vars, eterm, vcat(einf[1],einf[2]))
    push!(jterm, fterm)
  end
  return jterm
end

# :L5 == econdsof(:R4, cf, bf,[])[1]
# terms2julia(cf, :fact, [literalof(:L5,cf).body])
# => "!gt(i_C3, 0)"

function inputcond(ins, einf)
  icond = []
  for ix in 1:length(ins)
    iv = ins[ix]
    ev = einf[1][ix]
    if iv != ev
      push!(icond, :(==($iv,$ev)))
    end
  end
  return icond
end

function makeands(ts)
  if length(ts) == 0; return :[] end
  ands = ts[1]
  if length(ts) == 1; return ands end
  for ix in 2:length(ts)
    ands = :($ands && $(ts[ix]))
  end
  return ands
end

function isidaxiom(einf, fname)
  oform = einf[2][1]
  return isa(oform, Expr) && oform.args[1] == fname
end

function printbody(io,core, fname, iov, einfs, sdb)
  println(io)
  body = []
  for einf in einfs
    iconds = inputcond(iov[1],einf)
    econds = econdsof(einf[3], core, sdb, [])
    aconds = vcat(iconds, terms2julia(core, fname, econds, einf))
    ands = makeands(aconds)
    if isidaxiom(einf, fname); print(io,"# idaxiom: ") end
    if ands != []; print(io,"if $ands; ") end
    printlist(io,iov[2]);print(io,"=")

    printaline(io,core, fname, iov, einf, sdb)
    if ands != []; println(io, " end") end
  end
end

function printlist(io, syms::Array)
  for ix in 1:length(syms)
    print(io,syms[ix])
    if ix != length(syms); print(io,",") end
  end
end

function printargs(io,syms::Array)
  print(io,"(")
  printlist(io, syms)
  print(io,")")
end

function coder_julia(outfname, core, fname, iov, einfs, sdb)
  open(outfname, "w") do io
    println(io, "# $fname")
    print(io,"function $fname")
    printargs(io,iov[1])
    printbody(io,core, fname, iov, einfs, sdb)
    print(io,"return "); printlist(io,iov[2]);println(io)
    println(io,"end")
  end
end

#### combined
function lps(outfname, specp, fname, steplimit, contralimit)
  df,cf,bf,ins,ons=lps_prover(specp, steplimit, contralimit)
  invs = ins.args
  onvs = ons.args

  ves, eis=extract(cf, invs, onvs)
  #coder_julia(cf, fname, [invs,onvs], eis, bf)
  coder_julia(outfname, cf, fname, ves, eis, bf)

  return df,cf,bf,eis,ins,ons
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
module BenchFilterSum

import Random
using BenchmarkTools
using Transducers

const SUITE = BenchmarkGroup()

function naive_sum(xs, acc = false)
    for x in xs
        acc += x
    end
    return acc
end

Random.seed!(12345)
for n in [1000, 10000]
    s0 = SUITE[string(n)] = BenchmarkGroup()

    for (xslabel, xs, init) in [
        ("UnitRange", (x for x in 1:n if isodd(x)), 0),
        ("RandomFloats", (x for x in randn(n) if x > 0), 0.0),
    ]
        s1 = s0[xslabel] = BenchmarkGroup()

        s2 = s1["noinit"] = BenchmarkGroup()
        s2["naive"] = @benchmarkable naive_sum($xs)
        s2["base"] = @benchmarkable sum($xs)
        s2["xf"] = @benchmarkable sum($(eduction(xs)))

        s2 = s1["withinit"] = BenchmarkGroup()
        s2["naive"] = @benchmarkable naive_sum($xs, $init)
        s2["base"] = @benchmarkable foldl(+, $xs; init = $init)
        s2["xf"] = @benchmarkable foldl(+, $(eduction(xs)); init = $init)
    end
end

end  # module
BenchFilterSum.SUITE


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(x^3*sinh(a+b*x^2), x) == :(-1//2*b^-2*sinh(a+b*x^2)+(1/2)*b^-1*x^2*cosh(a+b*x^2))
@test integrate(x^2*sinh(a+b*x^2), x) == :((1/2)*x*b^-1*cosh(a+b*x^2)+-1//8*pi^(1/2)*b^-3//2*erf(x*b^(1/2))*exp(-1a)+-1//8*pi^(1/2)*b^-3//2*erfi(x*b^(1/2))*exp(a))
@test integrate(x*sinh(a+b*x^2), x) == :((1/2)*b^-1*cosh(a+b*x^2))
@test integrate(sinh(a+b*x^2), x) == :(-1//4*pi^(1/2)*b^-1//2*erf(x*b^(1/2))*exp(-1a)+1//4*pi^(1/2)*b^-1//2*erfi(x*b^(1/2))*exp(a))
@test integrate(x^-1*sinh(a+b*x^2), x) == :((1/2)*cosh(a)*sinhintegral(b*x^2)+(1/2)*coshintegral(b*x^2)*sinh(a))
@test integrate(x^-2*sinh(a+b*x^2), x) == :(-1*x^-1*sinh(a+b*x^2)+(1/2)*pi^(1/2)*b^(1/2)*erf(x*b^(1/2))*exp(-1a)+(1/2)*pi^(1/2)*b^(1/2)*erfi(x*b^(1/2))*exp(a))
@test integrate(x^-3*sinh(a+b*x^2), x) == :(-1//2*x^-2*sinh(a+b*x^2)+(1/2)*b*cosh(a)*coshintegral(b*x^2)+(1/2)*b*sinh(a)*sinhintegral(b*x^2))
@test integrate(x^3*sinh(a+b*x^2)^2, x) == :(-1//8*x^4+-1//8*b^-2*sinh(a+b*x^2)^2+1//4*b^-1*x^2*cosh(a+b*x^2)*sinh(a+b*x^2))
@test integrate(x^2*sinh(a+b*x^2)^2, x) == :(-1//6*x^3+1//8*x*b^-1*sinh(2a+2*b*x^2)+-1//64*2^(1/2)*pi^(1/2)*b^-3//2*erfi(x*2^(1/2)*b^(1/2))*exp(2a)+1//64*2^(1/2)*pi^(1/2)*b^-3//2*erf(x*2^(1/2)*b^(1/2))*exp(-2a))
@test integrate(x*sinh(a+b*x^2)^2, x) == :(-1//4*x^2+1//4*b^-1*cosh(a+b*x^2)*sinh(a+b*x^2))
@test integrate(sinh(a+b*x^2)^2, x) == :(-1//2*x+1//16*2^(1/2)*pi^(1/2)*b^-1//2*erf(x*2^(1/2)*b^(1/2))*exp(-2a)+1//16*2^(1/2)*pi^(1/2)*b^-1//2*erfi(x*2^(1/2)*b^(1/2))*exp(2a))
@test integrate(x^-1*sinh(a+b*x^2)^2, x) == :(-1//2*log(x)+1//4*cosh(2a)*coshintegral(2*b*x^2)+1//4*sinh(2a)*sinhintegral(2*b*x^2))
@test integrate(x^-2*sinh(a+b*x^2)^2, x) == :(-1*x^-1*sinh(a+b*x^2)^2+-1//4*2^(1/2)*pi^(1/2)*b^(1/2)*erf(x*2^(1/2)*b^(1/2))*exp(-2a)+1//4*2^(1/2)*pi^(1/2)*b^(1/2)*erfi(x*2^(1/2)*b^(1/2))*exp(2a))
@test integrate(x^-3*sinh(a+b*x^2)^2, x) == :(1//4*x^-2+-1//4*x^-2*cosh(2a+2*b*x^2)+(1/2)*b*cosh(2a)*sinhintegral(2*b*x^2)+(1/2)*b*coshintegral(2*b*x^2)*sinh(2a))
@test integrate(x^3*sinh(a+b*x^2)^3, x) == :(-1//18*b^-2*sinh(a+b*x^2)^3+1//3*b^-2*sinh(a+b*x^2)+-1//3*b^-1*x^2*cosh(a+b*x^2)+1//6*b^-1*x^2*sinh(a+b*x^2)^2*cosh(a+b*x^2))
@test integrate(x^2*sinh(a+b*x^2)^3, x) == :(-3//8*x*b^-1*cosh(a+b*x^2)+1//24*x*b^-1*cosh(3a+3*b*x^2)+3//32*pi^(1/2)*b^-3//2*erf(x*b^(1/2))*exp(-1a)+3//32*pi^(1/2)*b^-3//2*erfi(x*b^(1/2))*exp(a)+-1//288*3^(1/2)*pi^(1/2)*b^-3//2*erf(x*3^(1/2)*b^(1/2))*exp(-3a)+-1//288*3^(1/2)*pi^(1/2)*b^-3//2*erfi(x*3^(1/2)*b^(1/2))*exp(3a))
@test integrate(x*sinh(a+b*x^2)^3, x) == :(-1//2*b^-1*cosh(a+b*x^2)+1//6*b^-1*cosh(a+b*x^2)^3)
@test integrate(sinh(a+b*x^2)^3, x) == :(-3//16*pi^(1/2)*b^-1//2*erfi(x*b^(1/2))*exp(a)+3//16*pi^(1/2)*b^-1//2*erf(x*b^(1/2))*exp(-1a)+-1//48*3^(1/2)*pi^(1/2)*b^-1//2*erf(x*3^(1/2)*b^(1/2))*exp(-3a)+1//48*3^(1/2)*pi^(1/2)*b^-1//2*erfi(x*3^(1/2)*b^(1/2))*exp(3a))
@test integrate(x^-1*sinh(a+b*x^2)^3, x) == :(-3//8*cosh(a)*sinhintegral(b*x^2)+-3//8*coshintegral(b*x^2)*sinh(a)+1//8*cosh(3a)*sinhintegral(3*b*x^2)+1//8*coshintegral(3*b*x^2)*sinh(3a))
@test integrate(x^-2*sinh(a+b*x^2)^3, x) == :(-1*x^-1*sinh(a+b*x^2)^3+-3//8*pi^(1/2)*b^(1/2)*erf(x*b^(1/2))*exp(-1a)+-3//8*pi^(1/2)*b^(1/2)*erfi(x*b^(1/2))*exp(a)+1//8*3^(1/2)*pi^(1/2)*b^(1/2)*erf(x*3^(1/2)*b^(1/2))*exp(-3a)+1//8*3^(1/2)*pi^(1/2)*b^(1/2)*erfi(x*3^(1/2)*b^(1/2))*exp(3a))
@test integrate(x^-3*sinh(a+b*x^2)^3, x) == :(-1//8*x^-2*sinh(3a+3*b*x^2)+3//8*x^-2*sinh(a+b*x^2)+-3//8*b*cosh(a)*coshintegral(b*x^2)+-3//8*b*sinh(a)*sinhintegral(b*x^2)+3//8*b*cosh(3a)*coshintegral(3*b*x^2)+3//8*b*sinh(3a)*sinhintegral(3*b*x^2))
@test integrate(x*sinh(a+b*x^2)^7, x) == :((1/2)*b^-1*cosh(a+b*x^2)^3+-3//10*b^-1*cosh(a+b*x^2)^5+-1//2*b^-1*cosh(a+b*x^2)+1//14*b^-1*cosh(a+b*x^2)^7)
@test integrate((e*x)^m*sinh(a+b*x^2)^3, x) == :(-3//16*e^-1*(b*x^2)^(-1//2+-1//2*m)*(e*x)^(1+m)*SpecialFunctions.gamma(1/2+(1/2)*m,b*x^2)*exp(-1a)+3//16*e^-1*(e*x)^(1+m)*(-1*b*x^2)^(-1//2+-1//2*m)*SpecialFunctions.gamma(1/2+(1/2)*m,-1*b*x^2)*exp(a)+-1//16*3^(-1//2+-1//2*m)*e^-1*(e*x)^(1+m)*(-1*b*x^2)^(-1//2+-1//2*m)*SpecialFunctions.gamma(1/2+(1/2)*m,-3*b*x^2)*exp(3a)+1//16*3^(-1//2+-1//2*m)*e^-1*(b*x^2)^(-1//2+-1//2*m)*(e*x)^(1+m)*SpecialFunctions.gamma(1/2+(1/2)*m,3*b*x^2)*exp(-3a))
@test integrate((e*x)^m*sinh(a+b*x^2)^2, x) == :(-1//2*e^-1*(e*x)^(1+m)*(1+m)^-1+-1*2^(-7//2+-1//2*m)*e^-1*(b*x^2)^(-1//2+-1//2*m)*(e*x)^(1+m)*SpecialFunctions.gamma(1/2+(1/2)*m,2*b*x^2)*exp(-2a)+-1*2^(-7//2+-1//2*m)*e^-1*(e*x)^(1+m)*(-1*b*x^2)^(-1//2+-1//2*m)*SpecialFunctions.gamma(1/2+(1/2)*m,-2*b*x^2)*exp(2a))
@test integrate((e*x)^m*sinh(a+b*x^2), x) == :(-1//4*e^-1*(e*x)^(1+m)*(-1*b*x^2)^(-1//2+-1//2*m)*SpecialFunctions.gamma(1/2+(1/2)*m,-1*b*x^2)*exp(a)+1//4*e^-1*(b*x^2)^(-1//2+-1//2*m)*(e*x)^(1+m)*SpecialFunctions.gamma(1/2+(1/2)*m,b*x^2)*exp(-1a))
@test integrate(x^3*sinh(a+b*x^4), x) == :(1//4*b^-1*cosh(a+b*x^4))
@test integrate(x^2*sinh(a+b*x^-1), x) == :(1//3*x^3*sinh(a+b*x^-1)+-1//6*b^3*cosh(a)*coshintegral(b*x^-1)+-1//6*b^3*sinh(a)*sinhintegral(b*x^-1)+1//6*b*x^2*cosh(a+b*x^-1)+1//6*x*b^2*sinh(a+b*x^-1))
@test integrate(x*sinh(a+b*x^-1), x) == :((1/2)*x^2*sinh(a+b*x^-1)+(1/2)*b*x*cosh(a+b*x^-1)+-1//2*b^2*cosh(a)*sinhintegral(b*x^-1)+-1//2*b^2*coshintegral(b*x^-1)*sinh(a))
@test integrate(sinh(a+b*x^-1), x) == :(x*sinh(a+b*x^-1)+-1*b*cosh(a)*coshintegral(b*x^-1)+-1*b*sinh(a)*sinhintegral(b*x^-1))
@test integrate(x^-1*sinh(a+b*x^-1), x) == :(-1*cosh(a)*sinhintegral(b*x^-1)+-1*coshintegral(b*x^-1)*sinh(a))
@test integrate(x^-2*sinh(a+b*x^-1), x) == :(-1*b^-1*cosh(a+b*x^-1))
@test integrate(x^-3*sinh(a+b*x^-1), x) == :(b^-2*sinh(a+b*x^-1)+-1*b^-1*x^-1*cosh(a+b*x^-1))
@test integrate(x^-4*sinh(a+b*x^-1), x) == :(-2*b^-3*cosh(a+b*x^-1)+-1*b^-1*x^-2*cosh(a+b*x^-1)+2*b^-2*x^-1*sinh(a+b*x^-1))
@test integrate(x^-5*sinh(a+b*x^-1), x) == :(6*b^-4*sinh(a+b*x^-1)+-1*b^-1*x^-3*cosh(a+b*x^-1)+-6*b^-3*x^-1*cosh(a+b*x^-1)+3*b^-2*x^-2*sinh(a+b*x^-1))
@test integrate((e*x)^m*sinh(a+b*x^-1)^3, x) == :(3//8*b*(b*x^-1)^m*(e*x)^m*SpecialFunctions.gamma(-1+-1m,b*x^-1)*exp(-1a)+3//8*b*(e*x)^m*(-1*b*x^-1)^m*SpecialFunctions.gamma(-1+-1m,-1*b*x^-1)*exp(a)+-1//8*b*3^(1+m)*(b*x^-1)^m*(e*x)^m*SpecialFunctions.gamma(-1+-1m,3*b*x^-1)*exp(-3a)+-1//8*b*3^(1+m)*(e*x)^m*(-1*b*x^-1)^m*SpecialFunctions.gamma(-1+-1m,-3*b*x^-1)*exp(3a))
@test integrate((e*x)^m*sinh(a+b*x^-1)^2, x) == :(-1*x*(e*x)^m*(2+2m)^-1+b*2^(-1+m)*(b*x^-1)^m*(e*x)^m*SpecialFunctions.gamma(-1+-1m,2*b*x^-1)*exp(-2a)+-1*b*2^(-1+m)*(e*x)^m*(-1*b*x^-1)^m*SpecialFunctions.gamma(-1+-1m,-2*b*x^-1)*exp(2a))
@test integrate((e*x)^m*sinh(a+b*x^-1), x) == :(-1//2*b*(b*x^-1)^m*(e*x)^m*SpecialFunctions.gamma(-1+-1m,b*x^-1)*exp(-1a)+-1//2*b*(e*x)^m*(-1*b*x^-1)^m*SpecialFunctions.gamma(-1+-1m,-1*b*x^-1)*exp(a))
@test integrate(x^4*sinh(a+b*x^-2), x) == :(1//5*x^5*sinh(a+b*x^-2)+2//15*b*x^3*cosh(a+b*x^-2)+4//15*x*b^2*sinh(a+b*x^-2)+-2//15*pi^(1/2)*b^5//2*erf(b^(1/2)*x^-1)*exp(-1a)+-2//15*pi^(1/2)*b^5//2*erfi(b^(1/2)*x^-1)*exp(a))
@test integrate(x^3*sinh(a+b*x^-2), x) == :(1//4*x^4*sinh(a+b*x^-2)+-1//4*b^2*cosh(a)*sinhintegral(b*x^-2)+-1//4*b^2*coshintegral(b*x^-2)*sinh(a)+1//4*b*x^2*cosh(a+b*x^-2))
@test integrate(x^2*sinh(a+b*x^-2), x) == :(1//3*x^3*sinh(a+b*x^-2)+2//3*b*x*cosh(a+b*x^-2)+-1//3*pi^(1/2)*b^3//2*erfi(b^(1/2)*x^-1)*exp(a)+1//3*pi^(1/2)*b^3//2*erf(b^(1/2)*x^-1)*exp(-1a))
@test integrate(x*sinh(a+b*x^-2), x) == :((1/2)*x^2*sinh(a+b*x^-2)+-1//2*b*cosh(a)*coshintegral(b*x^-2)+-1//2*b*sinh(a)*sinhintegral(b*x^-2))
@test integrate(sinh(a+b*x^-2), x) == :(x*sinh(a+b*x^-2)+-1//2*pi^(1/2)*b^(1/2)*erf(b^(1/2)*x^-1)*exp(-1a)+-1//2*pi^(1/2)*b^(1/2)*erfi(b^(1/2)*x^-1)*exp(a))
@test integrate(x^-1*sinh(a+b*x^-2), x) == :(-1//2*cosh(a)*sinhintegral(b*x^-2)+-1//2*coshintegral(b*x^-2)*sinh(a))
@test integrate(x^-2*sinh(a+b*x^-2), x) == :(-1//4*pi^(1/2)*b^-1//2*erfi(b^(1/2)*x^-1)*exp(a)+1//4*pi^(1/2)*b^-1//2*erf(b^(1/2)*x^-1)*exp(-1a))
@test integrate(x^-3*sinh(a+b*x^-2), x) == :(-1//2*b^-1*cosh(a+b*x^-2))
@test integrate(x^-4*sinh(a+b*x^-2), x) == :(-1//2*b^-1*x^-1*cosh(a+b*x^-2)+1//8*pi^(1/2)*b^-3//2*erf(b^(1/2)*x^-1)*exp(-1a)+1//8*pi^(1/2)*b^-3//2*erfi(b^(1/2)*x^-1)*exp(a))
@test integrate(x^-5*sinh(a+b*x^-2), x) == :((1/2)*b^-2*sinh(a+b*x^-2)+-1//2*b^-1*x^-2*cosh(a+b*x^-2))
@test integrate(x^-6*sinh(a+b*x^-2), x) == :(-1//2*b^-1*x^-3*cosh(a+b*x^-2)+3//4*b^-2*x^-1*sinh(a+b*x^-2)+-3//16*pi^(1/2)*b^-5//2*erfi(b^(1/2)*x^-1)*exp(a)+3//16*pi^(1/2)*b^-5//2*erf(b^(1/2)*x^-1)*exp(-1a))
@test integrate(x^-7*sinh(a+b*x^-2), x) == :(-1*b^-3*cosh(a+b*x^-2)+b^-2*x^-2*sinh(a+b*x^-2)+-1//2*b^-1*x^-4*cosh(a+b*x^-2))
@test integrate((e*x)^m*sinh(a+b*x^-2)^3, x) == :(-3//16*x*(e*x)^m*(-1*b*x^-2)^(1/2+(1/2)*m)*SpecialFunctions.gamma(-1//2+-1//2*m,-1*b*x^-2)*exp(a)+3//16*x*(b*x^-2)^(1/2+(1/2)*m)*(e*x)^m*SpecialFunctions.gamma(-1//2+-1//2*m,b*x^-2)*exp(-1a)+-1//16*x*3^(1/2+(1/2)*m)*(b*x^-2)^(1/2+(1/2)*m)*(e*x)^m*SpecialFunctions.gamma(-1//2+-1//2*m,3*b*x^-2)*exp(-3a)+1//16*x*3^(1/2+(1/2)*m)*(e*x)^m*(-1*b*x^-2)^(1/2+(1/2)*m)*SpecialFunctions.gamma(-1//2+-1//2*m,-3*b*x^-2)*exp(3a))
@test integrate((e*x)^m*sinh(a+b*x^-2)^2, x) == :(-1*x*(e*x)^m*(2+2m)^-1+x*2^(-5//2+(1/2)*m)*(b*x^-2)^(1/2+(1/2)*m)*(e*x)^m*SpecialFunctions.gamma(-1//2+-1//2*m,2*b*x^-2)*exp(-2a)+x*2^(-5//2+(1/2)*m)*(e*x)^m*(-1*b*x^-2)^(1/2+(1/2)*m)*SpecialFunctions.gamma(-1//2+-1//2*m,-2*b*x^-2)*exp(2a))
@test integrate((e*x)^m*sinh(a+b*x^-2), x) == :(-1//4*x*(b*x^-2)^(1/2+(1/2)*m)*(e*x)^m*SpecialFunctions.gamma(-1//2+-1//2*m,b*x^-2)*exp(-1a)+1//4*x*(e*x)^m*(-1*b*x^-2)^(1/2+(1/2)*m)*SpecialFunctions.gamma(-1//2+-1//2*m,-1*b*x^-2)*exp(a))
@test integrate(x^-1//2*sinh(x^(1/2)), x) == :(2*cosh(x^(1/2)))
@test integrate(x^2*sinh(a+b*x^n), x) == :((1/2)*n^-1*x^3*(b*x^n)^(-3*n^-1)*SpecialFunctions.gamma(3*n^-1,b*x^n)*exp(-1a)+-1//2*n^-1*x^3*(-1*b*x^n)^(-3*n^-1)*SpecialFunctions.gamma(3*n^-1,-1*b*x^n)*exp(a))
@test integrate(x*sinh(a+b*x^n), x) == :((1/2)*n^-1*x^2*(b*x^n)^(-2*n^-1)*SpecialFunctions.gamma(2*n^-1,b*x^n)*exp(-1a)+-1//2*n^-1*x^2*(-1*b*x^n)^(-2*n^-1)*SpecialFunctions.gamma(2*n^-1,-1*b*x^n)*exp(a))
@test integrate(sinh(a+b*x^n), x) == :((1/2)*x*n^-1*(b*x^n)^(-1*n^-1)*SpecialFunctions.gamma(n^-1,b*x^n)*exp(-1a)+-1//2*x*n^-1*(-1*b*x^n)^(-1*n^-1)*SpecialFunctions.gamma(n^-1,-1*b*x^n)*exp(a))
@test integrate(x^-1*sinh(a+b*x^n), x) == :(n^-1*cosh(a)*sinhintegral(b*x^n)+n^-1*coshintegral(b*x^n)*sinh(a))
@test integrate(x^-2*sinh(a+b*x^n), x) == :((1/2)*n^-1*x^-1*(b*x^n)^(n^-1)*SpecialFunctions.gamma(-1*n^-1,b*x^n)*exp(-1a)+-1//2*n^-1*x^-1*(-1*b*x^n)^(n^-1)*SpecialFunctions.gamma(-1*n^-1,-1*b*x^n)*exp(a))
@test integrate(x^-3*sinh(a+b*x^n), x) == :((1/2)*n^-1*x^-2*(b*x^n)^(2*n^-1)*SpecialFunctions.gamma(-2*n^-1,b*x^n)*exp(-1a)+-1//2*n^-1*x^-2*(-1*b*x^n)^(2*n^-1)*SpecialFunctions.gamma(-2*n^-1,-1*b*x^n)*exp(a))
@test integrate(x^2*sinh(a+b*x^n)^2, x) == :(-1//6*x^3+-1*2^(-2+-3*n^-1)*n^-1*x^3*(b*x^n)^(-3*n^-1)*SpecialFunctions.gamma(3*n^-1,2*b*x^n)*exp(-2a)+-1*2^(-2+-3*n^-1)*n^-1*x^3*(-1*b*x^n)^(-3*n^-1)*SpecialFunctions.gamma(3*n^-1,-2*b*x^n)*exp(2a))
@test integrate(x*sinh(a+b*x^n)^2, x) == :(-1//4*x^2+-1*4^(-1+-1*n^-1)*n^-1*x^2*(b*x^n)^(-2*n^-1)*SpecialFunctions.gamma(2*n^-1,2*b*x^n)*exp(-2a)+-1*4^(-1+-1*n^-1)*n^-1*x^2*(-1*b*x^n)^(-2*n^-1)*SpecialFunctions.gamma(2*n^-1,-2*b*x^n)*exp(2a))
@test integrate(sinh(a+b*x^n)^2, x) == :(-1//2*x+-1*x*2^(-2+-1*n^-1)*n^-1*(b*x^n)^(-1*n^-1)*SpecialFunctions.gamma(n^-1,2*b*x^n)*exp(-2a)+-1*x*2^(-2+-1*n^-1)*n^-1*(-1*b*x^n)^(-1*n^-1)*SpecialFunctions.gamma(n^-1,-2*b*x^n)*exp(2a))
@test integrate(x^-1*sinh(a+b*x^n)^2, x) == :(-1//2*log(x)+(1/2)*n^-1*cosh(2a)*coshintegral(2*b*x^n)+(1/2)*n^-1*sinh(2a)*sinhintegral(2*b*x^n))
@test integrate(x^-2*sinh(a+b*x^n)^2, x) == :((1/2)*x^-1+-1*2^(-2+n^-1)*n^-1*x^-1*(b*x^n)^(n^-1)*SpecialFunctions.gamma(-1*n^-1,2*b*x^n)*exp(-2a)+-1*2^(-2+n^-1)*n^-1*x^-1*(-1*b*x^n)^(n^-1)*SpecialFunctions.gamma(-1*n^-1,-2*b*x^n)*exp(2a))
@test integrate(x^2*sinh(a+b*x^n)^3, x) == :(-3//8*n^-1*x^3*(b*x^n)^(-3*n^-1)*SpecialFunctions.gamma(3*n^-1,b*x^n)*exp(-1a)+3//8*n^-1*x^3*(-1*b*x^n)^(-3*n^-1)*SpecialFunctions.gamma(3*n^-1,-1*b*x^n)*exp(a)+-1//8*3^(-3*n^-1)*n^-1*x^3*(-1*b*x^n)^(-3*n^-1)*SpecialFunctions.gamma(3*n^-1,-3*b*x^n)*exp(3a)+1//8*3^(-3*n^-1)*n^-1*x^3*(b*x^n)^(-3*n^-1)*SpecialFunctions.gamma(3*n^-1,3*b*x^n)*exp(-3a))
@test integrate(x*sinh(a+b*x^n)^3, x) == :(-3//8*n^-1*x^2*(b*x^n)^(-2*n^-1)*SpecialFunctions.gamma(2*n^-1,b*x^n)*exp(-1a)+3//8*n^-1*x^2*(-1*b*x^n)^(-2*n^-1)*SpecialFunctions.gamma(2*n^-1,-1*b*x^n)*exp(a)+-1//8*9^(-1*n^-1)*n^-1*x^2*(-1*b*x^n)^(-2*n^-1)*SpecialFunctions.gamma(2*n^-1,-3*b*x^n)*exp(3a)+1//8*9^(-1*n^-1)*n^-1*x^2*(b*x^n)^(-2*n^-1)*SpecialFunctions.gamma(2*n^-1,3*b*x^n)*exp(-3a))
@test integrate(sinh(a+b*x^n)^3, x) == :(-3//8*x*n^-1*(b*x^n)^(-1*n^-1)*SpecialFunctions.gamma(n^-1,b*x^n)*exp(-1a)+3//8*x*n^-1*(-1*b*x^n)^(-1*n^-1)*SpecialFunctions.gamma(n^-1,-1*b*x^n)*exp(a)+-1//8*x*3^(-1*n^-1)*n^-1*(-1*b*x^n)^(-1*n^-1)*SpecialFunctions.gamma(n^-1,-3*b*x^n)*exp(3a)+1//8*x*3^(-1*n^-1)*n^-1*(b*x^n)^(-1*n^-1)*SpecialFunctions.gamma(n^-1,3*b*x^n)*exp(-3a))
@test integrate(x^-1*sinh(a+b*x^n)^3, x) == :(-3//4*n^-1*cosh(a)*sinhintegral(b*x^n)+-3//4*n^-1*coshintegral(b*x^n)*sinh(a)+1//4*n^-1*cosh(3a)*sinhintegral(3*b*x^n)+1//4*n^-1*coshintegral(3*b*x^n)*sinh(3a))
@test integrate(x^-2*sinh(a+b*x^n)^3, x) == :(-3//8*n^-1*x^-1*(b*x^n)^(n^-1)*SpecialFunctions.gamma(-1*n^-1,b*x^n)*exp(-1a)+3//8*n^-1*x^-1*(-1*b*x^n)^(n^-1)*SpecialFunctions.gamma(-1*n^-1,-1*b*x^n)*exp(a)+-1//8*3^(n^-1)*n^-1*x^-1*(-1*b*x^n)^(n^-1)*SpecialFunctions.gamma(-1*n^-1,-3*b*x^n)*exp(3a)+1//8*3^(n^-1)*n^-1*x^-1*(b*x^n)^(n^-1)*SpecialFunctions.gamma(-1*n^-1,3*b*x^n)*exp(-3a))
@test integrate((e*x)^m*sinh(a+b*x^n)^3, x) == :(-3//8*e^-1*n^-1*(b*x^n)^(-1*n^-1*(1+m))*(e*x)^(1+m)*SpecialFunctions.gamma(n^-1*(1+m),b*x^n)*exp(-1a)+3//8*e^-1*n^-1*(e*x)^(1+m)*(-1*b*x^n)^(-1*n^-1*(1+m))*SpecialFunctions.gamma(n^-1*(1+m),-1*b*x^n)*exp(a)+-1//8*3^(-1*n^-1*(1+m))*e^-1*n^-1*(e*x)^(1+m)*(-1*b*x^n)^(-1*n^-1*(1+m))*SpecialFunctions.gamma(n^-1*(1+m),-3*b*x^n)*exp(3a)+1//8*3^(-1*n^-1*(1+m))*e^-1*n^-1*(b*x^n)^(-1*n^-1*(1+m))*(e*x)^(1+m)*SpecialFunctions.gamma(n^-1*(1+m),3*b*x^n)*exp(-3a))
@test integrate((e*x)^m*sinh(a+b*x^n)^2, x) == :(-1//2*e^-1*(e*x)^(1+m)*(1+m)^-1+-1*2^(-1*n^-1*(1+m+2n))*e^-1*n^-1*(b*x^n)^(-1*n^-1*(1+m))*(e*x)^(1+m)*SpecialFunctions.gamma(n^-1*(1+m),2*b*x^n)*exp(-2a)+-1*2^(-1*n^-1*(1+m+2n))*e^-1*n^-1*(e*x)^(1+m)*(-1*b*x^n)^(-1*n^-1*(1+m))*SpecialFunctions.gamma(n^-1*(1+m),-2*b*x^n)*exp(2a))
@test integrate((e*x)^m*sinh(a+b*x^n), x) == :((1/2)*e^-1*n^-1*(b*x^n)^(-1*n^-1*(1+m))*(e*x)^(1+m)*SpecialFunctions.gamma(n^-1*(1+m),b*x^n)*exp(-1a)+-1//2*e^-1*n^-1*(e*x)^(1+m)*(-1*b*x^n)^(-1*n^-1*(1+m))*SpecialFunctions.gamma(n^-1*(1+m),-1*b*x^n)*exp(a))
@test integrate(x^(-1+-1n)*sinh(a+b*x^n), x) == :(-1*n^-1*x^(-1n)*sinh(a+b*x^n)+b*n^-1*cosh(a)*coshintegral(b*x^n)+b*n^-1*sinh(a)*sinhintegral(b*x^n))
@test integrate(x^(-1+-1n)*sinh(a+b*x^n)^2, x) == :((1/2)*n^-1*x^(-1n)+-1//2*n^-1*x^(-1n)*cosh(2a+2*b*x^n)+b*n^-1*cosh(2a)*sinhintegral(2*b*x^n)+b*n^-1*coshintegral(2*b*x^n)*sinh(2a))
@test integrate(x^(-1+-1n)*sinh(a+b*x^n)^3, x) == :(-1//4*n^-1*x^(-1n)*sinh(3a+3*b*x^n)+3//4*n^-1*x^(-1n)*sinh(a+b*x^n)+-3//4*b*n^-1*cosh(a)*coshintegral(b*x^n)+-3//4*b*n^-1*sinh(a)*sinhintegral(b*x^n)+3//4*b*n^-1*cosh(3a)*coshintegral(3*b*x^n)+3//4*b*n^-1*sinh(3a)*sinhintegral(3*b*x^n))
@test integrate(x^(-1+(1/2)*n)*sinh(a+b*x^n), x) == :((1/2)*pi^(1/2)*b^-1//2*n^-1*erfi(b^(1/2)*x^((1/2)*n))*exp(a)+-1//2*pi^(1/2)*b^-1//2*n^-1*erf(b^(1/2)*x^((1/2)*n))*exp(-1a))
@test integrate(x^2*sinh((a+b*x)^2), x) == :((1/2)*b^-3*(a+b*x)*cosh((a+b*x)^2)+-1*a*b^-3*cosh((a+b*x)^2)+-1//8*pi^(1/2)*b^-3*erf(a+b*x)+-1//8*pi^(1/2)*b^-3*erfi(a+b*x)+-1//4*pi^(1/2)*a^2*b^-3*erf(a+b*x)+1//4*pi^(1/2)*a^2*b^-3*erfi(a+b*x))
@test integrate(x*sinh((a+b*x)^2), x) == :((1/2)*b^-2*cosh((a+b*x)^2)+-1//4*a*pi^(1/2)*b^-2*erfi(a+b*x)+1//4*a*pi^(1/2)*b^-2*erf(a+b*x))
@test integrate(sinh((a+b*x)^2), x) == :(-1//4*pi^(1/2)*b^-1*erf(a+b*x)+1//4*pi^(1/2)*b^-1*erfi(a+b*x))
@test integrate(x^2*sinh(a+b*(c+d*x)^(1/2)), x) == :(-240*b^-6*d^-3*sinh(a+b*(c+d*x)^(1/2))+-1*b^-4*d^-3*(120c+120*d*x)*sinh(a+b*(c+d*x)^(1/2))+-10*b^-2*d^-3*(c+d*x)^2*sinh(a+b*(c+d*x)^(1/2))+-2*b^-2*c^2*d^-3*sinh(a+b*(c+d*x)^(1/2))+2*b^-1*d^-3*(c+d*x)^5//2*cosh(a+b*(c+d*x)^(1/2))+24*c*b^-4*d^-3*sinh(a+b*(c+d*x)^(1/2))+40*b^-3*d^-3*(c+d*x)^3//2*cosh(a+b*(c+d*x)^(1/2))+240*b^-5*d^-3*(c+d*x)^(1/2)*cosh(a+b*(c+d*x)^(1/2))+-24*c*b^-3*d^-3*(c+d*x)^(1/2)*cosh(a+b*(c+d*x)^(1/2))+-4*c*b^-1*d^-3*(c+d*x)^3//2*cosh(a+b*(c+d*x)^(1/2))+2*b^-1*c^2*d^-3*(c+d*x)^(1/2)*cosh(a+b*(c+d*x)^(1/2))+12*c*b^-2*d^-3*(c+d*x)*sinh(a+b*(c+d*x)^(1/2)))
@test integrate(x*sinh(a+b*(c+d*x)^(1/2)), x) == :(-12*b^-4*d^-2*sinh(a+b*(c+d*x)^(1/2))+-1*b^-2*d^-2*(6c+6*d*x)*sinh(a+b*(c+d*x)^(1/2))+2*c*b^-2*d^-2*sinh(a+b*(c+d*x)^(1/2))+2*b^-1*d^-2*(c+d*x)^3//2*cosh(a+b*(c+d*x)^(1/2))+12*b^-3*d^-2*(c+d*x)^(1/2)*cosh(a+b*(c+d*x)^(1/2))+-2*c*b^-1*d^-2*(c+d*x)^(1/2)*cosh(a+b*(c+d*x)^(1/2)))
@test integrate(sinh(a+b*(c+d*x)^(1/2)), x) == :(-2*b^-2*d^-1*sinh(a+b*(c+d*x)^(1/2))+2*b^-1*d^-1*(c+d*x)^(1/2)*cosh(a+b*(c+d*x)^(1/2)))
@test integrate(x^-1*sinh(a+b*(c+d*x)^(1/2)), x) == :(cosh(a+-1*b*c^(1/2))*sinhintegral(b*(c^(1/2)+(c+d*x)^(1/2)))+coshintegral(b*(c^(1/2)+(c+d*x)^(1/2)))*sinh(a+-1*b*c^(1/2))+coshintegral(b*(c^(1/2)+-1*(c+d*x)^(1/2)))*sinh(a+b*c^(1/2))+-1*cosh(a+b*c^(1/2))*sinhintegral(b*(c^(1/2)+-1*(c+d*x)^(1/2))))
@test integrate(x^-2*sinh(a+b*(c+d*x)^(1/2)), x) == :(-1*x^-1*sinh(a+b*(c+d*x)^(1/2))+(1/2)*b*d*c^-1//2*cosh(a+b*c^(1/2))*coshintegral(b*(c^(1/2)+-1*(c+d*x)^(1/2)))+-1//2*b*d*c^-1//2*cosh(a+-1*b*c^(1/2))*coshintegral(b*(c^(1/2)+(c+d*x)^(1/2)))+-1//2*b*d*c^-1//2*sinh(a+b*c^(1/2))*sinhintegral(b*(c^(1/2)+-1*(c+d*x)^(1/2)))+-1//2*b*d*c^-1//2*sinh(a+-1*b*c^(1/2))*sinhintegral(b*(c^(1/2)+(c+d*x)^(1/2))))
@test integrate(x^2*sinh(a+b*(c+d*x)^1//3), x) == :(120960*b^-9*d^-3*cosh(a+b*(c+d*x)^1//3)+-1*b^-6*d^-3*(20160c+20160*d*x)*sinh(a+b*(c+d*x)^1//3)+-120960*b^-8*d^-3*(c+d*x)^1//3*sinh(a+b*(c+d*x)^1//3)+-1008*b^-4*d^-3*(c+d*x)^5//3*sinh(a+b*(c+d*x)^1//3)+-24*b^-2*d^-3*(c+d*x)^7//3*sinh(a+b*(c+d*x)^1//3)+3*b^-1*d^-3*(c+d*x)^8//3*cosh(a+b*(c+d*x)^1//3)+6*b^-3*c^2*d^-3*cosh(a+b*(c+d*x)^1//3)+168*b^-3*d^-3*(c+d*x)^2*cosh(a+b*(c+d*x)^1//3)+720*c*b^-6*d^-3*sinh(a+b*(c+d*x)^1//3)+5040*b^-5*d^-3*(c+d*x)^4//3*cosh(a+b*(c+d*x)^1//3)+60480*b^-7*d^-3*(c+d*x)^2//3*cosh(a+b*(c+d*x)^1//3)+-720*c*b^-5*d^-3*(c+d*x)^1//3*cosh(a+b*(c+d*x)^1//3)+-120*c*b^-3*d^-3*(c+d*x)*cosh(a+b*(c+d*x)^1//3)+-6*c*b^-1*d^-3*(c+d*x)^5//3*cosh(a+b*(c+d*x)^1//3)+-6*b^-2*c^2*d^-3*(c+d*x)^1//3*sinh(a+b*(c+d*x)^1//3)+3*b^-1*c^2*d^-3*(c+d*x)^2//3*cosh(a+b*(c+d*x)^1//3)+30*c*b^-2*d^-3*(c+d*x)^4//3*sinh(a+b*(c+d*x)^1//3)+360*c*b^-4*d^-3*(c+d*x)^2//3*sinh(a+b*(c+d*x)^1//3))
@test integrate(x*sinh(a+b*(c+d*x)^1//3), x) == :(-360*b^-6*d^-2*sinh(a+b*(c+d*x)^1//3)+b^-3*d^-2*(60c+60*d*x)*cosh(a+b*(c+d*x)^1//3)+-180*b^-4*d^-2*(c+d*x)^2//3*sinh(a+b*(c+d*x)^1//3)+-15*b^-2*d^-2*(c+d*x)^4//3*sinh(a+b*(c+d*x)^1//3)+-6*c*b^-3*d^-2*cosh(a+b*(c+d*x)^1//3)+3*b^-1*d^-2*(c+d*x)^5//3*cosh(a+b*(c+d*x)^1//3)+360*b^-5*d^-2*(c+d*x)^1//3*cosh(a+b*(c+d*x)^1//3)+-3*c*b^-1*d^-2*(c+d*x)^2//3*cosh(a+b*(c+d*x)^1//3)+6*c*b^-2*d^-2*(c+d*x)^1//3*sinh(a+b*(c+d*x)^1//3))
@test integrate(sinh(a+b*(c+d*x)^1//3), x) == :(6*b^-3*d^-1*cosh(a+b*(c+d*x)^1//3)+-6*b^-2*d^-1*(c+d*x)^1//3*sinh(a+b*(c+d*x)^1//3)+3*b^-1*d^-1*(c+d*x)^2//3*cosh(a+b*(c+d*x)^1//3))
@test integrate(x^-1*sinh(a+b*(c+d*x)^1//3), x) == :(cosh(a+-1*b*(-1)^1//3*c^1//3)*sinhintegral(b*((c+d*x)^1//3+(-1)^1//3*c^1//3))+coshintegral(b*(c^1//3+-1*(c+d*x)^1//3))*sinh(a+b*c^1//3)+coshintegral(b*((c+d*x)^1//3+(-1)^1//3*c^1//3))*sinh(a+-1*b*(-1)^1//3*c^1//3)+coshintegral(-1*b*(-1*(c+d*x)^1//3+(-1)^2//3*c^1//3))*sinh(a+b*(-1)^2//3*c^1//3)+-1*cosh(a+b*c^1//3)*sinhintegral(b*(c^1//3+-1*(c+d*x)^1//3))+-1*cosh(a+b*(-1)^2//3*c^1//3)*sinhintegral(b*(-1*(c+d*x)^1//3+(-1)^2//3*c^1//3)))
@test integrate(x^-2*sinh(a+b*(c+d*x)^1//3), x) == :(-1*x^-1*sinh(a+b*(c+d*x)^1//3)+-1//3*b*d*c^-2//3*sinh(a+b*c^1//3)*sinhintegral(b*(c^1//3+-1*(c+d*x)^1//3))+1//3*b*d*c^-2//3*cosh(a+b*c^1//3)*coshintegral(b*(c^1//3+-1*(c+d*x)^1//3))+-1//3*b*d*(-1)^1//3*c^-2//3*cosh(a+-1*b*(-1)^1//3*c^1//3)*coshintegral(b*((c+d*x)^1//3+(-1)^1//3*c^1//3))+-1//3*b*d*(-1)^1//3*c^-2//3*sinh(a+-1*b*(-1)^1//3*c^1//3)*sinhintegral(b*((c+d*x)^1//3+(-1)^1//3*c^1//3))+-1//3*b*d*(-1)^2//3*c^-2//3*sinh(a+b*(-1)^2//3*c^1//3)*sinhintegral(b*(-1*(c+d*x)^1//3+(-1)^2//3*c^1//3))+1//3*b*d*(-1)^2//3*c^-2//3*cosh(a+b*(-1)^2//3*c^1//3)*coshintegral(-1*b*(-1*(c+d*x)^1//3+(-1)^2//3*c^1//3)))

max_stars_repo_path: test/6 Hyperbolic functions/6.1 Hyperbolic sine/6.1.3 (e x)^m (a+b sinh(c+d x^n))^p.jl


text:
if (Base.libllvm_version ≥ v"7" && VectorizationBase.AVX512F) || Base.libllvm_version ≥ v"9"
    function vfilter!(f::F, x::Vector{T}, y::AbstractArray{T}) where {F,T <: NativeTypes}
        W, Wshift = VectorizationBase.pick_vector_width_shift(T)
        N = length(y)
        Nrep = N >>> Wshift
        Nrem = N & (W - 1)
        j = 0
        GC.@preserve x y begin
            ptr_x = pointer(x)
            ptr_y = pointer(y)
            for _ ∈ 1:Nrep
                vy = vload(Vec{W,T}, ptr_y)
                mask = f(SVec(vy))
                SIMDPirates.compressstore!(gep(ptr_x, j), vy, mask)
                ptr_y = gepbyte(ptr_y, VectorizationBase.REGISTER_SIZE)
                j = vadd(j, count_ones(mask))
            end
            rem_mask = VectorizationBase.mask(T, Nrem)
            vy = vload(Vec{W,T}, ptr_y, rem_mask)
            mask = rem_mask & f(SVec(vy))
            SIMDPirates.compressstore!(gep(ptr_x, j), vy, mask)
            j = vadd(j, count_ones(mask))
            Base._deleteend!(x, N-j) # resize!(x, j)
        end
        x
    end
    vfilter!(f::F, x::Vector{T}) where {F, T<:NativeTypes} = vfilter!(f, x, x)
    vfilter(f::F, y::AbstractArray{T}) where {F, T<:NativeTypes} = vfilter!(f, Vector{T}(undef, length(y)), y)
end
vfilter(f::F, y) where {F} = filter(f, y)
vfilter!(f::F, y) where {F} = filter!(f, y)

"""
    vfilter(f, a::AbstractArray)

SIMD-vectorized `filter`, returning an array containing the elements of `a` for which `f` return `true`.
"""
vfilter

"""
    vfilter!(f, a::AbstractArray)

SIMD-vectorized `filter!`, removing the element of `a` for which `f` is false.
"""
vfilter!


text:
function _solve_tsp(name; log="off", kwargs...)
    if log == "off" || Sys.iswindows()
        status = run(`$(LKH.LKH_EXECUTABLE) $(name).par`, wait=false)
        while !success(status)
            # Wait until it finishes
            # In Windows, it stops with "Press any key to continue . . . "
        end
    else
        run(`$(LKH.LKH_EXECUTABLE) $(name).par`)
    end

    opt_tour, opt_len = read_output(name)

    cleanup(name)

    return opt_tour, opt_len
end

function solve_tsp(dist_mtx::Matrix{Int}; log="off", kwargs...)
    if issymmetric(dist_mtx)
        name, tsp_filepath = write_tsp(dist_mtx)
    else
        name, tsp_filepath = write_atsp(dist_mtx)
    end

    write_par(name; kwargs...)
    return _solve_tsp(name; log=log, kwargs...)
end

function solve_tsp(x::Vector{Float64}, y::Vector{Float64}; dist="EUC_2D", log="off", kwargs...)
    name, tsp_filepath = write_tsp(x, y; dist=dist)
    write_par(name; kwargs...)
    return _solve_tsp(name; log=log, kwargs...)
end

function solve_tsp(tsp_file::String; log="off", kwargs...)
    if !isfile(tsp_file)
        error("$(tsp_file) is not a file.")
    end

    name = randstring(10)
    filepath = name * ".tsp"
    cp(tsp_file, filepath)

    write_par(name; kwargs...)
    return _solve_tsp(name; log=log, kwargs...)
end


text:
module dendogram

import Term: fint, int, cint, truncate, loop_firstlast, highlight, textlen

import ..box: get_rrow, get_lrow, get_row, SQUARE
import ..layout: pad
import ..segment: Segment
import ..measure: Measure
import ..renderables: AbstractRenderable
import ..style: apply_style

import MyterialColors: yellow, blue_light, green_light, salmon_light, blue_grey_light

export Dendogram, link

PRETITLE_STYLE = blue_grey_light
LEAVES_STYLE = blue_grey_light
LINES_STYLE = blue_light * " dim bold"
CELLWIDTH = 10
SPACING = 1

# ---------------------------------------------------------------------------- #
#                                     LEAF                                     #
# ---------------------------------------------------------------------------- #

"""
    Leaf <: AbstractRenderable

The terminal element of a `Dendogram`.
"""
struct Leaf <: AbstractRenderable
    segments::Vector{Segment}
    measure::Measure
    text::String
    midpoint::Int
end

"""
    Leaf(leaf)

Construct a `Leaf` out of any object.
"""
function Leaf(leaf)
    # get a string representation of the appropriate length
    leaf = string(leaf)
    if textlen(leaf) > CELLWIDTH
        leaf = truncate(leaf, CELLWIDTH)
    else
        leaf = pad(leaf, CELLWIDTH+1, :center)
    end

    # create renderable Leaf
    midpoint = fint(textlen(leaf)/2)
    seg = Segment(" "*leaf*" ", LEAVES_STYLE)
    return Leaf([seg], Measure(seg), leaf, midpoint)
end

function Base.string(leaf::Leaf, isfirst::Bool, islast::Bool, spacing::Int)
    l = isfirst ? "" : " "^spacing
    r = islast ? "" : " "^spacing
    return l * leaf.text * r
end

# ---------------------------------------------------------------------------- #
#                                   DENDOGRAM                                  #
# ---------------------------------------------------------------------------- #

"""
    Dendogram <: AbstractRenderable

A Dendogram tree renderable.
"""
struct Dendogram <: AbstractRenderable
    segments::Vector{Segment}
    measure::Measure
    midpoint::Int  # width of 'center'
end

"""
    Dendogram(head, args::Vector; first_arg=nothing, pretitle=nothing)

Construct a single `Dendogram`.

Construct a dendogram with one head node (`head`) and any number of leaves.

`first_arg` is used to create dendograms for `Expr` objects and gives special
status to an expression's first argument which gets printed with the head.
`pretitle` is used to create dendograms for `Expr`, pretitle is insrted as
three `Segment`s (a string and an upward arrow) before the head of the dendogram.
"""
function Dendogram(head, args::Vector; first_arg=nothing, pretitle=nothing)
    # get leaves
    leaves = Leaf.(args)
    leaves_line = join(
        map(
            nl -> string(nl[3], nl[1], nl[2], SPACING),
            loop_firstlast(leaves)
            )
        )
    width = textlen(leaves_line)

    # get tree structure line
    if length(leaves) > 1
        widths = repeat([CELLWIDTH+1 + SPACING], length(leaves)-1)

        line = get_row(SQUARE, widths, :top)
        line = pad(replace_line_midpoint(line)[1], width, :center)
    else
        widths = [CELLWIDTH]
        line = pad(string(SQUARE.bottom.vertical), CELLWIDTH, :center)
    end

    # get title
    if isnothing(first_arg)
        title = pad(apply_style("$(head)", salmon_light), width, :center)

    else
        _title = ": [bold default underline $salmon_light]$first_arg[/bold default underline $salmon_light]"
        title = pad(apply_style("$(head)$_title", salmon_light * " dim"), width-4, :center) * " "^4
    end

    # put together
    segments = [
        Segment(title),
        Segment(line, LINES_STYLE),
        Segment(leaves_line, LEAVES_STYLE)
    ]

    # add 'pretitle' lines (for expressions only)
    if !isnothing(pretitle)
        prepend!(
            segments,
            [
                Segment(" " * pad(truncate(pretitle, CELLWIDTH), width-1, :center), PRETITLE_STYLE),
                Segment(" " * pad("⋀", width-1, :center), "green bold"),
                Segment(" " * pad(string(SQUARE.bottom.vertical), width-1, :center), "green dim"),
            ]
        )
    end

    return Dendogram(segments, Measure(segments), int(width/2))
end

Dendogram(head, args...; kwargs...) = Dendogram(head, collect(args); kwargs...)

"""
    Dendogram(e::Expr; pretitle=nothing)

Create a Dendogram representation for an `Expr`.
For expressions whose arguments are themselves `Expr` objects, dendograms
are created recursively. In the end all dendograms are linked using `link`
to create a single dendogram object (possibly nested).

`pretitle` is used to create dendograms for `Expr`, pretitle is insrted as
three `Segment`s (a string and an upward arrow) before the head of the dendogram.
"""
function Dendogram(e::Expr; pretitle=nothing)
    length(e.args) == 1 && return Dendogram(e.head, e.args)

    # if there's no more nested expressions, return a dendogram
    !any(isa.(e.args[2:end], Expr)) && return Dendogram(
                        e.head,
                        e.args[2:end];
                        first_arg=e.args[1],
                        pretitle=replace(string(e), ' '=>"")
                )

    # recursively get leaves, creating dendograms for sub expressions
    leaves = map(
        arg -> arg isa Expr ? Dendogram(arg; pretitle=replace(string(arg), ' '=>"")) : Leaf(arg),
        length(e.args) > 2 ? e.args[2:end] : e.args
    )

    # make dendogram by linking individual elements
    title = apply_style("$(e.head): [bold underline default $salmon_light]$(e.args[1])[/bold underline default $salmon_light]", salmon_light)
    return link(leaves...; title=title, shifttitle=true, pretitle=pretitle)
end

# ---------------------------------------------------------------------------- #
#                                     LINK                                     #
# ---------------------------------------------------------------------------- #
"""
    link(dendos...; title="", shifttitle=false, pretitle=nothing)::Dendogram

Link a variable number of `Dendogram` and `Leaf` objects in a new `Dendogram`.

Create a `Dendogram` whose leaves are other leaves and dendograms.
This is done carefully to ensure that all spacings are correct and text
is aligned as much as possible.

Annoyingly, the code is very similar to that of `Dendogram` but different enough
that it can't be refactored into single functions.
"""
function link(dendos...; title="", shifttitle=false, pretitle=nothing)::Dendogram
    # get the widths, spacing between branching points of the dendogram's line
    length(dendos) == 1 && return dendos[1]

    if length(dendos) > 2
        widths = collect(map(
            d -> d[1] == 1 ?
                d[2].measure.w -1  :
                d[2].measure.w + adjust_width(dendos[d[1]-1], d[2]) - 1,

            enumerate(dendos)
        ))[2:end]
    else
        d1, d2 = dendos
        widths = (d1.measure.w - d1.midpoint) + d2.midpoint
    end

    # get elements of linking line
    line = get_row(SQUARE, widths, :top)
    line, midpoint = replace_line_midpoint(line; widths = widths .+ 1)

    if shifttitle
        title = pad(apply_style(title, salmon_light * " dim"), textwidth(line) -5, :center) * " "^5
    else
        title = pad(apply_style(title, salmon_light), textwidth(line), :center)
    end
    space = " "^(dendos[1].midpoint)

    # ensure all elements have the right width
    width = sum(map(d -> d.measure.w, dendos)) - length(space)
    line = pad(line, width, :left)
    title = pad(title, width, :left)

    # create dendogram's segments
    segments::Vector{Segment} = [
        Segment(space * title),
        Segment(space * line, LINES_STYLE,),
        *(dendos...).segments...
    ]

    # add 'pretitle' lines (for expressions only)
    if !isnothing(pretitle)
        pretitle = truncate(pretitle, CELLWIDTH)
        l(txt) = fint(midpoint - textwidth(txt)/2) + length(space)
        r(txt) = cint((width - midpoint - textwidth(txt)/2))
        prepend!(
            segments,
            [
                Segment(" "^l(pretitle) * pretitle * " "^r(pretitle), PRETITLE_STYLE),
                Segment(" "^l("⋀") * "⋀" * " "^r("⋀"), "green bold"),
                Segment(" "^l(SQUARE.bottom.vertical) * SQUARE.bottom.vertical * " "^r(SQUARE.bottom.vertical), "green dim"),

            ]
        )
    end

    return Dendogram(segments, Measure(segments), midpoint+length(space))
end

# ----------------------------------- utils ---------------------------------- #

"""
    adjust_width(x, y)::Int

Width correction factor.

When creating a link between dendograms, the width of the spacing
between branches of the dendogram line needs to be adjusted
depending on the type (`Leaf` or `Dendogram`) of the line.
"""
function adjust_width(x, y)::Int
    _x, _y = typeof(x), typeof(y)
    _x == _y && return 0
    _x == Dendogram && return y.midpoint
    return -x.midpoint
end

"""
    replace_line_midpoint(line::String; widths=nothing)

Replace the mid character of a dendogram tree line with a vertical line for the title.

If the mid character is also the location of a branch (one of the entries in `widths`),
then use the appropriate double branching character.
"""
function replace_line_midpoint(line::String; widths=nothing)
    w1 = prevind(line, int(ncodeunits(line)/2)-1)
    w2 = nextind(line, int(ncodeunits(line)/2)+1)

    char=SQUARE.bottom.vertical
    if !isnothing(widths)
        if textwidth(line[1:w1]) in widths
            char = SQUARE.row.vertical
        end
    end
    line = line[1:w1] * char * line[w2:end]
    return line, length(line[1:w1])
end

end