text:
import DSP.Periodograms: compute_window, forward_plan, fft2oneortwosided!, nextfastfft, arraysplit
import DSP: fftouttype

function stft!(out,s::AbstractVector{T}, n::Int=length(s)>>3, noverlap::Int=n>>1;
              onesided::Bool=eltype(s)<:Real, nfft::Int=nextfastfft(n), fs::Real=1,
              window::Union{Function,AbstractVector,Nothing}=nothing) where T
    #onesided && T <: Complex && error("cannot compute one-sided FFT of a complex signal")

    win, norm2 = compute_window(window, n)
    sig_split = arraysplit(s, n, noverlap, nfft, win)
    nout = onesided ? (nfft >> 1)+1 : nfft
    #out = zeros(stfttype(T, psdonly), nout, length(sig_split))
    tmp = Vector{fftouttype(T)}(undef, T<:Real ? (nfft >> 1)+1 : nfft)
    r = fs*norm2

    plan = forward_plan(sig_split.buf, tmp)
    offset = 0
    for sig in sig_split
        mul!(tmp, plan, sig)
        fft2oneortwosided!(out, tmp, nfft, onesided, offset)
        offset += nout
    end
    out
end

# flat_t is a helper function to reshape a S x M matrix to a vector of vectors
function flat_t(t::AbstractArray{T,N}) where {T <: AbstractFloat, N}
    if size(t,1) < size(t,2)
        t = permutedims(t)
    end
    ta = Array{Vector{T}}(undef,size(t,2))
    for i in axes(t,2)
        ta[i] = view(t,:,i)
    end
    return ta
end

"""
    csm(t;n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum")

Calculate cross-spectral matrix from time series `t` which is `S x M` dimensional,
where `S` is the number of samples and `M`the number of microphones.
"""
function csm(t::AbstractArray{T};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum") where T <: AbstractFloat
    csm(flat_t(t);n=n,noverlap=noverlap,fs=fs,win=win,scaling=scaling)
end

function csm(t::Vector{Vector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum") where T <: AbstractFloat
    M = length(t)
    Nf = div(n,2)+1
    nout = div((length(t[1]) - n), n - noverlap)+1
    fc = [k*div(fs,n) for k in range(0,stop=noverlap)]
    C = Array{Complex{T}}(undef,M,M,Nf)
    S = DSP.stft.(t, n, noverlap; fs=fs, window=win, onesided=true)
    Sc = conj.(S)
    for m in 1:M
        for j in m:M
            C[j,m,:] .= dropdims(mean(LazyArray(@~ Sc[m].*S[j]),dims=2);dims=2)
        end
    end

    if scaling == "density"
        scale = 1/(fs*sum(abs2,win))
    elseif scaling == "spectrum"
        scale = 1/(sum(win).^2)
    end

    C .*= scale
    C[:,:,2:end-1] .*= 2

    for ω in 1:Nf
        C[:,:,ω] = Hermitian(C[:,:,ω],:L)
    end
    return FreqArray(C,fc)
end

"""
    csm_slow(t;n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum")

Calculate cross-spectral matrix from time series `t` which is `S x M` dimensional,
where `S` is the number of samples and `M`the number of microphones. This version is slower than
`csm` but allocates much less (in some cases), therefore `csm_slow` could be used for long time signals.
"""
function csm_slow(t::AbstractArray{T};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum") where T <: AbstractFloat
    csm_slow(flat_t(t);n=n,noverlap=noverlap,fs=fs,win=win,scaling=scaling)
end

function csm_slow(t::Vector{Vector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum") where T <: AbstractFloat
    M = length(t)
    Nf = div(n,2)+1
    nout = div((length(t[1]) - n), n - noverlap)+1
    weight = sum(abs2,win)
    fc = [k*div(fs,n) for k in range(0,stop=noverlap)]
    C = Array{Complex{T}}(undef,M,M,Nf)
    S = Array{Complex{T}}(undef,Nf,nout)
    Sc = similar(S)
    for m in 1:M
        stft!(Sc,t[m],n,noverlap; fs=fs, window=win, onesided=true)
        conj!(Sc)
        for j in m:M
            stft!(S,t[j],n,noverlap; fs=fs, window=win, onesided=true)
            C[j,m,:] .= dropdims(mean(LazyArray(@~ Sc.*S),dims=2);dims=2)
        end
    end

    if scaling == "density"
        scale = 1/(fs*sum(abs2,win))
    elseif scaling == "spectrum"
        scale = 1/(sum(win).^2)
    end

    C .*= scale
    C[:,:,2:end-1] .*= 2

    for ω in 1:Nf
        C[:,:,ω] = Hermitian(C[:,:,ω],:L)
    end
    return FreqArray(C,fc)
end


text:
# Consider using Unitful.jl.
export MaxwellUnit

struct MaxwellUnit
    # Quantities whose units are fundamental constants
    ε::Float  # permittivity unit
    μ::Float  # permeability unit
    c::Float  # speed unit
    η::Float  # impedance unit

    # General quantities (not specific to electromagnetism)
    L::Float  # length unit
    ω::Float  # angular frequency unit

    # Electromagnetic quantities, in order of frequency of appearance in Maxwell's equations
    E::Float  # E-field unit
    H::Float  # H-field unit

    J::Float  # electric current density unit
    M::Float  # magnetic current density unit

    S::Float  # Poynting vector unit
    P::Float  # power unit

    V::Float  # voltage unit
    I::Float  # electric current unit

    function MaxwellUnit(L₀::Real)
        # The length unit L₀ is set by users.
        L₀ > 0 || throw(ArgumentError("L₀ = $L₀ must be positive."));

        # The E-field unit is arbitrarily chosen as E₀ = 1 V/m.
        E₀ = 1.0  # E-field is measured in V/m

        # The following quantities are measured in their fundamental constants:
        # - Permittivity is measured in vacuum permittivity ε₀
        # - Permeability is measured in vacuum permeability μ₀
        # - Speed is measured in speed c₀ of light in vacuum
        # - Impedance is measured in vacuum impedance η₀ = √(μ₀/ε₀)

        # The units of all the other quantities are derived from the above units, such that
        # they are consistent with the above units in Maxwell's equations.  Among them, the
        # following units are independent of the length unit L₀:
        H₀ = E₀ / η₀  # H-field is measured in in E₀/η₀ (L₀-independent)
        S₀ = E₀ * H₀  # Poynting vector is measured in E₀H₀ (L₀-independent)

        # The following units are dependent on the length unit L₀:
        ω₀ = c₀ / L₀  # frequency is measured in c₀/L₀
        J₀ = H₀ / L₀  # electric current density is measured in H₀/L₀
        M₀ = E₀ / L₀  # magnetic current density is measured in E₀/L₀
        I₀ = J₀ * L₀^2  # electric current is measured in J₀L₀²
        V₀ = E₀ * L₀  # voltage is measured in E₀L₀
        P₀ = I₀ * V₀  # power is measured in I₀V₀ = S₀L₀²

        new(ε₀, μ₀, c₀, η₀, L₀, ω₀, E₀, H₀, J₀, M₀, S₀, P₀, V₀, I₀)
    end
end

ext: jl


text:
using AdventOfCodeSolutions
using Test

function input(puzzle::Puzzle{2020, 3, n}) where n
    io = openInput(puzzle)
    A = split(read(io, String), '\n', keepempty=false)
    return permutedims(hcat(map(collect, A)...))
end

function treesHitUsingSlope(A, slope)
    position = [1, 1]
    trees = 0
    while (position[1] < size(A, 1))
        position[1] = position[1] + slope[1]
        position[2] = mod1(position[2] + slope[2], size(A, 2))

        if (A[position[1], position[2]] == '#')
            trees += 1
        end
    end

    return trees
end

function solve(::Puzzle{2020, 3, 1}, A)
    return treesHitUsingSlope(A, [1, 3])
end

function solve(::Puzzle{2020, 3, 2}, A)
    slopes = [[1, 1], [1, 3], [1, 5], [1, 7], [2, 1]]
    return prod(slope -> treesHitUsingSlope(A, slope), slopes)
end

@testset "Day 3" begin
    @test(solve(Puzzle(3, 1), input(Puzzle(3, 2))) == 173)
    @test(solve(Puzzle(3, 2), input(Puzzle(3, 1))) == 4385176320)
end


text:
@testset "6.4.1 (c+d x)^m (a+b coth)^n" begin
    (a, b, c, d, e, f, m, x, ) = @variables a b c d e f m x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*coth(e+f*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(b*coth(e+f*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*coth(e+f*x)^n=#

    @test_int [x^3*coth(a + b*x), x, 6, -(x^4/4) + (x^3*log(1 - ℯ^(2*(a + b*x))))/b + (3*x^2*PolyLog(2, ℯ^(2*(a + b*x))))/(2*b^2) - (3*x*PolyLog(3, ℯ^(2*(a + b*x))))/(2*b^3) + (3*PolyLog(4, ℯ^(2*(a + b*x))))/(4*b^4)]
    @test_int [x^2*coth(a + b*x), x, 5, -(x^3/3) + (x^2*log(1 - ℯ^(2*(a + b*x))))/b + (x*PolyLog(2, ℯ^(2*(a + b*x))))/b^2 - PolyLog(3, ℯ^(2*(a + b*x)))/(2*b^3)]
    @test_int [x^1*coth(a + b*x), x, 4, -(x^2/2) + (x*log(1 - ℯ^(2*(a + b*x))))/b + PolyLog(2, ℯ^(2*(a + b*x)))/(2*b^2)]
    @test_int [coth(a + b*x)/x^1, x, 0, Unintegrable(coth(a + b*x)/x, x)]
    @test_int [coth(a + b*x)/x^2, x, 0, Unintegrable(coth(a + b*x)/x^2, x)]

    @test_int [x^3*coth(a + b*x)^2, x, 7, -(x^3/b) + x^4/4 - (x^3*coth(a + b*x))/b + (3*x^2*log(1 - ℯ^(2*(a + b*x))))/b^2 + (3*x*PolyLog(2, ℯ^(2*(a + b*x))))/b^3 - (3*PolyLog(3, ℯ^(2*(a + b*x))))/(2*b^4)]
    @test_int [x^2*coth(a + b*x)^2, x, 6, -(x^2/b) + x^3/3 - (x^2*coth(a + b*x))/b + (2*x*log(1 - ℯ^(2*(a + b*x))))/b^2 + PolyLog(2, ℯ^(2*(a + b*x)))/b^3]
    @test_int [x^1*coth(a + b*x)^2, x, 3, x^2/2 - (x*coth(a + b*x))/b + log(sinh(a + b*x))/b^2]
    @test_int [coth(a + b*x)^2/x^1, x, 0, Unintegrable(coth(a + b*x)^2/x, x)]
    @test_int [coth(a + b*x)^2/x^2, x, 0, Unintegrable(coth(a + b*x)^2/x^2, x)]

    @test_int [x^3*coth(a + b*x)^3, x, 13, -((3*x^2)/(2*b^2)) + x^3/(2*b) - x^4/4 - (3*x^2*coth(a + b*x))/(2*b^2) - (x^3*coth(a + b*x)^2)/(2*b) + (3*x*log(1 - ℯ^(2*(a + b*x))))/b^3 + (x^3*log(1 - ℯ^(2*(a + b*x))))/b + (3*PolyLog(2, ℯ^(2*(a + b*x))))/(2*b^4) + (3*x^2*PolyLog(2, ℯ^(2*(a + b*x))))/(2*b^2) - (3*x*PolyLog(3, ℯ^(2*(a + b*x))))/(2*b^3) + (3*PolyLog(4, ℯ^(2*(a + b*x))))/(4*b^4)]
    @test_int [x^2*coth(a + b*x)^3, x, 9, x^2/(2*b) - x^3/3 - (x*coth(a + b*x))/b^2 - (x^2*coth(a + b*x)^2)/(2*b) + (x^2*log(1 - ℯ^(2*(a + b*x))))/b + log(sinh(a + b*x))/b^3 + (x*PolyLog(2, ℯ^(2*(a + b*x))))/b^2 - PolyLog(3, ℯ^(2*(a + b*x)))/(2*b^3)]
    @test_int [x^1*coth(a + b*x)^3, x, 7, x/(2*b) - x^2/2 - coth(a + b*x)/(2*b^2) - (x*coth(a + b*x)^2)/(2*b) + (x*log(1 - ℯ^(2*(a + b*x))))/b + PolyLog(2, ℯ^(2*(a + b*x)))/(2*b^2)]
    @test_int [coth(a + b*x)^3/x^1, x, 0, Unintegrable(coth(a + b*x)^3/x, x)]
    @test_int [coth(a + b*x)^3/x^2, x, 0, Unintegrable(coth(a + b*x)^3/x^2, x)]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*x^m*coth(e+f*x)^(n/2)=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*coth(e+f*x))^n*with*a^2-b^2=0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+a*coth(e+f*x))^n=#

    #= ::Subsubsection:: =#
    #=n>0=#

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [(c + d*x)^3/(a + a*coth(e + f*x)), x, 5, (3*d^3*x)/(8*a*f^3) + (3*d*(c + d*x)^2)/(8*a*f^2) + (c + d*x)^3/(4*a*f) + (c + d*x)^4/(8*a*d) - (3*d^3)/(8*f^4*(a + a*coth(e + f*x))) - (3*d^2*(c + d*x))/(4*f^3*(a + a*coth(e + f*x))) - (3*d*(c + d*x)^2)/(4*f^2*(a + a*coth(e + f*x))) - (c + d*x)^3/(2*f*(a + a*coth(e + f*x)))]
    @test_int [(c + d*x)^2/(a + a*coth(e + f*x)), x, 4, (d^2*x)/(4*a*f^2) + (c + d*x)^2/(4*a*f) + (c + d*x)^3/(6*a*d) - d^2/(4*f^3*(a + a*coth(e + f*x))) - (d*(c + d*x))/(2*f^2*(a + a*coth(e + f*x))) - (c + d*x)^2/(2*f*(a + a*coth(e + f*x)))]
    @test_int [(c + d*x)^1/(a + a*coth(e + f*x)), x, 3, (d*x)/(4*a*f) + (c + d*x)^2/(4*a*d) - d/(4*f^2*(a + a*coth(e + f*x))) - (c + d*x)/(2*f*(a + a*coth(e + f*x)))]
    @test_int [1/((c + d*x)^1*(a + a*coth(e + f*x))), x, 7, -(cosh(2*e - (2*c*f)/d)*CoshIntegral((2*c*f)/d + 2*f*x))/(2*a*d) + log(c + d*x)/(2*a*d) + (CoshIntegral((2*c*f)/d + 2*f*x)*sinh(2*e - (2*c*f)/d))/(2*a*d) + (cosh(2*e - (2*c*f)/d)*SinhIntegral((2*c*f)/d + 2*f*x))/(2*a*d) - (sinh(2*e - (2*c*f)/d)*SinhIntegral((2*c*f)/d + 2*f*x))/(2*a*d)]
    @test_int [1/((c + d*x)^2*(a + a*coth(e + f*x))), x, 7, (f*cosh(2*e - (2*c*f)/d)*CoshIntegral((2*c*f)/d + 2*f*x))/(a*d^2) - 1/(d*(c + d*x)*(a + a*coth(e + f*x))) - (f*CoshIntegral((2*c*f)/d + 2*f*x)*sinh(2*e - (2*c*f)/d))/(a*d^2) - (f*cosh(2*e - (2*c*f)/d)*SinhIntegral((2*c*f)/d + 2*f*x))/(a*d^2) + (f*sinh(2*e - (2*c*f)/d)*SinhIntegral((2*c*f)/d + 2*f*x))/(a*d^2)]
    @test_int [1/((c + d*x)^3*(a + a*coth(e + f*x))), x, 8, -f/(2*a*d^2*(c + d*x)) - (f^2*cosh(2*e - (2*c*f)/d)*CoshIntegral((2*c*f)/d + 2*f*x))/(a*d^3) - 1/(2*d*(c + d*x)^2*(a + a*coth(e + f*x))) + f/(d^2*(c + d*x)*(a + a*coth(e + f*x))) + (f^2*CoshIntegral((2*c*f)/d + 2*f*x)*sinh(2*e - (2*c*f)/d))/(a*d^3) + (f^2*cosh(2*e - (2*c*f)/d)*SinhIntegral((2*c*f)/d + 2*f*x))/(a*d^3) - (f^2*sinh(2*e - (2*c*f)/d)*SinhIntegral((2*c*f)/d + 2*f*x))/(a*d^3)]

    @test_int [(c + d*x)^3/(a + a*coth(e + f*x))^2, x, 10, (-3*d^3*ℯ^(-4*e - 4*f*x))/(512*a^2*f^4) + (3*d^3*ℯ^(-2*e - 2*f*x))/(16*a^2*f^4) - (3*d^2*ℯ^(-4*e - 4*f*x)*(c + d*x))/(128*a^2*f^3) + (3*d^2*ℯ^(-2*e - 2*f*x)*(c + d*x))/(8*a^2*f^3) - (3*d*ℯ^(-4*e - 4*f*x)*(c + d*x)^2)/(64*a^2*f^2) + (3*d*ℯ^(-2*e - 2*f*x)*(c + d*x)^2)/(8*a^2*f^2) - (ℯ^(-4*e - 4*f*x)*(c + d*x)^3)/(16*a^2*f) + (ℯ^(-2*e - 2*f*x)*(c + d*x)^3)/(4*a^2*f) + (c + d*x)^4/(16*a^2*d)]
    @test_int [(c + d*x)^2/(a + a*coth(e + f*x))^2, x, 8, -(d^2*ℯ^(-4*e - 4*f*x))/(128*a^2*f^3) + (d^2*ℯ^(-2*e - 2*f*x))/(8*a^2*f^3) - (d*ℯ^(-4*e - 4*f*x)*(c + d*x))/(32*a^2*f^2) + (d*ℯ^(-2*e - 2*f*x)*(c + d*x))/(4*a^2*f^2) - (ℯ^(-4*e - 4*f*x)*(c + d*x)^2)/(16*a^2*f) + (ℯ^(-2*e - 2*f*x)*(c + d*x)^2)/(4*a^2*f) + (c + d*x)^3/(12*a^2*d)]
    @test_int [(c + d*x)^1/(a + a*coth(e + f*x))^2, x, 7, (3*d*x)/(16*a^2*f) - (d*x^2)/(8*a^2) + (x*(c + d*x))/(4*a^2) - d/(16*f^2*(a + a*coth(e + f*x))^2) - (c + d*x)/(4*f*(a + a*coth(e + f*x))^2) - (3*d)/(16*f^2*(a^2 + a^2*coth(e + f*x))) - (c + d*x)/(4*f*(a^2 + a^2*coth(e + f*x)))]
    @test_int [1/((c + d*x)^1*(a + a*coth(e + f*x))^2), x, 21, -(cosh(2*e - (2*c*f)/d)*CoshIntegral((2*c*f)/d + 2*f*x))/(2*a^2*d) + (cosh(4*e - (4*c*f)/d)*CoshIntegral((4*c*f)/d + 4*f*x))/(4*a^2*d) + log(c + d*x)/(4*a^2*d) - (CoshIntegral((4*c*f)/d + 4*f*x)*sinh(4*e - (4*c*f)/d))/(4*a^2*d) + (CoshIntegral((2*c*f)/d + 2*f*x)*sinh(2*e - (2*c*f)/d))/(2*a^2*d) + (cosh(2*e - (2*c*f)/d)*SinhIntegral((2*c*f)/d + 2*f*x))/(2*a^2*d) - (sinh(2*e - (2*c*f)/d)*SinhIntegral((2*c*f)/d + 2*f*x))/(2*a^2*d) - (cosh(4*e - (4*c*f)/d)*SinhIntegral((4*c*f)/d + 4*f*x))/(4*a^2*d) + (sinh(4*e - (4*c*f)/d)*SinhIntegral((4*c*f)/d + 4*f*x))/(4*a^2*d)]
    @test_int [1/((c + d*x)^2*(a + a*coth(e + f*x))^2), x, 24, -(1/(4*a^2*d*(c + d*x))) + cosh(2*e + 2*f*x)/(2*a^2*d*(c + d*x)) - cosh(2*e + 2*f*x)^2/(4*a^2*d*(c + d*x)) + (f*cosh(2*e - (2*c*f)/d)*CoshIntegral((2*c*f)/d + 2*f*x))/(a^2*d^2) - (f*cosh(4*e - (4*c*f)/d)*CoshIntegral((4*c*f)/d + 4*f*x))/(a^2*d^2) + (f*CoshIntegral((4*c*f)/d + 4*f*x)*sinh(4*e - (4*c*f)/d))/(a^2*d^2) - (f*CoshIntegral((2*c*f)/d + 2*f*x)*sinh(2*e - (2*c*f)/d))/(a^2*d^2) - sinh(2*e + 2*f*x)/(2*a^2*d*(c + d*x)) - sinh(2*e + 2*f*x)^2/(4*a^2*d*(c + d*x)) + sinh(4*e + 4*f*x)/(4*a^2*d*(c + d*x)) - (f*cosh(2*e - (2*c*f)/d)*SinhIntegral((2*c*f)/d + 2*f*x))/(a^2*d^2) + (f*sinh(2*e - (2*c*f)/d)*SinhIntegral((2*c*f)/d + 2*f*x))/(a^2*d^2) + (f*cosh(4*e - (4*c*f)/d)*SinhIntegral((4*c*f)/d + 4*f*x))/(a^2*d^2) - (f*sinh(4*e - (4*c*f)/d)*SinhIntegral((4*c*f)/d + 4*f*x))/(a^2*d^2)]

    @test_int [(c + d*x)^3/(a + a*coth(e + f*x))^3, x, 14, (d^3*ℯ^(-6*e - 6*f*x))/(1728*a^3*f^4) - (9*d^3*ℯ^(-4*e - 4*f*x))/(1024*a^3*f^4) + (9*d^3*ℯ^(-2*e - 2*f*x))/(64*a^3*f^4) + (d^2*ℯ^(-6*e - 6*f*x)*(c + d*x))/(288*a^3*f^3) - (9*d^2*ℯ^(-4*e - 4*f*x)*(c + d*x))/(256*a^3*f^3) + (9*d^2*ℯ^(-2*e - 2*f*x)*(c + d*x))/(32*a^3*f^3) + (d*ℯ^(-6*e - 6*f*x)*(c + d*x)^2)/(96*a^3*f^2) - (9*d*ℯ^(-4*e - 4*f*x)*(c + d*x)^2)/(128*a^3*f^2) + (9*d*ℯ^(-2*e - 2*f*x)*(c + d*x)^2)/(32*a^3*f^2) + (ℯ^(-6*e - 6*f*x)*(c + d*x)^3)/(48*a^3*f) - (3*ℯ^(-4*e - 4*f*x)*(c + d*x)^3)/(32*a^3*f) + (3*ℯ^(-2*e - 2*f*x)*(c + d*x)^3)/(16*a^3*f) + (c + d*x)^4/(32*a^3*d)]
    @test_int [(c + d*x)^2/(a + a*coth(e + f*x))^3, x, 11, (d^2*ℯ^(-6*e - 6*f*x))/(864*a^3*f^3) - (3*d^2*ℯ^(-4*e - 4*f*x))/(256*a^3*f^3) + (3*d^2*ℯ^(-2*e - 2*f*x))/(32*a^3*f^3) + (d*ℯ^(-6*e - 6*f*x)*(c + d*x))/(144*a^3*f^2) - (3*d*ℯ^(-4*e - 4*f*x)*(c + d*x))/(64*a^3*f^2) + (3*d*ℯ^(-2*e - 2*f*x)*(c + d*x))/(16*a^3*f^2) + (ℯ^(-6*e - 6*f*x)*(c + d*x)^2)/(48*a^3*f) - (3*ℯ^(-4*e - 4*f*x)*(c + d*x)^2)/(32*a^3*f) + (3*ℯ^(-2*e - 2*f*x)*(c + d*x)^2)/(16*a^3*f) + (c + d*x)^3/(24*a^3*d)]
    @test_int [(c + d*x)^1/(a + a*coth(e + f*x))^3, x, 11, (11*d*x)/(96*a^3*f) - (d*x^2)/(16*a^3) + (x*(c + d*x))/(8*a^3) - d/(36*f^2*(a + a*coth(e + f*x))^3) - (c + d*x)/(6*f*(a + a*coth(e + f*x))^3) - (5*d)/(96*a*f^2*(a + a*coth(e + f*x))^2) - (c + d*x)/(8*a*f*(a + a*coth(e + f*x))^2) - (11*d)/(96*f^2*(a^3 + a^3*coth(e + f*x))) - (c + d*x)/(8*f*(a^3 + a^3*coth(e + f*x)))]
    @test_int [1/((c + d*x)^1*(a + a*coth(e + f*x))^3), x, 53, (-3*cosh(2*e - (2*c*f)/d)*CoshIntegral((2*c*f)/d + 2*f*x))/(8*a^3*d) + (3*cosh(4*e - (4*c*f)/d)*CoshIntegral((4*c*f)/d + 4*f*x))/(8*a^3*d) - (cosh(6*e - (6*c*f)/d)*CoshIntegral((6*c*f)/d + 6*f*x))/(8*a^3*d) + log(c + d*x)/(8*a^3*d) + (CoshIntegral((6*c*f)/d + 6*f*x)*sinh(6*e - (6*c*f)/d))/(8*a^3*d) - (3*CoshIntegral((4*c*f)/d + 4*f*x)*sinh(4*e - (4*c*f)/d))/(8*a^3*d) + (3*CoshIntegral((2*c*f)/d + 2*f*x)*sinh(2*e - (2*c*f)/d))/(8*a^3*d) + (3*cosh(2*e - (2*c*f)/d)*SinhIntegral((2*c*f)/d + 2*f*x))/(8*a^3*d) - (3*sinh(2*e - (2*c*f)/d)*SinhIntegral((2*c*f)/d + 2*f*x))/(8*a^3*d) - (3*cosh(4*e - (4*c*f)/d)*SinhIntegral((4*c*f)/d + 4*f*x))/(8*a^3*d) + (3*sinh(4*e - (4*c*f)/d)*SinhIntegral((4*c*f)/d + 4*f*x))/(8*a^3*d) + (cosh(6*e - (6*c*f)/d)*SinhIntegral((6*c*f)/d + 6*f*x))/(8*a^3*d) - (sinh(6*e - (6*c*f)/d)*SinhIntegral((6*c*f)/d + 6*f*x))/(8*a^3*d)]
    @test_int [1/((c + d*x)^2*(a + a*coth(e + f*x))^3), x, 60, -(1/(8*a^3*d*(c + d*x))) + (9*cosh(2*e + 2*f*x))/(32*a^3*d*(c + d*x)) - (3*cosh(2*e + 2*f*x)^2)/(8*a^3*d*(c + d*x)) + cosh(2*e + 2*f*x)^3/(8*a^3*d*(c + d*x)) + (3*cosh(6*e + 6*f*x))/(32*a^3*d*(c + d*x)) + (3*f*cosh(2*e - (2*c*f)/d)*CoshIntegral((2*c*f)/d + 2*f*x))/(4*a^3*d^2) - (3*f*cosh(4*e - (4*c*f)/d)*CoshIntegral((4*c*f)/d + 4*f*x))/(2*a^3*d^2) + (3*f*cosh(6*e - (6*c*f)/d)*CoshIntegral((6*c*f)/d + 6*f*x))/(4*a^3*d^2) - (3*f*CoshIntegral((6*c*f)/d + 6*f*x)*sinh(6*e - (6*c*f)/d))/(4*a^3*d^2) + (3*f*CoshIntegral((4*c*f)/d + 4*f*x)*sinh(4*e - (4*c*f)/d))/(2*a^3*d^2) - (3*f*CoshIntegral((2*c*f)/d + 2*f*x)*sinh(2*e - (2*c*f)/d))/(4*a^3*d^2) - (15*sinh(2*e + 2*f*x))/(32*a^3*d*(c + d*x)) - (3*sinh(2*e + 2*f*x)^2)/(8*a^3*d*(c + d*x)) - sinh(2*e + 2*f*x)^3/(8*a^3*d*(c + d*x)) + (3*sinh(4*e + 4*f*x))/(8*a^3*d*(c + d*x)) - (3*sinh(6*e + 6*f*x))/(32*a^3*d*(c + d*x)) - (3*f*cosh(2*e - (2*c*f)/d)*SinhIntegral((2*c*f)/d + 2*f*x))/(4*a^3*d^2) + (3*f*sinh(2*e - (2*c*f)/d)*SinhIntegral((2*c*f)/d + 2*f*x))/(4*a^3*d^2) + (3*f*cosh(4*e - (4*c*f)/d)*SinhIntegral((4*c*f)/d + 4*f*x))/(2*a^3*d^2) - (3*f*sinh(4*e - (4*c*f)/d)*SinhIntegral((4*c*f)/d + 4*f*x))/(2*a^3*d^2) - (3*f*cosh(6*e - (6*c*f)/d)*SinhIntegral((6*c*f)/d + 6*f*x))/(4*a^3*d^2) + (3*f*sinh(6*e - (6*c*f)/d)*SinhIntegral((6*c*f)/d + 6*f*x))/(4*a^3*d^2)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+a*coth(e+f*x))^n*with*m*symbolic=#

    @test_int [(c + d*x)^m*(a + a*coth(e + f*x))^2, x, 0, Unintegrable((c + d*x)^m*(a + a*coth(e + f*x))^2, x)]
    @test_int [(c + d*x)^m*(a + a*coth(e + f*x))^1, x, 0, Unintegrable((c + d*x)^m*(a + a*coth(e + f*x)), x)]
    @test_int [(c + d*x)^m/(a + a*coth(e + f*x))^1, x, 2, (c + d*x)^(1 + m)/(2*a*d*(1 + m)) + (2^(-2 - m)*ℯ^(-2*e + (2*c*f)/d)*(c + d*x)^m*Gamma(1 + m, (2*f*(c + d*x))/d))/(a*f*((f*(c + d*x))/d)^m)]
    @test_int [(c + d*x)^m/(a + a*coth(e + f*x))^2, x, 4, (c + d*x)^(1 + m)/(4*a^2*d*(1 + m)) + (2^(-2 - m)*ℯ^(-2*e + (2*c*f)/d)*(c + d*x)^m*Gamma(1 + m, (2*f*(c + d*x))/d))/(a^2*f*((f*(c + d*x))/d)^m) - (4^(-2 - m)*ℯ^(-4*e + (4*c*f)/d)*(c + d*x)^m*Gamma(1 + m, (4*f*(c + d*x))/d))/(a^2*f*((f*(c + d*x))/d)^m)]
    @test_int [(c + d*x)^m/(a + a*coth(e + f*x))^3, x, 5, (c + d*x)^(1 + m)/(8*a^3*d*(1 + m)) + (3*2^(-4 - m)*ℯ^(-2*e + (2*c*f)/d)*(c + d*x)^m*Gamma(1 + m, (2*f*(c + d*x))/d))/(a^3*f*((f*(c + d*x))/d)^m) - (3*2^(-5 - 2*m)*ℯ^(-4*e + (4*c*f)/d)*(c + d*x)^m*Gamma(1 + m, (4*f*(c + d*x))/d))/(a^3*f*((f*(c + d*x))/d)^m) + (2^(-4 - m)*3^(-1 - m)*ℯ^(-6*e + (6*c*f)/d)*(c + d*x)^m*Gamma(1 + m, (6*f*(c + d*x))/d))/(a^3*f*((f*(c + d*x))/d)^m)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*coth(e+f*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*coth(e+f*x))^n=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [(c + d*x)^3*(a + b*coth(e + f*x)), x, 8, (a*(c + d*x)^4)/(4*d) - (b*(c + d*x)^4)/(4*d) + (b*(c + d*x)^3*log(1 - ℯ^(2*(e + f*x))))/f + (3*b*d*(c + d*x)^2*PolyLog(2, ℯ^(2*(e + f*x))))/(2*f^2) - (3*b*d^2*(c + d*x)*PolyLog(3, ℯ^(2*(e + f*x))))/(2*f^3) + (3*b*d^3*PolyLog(4, ℯ^(2*(e + f*x))))/(4*f^4)]
    @test_int [(c + d*x)^2*(a + b*coth(e + f*x)), x, 7, (a*(c + d*x)^3)/(3*d) - (b*(c + d*x)^3)/(3*d) + (b*(c + d*x)^2*log(1 - ℯ^(2*(e + f*x))))/f + (b*d*(c + d*x)*PolyLog(2, ℯ^(2*(e + f*x))))/f^2 - (b*d^2*PolyLog(3, ℯ^(2*(e + f*x))))/(2*f^3)]
    @test_int [(c + d*x)^1*(a + b*coth(e + f*x)), x, 6, (a*(c + d*x)^2)/(2*d) - (b*(c + d*x)^2)/(2*d) + (b*(c + d*x)*log(1 - ℯ^(2*(e + f*x))))/f + (b*d*PolyLog(2, ℯ^(2*(e + f*x))))/(2*f^2)]
    @test_int [(a + b*coth(e + f*x))/(c + d*x)^1, x, 0, Unintegrable((a + b*coth(e + f*x))/(c + d*x), x)]
    @test_int [(a + b*coth(e + f*x))/(c + d*x)^2, x, 0, Unintegrable((a + b*coth(e + f*x))/(c + d*x)^2, x)]

    @test_int [(c + d*x)^3*(a + b*coth(e + f*x))^2, x, 15, -((b^2*(c + d*x)^3)/f) + (a^2*(c + d*x)^4)/(4*d) - (a*b*(c + d*x)^4)/(2*d) + (b^2*(c + d*x)^4)/(4*d) - (b^2*(c + d*x)^3*coth(e + f*x))/f + (3*b^2*d*(c + d*x)^2*log(1 - ℯ^(2*(e + f*x))))/f^2 + (2*a*b*(c + d*x)^3*log(1 - ℯ^(2*(e + f*x))))/f + (3*b^2*d^2*(c + d*x)*PolyLog(2, ℯ^(2*(e + f*x))))/f^3 + (3*a*b*d*(c + d*x)^2*PolyLog(2, ℯ^(2*(e + f*x))))/f^2 - (3*b^2*d^3*PolyLog(3, ℯ^(2*(e + f*x))))/(2*f^4) - (3*a*b*d^2*(c + d*x)*PolyLog(3, ℯ^(2*(e + f*x))))/f^3 + (3*a*b*d^3*PolyLog(4, ℯ^(2*(e + f*x))))/(2*f^4)]
    @test_int [(c + d*x)^2*(a + b*coth(e + f*x))^2, x, 13, -((b^2*(c + d*x)^2)/f) + (a^2*(c + d*x)^3)/(3*d) - (2*a*b*(c + d*x)^3)/(3*d) + (b^2*(c + d*x)^3)/(3*d) - (b^2*(c + d*x)^2*coth(e + f*x))/f + (2*b^2*d*(c + d*x)*log(1 - ℯ^(2*(e + f*x))))/f^2 + (2*a*b*(c + d*x)^2*log(1 - ℯ^(2*(e + f*x))))/f + (b^2*d^2*PolyLog(2, ℯ^(2*(e + f*x))))/f^3 + (2*a*b*d*(c + d*x)*PolyLog(2, ℯ^(2*(e + f*x))))/f^2 - (a*b*d^2*PolyLog(3, ℯ^(2*(e + f*x))))/f^3]
    @test_int [(c + d*x)^1*(a + b*coth(e + f*x))^2, x, 9, b^2*c*x + (b^2*d*x^2)/2 + (a^2*(c + d*x)^2)/(2*d) - (a*b*(c + d*x)^2)/d - (b^2*(c + d*x)*coth(e + f*x))/f + (2*a*b*(c + d*x)*log(1 - ℯ^(2*(e + f*x))))/f + (b^2*d*log(sinh(e + f*x)))/f^2 + (a*b*d*PolyLog(2, ℯ^(2*(e + f*x))))/f^2]
    @test_int [(a + b*coth(e + f*x))^2/(c + d*x)^1, x, 0, Unintegrable((a + b*coth(e + f*x))^2/(c + d*x), x)]
    @test_int [(a + b*coth(e + f*x))^2/(c + d*x)^2, x, 0, Unintegrable((a + b*coth(e + f*x))^2/(c + d*x)^2, x)]

    @test_int [(c + d*x)^3*(a + b*coth(e + f*x))^3, x, 28, (-3*b^3*d*(c + d*x)^2)/(2*f^2) - (3*a*b^2*(c + d*x)^3)/f + (b^3*(c + d*x)^3)/(2*f) + (a^3*(c + d*x)^4)/(4*d) - (3*a^2*b*(c + d*x)^4)/(4*d) + (3*a*b^2*(c + d*x)^4)/(4*d) - (b^3*(c + d*x)^4)/(4*d) - (3*b^3*d*(c + d*x)^2*coth(e + f*x))/(2*f^2) - (3*a*b^2*(c + d*x)^3*coth(e + f*x))/f - (b^3*(c + d*x)^3*coth(e + f*x)^2)/(2*f) + (3*b^3*d^2*(c + d*x)*log(1 - ℯ^(2*(e + f*x))))/f^3 + (9*a*b^2*d*(c + d*x)^2*log(1 - ℯ^(2*(e + f*x))))/f^2 + (3*a^2*b*(c + d*x)^3*log(1 - ℯ^(2*(e + f*x))))/f + (b^3*(c + d*x)^3*log(1 - ℯ^(2*(e + f*x))))/f + (3*b^3*d^3*PolyLog(2, ℯ^(2*(e + f*x))))/(2*f^4) + (9*a*b^2*d^2*(c + d*x)*PolyLog(2, ℯ^(2*(e + f*x))))/f^3 + (9*a^2*b*d*(c + d*x)^2*PolyLog(2, ℯ^(2*(e + f*x))))/(2*f^2) + (3*b^3*d*(c + d*x)^2*PolyLog(2, ℯ^(2*(e + f*x))))/(2*f^2) - (9*a*b^2*d^3*PolyLog(3, ℯ^(2*(e + f*x))))/(2*f^4) - (9*a^2*b*d^2*(c + d*x)*PolyLog(3, ℯ^(2*(e + f*x))))/(2*f^3) - (3*b^3*d^2*(c + d*x)*PolyLog(3, ℯ^(2*(e + f*x))))/(2*f^3) + (9*a^2*b*d^3*PolyLog(4, ℯ^(2*(e + f*x))))/(4*f^4) + (3*b^3*d^3*PolyLog(4, ℯ^(2*(e + f*x))))/(4*f^4)]
    @test_int [(c + d*x)^2*(a + b*coth(e + f*x))^3, x, 22, (b^3*c*d*x)/f + (b^3*d^2*x^2)/(2*f) - (3*a*b^2*(c + d*x)^2)/f + (a^3*(c + d*x)^3)/(3*d) - (a^2*b*(c + d*x)^3)/d + (a*b^2*(c + d*x)^3)/d - (b^3*(c + d*x)^3)/(3*d) - (b^3*d*(c + d*x)*coth(e + f*x))/f^2 - (3*a*b^2*(c + d*x)^2*coth(e + f*x))/f - (b^3*(c + d*x)^2*coth(e + f*x)^2)/(2*f) + (6*a*b^2*d*(c + d*x)*log(1 - ℯ^(2*(e + f*x))))/f^2 + (3*a^2*b*(c + d*x)^2*log(1 - ℯ^(2*(e + f*x))))/f + (b^3*(c + d*x)^2*log(1 - ℯ^(2*(e + f*x))))/f + (b^3*d^2*log(sinh(e + f*x)))/f^3 + (3*a*b^2*d^2*PolyLog(2, ℯ^(2*(e + f*x))))/f^3 + (3*a^2*b*d*(c + d*x)*PolyLog(2, ℯ^(2*(e + f*x))))/f^2 + (b^3*d*(c + d*x)*PolyLog(2, ℯ^(2*(e + f*x))))/f^2 - (3*a^2*b*d^2*PolyLog(3, ℯ^(2*(e + f*x))))/(2*f^3) - (b^3*d^2*PolyLog(3, ℯ^(2*(e + f*x))))/(2*f^3)]
    @test_int [(c + d*x)^1*(a + b*coth(e + f*x))^3, x, 16, 3*a*b^2*c*x + (b^3*d*x)/(2*f) + (3*a*b^2*d*x^2)/2 + (a^3*(c + d*x)^2)/(2*d) - (3*a^2*b*(c + d*x)^2)/(2*d) - (b^3*(c + d*x)^2)/(2*d) - (b^3*d*coth(e + f*x))/(2*f^2) - (3*a*b^2*(c + d*x)*coth(e + f*x))/f - (b^3*(c + d*x)*coth(e + f*x)^2)/(2*f) + (3*a^2*b*(c + d*x)*log(1 - ℯ^(2*(e + f*x))))/f + (b^3*(c + d*x)*log(1 - ℯ^(2*(e + f*x))))/f + (3*a*b^2*d*log(sinh(e + f*x)))/f^2 + (3*a^2*b*d*PolyLog(2, ℯ^(2*(e + f*x))))/(2*f^2) + (b^3*d*PolyLog(2, ℯ^(2*(e + f*x))))/(2*f^2)]
    @test_int [(a + b*coth(e + f*x))^3/(c + d*x)^1, x, 0, Unintegrable((a + b*coth(e + f*x))^3/(c + d*x), x)]
    @test_int [(a + b*coth(e + f*x))^3/(c + d*x)^2, x, 0, Unintegrable((a + b*coth(e + f*x))^3/(c + d*x)^2, x)]

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [(c + d*x)^3/(a + b*coth(e + f*x)), x, 6, (c + d*x)^4/(4*(a + b)*d) - (b*(c + d*x)^3*log(1 - (a - b)/(ℯ^(2*(e + f*x))*(a + b))))/((a^2 - b^2)*f) + (3*b*d*(c + d*x)^2*PolyLog(2, (a - b)/(ℯ^(2*(e + f*x))*(a + b))))/(2*(a^2 - b^2)*f^2) + (3*b*d^2*(c + d*x)*PolyLog(3, (a - b)/(ℯ^(2*(e + f*x))*(a + b))))/(2*(a^2 - b^2)*f^3) + (3*b*d^3*PolyLog(4, (a - b)/(ℯ^(2*(e + f*x))*(a + b))))/(4*(a^2 - b^2)*f^4)]
    @test_int [(c + d*x)^2/(a + b*coth(e + f*x)), x, 5, (c + d*x)^3/(3*(a + b)*d) - (b*(c + d*x)^2*log(1 - (a - b)/(ℯ^(2*(e + f*x))*(a + b))))/((a^2 - b^2)*f) + (b*d*(c + d*x)*PolyLog(2, (a - b)/(ℯ^(2*(e + f*x))*(a + b))))/((a^2 - b^2)*f^2) + (b*d^2*PolyLog(3, (a - b)/(ℯ^(2*(e + f*x))*(a + b))))/(2*(a^2 - b^2)*f^3)]
    @test_int [(c + d*x)^1/(a + b*coth(e + f*x)), x, 4, (c + d*x)^2/(2*(a + b)*d) - (b*(c + d*x)*log(1 - (a - b)/(ℯ^(2*(e + f*x))*(a + b))))/((a^2 - b^2)*f) + (b*d*PolyLog(2, (a - b)/(ℯ^(2*(e + f*x))*(a + b))))/(2*(a^2 - b^2)*f^2)]
    @test_int [1/((c + d*x)^1*(a + b*coth(e + f*x))), x, 0, Unintegrable(1/((c + d*x)*(a + b*coth(e + f*x))), x)]
    @test_int [1/((c + d*x)^2*(a + b*coth(e + f*x))), x, 0, Unintegrable(1/((c + d*x)^2*(a + b*coth(e + f*x))), x)]

    @test_int [(c + d*x)^3/(a + b*coth(e + f*x))^2, x, 28, -((2*b^2*(c + d*x)^3)/((a^2 - b^2)^2*f)) + (2*b^2*(c + d*x)^3)/((a - b)*(a + b)^2*(a - b - (a + b)*ℯ^(2*e + 2*f*x))*f) + (c + d*x)^4/(4*(a - b)^2*d) + (3*b^2*d*(c + d*x)^2*log(1 - ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a^2 - b^2)^2*f^2) - (2*b*(c + d*x)^3*log(1 - ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a - b)^2*(a + b)*f) + (2*b^2*(c + d*x)^3*log(1 - ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a^2 - b^2)^2*f) + (3*b^2*d^2*(c + d*x)*PolyLog(2, ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a^2 - b^2)^2*f^3) - (3*b*d*(c + d*x)^2*PolyLog(2, ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a - b)^2*(a + b)*f^2) + (3*b^2*d*(c + d*x)^2*PolyLog(2, ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a^2 - b^2)^2*f^2) - (3*b^2*d^3*PolyLog(3, ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/(2*(a^2 - b^2)^2*f^4) + (3*b*d^2*(c + d*x)*PolyLog(3, ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a - b)^2*(a + b)*f^3) - (3*b^2*d^2*(c + d*x)*PolyLog(3, ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a^2 - b^2)^2*f^3) - (3*b*d^3*PolyLog(4, ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/(2*(a - b)^2*(a + b)*f^4) + (3*b^2*d^3*PolyLog(4, ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/(2*(a^2 - b^2)^2*f^4)]
    @test_int [(c + d*x)^2/(a + b*coth(e + f*x))^2, x, 24, -((2*b^2*(c + d*x)^2)/((a^2 - b^2)^2*f)) + (2*b^2*(c + d*x)^2)/((a - b)*(a + b)^2*(a - b - (a + b)*ℯ^(2*e + 2*f*x))*f) + (c + d*x)^3/(3*(a - b)^2*d) + (2*b^2*d*(c + d*x)*log(1 - ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a^2 - b^2)^2*f^2) - (2*b*(c + d*x)^2*log(1 - ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a - b)^2*(a + b)*f) + (2*b^2*(c + d*x)^2*log(1 - ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a^2 - b^2)^2*f) + (b^2*d^2*PolyLog(2, ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a^2 - b^2)^2*f^3) - (2*b*d*(c + d*x)*PolyLog(2, ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a - b)^2*(a + b)*f^2) + (2*b^2*d*(c + d*x)*PolyLog(2, ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a^2 - b^2)^2*f^2) + (b*d^2*PolyLog(3, ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a - b)^2*(a + b)*f^3) - (b^2*d^2*PolyLog(3, ((a + b)*ℯ^(2*e + 2*f*x))/(a - b)))/((a^2 - b^2)^2*f^3)]
    @test_int [(c + d*x)^1/(a + b*coth(e + f*x))^2, x, 5, -((c + d*x)^2/(2*(a^2 - b^2)*d)) + (b*d - 2*a*c*f - 2*a*d*f*x)^2/(4*a*(a - b)*(a + b)^2*d*f^2) + (b*(c + d*x))/((a^2 - b^2)*f*(a + b*coth(e + f*x))) + (b*(b*d - 2*a*c*f - 2*a*d*f*x)*log(1 - (a - b)/(ℯ^(2*(e + f*x))*(a + b))))/((a^2 - b^2)^2*f^2) + (a*b*d*PolyLog(2, (a - b)/(ℯ^(2*(e + f*x))*(a + b))))/((a^2 - b^2)^2*f^2)]
    @test_int [1/((c + d*x)^1*(a + b*coth(e + f*x))^2), x, 0, Unintegrable(1/((c + d*x)*(a + b*coth(e + f*x))^2), x)]
    @test_int [1/((c + d*x)^2*(a + b*coth(e + f*x))^2), x, 0, Unintegrable(1/((c + d*x)^2*(a + b*coth(e + f*x))^2), x)]
    end

max_stars_repo_path: test/6 Hyperbolic functions/6.4 Hyperbolic cotangent/6.4.1 (c+d x)^m (a+b coth)^n.jl


text:
function reverse_side(side::HorizonSide)
    return HorizonSide((Int(side) + 2) % 4)
end

function movements!(r::Robot, side::HorizonSide)
    while !isborder(r, side)
        move!(r, side)
    end
end

function put_markers_to_border!(r::Robot, side::HorizonSide)
    while !isborder(r, side)
        move!(r, side)
        putmarker!(r)
    end
end

function move_with_path!(r::Robot, side::HorizonSide)
    c = 0
    while !isborder(r, side)
        move!(r, side)
        c += 1
    end
    return c
end

function movements!(r::Robot, side::HorizonSide, n::Int)
    while n > 0
        move!(r, side)
        n -= 1
    end
end

function move_n_steps_putmarkers!(r::Robot, side::HorizonSide, n::Int)
    putmarker!(r)
    while n > 0
        move!(r, side)
        putmarker!(r)
        n -= 1
    end
end

function steps_to_border!(r::Robot, side::HorizonSide)
    n = 0
    while !isborder(r, side)
        move!(r, side)
        n += 1
    end
    return n
end

function left(side::HorizonSide)
    return(HorizonSide((Int(side) + 1) % 4))
end

function try_move!(r::Robot,side::HorizonSide)::Bool
    ort_side = left(side)
    n=0
    while isborder(r,side)
        if !isborder(r,ort_side)
            move!(r,ort_side)
            n += 1
        else
            break
        end
    end
    if isborder(r,side)
        movements!(r,reverse_side(ort_side),n)
        return false
    end
    move!(r,side)
    while (isborder(r,reverse_side(ort_side)) && n != 0)
        move!(r,side)
    end
    movements!(r,reverse_side(ort_side),n)
    return true
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
# TODO: clean this up. Only @knet8 and @knet8r are needed.

module LibKnet8

export libknet8, @knet8, @knet8r, gpu
using CUDA, Libdl, Pkg.Artifacts
const libknet8 = Libdl.find_library(["libknet8"], [artifact"libknet8"])

include("ops.jl")

macro cudacall(lib,fun,returntype,argtypes,argvalues,errmsg=true,notfound=:(error("Cannot find $lib")))
    lib = string(lib); fun = string(fun)
	if isa(argtypes,Expr); argtypes = argtypes.args; end
    if isa(argvalues,Expr); argvalues = argvalues.args; end
	if lib == "cudnn"
		path = CUDA.find_cuda_library(lib,tk,cudnnver)
	elseif lib == "knet8"
		path = libknet8
	else
		path = CUDA.find_cuda_library(lib,tk,tkver)
	end
    if path==nothing || path==""; return notfound; end
    fx = Expr(:call, :ccall, Expr(:tuple,fun,path), returntype, Expr(:tuple,argtypes...), argvalues...)
    r = gensym()
    fx = Expr(:block,esc(:($r=$fx)))
    if errmsg # || TIMER
        # if TIMER
        #     push!(fx.args, esc(:(ccall(("cudaDeviceSynchronize","libcudart"),UInt32,()))))
        # end
        if errmsg
            push!(fx.args, esc(:(if $r!=0; error($(@__MODULE__).getErrorString($lib,$fun,$r)); end)))
        else
            push!(fx.args, esc(r))
        end
        # if TIMER
        #     fx = :(@timeit to $fun ($fx))
        # end
    end
    esc(fx)
end

# macro cudnn(fun, argtypes, argvalues...); :(@cudacall("cudnn",$fun,UInt32,$argtypes,$argvalues)); end
# macro cuda(fun, argtypes, argvalues...); :(@cudacall("cuda",$fun,UInt32,$argtypes,$argvalues)); end
# macro cudart(fun, argtypes, argvalues...); :(@cudacall("cudart",$fun,UInt32,$argtypes,$argvalues)); end
# macro cudart1(fun, argtypes, argvalues...); :(@cudacall("cudart",$fun,UInt32,$argtypes,$argvalues,false,-1)); end # don't throw error
# macro cublas(fun, argtypes, argvalues...); :(@cudacall("cublas",$fun,UInt32,$argtypes,$argvalues)); end
# macro curand(fun, argtypes, argvalues...); :(@cudacall("curand",$fun,UInt32,$argtypes,$argvalues)); end
# macro nvml(fun, argtypes, argvalues...); :(@cudacall("nvml",$fun,UInt32,$argtypes,$argvalues)); end
macro knet8(fun, argtypes, argvalues...); :(@cudacall("knet8",$fun,Nothing,$argtypes,$argvalues,false)); end
macro knet8r(fun, returntype, argtypes, argvalues...); :(@cudacall("knet8",$fun,$returntype,$argtypes,$argvalues,false)); end # specify return type

function gpu(x...)
    @warn "gpu() is deprecated, please use CUDA.device instead" maxlog=1
    CUDA.functional() ? CUDA.device().handle : -1
end

# TODO: pick best gpu in init, deprecate gpu

# """

# `gpu()` returns the id of the active GPU device or -1 if none are
# active.

# `gpu(true)` resets all GPU devices and activates the one with the most
# available memory.

# `gpu(false)` resets and deactivates all GPU devices.

# `gpu(d::Int)` activates the GPU device `d` if `0 <= d < gpuCount()`,
# otherwise deactivates devices.

# `gpu(true/false)` resets all devices.  If there are any allocated
# KnetArrays their pointers will be left dangling.  Thus
# `gpu(true/false)` should only be used during startup.  If you want to
# suspend GPU use temporarily, use `gpu(-1)`.

# `gpu(d::Int)` does not reset the devices.  You can select a previous
# device and find allocated memory preserved.  However trying to operate
# on arrays of an inactive device will result in error.

# """
# function gpu end

# GPU = -1
# gpu() = GPU
# gpu(d::Integer) = (global GPU; GPU==d ? d : (try CUDA.device!(d); GPU=Int(d); catch; GPU=-1; end))
# gpu(b::Bool) = (global GPU=try b ? Int(CUDA.device().handle) : (CUDA.device_reset!();-1); catch; -1; end) # TODO: pick gpu with most memory

# gpuCount() = length(CUDA.devices())
# gpufree() = Mem.info()[1] + CUDA.pool[].cached_memory()

#=
let GPU=-1, GPUCNT=-1, CUBLAS=nothing, CUDNN=nothing, CURAND=nothing
    global gpu, gpuCount, cublashandle, cudnnhandle, curandGenerator

    gpu()=GPU

    function gpu(usegpu::Bool)
        global cudaRuntimeVersion, cudaDriverVersion, nvmlDriverVersion, nvmlVersion, nvmlfound, cudartfound
        ENV["CUDA_DEVICE_ORDER"]="PCI_BUS_ID" # From @ekinakyurek #541
        if !isdefined(@__MODULE__,:cudartfound)
            try #if (cudartfound = (Libdl.find_library(["libcudart"],[]) != ""))
                cudaRuntimeVersion = (p=Cint[0];@cudart(cudaRuntimeGetVersion,(Ptr{Cint},),p);Int(p[1]))
                cudaDriverVersion  = (p=Cint[0];@cudart(cudaDriverGetVersion, (Ptr{Cint},),p);Int(p[1]))
                cudartfound = true
            catch
                cudartfound = false
            end
        end
        if !isdefined(@__MODULE__,:nvmlfound)
            try #if (nvmlfound = (Libdl.find_library(["libnvidia-ml"],[]) != ""))
                # This code is only run once if successful, so nvmlInit here is ok
                @nvml(nvmlInit,())
                s = zeros(UInt8,80)
                @nvml(nvmlSystemGetDriverVersion,(Ptr{Cchar},Cuint),s,80)
                nvmlDriverVersion = unsafe_string(pointer(s))
                @nvml(nvmlSystemGetNVMLVersion,(Ptr{Cchar},Cuint),s,80)
                nvmlVersion = unsafe_string(pointer(s))
                # Let us keep nvml initialized for future ops such as meminfo
                # @nvml(nvmlShutdown,())
                nvmlfound = true
            catch
                nvmlfound = false
            end
        end
        if usegpu && gpuCount() > 0
            if GPU >= 0
                return
            elseif gpuCount() == 1
                pick = 0
            else # Pick based on memory usage
                pick = free = same = -1
                for i=0:gpuCount()-1
                    if nvmlfound
                        freemem = nvmlDeviceGetMemoryInfo(nvmlid(i))[2]
                    else
                        @cudart(cudaSetDevice, (Cint,), i)
                        freemem = cudaMemGetInfo()[1]
                        @cudart(cudaDeviceReset,()) # otherwise we leave a memory footprint
                    end
                    if freemem > free
                        pick = i
                        free = freemem
                        same = 1
                    elseif freemem == free
                        # pick one of equal devices randomly
                        rand(1:(same+=1)) == 1 && (pick = i)
                    end
                end
            end
            @cudart(cudaDeviceReset,())  # #541: clear memory
            gpu(pick)
        elseif gpuCount() > 0
            @cudart(cudaDeviceReset,())  # #541: clear memory
            gpu(-1)
        else
            gpu(-1)
        end
    end

    function gpu(i::Int)
        if GPU == i
            # nothing to do
        elseif 0 <= i < gpuCount()
            GPU = i
            @cudart(cudaSetDevice, (Cint,), i)

            # Initializing CUDAnative helps with stability problems for some devices
            # However cuda, nvml and cu use different device numbers! (i) is the cuda device number, CUDAdrv uses cu numbers
            # We find the equivalent cu number from pciBusId:
            # https://stackoverflow.com/questions/13781738/how-does-cuda-assign-device-ids-to-gpus
            # 20200108: CUDAdrv 5.0 initializes to currently active device, so this is no longer needed.
            # CUDAnative.initialize(CUDAdrv.CuDevice(cuid(i)))
            # 20200707: @maleadt added this:
            CUDA.device!(cuid(i))

            # Initialize curand to guard against gpu memory fillup before first dropout (#181)
            curandGenerator()

        else
            GPU = -1
            # @cudart(cudaDeviceReset,()) # may still go back and use arrays allocated in a previous gpu
        end
        return GPU
    end

    function gpuCount() # should not bomb when there is no gpu or nvidia libs
        if GPUCNT == -1
            GPUCNT = try
	        p=Cuint[0]
                # if nvmlfound
                # # We prefer nvml because cudart takes up memory even if we don't use a device
                #   @nvml(nvmlDeviceGetCount,(Ptr{Cuint},),p)
                # #541: nvml does not respect visible devices (@ekinakyurek)
                # also OSX does not have the nvidia-ml library!
                if cudartfound
                    @cudart1(cudaGetDeviceCount,(Ptr{Cuint},),p)
                end
	        Int(p[1])
            catch
	        0
            end
        end
        return GPUCNT
    end

    function cublashandle(dev=gpu())
        dev==-1 && return nothing
        i = dev+2
        if CUBLAS === nothing
            CUBLAS=Array{Any}(nothing,gpuCount()+1)
        end
        if CUBLAS[i] === nothing
            @assert dev == gpu()
            CUBLAS[i]=cublasCreate()
        end
        return CUBLAS[i]
    end

    function cudnnhandle(dev=gpu())
        dev==-1 && return nothing
        i = dev+2
        if CUDNN === nothing
            CUDNN=Array{Any}(nothing,gpuCount()+1)
        end
        if CUDNN[i] === nothing
            @assert dev == gpu()
            CUDNN[i]=cudnnCreate()
        end
        return CUDNN[i]
    end

    function curandGenerator(dev=gpu(); seed=nothing)
        dev==-1 && return nothing
        i = dev+2
        if CURAND === nothing
            CURAND=Array{Any}(nothing,gpuCount()+1)
        end
        if CURAND[i] === nothing || seed !== nothing
            @assert dev == gpu()
            if CURAND[i] !== nothing
                @curand(curandDestroyGenerator,(Cptr,),CURAND[i])
                CURAND[i] = nothing  # to prevent double free if curandCreateGenerator errors
            end
            CURAND[i]=curandCreateGenerator(seed)
        end
        return CURAND[i]
    end
end

# cudaGetDeviceCount is deprecated, use gpuCount instead:
cudaGetDeviceCount()=(p=Cuint[0];@cudart1(cudaGetDeviceCount,(Ptr{Cuint},),p);Int(p[1])) # will not bomb when there is no gpu
cudaGetDevice()=(d=Cint[-1];@cudart(cudaGetDevice,(Ptr{Cint},),d);Int(d[1]))
"Returns free,total memory."
cudaMemGetInfo()=(f=Csize_t[0];m=Csize_t[0]; @cudart(cudaMemGetInfo,(Ptr{Csize_t},Ptr{Csize_t}),f,m); (Int(f[1]),Int(m[1])))
cudaDeviceSynchronize()=@cudart(cudaDeviceSynchronize,())

# cuda, cudart, nvml each number devices differntly
# we are using the cudart ids, here is some converters:

function cudaDeviceGetPCIBusId(i=gpu())
    pciBusId = zeros(UInt8,16)
    @cudart(cudaDeviceGetPCIBusId, (Cptr, Cint, Cint), pciBusId, length(pciBusId), i)
    return unsafe_string(pointer(pciBusId))
end

function cuid(i=gpu())
    pci = cudaDeviceGetPCIBusId(i)
    id = Cint[0]
    @cuda(cuDeviceGetByPCIBusId, (Cptr, Cptr), id, pointer(pci))
    return Int(id[1])
end

function nvmlid(i=gpu())
    pci = cudaDeviceGetPCIBusId(i)
    hnd = Cptr[0]; idx = Cuint[0];
    @nvml(nvmlDeviceGetHandleByPciBusId, (Cptr, Cptr), pointer(pci), hnd)
    @nvml(nvmlDeviceGetIndex, (Cptr, Cptr), hnd[1], idx)
    return Int(idx[1])
end

"Returns total,free,used memory."
function nvmlDeviceGetMemoryInfo(i=nvmlid(gpu()))
    0 <= i || return nothing
    dev = Cptr[0]
    mem = Array{Culonglong}(undef,3)
    @nvml("nvmlDeviceGetHandleByIndex",(Cuint,Ptr{Cptr}),i,dev)
    @nvml("nvmlDeviceGetMemoryInfo",(Cptr,Ptr{Culonglong}),dev[1],mem)
    ntuple(i->Int(mem[i]),length(mem))
end

function gpumem(i=gpu())
    if i < 0 || i >= gpuCount()
        return (0,0,0)
    elseif nvmlfound
        return nvmlDeviceGetMemoryInfo(nvmlid(i))
    elseif cudartfound
        dev=gpu(); gpu(i)
        free,total = cudaMemGetInfo()
        gpu(dev)
        return (total,free,total-free)
    else
        error("gpumem cannot find nvml or cudart.")
    end
end

gpufree(i=gpu())=gpumem(i)[2]

function cublasCreate()
    handleP = Cptr[0]
    @cublas(cublasCreate_v2, (Ptr{Cptr},), handleP)
    handle = handleP[1]
    atexit(()->@cublas(cublasDestroy_v2, (Cptr,), handle))
    global cublasVersion = (p=Cint[0];@cublas(cublasGetVersion_v2,(Cptr,Ptr{Cint}),handle,p);Int(p[1]))
    return handle
end

function cudnnCreate()
    path = CUDA.find_cuda_library("cudnn",tk,cudnnver)
    if path==nothing; error("Cannot find cudnn"); end
    handleP = Cptr[0]
    @cudnn(cudnnCreate,(Ptr{Cptr},), handleP)
    handle = handleP[1]
    atexit(()->@cudnn(cudnnDestroy,(Cptr,), handle))
    global cudnnVersion = Int(eval(:(ccall(("cudnnGetVersion",$path),Csize_t,()))))
    return handle
end

function curandCreateGenerator(seed = nothing)
    CURAND_RNG_PSEUDO_DEFAULT = 100 # Default pseudorandom generator
    r = Cptr[0]
    @curand(curandCreateGenerator,(Cptr,Cint),r,CURAND_RNG_PSEUDO_DEFAULT)
    seed !== nothing && @curand(curandSetPseudoRandomGeneratorSeed,(Cptr,Culonglong),r[1],seed)
    # This call ensures memory buffer allocation. Without it we get "out of memory" later:
    @curand(curandGenerateUniform,(Cptr,Ptr{Cfloat},Csize_t),r[1],C_NULL,0)
    return r[1]
end

const curanderrors = Dict(
    0 => "No errors",
    100 => "Header file and linked library version do not match",
    101 => "Generator not initialized",
    102 => "Memory allocation failed",
    103 => "Generator is wrong type",
    104 => "Argument out of range",
    105 => "Length requested is not a multiple of dimension",
    106 => "GPU does not have double precision required by MRG32k3a",
    201 => "Kernel launch failure",
    202 => "Preexisting failure on library entry",
    203 => "Initialization of CUDA failed",
    204 => "Architecture mismatch, GPU does not support requested feature",
    999 => "Internal library error",
)

const cublaserrors = Dict(
    0 => "CUBLAS_STATUS_SUCCESS",
    1 => "CUBLAS_STATUS_NOT_INITIALIZED",
    3 => "CUBLAS_STATUS_ALLOC_FAILED",
    7 => "CUBLAS_STATUS_INVALID_VALUE",
    8 => "CUBLAS_STATUS_ARCH_MISMATCH",
    11 => "CUBLAS_STATUS_MAPPING_ERROR",
    13 => "CUBLAS_STATUS_EXECUTION_FAILED",
    14 => "CUBLAS_STATUS_INTERNAL_ERROR",
    15 => "CUBLAS_STATUS_NOT_SUPPORTED",
    16 => "CUBLAS_STATUS_LICENSE_ERROR",
)

function getErrorString(lib,fun,ret)
    if lib == "cudnn"
	path = CUDA.find_cuda_library(lib,tk,cudnnver)
    else
	path = CUDA.find_cuda_library(lib,tk,tkver)
    end
    if lib == "cudart" && path != nothing
        str = unsafe_string(@eval(ccall(("cudaGetErrorString",$path),Cstring,(UInt32,),$ret)))
    elseif lib == "cudnn" && path != nothing
        str = unsafe_string(@eval(ccall(("cudnnGetErrorString",$path),Cstring,(UInt32,),$ret)))
    elseif lib == "nvml" && path != nothing
        str = unsafe_string(@eval(ccall(("nvmlErrorString",$path),Cstring,(UInt32,),$ret)))
    elseif lib == "curand"
        str = get(curanderrors,ret,"Unknown $lib error in $fun")
    elseif lib == "cublas"
        str = get(cublaserrors,ret,"Unknown $lib error in $fun")
    else
        str = "Unknown $lib error in $fun"
    end
    string(fun, ": ", ret, ": ", str)
end

=#

#const tk = CUDA.find_toolkit()
#const tkver = isempty(tk) ? v"0" : CUDA.parse_toolkit_version(CUDA.find_cuda_binary("ptxas", tk)) # v"0" because "nothing" will crash precompilation
#const cudnnver = v"7" # cudnn version not read automatically. hardcoded to v7
# const to = TimerOutput()
#function getErrorString end

# TODO: deprecate timer?
# moved profiling option from Knet.jl to gpu.jl to make it self contained for testing
# const TIMER = haskey(ENV,"KNET_TIMER")

end # module


text:
macro load_polo_algorithm(algname)
    serial = Symbol(algname,:_s)
    consistent = Symbol(algname,:_mtc)
    inconsistent = Symbol(algname,:_mti)
    master = Symbol(algname,:_psm)
    worker = Symbol(algname,:_psw)
    scheduler = Symbol(algname,:_pss)
    quote
        global $serial = Libdl.dlsym(polo_lib, $(Meta.quot(serial)))
        global $consistent = Libdl.dlsym(polo_lib, $(Meta.quot(consistent)))
        global $inconsistent = Libdl.dlsym(polo_lib, $(Meta.quot(inconsistent)))
        global $master = Libdl.dlsym(polo_lib, $(Meta.quot(master)))
        global $worker = Libdl.dlsym(polo_lib, $(Meta.quot(worker)))
        global $scheduler = Libdl.dlsym(polo_lib, $(Meta.quot(scheduler)))
    end
end

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("POLO not installed properly, run Pkg.build(\"POLO\"), restart Julia and try again")
end
include(depsjl_path)

function __init__()
    check_deps()

    global polo_lib = Libdl.dlopen(libpolo)

    # POLO solvers
    @load_polo_algorithm(gradient)
    @load_polo_algorithm(momentum)
    @load_polo_algorithm(nesterov)
    @load_polo_algorithm(adagrad)
    @load_polo_algorithm(adadelta)
    @load_polo_algorithm(adam)
    @load_polo_algorithm(nadam)
    # Serial
    global proxgradient_s = Libdl.dlsym(polo_lib, :proxgradient_s)
    global delete_proxgradient_s = Libdl.dlsym(polo_lib, :delete_proxgradient_s)
    global run_serial = Libdl.dlsym(polo_lib, :run_serial)
    global getf_s = Libdl.dlsym(polo_lib, :getf_s)
    global getx_s = Libdl.dlsym(polo_lib, :getx_s)

    # Multithread
    global proxgradient_mt = Libdl.dlsym(polo_lib, :proxgradient_mt)
    global delete_proxgradient_mt = Libdl.dlsym(polo_lib, :delete_proxgradient_mt)
    global run_multithread = Libdl.dlsym(polo_lib, :run_multithread)
    global getf_mt = Libdl.dlsym(polo_lib, :getf_mt)
    global getx_mt = Libdl.dlsym(polo_lib, :getx_mt)

    # Parameter server
    # Master
    global proxgradient_master = Libdl.dlsym(polo_lib, :proxgradient_master)
    global delete_proxgradient_master = Libdl.dlsym(polo_lib, :delete_proxgradient_master)
    global run_master = Libdl.dlsym(polo_lib, :run_master)
    # Scheduler
    global proxgradient_scheduler = Libdl.dlsym(polo_lib, :proxgradient_scheduler)
    global delete_proxgradient_scheduler = Libdl.dlsym(polo_lib, :delete_proxgradient_scheduler)
    global run_scheduler = Libdl.dlsym(polo_lib, :run_scheduler)
    # Worker
    global proxgradient_worker = Libdl.dlsym(polo_lib, :proxgradient_worker)
    global delete_proxgradient_worker = Libdl.dlsym(polo_lib, :delete_proxgradient_worker)
    global run_worker = Libdl.dlsym(polo_lib, :run_worker)
end


text:
#=##############################################################################
# DESCRIPTION
    Vahana vehicle kinematic maneuver definition.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Feb 2021
  * License   : MIT
=###############################################################################

function generate_maneuver_vahana(; disp_plot=false, add_rotors=true, V0=0.0001)

    # NOTE: All of the following number are non-dimensional and scaled between
    # 0 and 1, with 0 and 1 beginning and end of simulation, or a scaling
    # parameter

    # End time of each stage
    #  Stage 1: [0, t1]  -> Take off
    #  Stage 2: [t1, t2] -> Transition
    #  Stage 3: [t2, t3] -> Cruise
    #  Stage 4: [t3, t4] -> Transition
    #  Stage 5: [t4, 1]  -> Landing
    t1, t2, t3, t4 = 0.2, 0.3, 0.5, 0.6

    # Target velocity at each stage (ratio of cruise velocity)
    V1 = 0.25
    V2 = 0.10
    V3 = 1.00
    V4 = V2
    V5 = 0.5*V1

    # Target RPM at each stage (ratio of hover RPM)
    RPM1 = 1.10
    RPM2 = 1.50
    RPM3 = 0.75
    RPM3_stacked = 2/600
    RPM4 = 1.00
    RPM5 = 0.90

    r_RPMh_stup = 2.00             # Ratio between stacked and main rotors in hover
    r_RPMh_tw = 0.75               # Ratio between tandem and main wing RPM in hover

    # NOTE: -x is in the direction of flight and +z is climb with ground at z=0

    ############################################################################
    # AIRCRAFT VELOCITY
    ############################################################################
    """
        Receives a nondimensional time between 0 and 1, and returns the
        vector of velocity of the aircraft at that instant.
    """
    function Vaircraft(t)

        # ------------ TAKE OFF ------------------------------------------------
        if t<t1

            # Weibull acceleration to target climb
            if t<t1/2
                val = t / (t1/2)
                Vz = V0 + V1*(1-exp(-(2*val)^5))
            # Weibull decceleration to hover
            else
                val = 1 - (t-t1/2) / (t1/2)
                Vz = V1*(1-exp(-(2*val)^5))
            end

            return [0, 0, Vz]

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t2

            # Weibull acceleration to cruise
            val = 1.5 * (t-t1) / (t2-t1) + 0.25
            Vx = V3*(1-exp(-(val)^5))

            # Slight drop during transition
            val = (1/0.8) * (t-t1) / (t2-t1)
            Vz = 0.5 * (5 * val^4 * exp(-val^5))
            Vz = 0.2*V3 * Vz

            return [-Vx, 0, -Vz]

        # ------------ CRUISE --------------------------------------------------
        elseif t<t3

            return [-V3, 0, 0]

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t4

            # Weibull decceleration to hover
            val = 1.5 * (1 - (t-t3) / (t4-t3))
            Vx = V3*(1-exp(-(val)^5))

            # Slight climb during transition
            val = 2.0 * (t-t3) / (t4-t3)
            Vz = 0.5 * (5 * val^4 * exp(-val^5))
            Vz = 0.2*V3 * Vz

            return [-Vx, 0, Vz]

        # ------------ LANDING -------------------------------------------------
        else

            # Weibull acceleration to target descend
            if (t-t4)<(1-t4)*0.35
                val = (t-t4) / ((1-t4)*0.35)
                Vz = V5*(1-exp(-(2*val)^5))
            # Weibull decceleration to hover
            else
                val = 1 - ((t-t4) - (1-t4)*0.35) / ((1-t4)*(1-0.35))
                Vz = V0 + V5*(1-exp(-(2*val)^10))
            end

            return [0, 0, -Vz]

        end
    end

    ############################################################################
    # AIRCRAFT ANGLES
    ############################################################################
    """
        Receives a nondimensional time between 0 and 1, and returns the angle
        (in degrees) of the aircraft.
        Returns: (angle_aircraft)
    """
    function angleaircraft(t)

        # ------------ TAKE OFF ------------------------------------------------
        if t<t1
            return [0, 0, 0]

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t2
            val = 1.5 * (t-t1) / (t2-t1)
            val3 = 1.5*val
            val3 = (val3)^(val3 < 1 ? 3 : 1.5)
            angle_aircraft = -15*(1.5/0.75 * val3^0.5 * exp(-(val3^1.5)))
            return [0, angle_aircraft,0 ]

        # ------------ CRUISE --------------------------------------------------
        elseif t<t3
            return [0, 0, 0]

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t4
            # Weibull decceleration to hover
            val = 1.5 * (1 - (t-t3) / (t4-t3))
            # Slight climb during transition
            val = 2.0 * (t-t3) / (t4-t3)
            angle_aircraft = 15 * 0.5 * (5 * val^4 * exp(-val^5))
            return [0, angle_aircraft, 0]

        # ------------ LANDING -------------------------------------------------
        else
            return [0, 0, 0]
        end
    end

    ############################################################################
    # TILT-WINGS ANGLES
    ############################################################################
    """
        Receives a nondimensional time between 0 and 1, and returns the angles
        (in degrees) of the main wing relative to the aircraft axis.
    """
    function angle_main(t)

        # ------------ TAKE OFF ------------------------------------------------
        if t<t1
            return [0, 90, 0]

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t2
            val = 1.5 * (t-t1) / (t2-t1)
            val2 = val - 0.00
            angle = 90 - 90*(1-exp(-(val2)^6))
            return [0, angle, 0]

        # ------------ CRUISE --------------------------------------------------
        elseif t<t3
            return [0, 0, 0]

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t4
            # Weibull decceleration to hover
            val = 1.5 * (1 - (t-t3) / (t4-t3))
            angle = 90 - 90*(1-exp(-(val)^5))
            return [0, angle, 0]

        # ------------ LANDING -------------------------------------------------
        else
            return [0, 90, 0]
        end
    end

    """
        Receives a nondimensional time between 0 and 1, and returns the angles
        (in degrees) of the tandem wing relative to the aircraft axis.
    """
    function angle_tandem(t)

        # ------------ TAKE OFF ------------------------------------------------
        if t<t1
            return [0, 90, 0]

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t2
            val = 2.0 * (t-t1) / (t2-t1)
            val1 = val + 0.25
            angle = 90 - 90*(1-exp(-(val1)^5))
            return [0, angle, 0]

        # ------------ CRUISE --------------------------------------------------
        elseif t<t3
            return [0, 0, 0]

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t4
            # Weibull decceleration to hover
            val = 1.5 * (1 - (t-t3) / (t4-t3))
            angle = 90 - 90*(1-exp(-(val)^5))
            return [0, angle, 0]

        # ------------ LANDING -------------------------------------------------
        else
            return [0, 90, 0]
        end
    end

    ############################################################################
    # ROTORS RPM
    ############################################################################
    """
        Receives a nondimensional time between 0 and 1, and returns the
        main wing RPM non-dimensionalized by the hover RPM.
    """
    function RPM_main(t)

        # ------------ TAKE OFF ------------------------------------------------
        if t<t1
            # Weibull acceleration to target climb
            if t<t1/2
                val = t / (t1/2)
                RPM_w = RPM1*(1.05-exp(-(5*val)^2))
            # Weibull decceleration to hover
            else
                val = 1 - (t-t1/2) / (t1/2)
                RPM_w = 1.0 + (RPM1-1.0)*(1-exp(-(2*val)^5))
            end

            return RPM_w

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t2
            # Increases RPM to sustain forward flight and hover
            if (t-t1)<(t2-t1)*0.25
                val = (t-t1)/((t2-t1)*0.25)
                RPM_w = 1.0 + (RPM2-1.0)*(1-exp(-(2*val)^5))
            # Decreases RPM as it tilts the wing out of hover support
            else
                val = 1 - ((t2-t)/(t2-t1)-0.25) / (1-0.25)
                RPM_w = RPM2 + (RPM3-RPM2)*(1-exp(-(1.1*val)^5))
            end

            return RPM_w

        # ------------ CRUISE --------------------------------------------------
        elseif t<t3
            RPM_w = RPM3
            return RPM_w

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t4
            # Weibull decceleration to hover
            val = 1.5 * (1 - (t-t3) / (t4-t3))
            RPM_w = RPM4 + (RPM3-RPM4)*(1-exp(-(val)^3))
            return RPM_w

        # ------------ LANDING -------------------------------------------------
        else
            # Weibull acceleration to target descend
            if (t-t4)<(1-t4)*0.55
                val = (t-t4) / ((1-t4)*0.55)
                RPM_w = RPM4 + (RPM5-RPM4)*(1-exp(-(3*val)^5))
            # Weibull decceleration to hover
            else
                val = ((t-t4) - (1-t4)*0.55) / ((1-t4)*(1-0.55))
                val = val + 0.60
                RPM_w = RPM5*(3.0*(val)^(3.0-1)*exp(-(val)^5.0))/1.0
            end

            return RPM_w

        end
    end

    """
        Receives a nondimensional time between 0 and 1, and returns the RPM of
        stacked upper rotors in main wing non-dimensionalized by the hover RPM.
    """
    function RPM_stacked_up(t)

        # ------------ TAKE OFF ------------------------------------------------
        if t<t1
            # Weibull acceleration to target climb
            if t<t1/5
                val = t / (t1/5)
                RPM_w = RPM1*(1.05-exp(-(5*val)^2))
            # Weibull decceleration to hover
            else
                val = 1 - (t-t1/5) / (1-t1/5)
                RPM_w = 1.0 + (RPM1-1.0)*(1-exp(-(2*val)^5))
            end

            return r_RPMh_stup*RPM_w

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t2
            # Increases RPM to sustain forward flight and hover
            if (t-t1)<(t2-t1)*0.65
                val = (t-t1)/((t2-t1)*0.65)
                RPM_tw = r_RPMh_tw*(1.0 + (RPM2-1.0)*(1-exp(-(2*val)^5)))
            # Decreases RPM as it tilts the wing out of hover support
            else
                val = 1 - (t2-t)/(t2-((t2-t1)*0.65 + t1))
                RPM_tw = r_RPMh_tw*(RPM2 + val*(RPM3_stacked-RPM2))
            end

            return RPM_tw

        # ------------ CRUISE --------------------------------------------------
        elseif t<t3
            return RPM3_stacked

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t4
            # Weibull decceleration to hover
            val = 1.5 * (1 - (t-t3) / (t4-t3))
            RPM_w = RPM4 + (RPM3_stacked-RPM4)*(1-exp(-(val)^3))
            return r_RPMh_stup*RPM_w

        # ------------ LANDING -------------------------------------------------
        else
            # Weibull acceleration to target descend
            if (t-t4)<(1-t4)*0.75
                val = (t-t4) / ((1-t4)*0.75)
                RPM_w = RPM4 + (RPM5-RPM4)*(1-exp(-(3*val)^5))
            # Weibull decceleration to hover
            else
                val = ((t-t4) - (1-t4)*0.75) / ((1-t4)*(1-0.75))
                val = val + 0.60
                RPM_w = RPM5*(3.0*(val)^(3.0-1)*exp(-(val)^5.0))/1.0
            end

            return r_RPMh_stup*RPM_w
        end
    end

    """
        Receives a nondimensional time between 0 and 1, and returns the RPM of
        stacked lower rotors in main wing non-dimensionalized by the hover RPM.
    """
    function RPM_stacked_low(t)

        # ------------ TAKE OFF ------------------------------------------------
        if t<t1
            return RPM_stacked_up(t)

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t2
            return RPM_stacked_up(t)

        # ------------ CRUISE --------------------------------------------------
        elseif t<t3
            return RPM_stacked_up(t)

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t4
            return RPM_stacked_up(t)

        # ------------ LANDING -------------------------------------------------
        else
            return RPM_stacked_up(t)
        end
    end

    """
        Receives a nondimensional time between 0 and 1, and returns the
        tandem wing RPM non-dimensionalized by the hover RPM.
    """
    function RPM_tandem(t)

        # ------------ TAKE OFF ------------------------------------------------
        if t<t1
            return r_RPMh_tw*RPM_main(t)

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t2
            # Increases RPM to sustain forward flight and hover
            if (t-t1)<(t2-t1)*0.65
                val = (t-t1)/((t2-t1)*0.65)
                RPM_tw = r_RPMh_tw*(1.0 + (RPM2-1.0)*(1-exp(-(2*val)^5)))
            # Decreases RPM as it tilts the wing out of hover support
            else
                val = 1 - (t2-t)/(t2-((t2-t1)*0.65 + t1))
                RPM_tw = r_RPMh_tw*(RPM2 + val*(RPM3-RPM2))
            end

            return RPM_tw

        # ------------ CRUISE --------------------------------------------------
        elseif t<t3
            return r_RPMh_tw*RPM_main(t)

        # ------------ TRANSITION ----------------------------------------------
        elseif t<t4
            # Weibull decceleration to hover
            val = 1.5 * (1 - (t-t3) / (t4-t3))
            RPM_tw = r_RPMh_tw * (RPM4 + (RPM3-RPM4)*(1-exp(-(0.8*val)^8)))
            return RPM_tw

        # ------------ LANDING -------------------------------------------------
        else
            return r_RPMh_tw*RPM_main(t)
        end
    end

    ############################################################################
    # MANEUVER OBJECT
    ############################################################################
    angle = (angle_main, angle_tandem)
    if add_rotors
        RPM = (RPM_main, RPM_stacked_up, RPM_stacked_low, RPM_tandem)
    else
        RPM = ()
    end
    Vvehicle = Vaircraft
    anglevehicle = angleaircraft

    maneuver = uns.KinematicManeuver(angle, RPM, Vvehicle, anglevehicle)

    if disp_plot
        uns.plot(maneuver)
    end

    return maneuver
end


text:
using OnlineStats, Plots, PlotThemes, Colors, Random, CSV
gr()
theme(:juno)
Random.seed!(123)

logo_colors = hcat(
    RGB(135/255, 222/255, 183/255),  # not in logo
    RGB(64/255, 191/255, 169/255),
    RGB(64/255, 169/255, 191/255),
    RGB(64/255, 124/255, 191/255),
    RGB(64/255, 79/255, 191/255),
)

# #-----------------------------------------------------------------------# readme plot
# @info "Plot 1"
# y = cumsum(randn(10^6)) + 100randn(10^6)
# o = Partition(KHist(10))
# fit!(o, y)
# p1 = plot(o)

# @info "Plot 2"
# p2 = plot(
#     plot(fit!(KHist(50), y[1:5_000])),
#     plot(fit!(KHist(50), y[1:200_000])),
#     plot(fit!(KHist(50), y[1:500_000])),
#     plot(fit!(KHist(50), y))
# )

# @info "Plot 3"
# p3 = plot(fit!(PlotNN(200), randn(10^5, 2)), colorbar=false)

# @info "Plot 4"
# df = CSV.read("/Users/joshday/datasets/diamonds.csv"; allowmissing=:none)
# o = Partition(CountMap(String), 75)
# p4 = plot(fit!(o, string.(df[:cut])), ylim = (0, df[1][end]))

# @info "Joining Plots"
# p = plot(p1, p2, p4, p3; xaxis=false, yaxis=false,xticks=[], yticks=[], size=(1000,800), legend=false)
# png("/Users/joshday/Desktop/readme.png")
# p

#-----------------------------------------------------------------------# Weight animation
# n = 100
# Random.seed!(123)
# y = .5cumsum(rand(n)) + 3randn(n)

# os = (
#     Mean(weight = EqualWeight()),
#     Mean(weight = ExponentialWeight(.1)),
#     Mean(weight = LearningRate(.6)),
#     Mean(weight = HarmonicWeight(10.)),
#     Mean(weight = McclainWeight(.1))
# )

# p = plot(zeros(0), zeros(0, 6), seriestype = [:scatter :line :line :line :line :line], w=3,
#     label = ["data" "EqualWeight" "ExponentialWeight(.1)" "LearningRate(.6)" "HarmonicWeight(10)" "McclainWeight(.1)"], ms=3,
#     linestyle=:auto,
#     grid=false, legend = :topleft, xlab = "Number of Observations", title="Means",
#     color = hcat(:black, logo_colors...), ylim = (-5, 30))

# a = @animate for yi in y
#     for o in os
#         fit!(o, yi)
#     end
#     vals = vcat(yi, map(x -> value(x)[1], os)...)
#     push!(p, nobs(os[1]), vals)
# end

# gif(a, @__DIR__() * "/weights.gif")

#-----------------------------------------------------------------------# readme animation
# Updated 3/31/18
# n = 100
# Random.seed!(123)
# y = .5cumsum(rand(n)) + 3randn(n)

# o1 = Mean()
# o2 = Mean(weight = x -> .1)
# p = plot(zeros(0), zeros(0, 3), seriestype = [:scatter :line :line], w=3,
#     label = ["data", "equally weighted", "exponentially weighted"],
#     grid=false, legend = :topleft, xlab = "Number of Observations",
#     title = "Means", ylim = (-5, 30), color = :viridis
# )

# a = @animate for yi in y
#     fit!(o1, yi)
#     fit!(o2, yi)
#     push!(p, nobs(o1), vcat(yi, value(o1), value(o2)))
# end

# gif(a, "/Users/joshday/Desktop/readme.gif")

#-----------------------------------------------------------------------# weight
p1 = plot(EqualWeight())
p2 = plot(ExponentialWeight(.1))
p4 = plot(LearningRate())
p5 = plot(HarmonicWeight())
p6 = plot(McclainWeight())

p = plot(p1, p2, p4, p5, p6, size = (1200, 800), font = font("Bangla MN", 12), linecolor=logo_colors[2])

png(p, @__DIR__() * "/weights.png")

#-----------------------------------------------------------------------# Logo
# o = fit!(Partition(KHist(8)), sin.(0:.001:2π) + .4randn(length(0:.001:2π)))
# p = plot(o, color = :viridis, grid=false, axis=false, legend=false, ms=4
# png(p, "/Users/joshday/Desktop/logo.png")


text:
include(joinpath(DFC.DEPS_DIR, "abinitflags.jl"))
const AbinitFlags = _ABINITFLAGS()


text:
"""
    mul_matrix_from_left(A, X)

```
     Λ  ┌───┐
(1)─1X2─1 A 2─(2)
     V  └─3─┘
          │
         (3)
```
"""
function mul_matrix_from_left(A::AbstractTensor3, X::AbstractMatrix)
    @tensoropt (p = 1, (l, r, x) = χ) mul_matrix_from_left[l, r, p] := X[l, x] * A[x, r, p]
end

"""
    mul_matrix_from_right(A, X)

```
    ┌───┐  Λ
(1)─1 A 2─1X2─(2)
    └─3─┘  V
      │
     (3)
```
"""
function mul_matrix_from_right(A::AbstractTensor3, X::AbstractMatrix)
    @tensoropt (p = 1, (l, r, x) = χ) mul_matrix_from_right[l, r, p] := A[l, x, p] * X[x, r]
end

"""
    mul_operator_onsite(A, O)

```
    ┌───┐
(1)─1 A 2─(2)
    └─3─┘
      │
     .2.
    ( O )
     `1'
      │
     (3)
```
"""
function mul_operator_onsite(A::AbstractTensor3, O::AbstractMatrix)
    @tensoropt ((p, x) = 1, (l, r) = χ) mul_operator_onsite[l, r, p] := A[l, r, x] * O[p, x]
end

"""
    mul_operator_with_left(A, O, AL)

```
 ┌────────┐  ┌───┐
┌1   AL   2──1 A 2─┐
│└────3───┘  └─3─┘ │
│     │        │   │
│    .3────────4.  │
│   (      O     ) │
│    `1────────2'  │
│     │        │   │
│┌────3───┐    │   │
└1conj(AL)2─┐  │  ┌┘
 └────────┘ │  │  │
(1)─────────┘ (3) └(2)
```
"""
function mul_operator_with_left(A::AbstractTensor3, O::AbstractTensor4, AL::AbstractTensor3)
    @tensoropt ((blp, p, klp, krp) = 1, (l, r, x, y) = χ) mul_operator_with_left[l, r, p] :=
        conj(AL)[y, l, blp] * AL[y, x, klp] * A[x, r, krp] * O[blp, p, klp, krp]
end

"""
    mul_operator_with_right(A, O, AR)

```
    ┌───┐  ┌────────┐
  ┌─1 A 2──1   AR   2┐
  │ └─3─┘  └───3────┘│
  │   │        │     │
  │  .3────────4.    │
  │ (     O      )   │
  │  `1────────2'    │
  │   │        │     │
  │   │    ┌───3────┐│
  └┐  │  ┌─1conj(AR)2┘
   │  │  │ └────────┘
(1)┘ (3) └───────────(2)
```
"""
function mul_operator_with_right(
    A::AbstractTensor3,
    O::AbstractTensor4,
    AR::AbstractTensor3,
)
    @tensoropt ((p, brp, klp, krp) = 1, (l, r, x, y) = χ) mul_operator_with_right[
        l,
        r,
        p,
    ] := A[l, x, klp] * AR[x, y, krp] * conj(AR)[r, y, brp] * O[p, brp, klp, krp]
end


text:
"""
# Description

Given an expression, `e`, convert it into an anonymous function over tuples. It
is assumed that `index` defines a mapping from symbols to indices that contains
an entry for every symbol that should be treated as a column name in `e`.

# Arguments

1. `e::Any`: An expression to be converted into an anonymous function
    expression.
2. `tuple_name::Any`: The local name of the tuple that the function will act
    upon. To maximize the readability of the resulting function, pass in a
    readable like `:t`. Use `safe_tuple_name` to find a name that is guaranteed
    to be safe to use.
3. `index::Dict{Symbol, Int}`: A mapping from column name symbosl to indices in
    the input tuple.
4. `passes::NamedTuple`: A specification of which rewrite passes should be
    performed on the input. There are three that can be activated:
    1. `locals`: A pass to unescape locals as in `a + \$b`.
    2. `lift`: A pass to lift all function calls. See docs/lifting.md for more
        details.
    3. `tvl`: A pass to rewrite all uses of `||` and `&&` to use three-valued
        logic.
    Currently `lift` is disabled as the functionality is not yet complete. The
    others are enabled by default, but can be turned off for easier debugging.

# Return Values

1. `e′::Any`: A new expression that defines an anonymous function in tuple
    form. See docs/forms.md for details.

# Examples

```
julia> tuple_form(:(a + b), :t, Dict(:a => 1, :b => 2))
:(t->begin
          #= REPL[21]:55 =#
          ($(Expr(:escape, :+)))(t[1], t[2])
      end)
```
"""
function tuple_form(
    @nospecialize(e::Any),
    column_names::NTuple{N, ColumnName},
    passes::NamedTuple = (locals=true, lift=false, tvl=true),
) where N
    tuple_name = safe_tuple_name(e)

    index = index_column_names(column_names)

    body = rewrite_column_names(e, tuple_name, index)

    if passes.locals
        body = unescape_locals(body)
    end

    if passes.lift
        body = lift_function_calls(body)
    end

    if passes.tvl
        body = tvl(body)
    end

    body = escape_symbols(body, (tuple_name, :uses_default_lifting))

    :($tuple_name -> $body)
end


text:
function PSO_velocity(x, v, pbest, gbest, C1, C2, ω)
    r1 = C1 * rand()
    r2 = C2 * rand()
    ω * v + r1 * (pbest - x) + r2 * (gbest - x)
end

function velocity(x, v, pbest, gbest, parameters)
    PSO_velocity(x, v, pbest, gbest, parameters.C1, parameters.C2, parameters.ω)
end


text:
export cache

## CachedOperator

type CachedOperator{T<:Number,DM<:AbstractMatrix,M<:Operator,DS,RS,BI} <: Operator{T}
    op::M
    data::DM
    datasize::Tuple{Int,Int}
    domainspace::DS
    rangespace::RS
    bandinds::BI
    padding::Bool   # records whether the operator should be padded for QR operations
                    #  note that BandedMatrix/AlmostBandedMatrix don't need this
                    #  as the padding is done in the bandinds
end

CachedOperator(op::Operator,data::AbstractMatrix,sz::Tuple{Int,Int},ds,rs,bi,pd=false) =
    CachedOperator{eltype(data),typeof(data),typeof(op),
                    typeof(ds),typeof(rs),typeof(bi)}(op,data,sz,ds,rs,bi,pd)

CachedOperator(op::Operator,data::AbstractMatrix,sz::Tuple{Int,Int},pd=false) =
    CachedOperator(op,data,sz,domainspace(op),rangespace(op),bandinds(op),pd)
CachedOperator(op::Operator,data::AbstractMatrix,padding=false) = CachedOperator(op,data,size(data),padding)
function CachedOperator(op::Operator;padding::Bool=false)
    if isbanded(op)
        l,u=bandwidths(op)
        padding && (u+=l)
        data=BandedMatrix(eltype(op),0,0,l,u)
        CachedOperator(op,data,size(data),domainspace(op),rangespace(op),(-l,u),padding)
    elseif israggedbelow(op)
        CachedOperator(op,RaggedMatrix(eltype(op),0,Int[]),padding)
    else
        CachedOperator(op,Array(eltype(op),0,0),padding)
    end
end

cache(O::Operator;kwds...) = CachedOperator(O;kwds...)

Base.convert{T}(::Type{Operator{T}},S::CachedOperator{T}) = S
Base.convert{T}(::Type{Operator{T}},S::CachedOperator) =
    CachedOperator(Operator{T}(S.op),AbstractMatrix{T}(S.data),
                    S.datasize,S.domainspace,S.rangespace,S.bandinds)

domainspace(C::CachedOperator) = C.domainspace
rangespace(C::CachedOperator) = C.rangespace
bandinds(C::CachedOperator) = C.bandinds
Base.stride(C::CachedOperator) = stride(C.op)

# when bandinds are infinite, use colstart/colstop from operator
# TODO: cache this information as well
for func in (:(ApproxFun.colstart),:(ApproxFun.colstop),
                :(ApproxFun.rowstart),:(ApproxFun.rowstop))
    @eval $func{T,DM,M,DS,RS}(D::CachedOperator{T,DM,M,DS,RS,Tuple{Infinity{Bool},Infinity{Bool}}},
                         k::Integer) = $func(D.op,k)
end

function Base.getindex(B::CachedOperator,k::Integer,j::Integer)
    resizedata!(B,k,j)
    B.data[k,j]
end

function Base.getindex(B::CachedOperator,k::Range,j::Range)
    if !isempty(k) && !isempty(j)
        resizedata!(B,maximum(k),maximum(j))
        B.data[k,j]
    else
        Array(eltype(B),length(k),length(j))
    end
end

function Base.getindex(B::CachedOperator,k::Integer)
    if size(B,1)==1
        B[1,k]
    elseif size(B,2)==1
        B[k,1]
    else
        error("Not implemented")
    end
 end

function resizedata!{T<:Number}(B::CachedOperator{T,Matrix{T}},n::Integer,m::Integer)
    if n > size(B,1) || m > size(B,2)
        throw(ArgumentError("Cannot resize beyound size of operator"))
    end

    # this does nothing if already in dimensions
    N,M=size(B.data)
    if n > N && m > M
        B.data = unsafe_resize!(B.data,n,m)
    elseif n > N
        B.data = unsafe_resize!(B.data,n,:)
    elseif m > M
        B.data = unsafe_resize!(B.data,:,m)
    end

    if n ≤ B.datasize[1] && m ≤ B.datasize[2]
        # do nothing
        B
    elseif n ≤ B.datasize[1]
        kr,jr=1:B.datasize[1],B.datasize[2]+1:m
        B.data[kr,jr] = B.op[kr,jr]
        B.datasize = (B.datasize[1],m)
        B
    elseif m ≤ B.datasize[2]
        kr,jr=B.datasize[1]+1:n,1:B.datasize[2]
        B.data[kr,jr] = B.op[kr,jr]
        B.datasize = (n,B.datasize[2])
        B
    else
        # resize rows then columns
        resizedata!(resizedata!(B,n,B.datasize[2]),n,m)
    end
end

function resizedata!{T<:Number}(B::CachedOperator{T,BandedMatrix{T}},n::Integer,::Colon)
    if n > size(B,1)
        throw(ArgumentError("Cannot resize beyound size of operator"))
    end

    if n > B.datasize[1]
        pad!(B.data,2n,:)

        kr=B.datasize[1]+1:n
        jr=max(B.datasize[1]+1-B.data.l,1):n+B.data.u
        BLAS.axpy!(1.0,view(B.op,kr,jr),view(B.data,kr,jr))

        B.datasize = (n,n+B.data.u)
    end

    B
end

function resizedata!{T<:Number}(B::CachedOperator{T,RaggedMatrix{T}},::Colon,n::Integer)
    if n > size(B,2)
        throw(ArgumentError("Cannot resize beyound size of operator"))
    end

    if n > B.datasize[2]
        resize!(B.data.cols,n+1)

        if B.padding
            # K is largest colstop.  We get previous largest by looking at precalulated
            # cols
            K = B.datasize[2]==0?0:B.data.cols[B.datasize[2]+1]-B.data.cols[B.datasize[2]]

            for j = B.datasize[2]+1:n-1
                K = max(K,colstop(B.op,j))
                B.data.cols[j+1] = B.data.cols[j] + K
            end
            K = max(K,colstop(B.op,n))
        else
            for j = B.datasize[2]+1:n-1
                B.data.cols[j+1] = B.data.cols[j] + colstop(B.op,j)
            end
            K = colstop(B.op,n)
        end

        B.data.cols[n+1] = B.data.cols[n] + K
        pad!(B.data.data,B.data.cols[n+1]-1)
        B.data.m = K

        jr=B.datasize[2]+1:n
        kr=1:K
        BLAS.axpy!(1.0,view(B.op,kr,jr),view(B.data,kr,jr))

        B.datasize = (K,n)

    end

    B
end

resizedata!{T<:Number}(B::CachedOperator{T,BandedMatrix{T}},n::Integer,m::Integer) =
    resizedata!(B,n,:)

resizedata!{T<:Number}(B::CachedOperator{T,RaggedMatrix{T}},n::Integer,m::Integer) =
    resizedata!(B,:,m)

resizedata!(B::CachedOperator,::Colon,m::Integer) = resizedata!(B,size(B,1),m)
resizedata!(B::CachedOperator,n::Integer,::Colon) = resizedata!(B,n,size(B,2))


text:
export mlstyle_add_deprecation_msg!, mlstyle_report_deprecation_msg!
const threadsafe_msgs = Vector{String}[String[] for _ in Base.Threads.nthreads()]

function mlstyle_add_deprecation_msg!(msg::String)
    push!(threadsafe_msgs[Base.Threads.threadid()], msg)
end

function mlstyle_report_deprecation_msg!(ln::LineNumberNode)
    msgs = threadsafe_msgs[Base.Threads.threadid()]
    isempty(msgs) || begin
        @warn "Deprecated use detected at $(ln.file):$(ln.line)."
        for msg in msgs
            @warn msg
        end
        empty!(msgs)
    end
end