text:
# Implementations of filt, filtfilt, and fftfilt

#
# filt and filt!
#

## PolynomialRatio
_zerosi(f::PolynomialRatio{:z,T}, x::AbstractArray{S}) where {T,S} =
    zeros(promote_type(T, S), max(-firstindex(f.a), -firstindex(f.b)))

"""
    filt!(out, f, x[, si])

Same as [`filt()`](@ref) but writes the result into the `out`
argument. Output array `out` may not be an alias of `x`, i.e. filtering may
not be done in place.
"""
filt!(out, f::PolynomialRatio{:z}, x::AbstractArray, si=_zerosi(f, x)) =
    filt!(out, coefb(f), coefa(f), x, si)

"""
    filt(f, x[, si])

Apply filter or filter coefficients `f` along the first dimension
of array `x`. If `f` is a filter coefficient object, `si`
is an optional array representing the initial filter state (defaults
to zeros). If `f` is a `PolynomialRatio`, `Biquad`, or
`SecondOrderSections`, filtering is implemented directly. If
`f` is a `ZeroPoleGain` object, it is first converted to a
`SecondOrderSections` object.  If `f` is a Vector, it is
interpreted as an FIR filter, and a naïve or FFT-based algorithm is
selected based on the data and filter length.
"""
filt(f::PolynomialRatio{:z}, x, si=_zerosi(f, x)) = filt(coefb(f), coefa(f), x, si)

## SecondOrderSections
_zerosi(f::SecondOrderSections{:z,T,G}, x::AbstractArray{S}) where {T,G,S} =
    zeros(promote_type(T, G, S), 2, length(f.biquads))

# filt! algorithm (no checking, returns si)
function _filt!(out::AbstractArray, si::AbstractArray{S,N}, f::SecondOrderSections{:z},
                x::AbstractArray, col::Int) where {S,N}
    g = f.g
    biquads = f.biquads
    n = length(biquads)
    @inbounds for i = 1:size(x, 1)
        yi = x[i, col]
        for fi = 1:n
            biquad = biquads[fi]
            xi = yi
            yi = si[1, fi] + biquad.b0*xi
            si[1, fi] = si[2, fi] + biquad.b1*xi - biquad.a1*yi
            si[2, fi] = biquad.b2*xi - biquad.a2*yi
        end
        out[i, col] = yi*g
    end
    si
end

function filt!(out::AbstractArray, f::SecondOrderSections{:z}, x::AbstractArray,
                    si::AbstractArray{S,N}=_zerosi(f, x)) where {S,N}
    biquads = f.biquads
    ncols = Base.trailingsize(x, 2)

    size(x) != size(out) && error("out size must match x")
    (size(si, 1) != 2 || size(si, 2) != length(biquads) || (N > 2 && Base.trailingsize(si, 3) != ncols)) &&
        error("si must be 2 x nbiquads or 2 x nbiquads x nsignals")

    initial_si = si
    g = f.g
    n = length(biquads)
    for col = 1:ncols
        _filt!(out, initial_si[:, :, N > 2 ? col : 1], f, x, col)
    end
    out
end

filt(f::SecondOrderSections{:z,T,G}, x::AbstractArray{S}, si=_zerosi(f, x)) where {T,G,S<:Number} =
    filt!(Array{promote_type(T, G, S)}(undef, size(x)), f, x, si)

## Biquad
_zerosi(f::Biquad{:z,T}, x::AbstractArray{S}) where {T,S} =
    zeros(promote_type(T, S), 2)

# filt! algorithm (no checking, returns si)
function _filt!(out::AbstractArray, si1::Number, si2::Number, f::Biquad{:z},
                x::AbstractArray, col::Int)
    @inbounds for i = 1:size(x, 1)
        xi = x[i, col]
        yi = si1 + f.b0*xi
        si1 = si2 + f.b1*xi - f.a1*yi
        si2 = f.b2*xi - f.a2*yi
        out[i, col] = yi
    end
    (si1, si2)
end

# filt! variant that preserves si
function filt!(out::AbstractArray, f::Biquad{:z}, x::AbstractArray,
                    si::AbstractArray{S,N}=_zerosi(f, x)) where {S,N}
    ncols = Base.trailingsize(x, 2)

    size(x) != size(out) && error("out size must match x")
    (size(si, 1) != 2 || (N > 1 && Base.trailingsize(si, 2) != ncols)) &&
        error("si must have two rows and 1 or nsignals columns")

    initial_si = si
    for col = 1:ncols
        _filt!(out, initial_si[1, N > 1 ? col : 1], initial_si[2, N > 1 ? col : 1], f, x, col)
    end
    out
end

filt(f::Biquad{:z,T}, x::AbstractArray{S}, si=_zerosi(f, x)) where {T,S<:Number} =
    filt!(Array{promote_type(T, S)}(undef, size(x)), f, x, si)

## For arbitrary filters, convert to SecondOrderSections
filt(f::FilterCoefficients{:z}, x) = filt(convert(SecondOrderSections, f), x)
filt!(out, f::FilterCoefficients{:z}, x) = filt!(out, convert(SecondOrderSections, f), x)

"""
    DF2TFilter(coef[, si])

Construct a stateful direct form II transposed filter with
coefficients `coef`. `si` is an optional array representing the
initial filter state (defaults to zeros). If `f` is a
`PolynomialRatio`, `Biquad`, or `SecondOrderSections`,
filtering is implemented directly. If `f` is a `ZeroPoleGain`
object, it is first converted to a `SecondOrderSections` object.
"""
struct DF2TFilter{T<:FilterCoefficients{:z},S<:Array}
    coef::T
    state::S

    function DF2TFilter{Ti,Si}(coef::PolynomialRatio{:z}, state::Vector) where {Ti,Si}
        length(state) == length(coef.a)-1 == length(coef.b)-1 ||
            throw(ArgumentError("length of state vector must match filter order"))
        new{Ti,Si}(coef, state)
    end
    function DF2TFilter{Ti,Si}(coef::SecondOrderSections{:z}, state::Matrix) where {Ti,Si}
        (size(state, 1) == 2 && size(state, 2) == length(coef.biquads)) ||
            throw(ArgumentError("state must be 2 x nbiquads"))
        new{Ti,Si}(coef, state)
    end
    function DF2TFilter{Ti,Si}(coef::Biquad{:z}, state::Vector) where {Ti,Si}
        length(state) == 2 || throw(ArgumentError("length of state must be 2"))
        new{Ti,Si}(coef, state)
    end
end

## PolynomialRatio
DF2TFilter(coef::PolynomialRatio{:z,T},
           state::Vector{S}=zeros(T, max(length(coef.a), length(coef.b))-1)) where {T,S} =
    DF2TFilter{PolynomialRatio{:z,T}, Vector{S}}(coef, state)

function filt!(out::AbstractVector, f::DF2TFilter{<:PolynomialRatio,<:Vector}, x::AbstractVector)
    length(x) != length(out) && throw(ArgumentError("out size must match x"))

    si = f.state
    b = coefb(f.coef)
    a = coefa(f.coef)
    n = length(b)
    if n == 1
        mul!(out, x, b[1])
    else
        @inbounds for i=1:length(x)
            xi = x[i]
            val = si[1] + b[1]*xi
            for j=2:n-1
                si[j-1] = si[j] + b[j]*xi - a[j]*val
            end
            si[n-1] = b[n]*xi - a[n]*val
            out[i] = val
        end
    end
    out
end

## SecondOrderSections
DF2TFilter(coef::SecondOrderSections{:z,T,G},
           state::Matrix{S}=zeros(promote_type(T, G), 2, length(coef.biquads))) where {T,G,S} =
    DF2TFilter{SecondOrderSections{:z,T,G}, Matrix{S}}(coef, state)

function filt!(out::AbstractVector, f::DF2TFilter{<:SecondOrderSections,<:Matrix}, x::AbstractVector)
    length(x) != length(out) && throw(ArgumentError("out size must match x"))
    _filt!(out, f.state, f.coef, x, 1)
    out
end

## Biquad
DF2TFilter(coef::Biquad{:z,T}, state::Vector{S}=zeros(T, 2)) where {T,S} =
    DF2TFilter{Biquad{:z,T}, Vector{S}}(coef, state)
function filt!(out::AbstractVector, f::DF2TFilter{<:Biquad,<:Vector}, x::AbstractVector)
    length(x) != length(out) && throw(ArgumentError("out size must match x"))
    si = f.state
    (si[1], si[2]) =_filt!(out, si[1], si[2], f.coef, x, 1)
    out
end

# Variant that allocates the output
filt(f::DF2TFilter{<:FilterCoefficients{:z},S}, x::AbstractVector) where {S<:Array} =
    filt!(Vector{eltype(S)}(undef, length(x)), f, x)

# Fall back to SecondOrderSections
DF2TFilter(coef::FilterCoefficients{:z}) = DF2TFilter(convert(SecondOrderSections, coef))

#
# filtfilt
#

# Extrapolate the beginning of a signal for use by filtfilt. This
# computes:
#
# [(2 * x[1]) .- x[pad_length+1:-1:2],
#  x,
#  (2 * x[end]) .- x[end-1:-1:end-pad_length]]
#
# in place in output. The istart and n parameters determine the portion
# of the input signal x to extrapolate.
function extrapolate_signal!(out, ostart, sig, istart, n, pad_length)
    length(out) >= n+2*pad_length || error("output is incorrectly sized")
    x = 2*sig[istart]
    for i = 1:pad_length
        out[ostart+i-1] = x - sig[istart+pad_length+1-i]
    end
    copyto!(out, ostart+pad_length, sig, istart, n)
    x = 2*sig[istart+n-1]
    for i = 1:pad_length
        out[ostart+n+pad_length+i-1] = x - sig[istart+n-1-i]
    end
    out
end

# Zero phase digital filtering by processing data in forward and reverse direction
function iir_filtfilt(b::AbstractVector, a::AbstractVector, x::AbstractArray)
    zi = filt_stepstate(b, a)
    pad_length = 3 * (max(length(a), length(b)) - 1)
    t = Base.promote_eltype(b, a, x)
    zitmp = similar(zi, t)
    extrapolated = Vector{t}(undef, size(x, 1)+pad_length*2)
    out = similar(x, t)

    istart = 1
    for i = 1:Base.trailingsize(x, 2)
        extrapolate_signal!(extrapolated, 1, x, istart, size(x, 1), pad_length)
        reverse!(filt!(extrapolated, b, a, extrapolated, mul!(zitmp, zi, extrapolated[1])))
        filt!(extrapolated, b, a, extrapolated, mul!(zitmp, zi, extrapolated[1]))
        for j = 1:size(x, 1)
            @inbounds out[j, i] = extrapolated[end-pad_length+1-j]
        end
        istart += size(x, 1)
    end

    out
end

"""
    filtfilt(coef, x)

Filter `x` in the forward and reverse directions using filter
coefficients `coef`. The initial state of the filter is computed so
that its response to a step function is steady state. Before
filtering, the data is extrapolated at both ends with an
odd-symmetric extension of length
`3*(max(length(b), length(a))-1)`.

Because `filtfilt` applies the given filter twice, the effective
filter order is twice the order of `coef`. The resulting signal has
zero phase distortion.
"""
function filtfilt(b::AbstractVector, x::AbstractArray)
    nb = length(b)
    # Only need as much padding as the order of the filter
    t = Base.promote_eltype(b, x)
    extrapolated = similar(x, t, size(x, 1)+2nb-2, Base.trailingsize(x, 2))

    # Convolve b with its reverse
    newb = filt(b, reverse(b))
    resize!(newb, 2nb-1)
    for i = 1:nb-1
        newb[nb+i] = newb[nb-i]
    end

    # Extrapolate each column
    for i = 1:size(extrapolated, 2)
        extrapolate_signal!(extrapolated, (i-1)*size(extrapolated, 1)+1, x, (i-1)*size(x, 1)+1, size(x, 1), nb-1)
    end

    # Filter
    out = filt(newb, extrapolated)

    # Drop garbage at start
    reshape(out[2nb-1:end, :], size(x))
end

# Choose whether to use FIR or iir_filtfilt depending on
# length of a
function filtfilt(b::AbstractVector, a::AbstractVector, x::AbstractArray)
    if length(a) == 1
        if a[1] != 1
            b /= a[1]
        end
        filtfilt(b, x)
    else
        iir_filtfilt(b, a, x)
    end
end

# Zero phase digital filtering for second order sections
function filtfilt(f::SecondOrderSections{:z,T,G}, x::AbstractArray{S}) where {T,G,S}
    zi = filt_stepstate(f)
    pad_length = 6 * length(f.biquads)
    t = Base.promote_type(T, G, S)
    zitmp = similar(zi, t)
    extrapolated = Vector{t}(undef, size(x, 1)+pad_length*2)
    out = similar(x, t)

    istart = 1
    for i = 1:Base.trailingsize(x, 2)
        extrapolate_signal!(extrapolated, 1, x, istart, size(x, 1), pad_length)
        reverse!(filt!(extrapolated, f, extrapolated, mul!(zitmp, zi, extrapolated[1])))
        filt!(extrapolated, f, extrapolated, mul!(zitmp, zi, extrapolated[1]))
        for j = 1:size(x, 1)
            @inbounds out[j, i] = extrapolated[end-pad_length+1-j]
        end
        istart += size(x, 1)
    end

    out
end

# Support for other filter types
filtfilt(f::FilterCoefficients{:z}, x) = filtfilt(convert(SecondOrderSections, f), x)
filtfilt(f::PolynomialRatio{:z}, x) = filtfilt(coefb(f), coefa(f), x)

## Initial filter state

# Compute an initial state for filt with coefficients (b,a) such that its
# response to a step function is steady state.
function filt_stepstate(b::Union{AbstractVector{T}, T}, a::Union{AbstractVector{T}, T}) where T<:Number
    scale_factor = a[1]
    if scale_factor != 1.0
        a = a ./ scale_factor
        b = b ./ scale_factor
    end

    bs = length(b)
    as = length(a)
    sz = max(bs, as)
    sz > 0 || error("a and b must have at least one element each")
    sz == 1 && return T[]

    # Pad the coefficients with zeros if needed
    bs<sz && (b = copyto!(zeros(eltype(b), sz), b))
    as<sz && (a = copyto!(zeros(eltype(a), sz), a))

    # construct the companion matrix A and vector B:
    A = [-a[2:end] [I; zeros(T, 1, sz-2)]]
    B = b[2:end] - a[2:end] * b[1]
    # Solve si = A*si + B
    # (I - A)*si = B
    scale_factor \ (I - A) \ B
 end

function filt_stepstate(f::SecondOrderSections{:z,T}) where T
    biquads = f.biquads
    si = Matrix{T}(undef, 2, length(biquads))
    y = one(T)
    for i = 1:length(biquads)
        biquad = biquads[i]

        # At steady state, we have:
        #  y = s1 + b0*x
        # s1 = s2 + b1*x - a1*y
        # s2 = b2*x - a2*y
        # where x is the input and y is the output. Solving these
        # equations yields the following.
        si[1, i] = (-(biquad.a1 + biquad.a2)*biquad.b0 + biquad.b1 + biquad.b2)/
                   (1 + biquad.a1 + biquad.a2)*y
        si[2, i] = (biquad.a1*biquad.b2 - biquad.a2*(biquad.b0 + biquad.b1) + biquad.b2)/
                   (1 + biquad.a1 + biquad.a2)*y
        y *= (biquad.b0 + biquad.b1 + biquad.b2)/(1 + biquad.a1 + biquad.a2)
    end
    si
end

"""
    tdfilt(h, x)

Apply filter or filter coefficients `h` along the first dimension
of array `x` using a naïve time-domain algorithm
"""
function tdfilt(h::AbstractVector, x::AbstractArray{T}) where T<:Real
    filt!(Array{T}(undef, size(x)), h, ones(eltype(h), 1), x)
end

"""
    tdfilt!(out, h, x)

Like `tdfilt`, but writes the result into array `out`. Output array `out` may
not be an alias of `x`, i.e. filtering may not be done in place.
"""
function tdfilt!(out::AbstractArray, h::AbstractVector, x::AbstractArray)
    filt!(out, h, ones(eltype(h), 1), x)
end

filt(h::AbstractArray, x::AbstractArray) =
    filt!(Array{eltype(x)}(undef, size(x)), h, x)

#
# fftfilt and filt
#

"""
    fftfilt(h, x)

Apply FIR filter taps `h` along the first dimension of array `x`
using an FFT-based overlap-save algorithm.
"""
function fftfilt(b::AbstractVector{T}, x::AbstractArray{T},
                 nfft::Integer=optimalfftfiltlength(length(b), length(x))) where T<:Real
    _fftfilt!(Array{T}(undef, size(x)), b, x, nfft)
end

"""
    fftfilt!(out, h, x)

Like `fftfilt` but writes result into out array.
"""
function fftfilt!(
    out::AbstractArray{<:Real},
    b::AbstractVector{<:Real},
    x::AbstractArray{<:Real},
    nfft::Integer=optimalfftfiltlength(length(b), length(x))
)
    size(out) == size(x) || throw(ArgumentError("out and x must be the same size"))
    _fftfilt!(out, b, x, nfft)
end

# Like fftfilt! but does not check if out and x are the same size
function _fftfilt!(
    out::AbstractArray{<:Real},
    b::AbstractVector{<:Real},
    x::AbstractArray{T},
    nfft::Integer
) where T<:Real
    nb = length(b)
    nx = size(x, 1)
    normfactor = 1/nfft

    L = min(nx, nfft - (nb - 1))
    tmp1 = Vector{T}(undef, nfft)
    tmp2 = Vector{Complex{T}}(undef, nfft >> 1 + 1)

    p1 = plan_rfft(tmp1)
    p2 = plan_brfft(tmp2, nfft)

    # FFT of filter
    filterft = similar(tmp2)
    tmp1[1:nb] .= b .* normfactor
    tmp1[nb+1:end] .= zero(T)
    mul!(filterft, p1, tmp1)

    # FFT of chunks
    for colstart = 0:nx:length(x)-1
        off = 1
        while off <= nx
            npadbefore = max(0, nb - off)
            xstart = off - nb + npadbefore + 1
            n = min(nfft - npadbefore, nx - xstart + 1)

            tmp1[1:npadbefore] .= zero(T)
            tmp1[npadbefore+n+1:end] .= zero(T)

            copyto!(tmp1, npadbefore+1, x, colstart+xstart, n)
            mul!(tmp2, p1, tmp1)
            broadcast!(*, tmp2, tmp2, filterft)
            mul!(tmp1, p2, tmp2)

            # Copy to output
            copyto!(out, colstart+off, tmp1, nb, min(L, nx - off + 1))

            off += L
        end
    end

    out
end

# Filter x using FIR filter b, heuristically choosing to perform convolution in
# the time domain using tdfilt or in the frequency domain using fftfilt
function filt(b::AbstractVector{T}, x::AbstractArray{T}) where T<:Number
    filt_choose_alg!(Array{T}(undef, size(x)), b, x)
end

# Like filt but mutates output array
function filt!(out::AbstractArray, b::AbstractVector, x::AbstractArray)
    size(out) == size(x) || throw(ArgumentError("out must be the same size as x"))
    filt_choose_alg!(out, b, x)
end

# Perform FIR filtering with either time domain or fft based algorithms, writing
# output to `out` argument. Does not check that x and out are the same size.
function filt_choose_alg!(
    out::AbstractArray{<:Real},
    b::AbstractVector{<:Real},
    x::AbstractArray{<:Real}
)
    nb = length(b)

    if nb > SMALL_FILT_CUTOFF
        nx = size(x, 1)
        nfft = optimalfftfiltlength(nb, nx)
        _fftfilt!(out, b, x, nfft)
    else
        tdfilt!(out, b, x)
    end
end

function filt_choose_alg!(out::AbstractArray, b::AbstractArray, x::AbstractArray)
    tdfilt!(out, b, x)
end


text:
import UncertainData:
    AbstractUncertainIndexValueDataset,
    ConstrainedIndexValueResampling,
    SequentialResampling,
    SequentialInterpolatedResampling,
    constrain,
    resample

function causality(source::AbstractUncertainIndexValueDataset,
    target::AbstractUncertainIndexValueDataset,
    test::CausalityTest,
    resampling::ConstrainedIndexValueResampling{N_VARIABLES, N_DATASETS}) where {N_VARIABLES, N_DATASETS}
    @warn "Performing causality test on uncertain index-value datasets without explicitly providing sequential sampling constraints. Results may be nonsensical, because the values are not sampled according to increasing indices."

    N_DATASETS != 2 ? error("Too many constraints ($N_DATASETS) provided. This `causality` method only takes `source` and `target` as data inputs. Thus, you should provide an `ConstrainedIndexValueResampling` containing two (`index_constraints::CONSTRAINT_TYPES, value_constraints::CONSTRAINT_TYPES)`-tuples. For example, `ConstrainedIndexValueResampling((TruncateStd(0.5), TruncateQuantiles(0.3, 0.7)), (TruncateStd(0.7), TruncateStd(1.1)))` would work.") : nothing

    # Take the (index_constraints, values_constraints)-tuple for the source and target
    # and represent as separate `IndexValueResampling` instances.
    cs_src = ConstrainedIndexValueResampling(resampling[1])
    cs_targ = ConstrainedIndexValueResampling(resampling[2])

    constrained_source = constrain(source, cs_src)
    constrained_target = constrain(target, cs_targ)

    [causality(resample(constrained_source.values), resample(constrained_target.values), test) for i = 1:resampling.n]
end

function causality(source::AbstractUncertainIndexValueDataset,
    target::AbstractUncertainIndexValueDataset,
    test::CausalityTest,
    constrained_resampling::ConstrainedIndexValueResampling{N_VARIABLES, N_DATASETS},
    sequential_resampling::SequentialInterpolatedResampling{S, G}) where {N_VARIABLES, N_DATASETS, S, G}

    N_DATASETS != 2 ? error("Too many constraints ($N_DATASETS) provided. This `causality` method only takes `source` and `target` as data inputs. Thus, you should provide an `ConstrainedIndexValueResampling` containing two (`index_constraints::CONSTRAINT_TYPES, value_constraints::CONSTRAINT_TYPES)`-tuples. For example, `ConstrainedIndexValueResampling((TruncateStd(0.5), TruncateQuantiles(0.3, 0.7)), (TruncateStd(0.7), TruncateStd(1.1)))` would work.") : nothing

    # Take the (index_constraints, values_constraints)-tuple for the source and target
    # and represent as separate `IndexValueResampling` instances.
    cs_src = ConstrainedIndexValueResampling(constrained_resampling[1])
    cs_targ = ConstrainedIndexValueResampling(constrained_resampling[2])

    constrained_source = constrain(source, cs_src)
    constrained_target = constrain(target, cs_targ)

    [causality(resample(constrained_source, sequential_resampling)[2],
               resample(constrained_target, sequential_resampling)[2],
               test) for i = 1:constrained_resampling.n]
end

function causality(source::AbstractUncertainIndexValueDataset,
    target::AbstractUncertainIndexValueDataset,
    test::CausalityTest,
    constrained_resampling::ConstrainedIndexValueResampling{N_VARIABLES, N_DATASETS},
    sequential_resampling::SequentialResampling{S}) where {N_VARIABLES, N_DATASETS, S}

    N_DATASETS != 2 ? error("Too many constraints ($N_DATASETS) provided. This `causality` method only takes `source` and `target` as data inputs. Thus, you should provide an `ConstrainedIndexValueResampling` containing two (`index_constraints::CONSTRAINT_TYPES, value_constraints::CONSTRAINT_TYPES)`-tuples. For example, `ConstrainedIndexValueResampling((TruncateStd(0.5), TruncateQuantiles(0.3, 0.7)), (TruncateStd(0.7), TruncateStd(1.1)))` would work.") : nothing

    # Take the (index_constraints, values_constraints)-tuple for the source and target
    # and represent as separate `IndexValueResampling` instances.
    cs_src = ConstrainedIndexValueResampling(constrained_resampling[1])
    cs_targ = ConstrainedIndexValueResampling(constrained_resampling[2])

    constrained_source = constrain(source, cs_src)
    constrained_target = constrain(target, cs_targ)

    [causality(resample(constrained_source, sequential_resampling)[2],
               resample(constrained_target, sequential_resampling)[2],
               test) for i = 1:constrained_resampling.n]
end

ext: jl


text:
int_rules_6_1_12 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 6.1.12*(e*x)^m*(a+b*sinh(c+d*x^n))^p =#
    @apply_utils Antiderivative(sinh(~(c') + ~(d') * (~x) ^ ~n), ~x) => (1 / 2) * Antiderivative(ℯ ^ (~c + ~d * (~x) ^ ~n), ~x) - (1 / 2) * Antiderivative(ℯ ^ (-(~c) - ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d], ~x) && IGtQ(~n, 1)
    @apply_utils Antiderivative(cosh(~(c') + ~(d') * (~x) ^ ~n), ~x) => (1 / 2) * Antiderivative(ℯ ^ (~c + ~d * (~x) ^ ~n), ~x) + (1 / 2) * Antiderivative(ℯ ^ (-(~c) - ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d], ~x) && IGtQ(~n, 1)
    @apply_utils Antiderivative((~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~a + ~b * sinh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (IGtQ(~n, 1) && IGtQ(~p, 1))
    @apply_utils Antiderivative((~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~a + ~b * cosh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (IGtQ(~n, 1) && IGtQ(~p, 1))
    @apply_utils Antiderivative((~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => -(Subst(Antiderivative((~a + ~b * sinh(~c + ~d * (~x) ^ -(~n))) ^ ~p / (~x) ^ 2, ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (ILtQ(~n, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => -(Subst(Antiderivative((~a + ~b * cosh(~c + ~d * (~x) ^ -(~n))) ^ ~p / (~x) ^ 2, ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (ILtQ(~n, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Module([k = Denominator(~n)], k * Subst(Antiderivative((~x) ^ (k - 1) * (~a + ~b * sinh(~c + ~d * (~x) ^ (k * ~n))) ^ ~p, ~x), ~x, (~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (FractionQ(~n) && IntegerQ(~p))
    @apply_utils Antiderivative((~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Module([k = Denominator(~n)], k * Subst(Antiderivative((~x) ^ (k - 1) * (~a + ~b * cosh(~c + ~d * (~x) ^ (k * ~n))) ^ ~p, ~x), ~x, (~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (FractionQ(~n) && IntegerQ(~p))
    @apply_utils Antiderivative(sinh(~(c') + ~(d') * (~x) ^ ~n), ~x) => (1 / 2) * Antiderivative(ℯ ^ (~c + ~d * (~x) ^ ~n), ~x) - (1 / 2) * Antiderivative(ℯ ^ (-(~c) - ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(cosh(~(c') + ~(d') * (~x) ^ ~n), ~x) => (1 / 2) * Antiderivative(ℯ ^ (~c + ~d * (~x) ^ ~n), ~x) + (1 / 2) * Antiderivative(ℯ ^ (-(~c) - ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~n], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~a + ~b * sinh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~a + ~b * cosh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * sinh(~(c') + ~(d') * (~u) ^ ~n)) ^ ~(p'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~a + ~b * sinh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (IntegerQ(~p) && (LinearQ(~u, ~x) && NeQ(~u, ~x)))
    @apply_utils Antiderivative((~(a') + ~(b') * cosh(~(c') + ~(d') * (~u) ^ ~n)) ^ ~(p'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~a + ~b * cosh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (IntegerQ(~p) && (LinearQ(~u, ~x) && NeQ(~u, ~x)))
    @apply_utils Antiderivative((~(a') + ~(b') * sinh(~(c') + ~(d') * (~u) ^ ~n)) ^ ~p, ~x) => Unintegrable((~a + ~b * sinh(~c + ~d * (~u) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && LinearQ(~u, ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * cosh(~(c') + ~(d') * (~u) ^ ~n)) ^ ~p, ~x) => Unintegrable((~a + ~b * cosh(~c + ~d * (~u) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && LinearQ(~u, ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * sinh(~u)) ^ ~(p'), ~x) => Antiderivative((~a + ~b * sinh(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~(a') + ~(b') * cosh(~u)) ^ ~(p'), ~x) => Antiderivative((~a + ~b * cosh(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative(sinh(~(d') * (~x) ^ ~n) / ~x, ~x) => SinhIntegral(~d * (~x) ^ ~n) / ~n <-- FreeQ([~d, ~n], ~x)
    @apply_utils Antiderivative(cosh(~(d') * (~x) ^ ~n) / ~x, ~x) => CoshIntegral(~d * (~x) ^ ~n) / ~n <-- FreeQ([~d, ~n], ~x)
    @apply_utils Antiderivative(sinh(~c + ~(d') * (~x) ^ ~n) / ~x, ~x) => sinh(~c) * Antiderivative(cosh(~d * (~x) ^ ~n) / ~x, ~x) + cosh(~c) * Antiderivative(sinh(~d * (~x) ^ ~n) / ~x, ~x) <-- FreeQ([~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(cosh(~c + ~(d') * (~x) ^ ~n) / ~x, ~x) => cosh(~c) * Antiderivative(cosh(~d * (~x) ^ ~n) / ~x, ~x) + sinh(~c) * Antiderivative(sinh(~d * (~x) ^ ~n) / ~x, ~x) <-- FreeQ([~c, ~d, ~n], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * (~a + ~b * sinh(~c + ~d * ~x)) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (IntegerQ(Simplify((~m + 1) / ~n)) && (EqQ(~p, 1) || (EqQ(~m, ~n - 1) || IntegerQ(~p) && GtQ(Simplify((~m + 1) / ~n), 0))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * (~a + ~b * cosh(~c + ~d * ~x)) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (IntegerQ(Simplify((~m + 1) / ~n)) && (EqQ(~p, 1) || (EqQ(~m, ~n - 1) || IntegerQ(~p) && GtQ(Simplify((~m + 1) / ~n), 0))))
    @apply_utils Antiderivative((~e * ~x) ^ ~m * (~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * sinh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && IntegerQ(Simplify((~m + 1) / ~n))
    @apply_utils Antiderivative((~e * ~x) ^ ~m * (~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * cosh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && IntegerQ(Simplify((~m + 1) / ~n))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * sinh(~(c') + ~(d') * (~x) ^ ~n), ~x) => ((~e) ^ (~n - 1) * (~e * ~x) ^ ((~m - ~n) + 1) * cosh(~c + ~d * (~x) ^ ~n)) / (~d * ~n) - (((~e) ^ ~n * ((~m - ~n) + 1)) / (~d * ~n)) * Antiderivative((~e * ~x) ^ (~m - ~n) * cosh(~c + ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e], ~x) && (IGtQ(~n, 0) && LtQ(0, ~n, ~m + 1))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * cosh(~(c') + ~(d') * (~x) ^ ~n), ~x) => ((~e) ^ (~n - 1) * (~e * ~x) ^ ((~m - ~n) + 1) * sinh(~c + ~d * (~x) ^ ~n)) / (~d * ~n) - (((~e) ^ ~n * ((~m - ~n) + 1)) / (~d * ~n)) * Antiderivative((~e * ~x) ^ (~m - ~n) * sinh(~c + ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e], ~x) && (IGtQ(~n, 0) && LtQ(0, ~n, ~m + 1))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * sinh(~(c') + ~(d') * (~x) ^ ~n), ~x) => ((~e * ~x) ^ (~m + 1) * sinh(~c + ~d * (~x) ^ ~n)) / (~e * (~m + 1)) - ((~d * ~n) / ((~e) ^ ~n * (~m + 1))) * Antiderivative((~e * ~x) ^ (~m + ~n) * cosh(~c + ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e], ~x) && (IGtQ(~n, 0) && LtQ(~m, -1))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * cosh(~(c') + ~(d') * (~x) ^ ~n), ~x) => ((~e * ~x) ^ (~m + 1) * cosh(~c + ~d * (~x) ^ ~n)) / (~e * (~m + 1)) - ((~d * ~n) / ((~e) ^ ~n * (~m + 1))) * Antiderivative((~e * ~x) ^ (~m + ~n) * sinh(~c + ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e], ~x) && (IGtQ(~n, 0) && LtQ(~m, -1))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * sinh(~(c') + ~(d') * (~x) ^ ~n), ~x) => (1 / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (~c + ~d * (~x) ^ ~n), ~x) - (1 / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (-(~c) - ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e, ~m], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * cosh(~(c') + ~(d') * (~x) ^ ~n), ~x) => (1 / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (~c + ~d * (~x) ^ ~n), ~x) + (1 / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (-(~c) - ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e, ~m], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~x) ^ ~(m') * sinh(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => -(sinh(~a + ~b * (~x) ^ ~n) ^ ~p) / ((~n - 1) * (~x) ^ (~n - 1)) + ((~b * ~n * ~p) / (~n - 1)) * Antiderivative(sinh(~a + ~b * (~x) ^ ~n) ^ (~p - 1) * cosh(~a + ~b * (~x) ^ ~n), ~x) <-- FreeQ([~a, ~b], ~x) && (IntegersQ(~n, ~p) && (EqQ(~m + ~n, 0) && (GtQ(~p, 1) && NeQ(~n, 1))))
    @apply_utils Antiderivative((~x) ^ ~(m') * cosh(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => -(cosh(~a + ~b * (~x) ^ ~n) ^ ~p) / ((~n - 1) * (~x) ^ (~n - 1)) + ((~b * ~n * ~p) / (~n - 1)) * Antiderivative(cosh(~a + ~b * (~x) ^ ~n) ^ (~p - 1) * sinh(~a + ~b * (~x) ^ ~n), ~x) <-- FreeQ([~a, ~b], ~x) && (IntegersQ(~n, ~p) && (EqQ(~m + ~n, 0) && (GtQ(~p, 1) && NeQ(~n, 1))))
    @apply_utils Antiderivative((~x) ^ ~(m') * sinh(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => ((-(~n) * sinh(~a + ~b * (~x) ^ ~n) ^ ~p) / ((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2) + ((~x) ^ ~n * cosh(~a + ~b * (~x) ^ ~n) * sinh(~a + ~b * (~x) ^ ~n) ^ (~p - 1)) / (~b * ~n * ~p)) - ((~p - 1) / ~p) * Antiderivative((~x) ^ ~m * sinh(~a + ~b * (~x) ^ ~n) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (EqQ((~m - 2 * ~n) + 1) && GtQ(~p, 1))
    @apply_utils Antiderivative((~x) ^ ~(m') * cosh(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (-(~n) * cosh(~a + ~b * (~x) ^ ~n) ^ ~p) / ((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2) + ((~x) ^ ~n * sinh(~a + ~b * (~x) ^ ~n) * cosh(~a + ~b * (~x) ^ ~n) ^ (~p - 1)) / (~b * ~n * ~p) + ((~p - 1) / ~p) * Antiderivative((~x) ^ ~m * cosh(~a + ~b * (~x) ^ ~n) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (EqQ((~m - 2 * ~n) + 1) && GtQ(~p, 1))
    @apply_utils Antiderivative((~x) ^ ~(m') * sinh(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (((-(((~m - ~n) + 1)) * (~x) ^ ((~m - 2 * ~n) + 1) * sinh(~a + ~b * (~x) ^ ~n) ^ ~p) / ((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2) + ((~x) ^ ((~m - ~n) + 1) * cosh(~a + ~b * (~x) ^ ~n) * sinh(~a + ~b * (~x) ^ ~n) ^ (~p - 1)) / (~b * ~n * ~p)) - ((~p - 1) / ~p) * Antiderivative((~x) ^ ~m * sinh(~a + ~b * (~x) ^ ~n) ^ (~p - 2), ~x)) + ((((~m - ~n) + 1) * ((~m - 2 * ~n) + 1)) / ((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2)) * Antiderivative((~x) ^ (~m - 2 * ~n) * sinh(~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b], ~x) && (IntegersQ(~m, ~n) && (GtQ(~p, 1) && LtQ(0, 2 * ~n, ~m + 1)))
    @apply_utils Antiderivative((~x) ^ ~(m') * cosh(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (-(((~m - ~n) + 1)) * (~x) ^ ((~m - 2 * ~n) + 1) * cosh(~a + ~b * (~x) ^ ~n) ^ ~p) / ((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2) + ((~x) ^ ((~m - ~n) + 1) * sinh(~a + ~b * (~x) ^ ~n) * cosh(~a + ~b * (~x) ^ ~n) ^ (~p - 1)) / (~b * ~n * ~p) + ((~p - 1) / ~p) * Antiderivative((~x) ^ ~m * cosh(~a + ~b * (~x) ^ ~n) ^ (~p - 2), ~x) + ((((~m - ~n) + 1) * ((~m - 2 * ~n) + 1)) / ((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2)) * Antiderivative((~x) ^ (~m - 2 * ~n) * cosh(~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b], ~x) && (IntegersQ(~m, ~n) && (GtQ(~p, 1) && LtQ(0, 2 * ~n, ~m + 1)))
    @apply_utils Antiderivative((~x) ^ ~(m') * sinh(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (((~x) ^ (~m + 1) * sinh(~a + ~b * (~x) ^ ~n) ^ ~p) / (~m + 1) - (~b * ~n * ~p * (~x) ^ (~m + ~n + 1) * cosh(~a + ~b * (~x) ^ ~n) * sinh(~a + ~b * (~x) ^ ~n) ^ (~p - 1)) / ((~m + 1) * (~m + ~n + 1))) + (((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2) / ((~m + 1) * (~m + ~n + 1))) * Antiderivative((~x) ^ (~m + 2 * ~n) * sinh(~a + ~b * (~x) ^ ~n) ^ ~p, ~x) + (((~b) ^ 2 * (~n) ^ 2 * ~p * (~p - 1)) / ((~m + 1) * (~m + ~n + 1))) * Antiderivative((~x) ^ (~m + 2 * ~n) * sinh(~a + ~b * (~x) ^ ~n) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b], ~x) && (IntegersQ(~m, ~n) && (GtQ(~p, 1) && (LtQ(0, 2 * ~n, 1 - ~m) && NeQ(~m + ~n + 1, 0))))
    @apply_utils Antiderivative((~x) ^ ~(m') * cosh(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => ((((~x) ^ (~m + 1) * cosh(~a + ~b * (~x) ^ ~n) ^ ~p) / (~m + 1) - (~b * ~n * ~p * (~x) ^ (~m + ~n + 1) * sinh(~a + ~b * (~x) ^ ~n) * cosh(~a + ~b * (~x) ^ ~n) ^ (~p - 1)) / ((~m + 1) * (~m + ~n + 1))) + (((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2) / ((~m + 1) * (~m + ~n + 1))) * Antiderivative((~x) ^ (~m + 2 * ~n) * cosh(~a + ~b * (~x) ^ ~n) ^ ~p, ~x)) - (((~b) ^ 2 * (~n) ^ 2 * ~p * (~p - 1)) / ((~m + 1) * (~m + ~n + 1))) * Antiderivative((~x) ^ (~m + 2 * ~n) * cosh(~a + ~b * (~x) ^ ~n) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b], ~x) && (IntegersQ(~m, ~n) && (GtQ(~p, 1) && (LtQ(0, 2 * ~n, 1 - ~m) && NeQ(~m + ~n + 1, 0))))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * (~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => With([k = Denominator(~m)], (k / ~e) * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) * (~a + ~b * sinh(~c + (~d * (~x) ^ (k * ~n)) / (~e) ^ ~n)) ^ ~p, ~x), ~x, (~e * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IntegerQ(~p) && (IGtQ(~n, 0) && FractionQ(~m)))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * (~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => With([k = Denominator(~m)], (k / ~e) * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) * (~a + ~b * cosh(~c + (~d * (~x) ^ (k * ~n)) / (~e) ^ ~n)) ^ ~p, ~x), ~x, (~e * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IntegerQ(~p) && (IGtQ(~n, 0) && FractionQ(~m)))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~e * ~x) ^ ~m, (~a + ~b * sinh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IGtQ(~p, 1) && IGtQ(~n, 0))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~e * ~x) ^ ~m, (~a + ~b * cosh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IGtQ(~p, 1) && IGtQ(~n, 0))
    @apply_utils Antiderivative((~x) ^ ~(m') * sinh(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (((~x) ^ ~n * cosh(~a + ~b * (~x) ^ ~n) * sinh(~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) - (~n * sinh(~a + ~b * (~x) ^ ~n) ^ (~p + 2)) / ((~b) ^ 2 * (~n) ^ 2 * (~p + 1) * (~p + 2))) - ((~p + 2) / (~p + 1)) * Antiderivative((~x) ^ ~m * sinh(~a + ~b * (~x) ^ ~n) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (EqQ((~m - 2 * ~n) + 1, 0) && (LtQ(~p, -1) && NeQ(~p, -2)))
    @apply_utils Antiderivative((~x) ^ ~(m') * cosh(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (-((~x) ^ ~n) * sinh(~a + ~b * (~x) ^ ~n) * cosh(~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) + (~n * cosh(~a + ~b * (~x) ^ ~n) ^ (~p + 2)) / ((~b) ^ 2 * (~n) ^ 2 * (~p + 1) * (~p + 2)) + ((~p + 2) / (~p + 1)) * Antiderivative((~x) ^ ~m * cosh(~a + ~b * (~x) ^ ~n) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (EqQ((~m - 2 * ~n) + 1, 0) && (LtQ(~p, -1) && NeQ(~p, -2)))
    @apply_utils Antiderivative((~x) ^ ~(m') * sinh(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => ((((~x) ^ ((~m - ~n) + 1) * cosh(~a + ~b * (~x) ^ ~n) * sinh(~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) - (((~m - ~n) + 1) * (~x) ^ ((~m - 2 * ~n) + 1) * sinh(~a + ~b * (~x) ^ ~n) ^ (~p + 2)) / ((~b) ^ 2 * (~n) ^ 2 * (~p + 1) * (~p + 2))) - ((~p + 2) / (~p + 1)) * Antiderivative((~x) ^ ~m * sinh(~a + ~b * (~x) ^ ~n) ^ (~p + 2), ~x)) + ((((~m - ~n) + 1) * ((~m - 2 * ~n) + 1)) / ((~b) ^ 2 * (~n) ^ 2 * (~p + 1) * (~p + 2))) * Antiderivative((~x) ^ (~m - 2 * ~n) * sinh(~a + ~b * (~x) ^ ~n) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b], ~x) && (IntegersQ(~m, ~n) && (LtQ(~p, -1) && (NeQ(~p, -2) && LtQ(0, 2 * ~n, ~m + 1))))
    @apply_utils Antiderivative((~x) ^ ~(m') * cosh(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => ((-((~x) ^ ((~m - ~n) + 1)) * sinh(~a + ~b * (~x) ^ ~n) * cosh(~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) + (((~m - ~n) + 1) * (~x) ^ ((~m - 2 * ~n) + 1) * cosh(~a + ~b * (~x) ^ ~n) ^ (~p + 2)) / ((~b) ^ 2 * (~n) ^ 2 * (~p + 1) * (~p + 2)) + ((~p + 2) / (~p + 1)) * Antiderivative((~x) ^ ~m * cosh(~a + ~b * (~x) ^ ~n) ^ (~p + 2), ~x)) - ((((~m - ~n) + 1) * ((~m - 2 * ~n) + 1)) / ((~b) ^ 2 * (~n) ^ 2 * (~p + 1) * (~p + 2))) * Antiderivative((~x) ^ (~m - 2 * ~n) * cosh(~a + ~b * (~x) ^ ~n) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b], ~x) && (IntegersQ(~m, ~n) && (LtQ(~p, -1) && (NeQ(~p, -2) && LtQ(0, 2 * ~n, ~m + 1))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => -(Subst(Antiderivative((~a + ~b * sinh(~c + ~d * (~x) ^ -(~n))) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (IntegerQ(~p) && (ILtQ(~n, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => -(Subst(Antiderivative((~a + ~b * cosh(~c + ~d * (~x) ^ -(~n))) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (IntegerQ(~p) && (ILtQ(~n, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * (~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => With([k = Denominator(~m)], (-k / ~e) * Subst(Antiderivative((~a + ~b * sinh(~c + ~d / ((~e) ^ ~n * (~x) ^ (k * ~n)))) ^ ~p / (~x) ^ (k * (~m + 1) + 1), ~x), ~x, 1 / (~e * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IntegerQ(~p) && (ILtQ(~n, 0) && FractionQ(~m)))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * (~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => With([k = Denominator(~m)], (-k / ~e) * Subst(Antiderivative((~a + ~b * cosh(~c + ~d / ((~e) ^ ~n * (~x) ^ (k * ~n)))) ^ ~p / (~x) ^ (k * (~m + 1) + 1), ~x), ~x, 1 / (~e * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IntegerQ(~p) && (ILtQ(~n, 0) && FractionQ(~m)))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * (~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => -((~e * ~x) ^ ~m) * ((~x) ^ -1) ^ ~m * Subst(Antiderivative((~a + ~b * sinh(~c + ~d * (~x) ^ -(~n))) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IntegerQ(~p) && (ILtQ(~n, 0) && Not(RationalQ(~m))))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * (~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => -((~e * ~x) ^ ~m) * ((~x) ^ -1) ^ ~m * Subst(Antiderivative((~a + ~b * cosh(~c + ~d * (~x) ^ -(~n))) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IntegerQ(~p) && (ILtQ(~n, 0) && Not(RationalQ(~m))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Module([k = Denominator(~n)], k * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) * (~a + ~b * sinh(~c + ~d * (~x) ^ (k * ~n))) ^ ~p, ~x), ~x, (~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (IntegerQ(~p) && FractionQ(~n))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Module([k = Denominator(~n)], k * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) * (~a + ~b * cosh(~c + ~d * (~x) ^ (k * ~n))) ^ ~p, ~x), ~x, (~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (IntegerQ(~p) && FractionQ(~n))
    @apply_utils Antiderivative((~e * ~x) ^ ~m * (~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * sinh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IntegerQ(~p) && FractionQ(~n))
    @apply_utils Antiderivative((~e * ~x) ^ ~m * (~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * cosh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IntegerQ(~p) && FractionQ(~n))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / (~m + 1)) * Subst(Antiderivative((~a + ~b * sinh(~c + ~d * (~x) ^ Simplify(~n / (~m + 1)))) ^ ~p, ~x), ~x, (~x) ^ (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (IntegerQ(~p) && (NeQ(~m, -1) && (IGtQ(Simplify(~n / (~m + 1)), 0) && Not(IntegerQ(~n)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / (~m + 1)) * Subst(Antiderivative((~a + ~b * cosh(~c + ~d * (~x) ^ Simplify(~n / (~m + 1)))) ^ ~p, ~x), ~x, (~x) ^ (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (IntegerQ(~p) && (NeQ(~m, -1) && (IGtQ(Simplify(~n / (~m + 1)), 0) && Not(IntegerQ(~n)))))
    @apply_utils Antiderivative((~e * ~x) ^ ~m * (~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * sinh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && (IntegerQ(~p) && (NeQ(~m, -1) && (IGtQ(Simplify(~n / (~m + 1)), 0) && Not(IntegerQ(~n)))))
    @apply_utils Antiderivative((~e * ~x) ^ ~m * (~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * cosh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && (IntegerQ(~p) && (NeQ(~m, -1) && (IGtQ(Simplify(~n / (~m + 1)), 0) && Not(IntegerQ(~n)))))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * sinh(~(c') + ~(d') * (~x) ^ ~n), ~x) => (1 / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (~c + ~d * (~x) ^ ~n), ~x) - (1 / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (-(~c) - ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * cosh(~(c') + ~(d') * (~x) ^ ~n), ~x) => (1 / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (~c + ~d * (~x) ^ ~n), ~x) + (1 / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (-(~c) - ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * sinh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~e * ~x) ^ ~m, (~a + ~b * sinh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * cosh(~(c') + ~(d') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~e * ~x) ^ ~m, (~a + ~b * cosh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * sinh(~(c') + ~(d') * (~u) ^ ~n)) ^ ~(p'), ~x) => (1 / Coefficient(~u, ~x, 1) ^ (~m + 1)) * Subst(Antiderivative((~x - Coefficient(~u, ~x, 0)) ^ ~m * (~a + ~b * sinh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (LinearQ(~u, ~x) && (NeQ(~u, ~x) && IntegerQ(~m)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * cosh(~(c') + ~(d') * (~u) ^ ~n)) ^ ~(p'), ~x) => (1 / Coefficient(~u, ~x, 1) ^ (~m + 1)) * Subst(Antiderivative((~x - Coefficient(~u, ~x, 0)) ^ ~m * (~a + ~b * cosh(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (LinearQ(~u, ~x) && (NeQ(~u, ~x) && IntegerQ(~m)))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * sinh(~(c') + ~(d') * (~u) ^ ~n)) ^ ~(p'), ~x) => Unintegrable((~e * ~x) ^ ~m * (~a + ~b * sinh(~c + ~d * (~u) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && LinearQ(~u, ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * cosh(~(c') + ~(d') * (~u) ^ ~n)) ^ ~(p'), ~x) => Unintegrable((~e * ~x) ^ ~m * (~a + ~b * cosh(~c + ~d * (~u) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && LinearQ(~u, ~x)
    @apply_utils Antiderivative((~e * ~x) ^ ~(m') * (~(a') + ~(b') * sinh(~u)) ^ ~(p'), ~x) => Antiderivative((~e * ~x) ^ ~m * (~a + ~b * sinh(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~e, ~m, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~e * ~x) ^ ~(m') * (~(a') + ~(b') * cosh(~u)) ^ ~(p'), ~x) => Antiderivative((~e * ~x) ^ ~m * (~a + ~b * cosh(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~e, ~m, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~x) ^ ~(m') * sinh(~(a') + ~(b') * (~x) ^ ~n) ^ ~(p') * cosh(~(a') + ~(b') * (~x) ^ ~(n')), ~x) => sinh(~a + ~b * (~x) ^ ~n) ^ (~p + 1) / (~b * ~n * (~p + 1)) <-- FreeQ([~a, ~b, ~m, ~n, ~p], ~x) && (EqQ(~m, ~n - 1) && NeQ(~p, -1))
    @apply_utils Antiderivative((~x) ^ ~(m') * cosh(~(a') + ~(b') * (~x) ^ ~n) ^ ~(p') * sinh(~(a') + ~(b') * (~x) ^ ~(n')), ~x) => cosh(~a + ~b * (~x) ^ ~n) ^ (~p + 1) / (~b * ~n * (~p + 1)) <-- FreeQ([~a, ~b, ~m, ~n, ~p], ~x) && (EqQ(~m, ~n - 1) && NeQ(~p, -1))
    @apply_utils Antiderivative((~x) ^ ~(m') * sinh(~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(p') * cosh(~(a') + ~(b') * (~x) ^ ~(n')), ~x) => ((~x) ^ ((~m - ~n) + 1) * sinh(~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) - (((~m - ~n) + 1) / (~b * ~n * (~p + 1))) * Antiderivative((~x) ^ (~m - ~n) * sinh(~a + ~b * (~x) ^ ~n) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (LtQ(0, ~n, ~m + 1) && NeQ(~p, -1))
    @apply_utils Antiderivative((~x) ^ ~(m') * cosh(~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(p') * sinh(~(a') + ~(b') * (~x) ^ ~(n')), ~x) => ((~x) ^ ((~m - ~n) + 1) * cosh(~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) - (((~m - ~n) + 1) / (~b * ~n * (~p + 1))) * Antiderivative((~x) ^ (~m - ~n) * cosh(~a + ~b * (~x) ^ ~n) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (LtQ(0, ~n, ~m + 1) && NeQ(~p, -1))
end

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


text:
abstract type AbstractAtom end

mutable struct Atom <: AbstractAtom
    symbol::String
    position::Array{Float64, 1}

    Atom() = new(" ", [0.0, 0.0, 0.0])
    function Atom(sym::String, pos::Array{T, 1}) where T<:Real
        if length(pos) != 3
            throw(BoundsError(pos, length(pos)))
        else
            new(sym, Array{Float64, 1}(pos))
        end
    end
end

function (atom::AbstractAtom)()
    return atom.symbol, atom.position
end

function Base.:(==)(atom1::AbstractAtom, atom2::AbstractAtom)
    atom1.symbol == atom2.symbol || return false
    atom1.position == atom2.position || return false
    return true
end

function Base.show(io::IO, atom::AbstractAtom)
    print(io, atom.symbol, "\t", atom.position)
end


text:
using LinearAlgebra
using Printf

function _polynomial(a, x)
    accumulator = a[end] * one(x)
    for i in length(a)-1:-1:1
        accumulator = accumulator * x + a[i]
    end
    accumulator
end

# definition of polynom
struct Polynom{C}
    coefficients::C
    Polynom(coefficients::CC) where CC = coefficients[end] == 0 ? throw(ArgumentError("Coefficient of the highest exponent cannot be zero.")) : new{CC}(coefficients)
end

# based on https://github.com/JuliaMath/Polynomials.jl
function from_roots(roots::AbstractVector{T}; aₙ = one(T)) where {T}
    n = length(roots)
    c = zeros(T, n+1)
    c[1] = one(T)
    for j = 1:n
        for i = j:-1:1
            c[i+1] = c[i+1]-roots[j]*c[i]
        end
    end
    return Polynom(aₙ.*reverse(c))
end

(p::Polynom)(x) = _polynomial(p.coefficients, x)
degree(p::Polynom) = length(p.coefficients) - 1

function _derivativeof(p::Polynom)
    n = degree(p)
    n > 1 ? Polynom([(i - 1)*p.coefficients[i] for i in 2:n+1]) : error("Low degree of a polynomial.")
end
LinearAlgebra.adjoint(p::Polynom) = _derivativeof(p)

function Base.show(io::IO, p::Polynom)
    n = degree(p)
    a = reverse(p.coefficients)
    for (i, c) in enumerate(a[1:end-1])
        if (c != 0)
            c < 0 && print(io, " - ")
            c > 0 && i > 1 && print(io, " + ")
            print(io, "$(abs(c))x^$(n - i + 1)")
        end
    end
    c = a[end]
    c > 0 && print(io, " + $(c)")
    c < 0 && print(io, " - $(abs(c))")
end

# default optimization parameters
atol = 1e-12
maxiter = 100
stepsize = 0.95

# definition of optimization methods
abstract type RootFindingMethod end
struct Newton <: RootFindingMethod end
struct Secant <: RootFindingMethod end
struct Bisection <: RootFindingMethod end

init!(::Bisection, p, a, b) = sign(p(a)) != sign(p(b)) ? (a, b) : throw(ArgumentError("Signs at both ends are the same."))
init!(::RootFindingMethod, p, a, b) = (a, b)

function step!(::Newton, poly::Polynom, xᵢ, step_size)
    _, x̃ = xᵢ
    dp = p'
    x̃, x̃ - step_size*p(x̃)/dp(x̃)
end

function step!(::Secant, poly::Polynom, xᵢ, step_size)
    x, x̃ = xᵢ
    dpx = (p(x) - p(x̃))/(x - x̃)
    x̃, x̃ - stepsize*p(x̃)/dpx
end

function step!(::Bisection, poly::Polynom, xᵢ, step_size)
    x, x̃ = xᵢ
    midpoint = (x + x̃)/2
    if sign(p(midpoint)) == sign(p(x̃))
        x̃ = midpoint
    else
        x = midpoint
    end
    x, x̃
end

function find_root(p::Polynom, rfm=Newton, a=-5.0, b=5.0, max_iter=100, step_size=0.95, tol=1e-12)
    x, x̃ = init!(rfm, p, a, b)
    for i in 1:maxiter
        x, x̃ = step!(rfm, p, (x, x̃), step_size)
        val = p(x̃)
        @printf "x = %.5f | x̃ = %.5f | p(x̃) = %g\n" x x̃ val
        abs(val) < atol && return x̃
    end
    println("Method did not converge in $(max_iter) iterations to a root within $(tol) tolerance.")
    return x̃
end

# test code
poly = Polynom(rand(4))
p = from_roots([-3, -2, -1, 0, 1, 2, 3])
dp = p'
p(3.0), dp(3.0)

x₀ = find_root(p, Bisection(), -5.0, 5.0)
x₀ = find_root(p, Newton(), -5.0, 5.0)
x₀ = find_root(p, Secant(), -5.0, 5.0)


text:
# ===============================
# Test default values and aliases
# ===============================
o = Reachability.Options(:key1 => 1, :key3 => 1, :alias4 => "a")
dict1 = o.dict
dict2 = Dict{Symbol, Any}()

# do not modify dict1 for existing key
Reachability.check_aliases_and_add_default_value!(dict1, dict2, [:key1], 2)
@test haskey(dict1, :key1) && dict1[:key1] == 1
@test haskey(dict2, :key1) && dict2[:key1] == 1

# use default value
Reachability.check_aliases_and_add_default_value!(dict1, dict2, [:key2], 2, true)
@test haskey(dict1, :key2) && dict1[:key2] == 2
@test haskey(dict2, :key2) && dict2[:key2] == 2

# an unused alias has no effect
Reachability.check_aliases_and_add_default_value!(dict1, dict2, [:key3, :alias3], 2)
@test dict1[:key3] == 1
@test !haskey(dict1, :alias3)
@test haskey(dict2, :key3) && dict2[:key3] == 1
@test !haskey(dict2, :alias3)

# check alias detection
Reachability.check_aliases_and_add_default_value!(dict1, dict2, [:key4, :alias4], "b")
@test !haskey(dict1, :key4)
@test haskey(dict1, :alias4) && dict1[:alias4] == "a"
@test haskey(dict2, :key4) && dict2[:key4] == "a"
@test !haskey(dict2, :alias4)

# check alias detection together with default values
Reachability.check_aliases_and_add_default_value!(dict1, dict2, [:key5, :alias5], "b", true)
@test haskey(dict1, :key5) && dict1[:key5] == "b"
@test !haskey(dict1, :alias5)
@test haskey(dict2, :key5) && dict2[:key5] == "b"
@test !haskey(dict2, :alias5)


text:
using StatsPlots, SyntheticDatasets

blobs = SyntheticDatasets.make_blobs(   n_samples = 1000,
                                        n_features = 2,
                                        centers = [-1 1; -0.5 0.5],
                                        cluster_std = 0.25,
                                        center_box = (-2.0, 2.0),
                                        shuffle = true,
                                        random_state = nothing);

@df blobs scatter(:feature_1, :feature_2, group = :label, title = "Blobs")

gauss = SyntheticDatasets.make_gaussian_quantiles(  mean = [10,1],
                                                    cov = 2.0,
                                                    n_samples = 1000,
                                                    n_features = 2,
                                                    n_classes = 3,
                                                    shuffle = true,
                                                    random_state = 2);

@df gauss scatter(:feature_1, :feature_2, group = :label, title = "Gaussian Quantiles")

spirals = SyntheticDatasets.make_twospirals(n_samples = 2000,
                                            start_degrees = 90,
                                            total_degrees = 570,
                                            noise =0.1);

@df spirals scatter(:feature_1, :feature_2, group = :label, title = "Two Spirals")

kernel = SyntheticDatasets.make_halfkernel( n_samples = 1000,
                                            minx = -20,
                                            r1 = 20,
                                            r2 = 35,
                                            noise = 3.0,
                                            ratio = 0.6);

@df kernel scatter(:feature_1, :feature_2, group = :label, title = "Half Kernel")


text:
@testset "94.binary-tree-inorder-traversal.jl" begin
    @test inorder_traversal(TreeNode{Int}([1])) == [1]
    @test inorder_traversal(nothing) == Int[]
    @test inorder_traversal(TreeNode{Int}([1, nothing, 2, 3])) == [1, 3, 2]
    @test inorder_traversal(TreeNode{Int}([1, 2, 3, 4, 5, 6])) == [4, 2, 5, 1, 6, 3]
    @test inorder_traversal(TreeNode{Int}([3, 9, 4, nothing, nothing, 5, 7])) ==
        [9, 3, 5, 4, 7]
end


text:
using DiscreteAxis
using StaticArrays#
using Rotations
using PDEUtils

abstract type Field{T, N} <: FieldVector{N, Array{T, N}} end
abstract type CompositeField{T, J, N} <: FieldVector{J, Field{T,N}} end

const ε₀ = 8.85418782*10^-12
const μ₀ = 4*pi*10^-7
const c₀ = 1/√(ε₀*μ
const η₀ = √(μ₀/ε₀)
const Sc = 1/√2

struct VecField3{T<:Number, N} <: FieldVector{3, Array{T, N}}
    x::Array{T,N}
    y::Array{T,N}
    z::Array{T,N}
    VecField3(x::Array{T,N}, y::Array{T,N}, z::Array{T,N}) where {T,N} = new{T,N}(x,y,z)
    VecField3(s::Space2D{T}) where {T<:Number} =  new{T,2}(zeros(T, s.x.N, s.y.N),
                                                         zeros(T, s.x.N, s.y.N),
                                                         zeros(T, s.x.N, s.y.N))

    function VecField3(s::Space2D{T}, nLayer::Int, dim::Int) where {T<:Number}
        if dim == 1
            return new{T,2}(zeros(T, nLayer, s.y.N),
                            zeros(T, nLayer, s.y.N),
                            zeros(T, nLayer, s.y.N))
        elseif dim == 2
            return new{T,2}(zeros(T, s.x.N, nLayer),
                          zeros(T, s.x.N, nLayer),
                          zeros(T, s.x.N, nLayer))
        else
            throw("Only Integer dimension of 1-3 is supported! Check the dim argument.\n")
        end
    end
end

Base.size(field::Field{T,N}) where {T,N} = (size(field[1]), N)

colons(n::Int) = fill(:, n)

function flatten(F::Field{T,N}) where {T,N}
    flatfield = zeros(T, size(F)...)
    for (i,component) in enumerate(Field)
        flatfield[i, colons(N)...] .= component
    end
    return flatfield
end
Base.eltype(A::VecField3) = eltype(A[1])

struct PMLAux{T<:Number,N}
    max::Vector{VecField3{T}}
    min::Vector{VecField3{T}}
    function PMLAux(s::Space2D{T}, nLayer) where {T<:Number}
        max = [VecField3(s, nLayer, 1), VecField3(s, nLayer, 2)]
        min = [VecField3(s, nLayer, 1), VecField3(s, nLayer, 2)]
        new{T,2}(max,min)
    end
    PMLAux(max::Vector{VecField3{T,N}}, min::Vector{VecField3{T,N}}) where {T,N} = new{T,N}(max,min)
end

struct PML{T<:Number,M}
    E::PMLAux{T,M}
    H::PMLAux{T,M}
    N::Int64
    function PML(s::Space2D{T}, nLayer) where {T<:Number}
        E, H = (PMLAux(s, nLayer+1), PMLAux(s,nLayer+1))
        new{T,2}(E,H, nLayer+1)
    end
    PML(E::PMLAux{T,M}, H::PMLAux{T,M},N::Int64) where {T,M} = new{T,M}(E,H,N)
end

struct EMField{T<:Number,N} <:FieldVector{2, VecField3{T}}
    E::VecField3{T,N}
    H::VecField3{T,N}
    EMField(E::VecField3{T,N}, H::VecField3{T,N}) where {T,N} = new{T,N}(E,H)
    EMField(space::Space2D{T}) where {T<:Number}= new{T,2}(VecField3(space), VecField3(space))
    EMField(F::Tuple{Tuple{Array{T,N},Array{T,N},Array{T,N}},Tuple{Array{T,N},Array{T,N},Array{T,N}}}) where {T,N} = new{T,N}(VecField3(F[1]...), VecField3(F[2]...))
end

    function split_array(A::AbstractArray{T}) where {T <: AbstractVector}
        out = [zeros(eltype(first(A)), size(A)...) for i in 1:length(first(A))]
        for I in CartesianIndices(A)
            for (target, n) in enumerate(out)
                target[I] = A[I][n]
            end
        end
        return Tuple(out)
    end

    function vec_array(A::NTuple{N, T}) where {T <: AbstractArray, N}
        @assert all([size(A[i]) == size(A[1]) for i in 2:N])
        @assert all([eltype(A[i]) == eltype(A[1]) for i in 2:N])
        out = Array{SVector{N, T}, ndims(A[1])}(undef, size(A[1])...)
        for I in CartesianIndices{A}
            out[I] = SVector([A[i][I] for i in 1:N])
        end
        return out
    end

    power_density(F::EMField) = @. (F.E.x^2 + F.E.y^2 + F.E.z^2)/(2*η₀)
    function pontying(F::EMField{T}) where T
        out = Array{SVector{3,T}}(undef, size(F)...)
        for I in CartesianIndices(F.E)
            out[I] = SVector([F.E[2][I]*F.H[3][I] - F.E[3][I]*F.H[2][I],
                 F.E[3][I]*F.H[1][I] - F.E[1][I]*F.H[3][I],
                 F.E[1][I]*F.H[2][I] - F.E[2][I]*F.H[1][I]])

        end
        return out
    end
function deep_times(F::VecArray, k::Number)
    Fx = F[1].*k
    Fy = F[2].*k
    Fz = F[3].*k
    return (Fx,Fy,Fz)
end

function Base.getproperty(F::EHTuple{T}, s::Symbol)
    if symbol == :E
        return F[1]
    elseif symbol == :H
        return F[2]
    elseif symbol == :D
        return deep_times(F[1], ε₀)
    elseif symbol == :W
        return power_density(F)
    elseif symbol == :S
        return pontying(F)
    elseif
        throw(ArgumentError("Symbol $s is not a property of F"))
    end
end

function Base.getproperty(F::VecArray, s::Symbol)
    if symbol == :x
        return F[1]
    elseif symbol == :y
        return F[2]
    elseif symbol == :z
        return F[3]
    elseif symbol == :

    else
        throw(ArgumentError("Symbol $s is not a property of F"))
    end
end

Base.getindex(F::EHTuple{T,2}, i::Int, j::Int) where T = (SVector(F.E.x[i,j] F.E.y[i,j] F.E.z[i,j]), SVector(F.H.x[i,j] F.H.y[i,j] F.H.z[i,j]))

Base.getindex(F::EHTuple{T,3}, i::Int, j::Int, k::Int) where T = (SVector(F.E.x[i,j,k] F.E.y[i,j,k] F.E.z[i,j,k]), SVector(F.H.x[i,j,k] F.H.y[i,j,k] F.H.z[i,j,k]))
Base.getindex(F::EHTuple, I::CartesianIndex) = (SVector(F.E.x[I] F.E.y[I] F.E.z[I]), SVector(F.H.x[I] F.H.y[I] F.H.z[I]))
Base.getindex(F::EHTuple, I::CartesianIndex, i::Int) = SVector(F.E.x[I] F.E.y[I] F.E.z[I])
power_density(x::AbstractVector) = sum(x.^2)/2*η₀

Base.CartesianIndices(F::EMField) = CartesianIndices(F.H.x)
Base.size(F::EMField) = size(F.E.x)
με_default(r::AbstractVector) = one(Float64)
σ_default(r::AbstractVector) = zero(Float64)

struct Medium{T<:Number, N} <: Field{T,4}
    ε::Array{T, N}
    μ::Array{T, N}
    σ::Array{T, N}
    σm::Array{T, N}
    function Medium(s::Space2D{T}; θ=0.0, rotcentre = [0.0,0.0], ε = με_default, μ = με_default, σ = σ_default, σm = σ_default) where {T<:Number}
        R = [cos(θ) -sin(θ); sin(θ) cos(θ)] #generate rotation matrix
        X, Y= (s.x.pts .- rotcentre[1], s.y.pts .- rotcentre[2])
        r = [SVector{2,T}(R*[x,y]) for x in X, y in Y]
        eps = ε₀.*ε.(r)
        mu = μ₀.*μ.(r)
        sig = σ.(r)
        sigm = σm.(r)
        new{eltype(eps),2}(eps, mu, sig, sigm)
    end
    function Medium(s::Space3D{T}; azimuth = 0.0, elevation = 0.0, rotcentre = [0.0,0.0,0.0], ε = με_default, μ = με_default, σ = σ_default, σm = σ_default) where {T<:Number}
        R = RotZXY(azimuth, elevation, 0.0) #generate rotation matrix
        X, Y, Z = (s.x.pts .- rotcentre[1], s.y.pts .- rotcentre[2], s.z.pts .- rotcentre[3])
        r′ = [SVector{3,T}(R*[x,y,z] for x in X, y in Y, z in Z]
        eps = ε₀.*ε.(r′)
        mu = μ₀.*μ.(r′)
        sig = σ.(r′)
        sigm = σm.(r′)
        new{eltype(eps),3}(eps, mu, sig, sigm)
    end
end

struct Coefficients{T,N} <: FieldVector{11, UnionAll}
    E::NTuple{6, Array{T, N}}
    H::NTuple{6, Array{T, N}}
    s::Space2D{T}
    t::LinearAxis{T}
    δ::NTuple{2, Matrix{T}}
    bu::NTuple{2,Vector{T}}
    bl::NTuple{2,Vector{T}}
    cu::NTuple{2,Vector{T}}
    cl::NTuple{2,Vector{T}}
    Km::NTuple{2,Vector{T}}
    Ke::NTuple{2,Vector{T}}
end

    # This function assumes that Δx Δy and Δz are the same - Not too keen on the idea of implementing a non uniform grid right now

function Coefficients(s::Space2D, t, m::Medium{T,2}, f₀, d::Int) where T
    δtmp = (δ⁺(d+2), δ⁻(d+2))
    δ = map(D -> D./s.x.Δ, δtmp)
    # Calculate Coefficients for the layer
    layer = (0:d)

    kmax = 15
    M = 4
    Ma = 1
    #https://www.degruyter.com/downloadpdf/j/jee.2017.68.issue-1/jee-2017-0006/jee-2017-0006.pdf !!!!!!
    σemax = (M+1)*0.8/(η₀*s.x.Δ)
    σmmax = (M+1)*0.8*η₀/s.x.Δ
    aemax = 2π*f₀*ε₀/10
    ammax = 2π*f₀*μ₀/10
    for (i, z) in enumerate(PML)
        kmu[i] = 1+(kmax-1)*((z+0.5)/d)^M
        kml[i] = 1+(kmax-1)*(z/d)^M
        keu[i] = 1+(kmax-1)*(z/d)^M
        kel[i] = 1+(kmax-1)*((z+0.5)/d)^M
        σmu[i] = σmmax*((z+0.5)/d)^M
        σml[i] = σmmax*(z/d)^M
        σeu[i] = σemax*(z/d)^M
        σel[i] = σemax*((z+0.5)/d)^M
        αmu[i] = ammax*(((d - z+0.5)/d)^Ma)
        αml[i] = ammax*(((d - z)/d)^Ma)
        αeu[i] = aemax*(((d - z)/d)^Ma)
        αel[i] = aemax*(((d - z+0.5)/d)^Ma)
    end

    # Calculate auxillary Coefficients for the layer
    beu = exp.(-(σeu./keu .+ αeu).*t.Δ/ε₀)
    ceu = (beu.-1.0).*σeu ./ (σeu.*keu + keu.^2 .*αeu)
    bel = exp.(-(σel./kel .+ αel).*t.Δ/ε₀)
    cel = (bel.-1.0).*σel ./ (σel.*kel + kel.^2 .*αel)

    bmu = exp.(-(σmu./kmu .+ αmu).*t.Δ/μ₀)
    cmu = (bmu.-1.0).*σmu ./ (σmu.*kmu + kmu.^2 .*αmu)
    bml = exp.(-(σml./kml .+ αml).*t.Δ/μ₀)

    cml = (bml.-1.0).*σml ./ (σml.*kml + kml.^2 .*αml)

    tmp = 1.0 .+ m.σm.*t.Δ ./(2.0.*m.μ)
    #Δt ⪕ √(μ₀ε₀)
    Chxh = (2.0 .- tmp)./tmp #source chapter 9 page 8 https://www.eecs.wsu.edu/~schneidj/ufdtd/, with modifications to support inhomogeneous media
    Chyh = Chzh = Chxh # Start here if you want to implement anisotropy

    Chxe = t.Δ ./ (tmp .*m.μ .* s.x.Δ )
    Chye = Chze = Chxe

    tmp = 1.0 .+ m.σ ./(2.0.*m.ε)
    Cexe = (2.0.-tmp)./tmp
    Ceye = Ceze = Cexe # these are equal in the anisotropic uniform grid case, this kind of assignment also only allocates 1 matrix of memory, meaning the constants are passed as pointers to the same place in memory

    Cexh = t.Δ ./ (tmp .* m.ε.*s.x.Δ)
    Ceyh = Cezh = Cexh

    # put all the coeffs in to a nice neat wrapper
    H = (Chxh, Chyh, Chzh, Chxe, Chye, Chze)
    E = (Cexe, Ceye, Ceze, Cexh, Ceyh, Cezh)

    Ke = Vector[ones(s.x.N),ones(s.y.N)]
    Km = Vector[ones(s.x.N),ones(s.y.N)]
    for i in 1:2
        Km[i][1:(d+1)] .= 1 ./reverse(kml)
        Km[i][(end-d):end] .= 1 ./kmu
        Ke[i][1:d+1] .= 1 ./reverse(kel)
        Ke[i][(end-d):end] .= 1 ./keu
    end
    rbel = reverse(bel)
    rcel = reverse(cel)
    rbml = reverse(bml)
    rcml = reverse(cml)
    Coefficients{T,2}(E, H,s,t, δ, (beu,bmu), (rbel,rbml), (ceu,cmu), (rcel,  rcml), Tuple(Km), Tuple(Ke))
end

    function Coefficients(s, t, m::Medium{T,3}, f₀, d::Int) where {T}
        δtmp =(δ⁺(d+2), δ⁻(d+2))
        δ = map(D -> D./s.x.Δ, δtmp)
        # Calculate Coefficients for the PML
        PML = (0:d)

        kmax = 15
        M = 4
        Ma = 1
        #https://www.degruyter.com/downloadpdf/j/jee.2017.68.issue-1/jee-2017-0006/jee-2017-0006.pdf !!!!!!
        σemax = (M+1)*0.8/(η₀*s.x.Δ)
        σmmax = (M+1)*0.8*η₀/s.x.Δ
        aemax = 2π*f₀*ε₀/10
        ammax = 2π*f₀*μ₀/10

        for (i, z) in enumerate(PML)
            kmu[i] = 1+(kmax-1)*((z+0.5)/d)^M
            kml[i] = 1+(kmax-1)*(z/d)^M
            keu[i] = 1+(kmax-1)*(z/d)^M
            kel[i] = 1+(kmax-1)*((z+0.5)/d)^M
            σmu[i] = σmmax*((z+0.5)/d)^M
            σml[i] = σmmax*(z/d)^M
            σeu[i] = σemax*(z/d)^M
            σel[i] = σemax*((z+0.5)/d)^M
            αmu[i] = ammax*(((d - z+0.5)/d)^Ma)
            αml[i] = ammax*(((d - z)/d)^Ma)
            αeu[i] = aemax*(((d - z)/d)^Ma)
            αel[i] = aemax*(((d - z+0.5)/d)^Ma)
        end

        # Calculate auxillary Coefficients for the PML
        beu = exp.(-(σeu./keu .+ αeu).*t.Δ/ε₀)
        ceu = (beu.-1.0).*σeu ./ (σeu.*keu + keu.^2 .*αeu)
        bel = exp.(-(σel./kel .+ αel).*t.Δ/ε₀)
        cel = (bel.-1.0).*σel ./ (σel.*kel + kel.^2 .*αel)

        bmu = exp.(-(σmu./kmu .+ αmu).*t.Δ/μ₀)
        cmu = (bmu.-1.0).*σmu ./ (σmu.*kmu + kmu.^2 .*αmu)
        bml = exp.(-(σml./kml .+ αml).*t.Δ/μ₀)

        cml = (bml.-1.0).*σml ./ (σml.*kml + kml.^2 .*αml)

        tmp = 1.0 .+ m.σm.*t.Δ ./(2.0.*m.μ)
        #Δt ⪕ √(μ₀ε₀)
        Chxh = (2.0 .- tmp)./tmp #source chapter 9 page 8 https://www.eecs.wsu.edu/~schneidj/ufdtd/, with modifications to support inhomogeneous media
        Chyh = Chzh = Chxh # Start here if you want to implement anisotropy

        Chxe = t.Δ ./ (tmp .*m.μ .* s.x.Δ )
        Chye = Chze = Chxe

        tmp = 1.0 .+ m.σ ./(2.0.*m.ε)
        Cexe = (2.0.-tmp)./tmp
        Ceye = Ceze = Cexe # these are equal in the anisotropic uniform grid case, this kind of assignment also only allocates 1 matrix of memory, meaning the constants are passed as pointers to the same place in memory

        Cexh = t.Δ ./ (tmp .* m.ε.*s.x.Δ)
        Ceyh = Cezh = Cexh

        # put all the coeffs in to a nice neat wrapper
        H = (Chxh, Chyh, Chzh, Chxe, Chye, Chze)
        E = (Cexe, Ceye, Ceze, Cexh, Ceyh, Cezh)

        Ke = Vector[ones(s.x.N),ones(s.y.N),ones(s.z.N)]
        Km = Vector[ones(s.x.N),ones(s.y.N),ones(s.z.N)]
        for i in 1:3
            Km[i][1:(d+1)] .= 1 ./reverse(kml)
            Km[i][(end-d):end] .= 1 ./kmu
            Ke[i][1:d+1] .= 1 ./reverse(kel)
            Ke[i][(end-d):end] .= 1 ./keu
        end
        rbel = reverse(bel)
        rcel = reverse(cel)
        rbml = reverse(bml)
        rcml = reverse(cml)
        Coefficients{T}(E, H,s,t, δ, (beu,bmu), (rbel,rbml), (ceu,cmu), (rcel,  rcml), Tuple(Km), Tuple(Ke))
    end


text:
using SymbolicTensors2
using Test

@testset "SymbolicTensors2.jl" begin
    # Write your tests here.
end


text:
export dic

#################### Posterior Statistics ####################

"""
    dic(chain::Chains, logpdf::Function) -> (DIC, pD)

Compute the deviance information criterion (DIC).
(Smaller is better)

Note: DIC assumes that the posterior distribution is approx. multivariate Gaussian and tends to select overfitted models.

## Returns:
* `DIC`: The calculated deviance information criterion
* `pD`: The effective number of parameters

## Usage:

```
chn ... # sampling results
lpfun = function f(chain::Chains) # function to compute the logpdf values
    niter, nparams, nchains = size(chain)
    lp = zeros(niter + nchains) # resulting logpdf values
    for i = 1:nparams
        lp += map(p -> logpdf( ... , x), Array(chain[:,i,:]))
    end
    return lp
end

DIC, pD = dic(chn, lpfun)
```

"""
function dic(chain::Chains, logpdf::Function)

    # expectation of each parameter
    Eθ = reshape(mean(Array(chain), dims = [1,3]), 1,:,1)
    Echain = Chains(Eθ)
    EθD = -2*mean(logpdf(Echain))

    D = -2*logpdf(chain)
    ED = mean(D)

    pD = 2*(ED - EθD)

    DIC = EθD + pD

    return DIC, pD
end


text:
include("../aimajulia.jl");

using Test;

using Main.aimajulia;

#The following Agent tests are from the aima-python doctest

RVA = ReflexVacuumAgent();

@test execute(RVA.program, (aimajulia.loc_A, "Clean")) == "Right";

@test execute(RVA.program, (aimajulia.loc_B, "Clean")) == "Left";

@test execute(RVA.program, (aimajulia.loc_A, "Dirty")) == "Suck";

@test execute(RVA.program, (aimajulia.loc_B, "Dirty")) == "Suck";

TVE = TrivialVacuumEnvironment();

@test add_object(TVE, ModelBasedVacuumAgent()) == nothing;

@test run(TVE, steps=5) == nothing;

#=

    The following tests may fail sometimes because the tests check for the expected bounds.

    However, the results of tests that lie outside of expected bounds does not imply something is wrong.

=#

function colorize_testv_doctest_results(result::Bool)
    if (result)
        printstyled("Test Passed\n", color=:green, bold=true);
    else
        printstyled("Test Failed\n", color=:red, bold=true);
    end
end

envs = [TrivialVacuumEnvironment() for i in range(0, stop=(100 - 1))];

mbva_result = test_agent(ModelBasedVacuumAgent, 4, deepcopy(envs));
colorize_testv_doctest_results(7 < mbva_result < 11);
println("Expression: 7 < test_agent(ModelBasedVacuumAgent, 4, deepcopy(envs)) < 11");
println("Evaluated: 7 < ", mbva_result, " < 11");

refva_result = test_agent(ReflexVacuumAgent, 4, deepcopy(envs));
colorize_testv_doctest_results(5 < refva_result < 9);
println("Expression: 5 < test_agent(ReflexVacuumAgent, 4, deepcopy(envs)) < 9");
println("Evaluated: 5 < ", refva_result, " < 9");

tdva_result = test_agent(TableDrivenVacuumAgent, 4, deepcopy(envs));
colorize_testv_doctest_results(2 < tdva_result < 6);
println("Expression: 2 < test_agent(TableDrivenVacuumAgent, 4, deepcopy(envs)) < 6");
println("Evaluated: 2 < ", tdva_result, " < 6");

randva_result = test_agent(RandomVacuumAgent, 4, deepcopy(envs));
colorize_testv_doctest_results(0.5 < randva_result < 3);
println("Expression: 0.5 < test_agent(RandomVacuumAgent, 4, deepcopy(envs)) < 3");
println("Evaluated: 0.5 < ", randva_result, " < 3");

max_stars_repo_name: mikhail-j/aima-julia-aimacode-fork

max_issues_repo_name: mikhail-j/aima-julia-aimacode-fork

max_forks_repo_name: mikhail-j/aima-julia-aimacode-fork


text:
using Compat.Test

@testset "Core" begin
    include("Core.jl")
end

@testset "RegionIter" begin
    include("RegionIter.jl")
end

@testset "kron macros" begin
    include("KronMacro.jl")
end

@testset "Fused Hamiltonian" begin
    include("FusedHamiltonian.jl")
end


text:
module ReachabilityAnalysis

# ===========================================================
# Dependencies & User API
# ===========================================================
include("Initialization/init.jl")
include("Initialization/exports.jl")
include("Continuous/solve.jl")

# ===========================================================
# Structures to represent solutions of reachability problems
# ===========================================================
include("Flowpipes/setops.jl")
include("ReachSets/reachsets.jl")
include("Flowpipes/flowpipes.jl")
include("Flowpipes/arrayops.jl")
include("Flowpipes/clustering.jl")
include("Flowpipes/solutions.jl")

# ===========================================================
# Shared functionality for continuous systems
# ===========================================================
include("Continuous/fields.jl")
include("Continuous/normalization.jl")
include("Continuous/homogeneization.jl")
include("Continuous/linearization.jl")
include("Discretization/exponentiation.jl")
include("Discretization/discretization.jl")

# ===========================================================
# Reachability solver algorithms
# ===========================================================

# Linear
include("Algorithms/A20/A20.jl")
include("Algorithms/LGG09/LGG09.jl")
include("Algorithms/GLGM06/GLGM06.jl")
include("Algorithms/INT/INT.jl")
include("Algorithms/BOX/BOX.jl")
include("Algorithms/BFFPSV18/BFFPSV18.jl")
include("Algorithms/ASB07/ASB07.jl")
include("Algorithms/ORBIT/ORBIT.jl")
include("Algorithms/VREP/VREP.jl")

# Nonlinear
include("Algorithms/CARLIN/CARLIN.jl")
include("Algorithms/TMJets/TMJets21b/TMJets21b.jl")
include("Algorithms/TMJets/common.jl")
include("Algorithms/TMJets/TMJets21a/TMJets21a.jl")
include("Algorithms/TMJets/TMJets20/TMJets20.jl")
include("Algorithms/QINT/QINT.jl")

# ===========================================================
# Discrete post-operators for hybrid systems
# ===========================================================
include("Hybrid/constructors.jl")
include("Hybrid/waiting_list.jl")
include("Hybrid/transitions.jl")
include("Hybrid/time_triggered.jl")
include("Hybrid/solve.jl")

# visualization
include("Flowpipes/recipes.jl")

end # module


text:
abstract type ObjWrapper end
_manifold(ow::ObjWrapper) = ow.manifold
########
## C⁰ ##
########
struct NonDiffed{Tobj, Tman<:Manifold} <: ObjWrapper
    obj::Tobj
    manifold::Tman
    NonDiffed(f, manifold=Euclidean(0); infer::Bool=false) = new{infer ? typeof(f) : Any, typeof(manifold)}(f, manifold)
end
(od::NonDiffed)(x) = od.obj(x)
(od::NonDiffed)(x, arg, args...) = throw(ArgumentError("NonDiffed cannot be called with more than one arguement."))

########
## C¹ ##
########
struct OnceDiffed{Tobj, Tman<:Manifold, TAD} <: ObjWrapper
    obj::Tobj
    manifold::Tman
    source::TAD
    function OnceDiffed(f, manifold=Euclidean(0); infer=false)
        source = nothing # no AD gradient
        Tsource = Nothing # no AD gradient
        inferred = infer ? typeof(f) : Any
        return new{inferred, typeof(manifold), Nothing}(f, manifold, source)
    end
end
derivative_source(od::OnceDiffed{<:Any, <:Any, Nothing}) = "user supplied"
(od::OnceDiffed)(x) = od.obj(nothing, x)
(od::OnceDiffed)(x, ∇f) = od.obj(∇f, x)
(od::OnceDiffed)(x, F, J) = od.obj(J, F, x)
(od::OnceDiffed)(x, y, z, args...) = throw(ArgumentError("OnceDiffed cannot be called with more than three arguements."))

########
## C² ##
########
struct TwiceDiffed{Tobj, Tman<:Manifold, TAD} <: ObjWrapper
    obj::Tobj
    manifold::Tman
    source::TAD
    function TwiceDiffed(f, manifold=Euclidean(0); infer=false)
        source = nothing # no AD gradient
        Tsource = Nothing
        inferred = infer ? typeof(f) : Any
        return new{inferred, typeof(manifold), Tsource}(f, manifold, source)
    end
end
(td::TwiceDiffed)(x) = td.obj(nothing, nothing, x)
(td::TwiceDiffed)(x, ∇f) = td.obj(nothing, ∇f, x)
(td::TwiceDiffed)(x, ∇f, ∇²f) = td.obj(∇²f, ∇f, x)

struct LineObjective!{T1, T2, T3, T4, T5}
    obj::T1
    ∇fz::T2
    z::T3
    x::T3
    d::T4
    φ0::T5
    dφ0::T5
end
(le::LineObjective!)(λ)=(le.obj(retract!(_manifold(le.obj), le.z, le.x, le.d, λ)))
function (le::LineObjective!)(λ, calc_grad::Bool)
    f, g = (le.obj(retract!(_manifold(le.obj), le.z, le.x, le.d, λ), le.∇fz))
    f, dot(g, le.d)
end
struct LineObjective{T1, T2, T3, T4}
    obj::T1
    x::T2
    d::T3
    φ0::T4
    dφ0::T4
end
(le::LineObjective)(λ)=(le.obj(retract(_manifold(le.obj), le.x, le.d, λ)))
function (le::LineObjective)(λ, calc_grad::Bool)
    f, g = le.obj(retract(_manifold(le.obj), le.x, le.d, λ), true)
    f, dot(g, le.d)
end

struct MeritObjective{T1, T2, T3, T4}
  F::T1
  Fx::T2
  Jx::T3
  d::T4
end
_manifold(mo::MeritObjective) = mo.F.manifold
function (mo::MeritObjective)(x)
  Fx = mo.F(x, mo.Fx, nothing) # FIXME FIXME should not have to include J here
  (norm(Fx)^2)/2
end
function (mo::MeritObjective)(x, calc_grad::Bool)
  Fx, Jx = mo.F(x, mo.Fx, mo.Jx) # FIXME FIXME should not have to include J here
  (norm(Fx)^2)/2, dot(Fx, Jx*mo.d)
end
# (_y,_x)->norm(F(_x, Fx, Jx)[1])^2)


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb03wd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 20
    PMAX = 20
    LDA1 = NMAX
    LDA2 = NMAX
    LDTAU = NMAX-1
    LDZ1 = NMAX
    LDZ2 = NMAX
    LDZTA = NMAX
    LDWORK = max( NMAX, NMAX + PMAX - 2 )
    ZERO = 0.0e0
    ONE = 1.0e0
    A = Array{Float64,3}(undef, LDA1,LDA2,PMAX)
    TAU = Array{Float64,2}(undef, LDTAU,PMAX)
    Z = Array{Float64,3}(undef, LDZ1,LDZ2,PMAX)
    WR = Array{Float64,1}(undef, NMAX)
    WI = Array{Float64,1}(undef, NMAX)
    AS = Array{Float64,3}(undef, LDA1,LDA2,PMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    ZTA = Array{Float64,2}(undef, LDZTA,NMAX)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    P = parse(BlasInt, vs[2])
    ILO = parse(BlasInt, vs[3])
    IHI = parse(BlasInt, vs[4])
    ILOZ = parse(BlasInt, vs[5])
    IHIZ = parse(BlasInt, vs[6])
    JOB =  vs[7][1]
    COMPZ =  vs[8][1]
    if ( N<0 || N>min( LDA1, LDA2 ) )
        @error "Illegal N=$N"
    end
    if ( P<=0 || P>PMAX )
        @error "Illegal P=$P"
    end

    _isz,_jsz,_ksz = (N,N,P)
    for k in 1:_ksz
      vs = String[]
      while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
      end
      for i in 1:_isz
        _i0 = (i-1)*_jsz
        A[i,1:_jsz,k] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
      end
    end
    close(f)
# interp call 1

    #FOREIGN.dlacpy!( 'F', N, N, A(1,1,K), LDA1, AS(1,1,K), LDA1 )
    AS[1:N,1:N,1:P] .= A[1:N,1:N,1:P]

# interp call 2

    INFO = SLICOT.mb03vd!(N, P, ILO, IHI, A, TAU)

    if ( INFO!=0 )
        @warn "mb03vd returns info=$INFO"
        return
    end
    if ( LSAME( COMPZ, 'V' ) )
# interp call 3

    #FOREIGN.dlacpy!( 'L', N, N, A(1,1,K), LDA1, Z(1,1,K), LDZ1 )
        for k in 1:P
            Z[1:N,1:N,k] .= tril(A[1:N,1:N,k])
        end

# interp call 4

        INFO = SLICOT.mb03vy!(N, P, ILO, IHI, Z, TAU)
        @test INFO == 0
    INFO == 0 || return

        if ( INFO > 0 )
            @error "mb03vy returns info=$INFO"
        end
        # deviate from source here; continue even if COMPZ != 'V'
    end
# interp call 5

    INFO = SLICOT.mb03wd!(JOB, COMPZ, N, P, ILO, IHI, ILOZ, IHIZ, A, Z, WR, WI, LDWORK)

    if ( INFO > 0 )
        @warn "mb03wd returns info=$INFO"
        println(io, "valid eigvals:")
        imin = max(ILO, INFO+1)
        show(io, "text/plain", WR[imin:N] + im * WI[imin:N])
        println(io)
    end # if

# interp call 6

    INFO = SLICOT.mb03wx!(ILO-1, P, A, WR, WI)

    if  IHI<N
# interp call 7

        INFO = SLICOT.mb03wx!(N-IHI, P, view(A,IHI+1:N, IHI+1:N, 1:P))

    end
# unable to translate write statement:
#   in do block [('40', 'I', 'N')]
#   write WR(I), WI(I)
# interp output 1
    println(io, "A:")
    _nc = N
    _nr = N
    _nk = P
    show(io, "text/plain", A[1:_nr,1:_nc,1:_nk])
    println(io)

# interp output 2
    println(io, "Z:")
    _nc = N
    _nr = N
    _nk = P
    show(io, "text/plain", Z[1:_nr,1:_nc,1:_nk])
    println(io)

    SSQ = ZERO
    for K in 1:P
        KP1 = mod(K,P) + 1
# interp call 8

    # FOREIGN.dgemm!( 'T', 'N', N, N, N, ONE, Z(1,1,K), LDZ1, AS(1,1,K), LDA1, ZERO, ZTA, LDZTA )
        zta = Z[1:N,1:N,K]' * AS[1:N,1:N,K]
# interp call 9

    # FOREIGN.dgemm!( 'N', 'N', N, N, N, ONE, ZTA, LDZTA, Z(1,1,KP1), LDZ1, -ONE, A(1,1,K), LDA1 )
        A[1:N,1:N,K] .-= zta * Z[1:N,1:N,KP1]
        # SSQ = DLAPY2( SSQ, DLANGE( 'Frobenius', N, N, A(1,1,K), LDA1, DWORK ) )
        SSQ = hypot(SSQ, norm(A[1:N,1:N,K]))
    end
    println(io, "decomposition residual: ",SSQ)
end # run_mb03wd()


text:
module XAM

export
    SAM,
    BAM,
    header,
    eachoverlap,
    isfilled,
    seqname,
    hasseqname,
    leftposition,
    hasleftposition,
    rightposition,
    hasrightposition,
    sequence,
    hassequence

import BioGenerics: BioGenerics, distance, header, isfilled, seqname, hasseqname, sequence, hassequence, leftposition, rightposition, hasleftposition, hasrightposition
import GenomicFeatures: eachoverlap

include("sam/sam.jl")
include("bam/bam.jl")

end # module


text:
abstract type HuffmanTree end

struct HuffmanLeaf <: HuffmanTree
    ch::Char
    freq::Int
end

struct HuffmanNode <: HuffmanTree
    freq::Int
    left::HuffmanTree
    right::HuffmanTree
end

function makefreqdict(s::String)
    d = Dict{Char, Int}()
    for c in s
        if !haskey(d, c)
            d[c] = 1
        else
            d[c] += 1
        end
    end
    d
end

function huffmantree(ftable::Dict)
    trees::Vector{HuffmanTree} = [HuffmanLeaf(ch, fq) for (ch, fq) in ftable]
    while length(trees) > 1
        sort!(trees, lt = (x, y) -> x.freq < y.freq, rev = true)
        least = pop!(trees)
        nextleast = pop!(trees)
        push!(trees, HuffmanNode(least.freq + nextleast.freq, least, nextleast))
    end
    trees[1]
end

printencoding(lf::HuffmanLeaf, code) = println(lf.ch == ' ' ? "space" : lf.ch, "\t", lf.freq, "\t", code)

function printencoding(nd::HuffmanNode, code)
    code *= '0'
    printencoding(nd.left, code)
    code = code[1:end-1]

    code *= '1'
    printencoding(nd.right, code)
    code = code[1:end-1]
end

const msg = "this is an example for huffman encoding"

println("Char\tFreq\tHuffman code")

printencoding(huffmantree(makefreqdict(msg)), "")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# read in arc data

type Arc
    initNode::Int
    termNode::Int
    capacity::Float64
    freeflowtime::Float64
    flow::Float64
end

Arc(initNode::Int, termNode::Int, capacity::Float64, freeflowtime::Float64) =
    Arc(initNode, termNode, capacity, freeflowtime, 0.)

function arcData(arc_file)
    arcs = Dict{(Int, Int), Arc}()
    file = open(arc_file)
    inHeader=true
    for line in eachline(file)
        if inHeader
            inHeader = !contains(line, "Init node")
            continue
        end
        vals = split(line, )
        arcs[(int(vals[1]), int(vals[2]))] = Arc(int(vals[1]), int(vals[2]), float(vals[3]), float(vals[5]))
    end
    close(file)
    return arcs
end

# add flow data to arcs
function observFlow(arc_file, tapFlowDic)
    arcs = arcData(arc_file)
    ix = 0
    for key in keys(arcs)
        arcs[key].flow = tapFlowDic[key]
    end
    return arcs
end

# read in initial demand data
srand(123)
function iniDemand(trip_file, flag=0)
    file = open(trip_file)
    demands = Dict{(Int64,Int64), Float64}()
    s = 0
    for line in eachline(file)
        if contains(line, "Origin")
            s = int(split(line)[2])
        else
            pairs = split(line, ";")
            for pair in pairs
                if !contains(pair, "\n")
                    pair_vals = split(pair, ":")
                    t, demand = int(pair_vals[1]), float(pair_vals[2])
                    demands[(s,t)] = demand
                    if flag == 1
                        # perturb the ground truth demands slightly
                        # with perturbation factor uniformly distributed on [.9, 1.1)
                        pert_fac = 1 + 0.1 * (1 - 2 * rand())
                        demands[(s,t)] = demand * pert_fac
                    end
                end
            end
        end
    end
    close(file)
    return demands
end

function demandsDicToVec(demandsDic)
    demandsVec = zeros(length(odPairLabel_))
    for i = 1:length(demandsVec)
        demandsVec[i] = demandsDic[(odPairLabel_["$i"][1], odPairLabel_["$i"][2])]
    end
    return demandsVec
end

function demandsVecToDic(demandsVec)
    demandsDic = Dict{(Int64,Int64), Float64}()
    for i = 1:numNodes
        demandsDic[(i, i)] = 0
    end
    for i = 1:length(demandsVec)
        demandsDic[(odPairLabel_["$i"][1], odPairLabel_["$i"][2])] = demandsVec[i]
    end
    return demandsDic
end

# obtain important parameters of the network

include("load_network_uni-class.jl")

function paraNetwork(nameNetwork)
    ta_data = load_ta_network(nameNetwork)
    numNodes = maximum(map(pair->pair[1], keys(demandsDict[0])))
    start_node = ta_data.start_node
    capacity = ta_data.capacity
    free_flow_time = ta_data.free_flow_time
    numLinks = size(start_node)[1]
    numODpairs = numNodes * (numNodes - 1)
    return numNodes, numLinks, numODpairs, capacity, free_flow_time, ta_data
end

# obtain further info of the network

using JSON

function furInfo()

    #get number of routes
    numRoutes = readall("../temp_files/numRoutes_Sioux_simp.json")
    numRoutes = JSON.parse(numRoutes)

    #load OD pair-route incidence
    odPairRoute = readall("../temp_files/od_pair_route_incidence_Sioux_simp.json")
    odPairRoute = JSON.parse(odPairRoute)

    #load link-route incidence
    linkRoute = readall("../temp_files/link_route_incidence_Sioux_simp.json")
    linkRoute = JSON.parse(linkRoute)

    link_label_dict = readall("../temp_files/link_label_dict_Sioux_simp.json")
    link_label_dict = JSON.parse(link_label_dict)

    return numRoutes, odPairRoute, linkRoute, link_label_dict
end

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
using ITensors
using ITensors.NDTensors
using Test

@testset "ITensors priming and tagging" begin
  i = Index(2)
  T1 = randomITensor(i'', i')
  T2 = randomITensor(i', i)

  @test inds(@inferred(adjoint(T1))) == (i''', i'')
  @test inds(@inferred(prime(T1, 2))) == (i'''', i''')
  @test inds(@inferred(addtags(T1, "x"))) == (addtags(i, "x")'', addtags(i, "x")')
  @test inds(@inferred(T1 * T2)) == (i'', i)

  @test @inferred(order(T1)) == 2
  @test @inferred(ndims(T1)) == 2
  @test @inferred(dim(T1)) == 4
  @test @inferred(maxdim(T1)) == 2
end

@testset "NDTensors Dense contract" begin
  i = Index(2)
  T1 = randomTensor((i'', i'))
  T2 = randomTensor((i', i))
  R = randomTensor((i'', i))

  labelsT1 = (1, -1)
  labelsT2 = (-1, 2)
  labelsR = (1, 2)

  @test @inferred(NDTensors.contraction_output(T1, labelsT1, T2, labelsT2, labelsR)) isa
    DenseTensor{
    Float64,2,Tuple{Index{Int64},Index{Int64}},Dense{Float64,Vector{Float64}}
  }
  @test @inferred(NDTensors.contract(T1, labelsT1, T2, labelsT2, labelsR)) isa DenseTensor{
    Float64,2,Tuple{Index{Int64},Index{Int64}},Dense{Float64,Vector{Float64}}
  }
  @test @inferred(NDTensors.contract!!(R, labelsR, T1, labelsT1, T2, labelsT2)) isa
    DenseTensor{
    Float64,2,Tuple{Index{Int64},Index{Int64}},Dense{Float64,Vector{Float64}}
  }

  A = Base.ReshapedArray(randn(4), (2, 2), ())
  B = Base.ReshapedArray(randn(4), (2, 2), ())
  @inferred NDTensors._contract_scalar_perm!(B, A, (2, 1), 1.0, 1.0)
  @inferred NDTensors._contract_scalar_perm!(B, A, (2, 1), 1.0, 1)
end

@testset "NDTensors BlockSparse contract" begin
  i = Index([QN(0) => 2, QN(1) => 2])
  IT1 = randomITensor(i'', dag(i)')
  IT2 = randomITensor(i', dag(i))
  IR = randomITensor(i'', dag(i))
  T1, T2, R = Tensor.((IT1, IT2, IR))

  labelsT1 = (1, -1)
  labelsT2 = (-1, 2)
  labelsR = (1, 2)

  @test @inferred(NDTensors.contraction_output(T1, labelsT1, T2, labelsT2, labelsR)) isa
    Tuple{BlockSparseTensor,Vector{Tuple{Block{2},Block{2},Block{2}}}}
  @test @inferred(NDTensors.contract(T1, labelsT1, T2, labelsT2, labelsR)) isa
    BlockSparseTensor
  # TODO: this function doesn't exist yet
  #@test @inferred(NDTensors.contract!!(R, labelsR, T1, labelsT1, T2, labelsT2)) isa BlockSparseTensor

  b = Block(1, 1)
  B1 = T1[b]
  B2 = T2[b]
  BR = R[b]
  @test @inferred(
    NDTensors.contract!(BR, labelsR, B1, labelsT1, B2, labelsT2, 1.0, 0.0)
  ) isa DenseTensor
end

@testset "dmrg" begin
  N = 10
  sites = siteinds("S=1", N)
  opsum = OpSum()
  for j in 1:(N - 1)
    opsum += "Sz", j, "Sz", j + 1
    opsum += 0.5, "S+", j, "S-", j + 1
    opsum += 0.5, "S-", j, "S+", j + 1
  end
  H = MPO(opsum, sites)
  psi0 = randomMPS(sites; linkdims=10)
  sweeps = Sweeps(5)
  setmaxdim!(sweeps, 10, 20, 100, 100, 200)
  setcutoff!(sweeps, 1E-11)
  @test @inferred(Tuple{Any,MPS}, dmrg(H, psi0, sweeps; outputlevel=0)) isa
    Tuple{Float64,MPS}
end


text:
"""
    ia_model(d::StudentHousingData)

For debugging purposes get ia model without column generation framework.
"""
function ia_model(d::StudentHousingData)
    nhouses = length(d.houses)
    npatterns = length(d.patterns_allow)
    m = Model(solver = GurobiSolver(OutputFlag=0))
    @variables(m, begin
        assignment[1:nhouses, 1:npatterns] >= 0, Int
        invest[1:nhouses], Bin
    end)
    @constraints(m, begin
        [i = 1:nhouses], sum(assignment[i, p] * beds_needed(p, d) for p = 1:npatterns) <= beds_avail(i, d.houses) * invest[i]
        [p = 1:npatterns], sum(assignment[:, p]) <= d.demands[p, 1, 1]
        [i = 1:nhouses, p = 1:npatterns], assignment[i, p] <= house_fits_pattern(i, p, d)
        sum(invest[i] * maintenance(i, d.houses) for i = 1:nhouses) <= d.budget
    end)
    @objective(m, Max, sum(sum(assignment)))
    m
end

"""
    solve_ia_generation(d::StudentHousingData)

Solve problem of just investing and assigning with column generation.
"""
function solve_ia_generation(d::StudentHousingData)

    nhouses = length(d.houses)
    npatterns = length(d.patterns_allow)

    # Create an initial set of columns by randomly picking the first entity
    # allowed for any house
    V = zeros(npatterns, nhouses)
    for i = 1:nhouses
        # If there is a house that doesn't match any characteristics, skip it
        if isempty(house_allowedby(i,d))
            continue
        end
        fsb_pattern = house_allowedby(i,d)[1]
        V[fsb_pattern, i] = 1.0
    end

    # Our cost vector
    c = ones(npatterns)

    # =========================================================================
    # Master problem
    # =========================================================================
    m = Model(solver = GurobiSolver(OutputFlag=0))
    @variables(m, begin
        0 <= λ[1:nhouses] <= 1
        0 <= invest[1:nhouses] <= 1
    end)
    @objective(m, Max, sum(λ[i] * dot(c, V[:, i]) for i in 1:nhouses))
    @constraints(m, begin
        capsupply[p = 1:npatterns], sum(λ[i] * V[p, i] for i in 1:nhouses) <= d.demands[p, 1, 1]
        convexity[i = 1:nhouses], sum(λ[i]) <= 1
        ifinvested[i = 1:nhouses], λ[i] <= invest[i] # V is optional
        budget, sum(invest[i] * maintenance(i, d.houses) for i = 1:nhouses) <= d.budget
    end)

    # =========================================================================
    # The knapsack subproblem for each house
    # =========================================================================
    function get_kp(i::Int, π::Vector{Float64})
        sp = Model(solver = GurobiSolver(OutputFlag=0))
        @variable(sp, v[p = 1:npatterns] >= 0, Int)
        @constraints(sp, begin
            sum(v[p] * beds_needed(p, d) for p = 1:npatterns) <= beds_avail(i, d.houses)
            [p = 1:npatterns], v[p] <= house_fits_pattern(i, p, d)
        end)
        @objective(sp, Max, dot((c - π), v))
        sp
    end

    # Cache information about columns we generate
    V_generated = Vector{Float64}[]
    house_choice = Int[]
    λ_generated = Variable[]

    best_sp  = Model()
    best_i = 0
    best_rc = Inf
    iter = 1

    while best_rc > 0
        # Solve the restricted master problem
        solve(m)
        π = getdual(m[:capsupply])
        w = getdual(m[:convexity])
        ρ = getdual(m[:ifinvested])

        # Solve a subproblem for each house, get the column with the best reduced
        # cost over all houses
        best_rc = -Inf
        for i = 1:nhouses
            sp = get_kp(i, π)
            @assert solve(sp) == :Optimal
            col_rc = getobjectivevalue(sp) - w[i] - ρ[i]
            if col_rc > best_rc
                best_rc = col_rc
                best_sp = sp
                best_i = i
            end
        end

        best_v = getvalue(best_sp[:v])
        push!(V_generated, best_v)
        push!(house_choice, best_i)

        # Add the best column
        push!(λ_generated, @variable(m,
                lowerbound = 0.0,
                upperbound = 1.0,
                basename="_λnew_$iter",
                objective = dot(c, best_v),
                inconstraints = [capsupply; convexity[best_i]; ifinvested[best_i]],
                coefficients = [best_v; 1.0; 1.0]
              )
        )
        iter += 1
    end

    # Return the master problem and information about generated columns
    m, V_generated, λ_generated, house_choice
end


text:
using LinearAlgebra

import Base: getindex, size
import LinearAlgebra: reflectorApply!

struct QR2{T,S<:AbstractMatrix{T},V<:AbstractVector{T}} <: Factorization{T}
    factors::S
    τ::V
    QR2{T,S,V}(factors::AbstractMatrix{T}, τ::AbstractVector{T}) where {T,S<:AbstractMatrix,V<:AbstractVector} =
        new(factors, τ)
end
QR2(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = QR2{T,typeof(factors)}(factors, τ)

size(F::QR2, i::Integer...) = size(F.factors, i...)

# Similar to the definition in base but applies the reflector from the right
@inline function reflectorApply!(A::StridedMatrix, x::AbstractVector, τ::Number) # apply conjugate transpose reflector from right.
    m, n = size(A)
    if length(x) != n
        throw(DimensionMismatch("reflector must have same length as second dimension of matrix"))
    end
    @inbounds begin
        for i in 1:m
            Aiv = A[i, 1]
            for j in 2:n
                Aiv += A[i, j]*x[j]
            end
            Aiv = Aiv*τ
            A[i, 1] -= Aiv
            for j in 2:n
                A[i, j] -= Aiv*x[j]'
            end
        end
    end
    return A
end

# FixMe! Consider how to represent Q

# immutable Q{T,S<:QR2} <: AbstractMatrix{T}
#     data::S
# end

# size(A::Q) = size(A.data)
# size(A::Q, i::Integer) = size(A.data, i)

# getindex{T}(A::QR2{T}, ::Type{Tuple{:Q}}) = Q{T,typeof(A)}(A)

function getindex(A::QR2{T}, ::Type{Tuple{:R}}) where T
    m, n = size(A)
    if m >= n
        UpperTriangular(view(A.factors, 1:n, 1:n))
    else
        throw(ArgumentError("R matrix is trapezoid and cannot be extracted with indexing"))
    end
end

function getindex(A::QR2{T}, ::Type{Tuple{:QBlocked}}) where T
    m, n = size(A)
    mmn = min(m,n)
    F = A.factors
    τ = A.τ
    Tmat = zeros(T, mmn, mmn)
    for j = 1:mmn
        for i = 1:j - 1
            Tmat[i,j] = τ[i]*(F[j,i] + dot(view(F, j + 1:m, i), view(F, j + 1:m, j)))
        end
    end
    LinearAlgebra.inv!(LinearAlgebra.UnitUpperTriangular(Tmat))
    for j = 1:min(m,n)
        Tmat[j,j] = τ[j]
        for i = 1:j - 1
            Tmat[i,j] = Tmat[i,j]*τ[j]
        end
    end
    HouseholderBlock{T,typeof(F),Matrix{T}}(F, UpperTriangular(Tmat))
end

# qrUnblocked!(A::StridedMatrix) = LinearAlgebra.qrfactUnblocked!(A)
function qrUnblocked!(A::StridedMatrix{T},
                      τ::StridedVector{T} = fill(zero(T), min(size(A)...))) where {T}

    m, n = size(A)

    # Make a vector view of the first column
    a = view(A, :, 1)
    # Construct the elementary reflector of the first column in-place
    τ1 = LinearAlgebra.reflector!(a)
    # Store reflector loading in τ
    τ[1] = τ1

    # Apply the reflector to the trailing columns if any
    if n > 1
        reflectorApply!(a, τ1, view(A, :, 2:n))
    end

    # Call the function recursively on the submatrix excluding first row and column or return
    if m > 1 && n > 1
        qrUnblocked!(view(A, 2:m, 2:n), view(τ, 2:min(m, n)))
    end

    return QR2{T,typeof(A),typeof(τ)}(A, τ)
end

function qrBlocked!(A::StridedMatrix{T},
                    blocksize::Integer = 12,
                    τ::StridedVector{T} = fill(zero(T), min(size(A)...)),
                    work = Matrix{T}(undef, blocksize, size(A, 2))) where T

    m, n = size(A)

    # Make a vector view of the first column block
    A1 = view(A, :, 1:min(n, blocksize))
    # Construct QR factorization of first column block in-place
    F = qrUnblocked!(A1, view(τ, 1:min(m, n, blocksize)))

    # Apply the QR factorization to the trailing columns (if any) and recurse
    if n > blocksize
        # Make a view of the trailing columns
        A2 = view(A, :, blocksize + 1:n)
        # Apply Q' to the trailing columns
        lmul!(F[Tuple{:QBlocked}]', A2, view(work, 1:min(m, blocksize), 1:(n - blocksize)))
    end

    # Compute QR factorization of trailing block
    if m > blocksize && n > blocksize
        qrBlocked!(view(A, (blocksize + 1):m, (blocksize + 1):n),
                   blocksize,
                   view(τ, (blocksize + 1):min(m, n)),
                   work)
    end

    return QR2{T,typeof(A),typeof(τ)}(A, τ)
end


text:
"""An example of collating all model parameters from specification files."""

using Agtor
using CSV

# Data location
zone_spec_dir = "test/data/zones/"

# Load a Zone specification
zone_params = load_spec(zone_spec_dir)
tgt_zone = zone_params[:Zone_1]

@info tgt_zone

# Collate all the non-constant variables in the model
agparams = collect_agparams(tgt_zone)

@info agparams

CSV.write("collated_params.csv", agparams)

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl