text:
@testset "6.4.2 Hyperbolic cotangent functions" begin
    (a, b, c, d, e, m, n, p, x, ) = @variables a b c d e m n p x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*Involving*Hyperbolic*Cotangents=#

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*coth(c+d*x))^(n/2)=#

    @test_int [(b*coth(c + d*x))^(7/2),x, 7, (b^(7/2)*atan(sqrt(b*coth(c + d*x))/sqrt(b)))/d + (b^(7/2)*atanh(sqrt(b*coth(c + d*x))/sqrt(b)))/d - (2*b^3*sqrt(b*coth(c + d*x)))/d - (2*b*(b*coth(c + d*x))^(5/2))/(5*d)]
    @test_int [(b*coth(c + d*x))^(5/2),x, 6, -((b^(5/2)*atan(sqrt(b*coth(c + d*x))/sqrt(b)))/d) + (b^(5/2)*atanh(sqrt(b*coth(c + d*x))/sqrt(b)))/d - (2*b*(b*coth(c + d*x))^(3/2))/(3*d)]
    @test_int [(b*coth(c + d*x))^(3/2),x, 6, (b^(3/2)*atan(sqrt(b*coth(c + d*x))/sqrt(b)))/d + (b^(3/2)*atanh(sqrt(b*coth(c + d*x))/sqrt(b)))/d - (2*b*sqrt(b*coth(c + d*x)))/d]
    @test_int [(b*coth(c + d*x))^(1/2), x, 5, -((sqrt(b)*atan(sqrt(b*coth(c + d*x))/sqrt(b)))/d) + (sqrt(b)*atanh(sqrt(b*coth(c + d*x))/sqrt(b)))/d]
    @test_int [1/(b*coth(c + d*x))^(1/2), x, 5, atan(sqrt(b*coth(c + d*x))/sqrt(b))/(sqrt(b)*d) + atanh(sqrt(b*coth(c + d*x))/sqrt(b))/(sqrt(b)*d)]
    @test_int [1/(b*coth(c + d*x))^(3/2), x, 6, -(atan(sqrt(b*coth(c + d*x))/sqrt(b))/(b^(3/2)*d)) + atanh(sqrt(b*coth(c + d*x))/sqrt(b))/(b^(3/2)*d) - 2/(b*d*sqrt(b*coth(c + d*x)))]
    @test_int [1/(b*coth(c + d*x))^(5/2), x, 6, atan(sqrt(b*coth(c + d*x))/sqrt(b))/(b^(5/2)*d) + atanh(sqrt(b*coth(c + d*x))/sqrt(b))/(b^(5/2)*d) - 2/(3*b*d*(b*coth(c + d*x))^(3/2))]
    @test_int [1/(b*coth(c + d*x))^(7/2), x, 7, -(atan(sqrt(b*coth(c + d*x))/sqrt(b))/(b^(7/2)*d)) + atanh(sqrt(b*coth(c + d*x))/sqrt(b))/(b^(7/2)*d) - 2/(5*b*d*(b*coth(c + d*x))^(5/2)) - 2/(b^3*d*sqrt(b*coth(c + d*x)))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*coth(c+d*x))^(n/3)=#

    @test_int [(b*coth(c + d*x))^(4/3),x, 13, -((sqrt(3)*b^(4/3)*atan((1 - (2*(b*coth(c + d*x))^(1/3))/b^(1/3))/sqrt(3)))/(2*d)) + (sqrt(3)*b^(4/3)*atan((1 + (2*(b*coth(c + d*x))^(1/3))/b^(1/3))/sqrt(3)))/(2*d) + (b^(4/3)*atanh((b*coth(c + d*x))^(1/3)/b^(1/3)))/d - (3*b*(b*coth(c + d*x))^(1/3))/d - (b^(4/3)*log(b^(2/3) - b^(1/3)*(b*coth(c + d*x))^(1/3) + (b*coth(c + d*x))^(2/3)))/(4*d) + (b^(4/3)*log(b^(2/3) + b^(1/3)*(b*coth(c + d*x))^(1/3) + (b*coth(c + d*x))^(2/3)))/(4*d)]
    @test_int [(b*coth(c + d*x))^(2/3),x, 12, (sqrt(3)*b^(2/3)*atan((1 - (2*(b*coth(c + d*x))^(1/3))/b^(1/3))/sqrt(3)))/(2*d) - (sqrt(3)*b^(2/3)*atan((1 + (2*(b*coth(c + d*x))^(1/3))/b^(1/3))/sqrt(3)))/(2*d) + (b^(2/3)*atanh((b*coth(c + d*x))^(1/3)/b^(1/3)))/d - (b^(2/3)*log(b^(2/3) - b^(1/3)*(b*coth(c + d*x))^(1/3) + (b*coth(c + d*x))^(2/3)))/(4*d) + (b^(2/3)*log(b^(2/3) + b^(1/3)*(b*coth(c + d*x))^(1/3) + (b*coth(c + d*x))^(2/3)))/(4*d)]
    @test_int [(b*coth(c + d*x))^(1/3), x, 9, -((sqrt(3)*b^(1/3)*atan((b^(2/3) + 2*(b*coth(c + d*x))^(2/3))/(sqrt(3)*b^(2/3))))/(2*d)) - (b^(1/3)*log(b^(2/3) - (b*coth(c + d*x))^(2/3)))/(2*d) + (b^(1/3)*log(b^(4/3) + b^(2/3)*(b*coth(c + d*x))^(2/3) + (b*coth(c + d*x))^(4/3)))/(4*d)]
    @test_int [1/(b*coth(c + d*x))^(1/3), x, 9, (sqrt(3)*atan((b^(2/3) + 2*(b*coth(c + d*x))^(2/3))/(sqrt(3)*b^(2/3))))/(2*b^(1/3)*d) - log(b^(2/3) - (b*coth(c + d*x))^(2/3))/(2*b^(1/3)*d) + log(b^(4/3) + b^(2/3)*(b*coth(c + d*x))^(2/3) + (b*coth(c + d*x))^(4/3))/(4*b^(1/3)*d)]
    @test_int [1/(b*coth(c + d*x))^(2/3), x, 12, -((sqrt(3)*atan((1 - (2*(b*coth(c + d*x))^(1/3))/b^(1/3))/sqrt(3)))/(2*b^(2/3)*d)) + (sqrt(3)*atan((1 + (2*(b*coth(c + d*x))^(1/3))/b^(1/3))/sqrt(3)))/(2*b^(2/3)*d) + atanh((b*coth(c + d*x))^(1/3)/b^(1/3))/(b^(2/3)*d) - log(b^(2/3) - b^(1/3)*(b*coth(c + d*x))^(1/3) + (b*coth(c + d*x))^(2/3))/(4*b^(2/3)*d) + log(b^(2/3) + b^(1/3)*(b*coth(c + d*x))^(1/3) + (b*coth(c + d*x))^(2/3))/(4*b^(2/3)*d)]
    @test_int [1/(b*coth(c + d*x))^(4/3), x, 13, (sqrt(3)*atan((1 - (2*(b*coth(c + d*x))^(1/3))/b^(1/3))/sqrt(3)))/(2*b^(4/3)*d) - (sqrt(3)*atan((1 + (2*(b*coth(c + d*x))^(1/3))/b^(1/3))/sqrt(3)))/(2*b^(4/3)*d) + atanh((b*coth(c + d*x))^(1/3)/b^(1/3))/(b^(4/3)*d) - 3/(b*d*(b*coth(c + d*x))^(1/3)) - log(b^(2/3) - b^(1/3)*(b*coth(c + d*x))^(1/3) + (b*coth(c + d*x))^(2/3))/(4*b^(4/3)*d) + log(b^(2/3) + b^(1/3)*(b*coth(c + d*x))^(1/3) + (b*coth(c + d*x))^(2/3))/(4*b^(4/3)*d)]

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

    @test_int [coth(a + b*x)^n, x, 2, (coth(a + b*x)^(1 + n)*HypergeometricFunctions._₂F₁(1, (1 + n)/2, (3 + n)/2, coth(a + b*x)^2))/(b*(1 + n))]
    @test_int [(b*coth(c + d*x))^n,x, 2, ((b*coth(c + d*x))^(1 + n)*HypergeometricFunctions._₂F₁(1, (1 + n)/2, (3 + n)/2, coth(c + d*x)^2))/(b*d*(1 + n))]

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

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

    @test_int [(b*coth(c + d*x)^2)^n,x, 3, (coth(c + d*x)*(b*coth(c + d*x)^2)^n*HypergeometricFunctions._₂F₁(1, (1/2)*(1 + 2*n), (1/2)*(3 + 2*n), coth(c + d*x)^2))/(d*(1 + 2*n))]

    @test_int [(b*coth(c + d*x)^2)^(3/2),x, 3, -((b*coth(c + d*x)*sqrt(b*coth(c + d*x)^2))/(2*d)) + (b*sqrt(b*coth(c + d*x)^2)*log(sinh(c + d*x))*tanh(c + d*x))/d]
    @test_int [(b*coth(c + d*x)^2)^(1/2), x, 2, (sqrt(b*coth(c + d*x)^2)*log(sinh(c + d*x))*tanh(c + d*x))/d]
    @test_int [1/(b*coth(c + d*x)^2)^(1/2), x, 2, (coth(c + d*x)*log(cosh(c + d*x)))/(d*sqrt(b*coth(c + d*x)^2))]
    @test_int [1/(b*coth(c + d*x)^2)^(3/2), x, 3, (coth(c + d*x)*log(cosh(c + d*x)))/(b*d*sqrt(b*coth(c + d*x)^2)) - tanh(c + d*x)/(2*b*d*sqrt(b*coth(c + d*x)^2))]

    @test_int [(b*coth(c + d*x)^2)^(4/3),x, 14, (sqrt(3)*b*atan((1 - 2*coth(c + d*x)^(1/3))/sqrt(3))*(b*coth(c + d*x)^2)^(1/3))/(2*d*coth(c + d*x)^(2/3)) - (sqrt(3)*b*atan((1 + 2*coth(c + d*x)^(1/3))/sqrt(3))*(b*coth(c + d*x)^2)^(1/3))/(2*d*coth(c + d*x)^(2/3)) + (b*atanh(coth(c + d*x)^(1/3))*(b*coth(c + d*x)^2)^(1/3))/(d*coth(c + d*x)^(2/3)) - (3*b*coth(c + d*x)*(b*coth(c + d*x)^2)^(1/3))/(5*d) - (b*(b*coth(c + d*x)^2)^(1/3)*log(1 - coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*coth(c + d*x)^(2/3)) + (b*(b*coth(c + d*x)^2)^(1/3)*log(1 + coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*coth(c + d*x)^(2/3))]
    @test_int [(b*coth(c + d*x)^2)^(2/3),x, 14, -((sqrt(3)*atan((1 - 2*coth(c + d*x)^(1/3))/sqrt(3))*(b*coth(c + d*x)^2)^(2/3))/(2*d*coth(c + d*x)^(4/3))) + (sqrt(3)*atan((1 + 2*coth(c + d*x)^(1/3))/sqrt(3))*(b*coth(c + d*x)^2)^(2/3))/(2*d*coth(c + d*x)^(4/3)) + (atanh(coth(c + d*x)^(1/3))*(b*coth(c + d*x)^2)^(2/3))/(d*coth(c + d*x)^(4/3)) - ((b*coth(c + d*x)^2)^(2/3)*log(1 - coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*coth(c + d*x)^(4/3)) + ((b*coth(c + d*x)^2)^(2/3)*log(1 + coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*coth(c + d*x)^(4/3)) - (3*(b*coth(c + d*x)^2)^(2/3)*tanh(c + d*x))/d]
    @test_int [(b*coth(c + d*x)^2)^(1/3), x, 13, (sqrt(3)*atan((1 - 2*coth(c + d*x)^(1/3))/sqrt(3))*(b*coth(c + d*x)^2)^(1/3))/(2*d*coth(c + d*x)^(2/3)) - (sqrt(3)*atan((1 + 2*coth(c + d*x)^(1/3))/sqrt(3))*(b*coth(c + d*x)^2)^(1/3))/(2*d*coth(c + d*x)^(2/3)) + (atanh(coth(c + d*x)^(1/3))*(b*coth(c + d*x)^2)^(1/3))/(d*coth(c + d*x)^(2/3)) - ((b*coth(c + d*x)^2)^(1/3)*log(1 - coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*coth(c + d*x)^(2/3)) + ((b*coth(c + d*x)^2)^(1/3)*log(1 + coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*coth(c + d*x)^(2/3))]
    @test_int [1/(b*coth(c + d*x)^2)^(1/3), x, 13, -((sqrt(3)*atan((1 - 2*coth(c + d*x)^(1/3))/sqrt(3))*coth(c + d*x)^(2/3))/(2*d*(b*coth(c + d*x)^2)^(1/3))) + (sqrt(3)*atan((1 + 2*coth(c + d*x)^(1/3))/sqrt(3))*coth(c + d*x)^(2/3))/(2*d*(b*coth(c + d*x)^2)^(1/3)) + (atanh(coth(c + d*x)^(1/3))*coth(c + d*x)^(2/3))/(d*(b*coth(c + d*x)^2)^(1/3)) - (coth(c + d*x)^(2/3)*log(1 - coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*(b*coth(c + d*x)^2)^(1/3)) + (coth(c + d*x)^(2/3)*log(1 + coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*(b*coth(c + d*x)^2)^(1/3))]
    @test_int [1/(b*coth(c + d*x)^2)^(2/3), x, 14, -((3*coth(c + d*x))/(d*(b*coth(c + d*x)^2)^(2/3))) + (sqrt(3)*atan((1 - 2*coth(c + d*x)^(1/3))/sqrt(3))*coth(c + d*x)^(4/3))/(2*d*(b*coth(c + d*x)^2)^(2/3)) - (sqrt(3)*atan((1 + 2*coth(c + d*x)^(1/3))/sqrt(3))*coth(c + d*x)^(4/3))/(2*d*(b*coth(c + d*x)^2)^(2/3)) + (atanh(coth(c + d*x)^(1/3))*coth(c + d*x)^(4/3))/(d*(b*coth(c + d*x)^2)^(2/3)) - (coth(c + d*x)^(4/3)*log(1 - coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*(b*coth(c + d*x)^2)^(2/3)) + (coth(c + d*x)^(4/3)*log(1 + coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*(b*coth(c + d*x)^2)^(2/3))]
    @test_int [1/(b*coth(c + d*x)^2)^(4/3), x, 14, -((sqrt(3)*atan((1 - 2*coth(c + d*x)^(1/3))/sqrt(3))*coth(c + d*x)^(2/3))/(2*b*d*(b*coth(c + d*x)^2)^(1/3))) + (sqrt(3)*atan((1 + 2*coth(c + d*x)^(1/3))/sqrt(3))*coth(c + d*x)^(2/3))/(2*b*d*(b*coth(c + d*x)^2)^(1/3)) + (atanh(coth(c + d*x)^(1/3))*coth(c + d*x)^(2/3))/(b*d*(b*coth(c + d*x)^2)^(1/3)) - (coth(c + d*x)^(2/3)*log(1 - coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*b*d*(b*coth(c + d*x)^2)^(1/3)) + (coth(c + d*x)^(2/3)*log(1 + coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*b*d*(b*coth(c + d*x)^2)^(1/3)) - (3*tanh(c + d*x))/(5*b*d*(b*coth(c + d*x)^2)^(1/3))]

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

    @test_int [(b*coth(c + d*x)^3)^n,x, 3, (coth(c + d*x)*(b*coth(c + d*x)^3)^n*HypergeometricFunctions._₂F₁(1, (1/2)*(1 + 3*n), (3*(1 + n))/2, coth(c + d*x)^2))/(d*(1 + 3*n))]

    @test_int [(b*coth(c + d*x)^3)^(3/2),x, 8, -((2*b*sqrt(b*coth(c + d*x)^3))/(3*d)) - (b*atan(sqrt(coth(c + d*x)))*sqrt(b*coth(c + d*x)^3))/(d*coth(c + d*x)^(3/2)) + (b*atanh(sqrt(coth(c + d*x)))*sqrt(b*coth(c + d*x)^3))/(d*coth(c + d*x)^(3/2)) - (2*b*coth(c + d*x)^2*sqrt(b*coth(c + d*x)^3))/(7*d)]
    @test_int [(b*coth(c + d*x)^3)^(1/2), x, 7, (atan(sqrt(coth(c + d*x)))*sqrt(b*coth(c + d*x)^3))/(d*coth(c + d*x)^(3/2)) + (atanh(sqrt(coth(c + d*x)))*sqrt(b*coth(c + d*x)^3))/(d*coth(c + d*x)^(3/2)) - (2*sqrt(b*coth(c + d*x)^3)*tanh(c + d*x))/d]
    @test_int [1/(b*coth(c + d*x)^3)^(1/2), x, 7, -((2*coth(c + d*x))/(d*sqrt(b*coth(c + d*x)^3))) - (atan(sqrt(coth(c + d*x)))*coth(c + d*x)^(3/2))/(d*sqrt(b*coth(c + d*x)^3)) + (atanh(sqrt(coth(c + d*x)))*coth(c + d*x)^(3/2))/(d*sqrt(b*coth(c + d*x)^3))]
    @test_int [1/(b*coth(c + d*x)^3)^(3/2), x, 8, -(2/(3*b*d*sqrt(b*coth(c + d*x)^3))) + (atan(sqrt(coth(c + d*x)))*coth(c + d*x)^(3/2))/(b*d*sqrt(b*coth(c + d*x)^3)) + (atanh(sqrt(coth(c + d*x)))*coth(c + d*x)^(3/2))/(b*d*sqrt(b*coth(c + d*x)^3)) - (2*tanh(c + d*x)^2)/(7*b*d*sqrt(b*coth(c + d*x)^3))]

    @test_int [(b*coth(c + d*x)^3)^(4/3),x, 4, -((b*(b*coth(c + d*x)^3)^(1/3))/d) - (b*coth(c + d*x)^2*(b*coth(c + d*x)^3)^(1/3))/(3*d) + b*x*(b*coth(c + d*x)^3)^(1/3)*tanh(c + d*x)]
    @test_int [(b*coth(c + d*x)^3)^(2/3),x, 3, -(((b*coth(c + d*x)^3)^(2/3)*tanh(c + d*x))/d) + x*(b*coth(c + d*x)^3)^(2/3)*tanh(c + d*x)^2]
    @test_int [(b*coth(c + d*x)^3)^(1/3), x, 2, ((b*coth(c + d*x)^3)^(1/3)*log(sinh(c + d*x))*tanh(c + d*x))/d]
    @test_int [1/(b*coth(c + d*x)^3)^(1/3), x, 2, (coth(c + d*x)*log(cosh(c + d*x)))/(d*(b*coth(c + d*x)^3)^(1/3))]
    @test_int [1/(b*coth(c + d*x)^3)^(2/3), x, 3, -(coth(c + d*x)/(d*(b*coth(c + d*x)^3)^(2/3))) + (x*coth(c + d*x)^2)/(b*coth(c + d*x)^3)^(2/3)]
    @test_int [1/(b*coth(c + d*x)^3)^(4/3), x, 4, -(1/(b*d*(b*coth(c + d*x)^3)^(1/3))) + (x*coth(c + d*x))/(b*(b*coth(c + d*x)^3)^(1/3)) - tanh(c + d*x)^2/(3*b*d*(b*coth(c + d*x)^3)^(1/3))]

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

    @test_int [(b*coth(c + d*x)^4)^n,x, 3, (coth(c + d*x)*(b*coth(c + d*x)^4)^n*HypergeometricFunctions._₂F₁(1, (1/2)*(1 + 4*n), (1/2)*(3 + 4*n), coth(c + d*x)^2))/(d*(1 + 4*n))]

    @test_int [(b*coth(c + d*x)^4)^(3/2),x, 5, -((b*coth(c + d*x)*sqrt(b*coth(c + d*x)^4))/(3*d)) - (b*coth(c + d*x)^3*sqrt(b*coth(c + d*x)^4))/(5*d) - (b*sqrt(b*coth(c + d*x)^4)*tanh(c + d*x))/d + b*x*sqrt(b*coth(c + d*x)^4)*tanh(c + d*x)^2]
    @test_int [(b*coth(c + d*x)^4)^(1/2), x, 3, -((sqrt(b*coth(c + d*x)^4)*tanh(c + d*x))/d) + x*sqrt(b*coth(c + d*x)^4)*tanh(c + d*x)^2]
    @test_int [1/(b*coth(c + d*x)^4)^(1/2), x, 3, -(coth(c + d*x)/(d*sqrt(b*coth(c + d*x)^4))) + (x*coth(c + d*x)^2)/sqrt(b*coth(c + d*x)^4)]
    @test_int [1/(b*coth(c + d*x)^4)^(3/2), x, 5, -(coth(c + d*x)/(b*d*sqrt(b*coth(c + d*x)^4))) + (x*coth(c + d*x)^2)/(b*sqrt(b*coth(c + d*x)^4)) - tanh(c + d*x)/(3*b*d*sqrt(b*coth(c + d*x)^4)) - tanh(c + d*x)^3/(5*b*d*sqrt(b*coth(c + d*x)^4))]

    @test_int [(b*coth(c + d*x)^4)^(4/3),x, 16, -((sqrt(3)*b*atan((1 - 2*coth(c + d*x)^(1/3))/sqrt(3))*(b*coth(c + d*x)^4)^(1/3))/(2*d*coth(c + d*x)^(4/3))) + (sqrt(3)*b*atan((1 + 2*coth(c + d*x)^(1/3))/sqrt(3))*(b*coth(c + d*x)^4)^(1/3))/(2*d*coth(c + d*x)^(4/3)) + (b*atanh(coth(c + d*x)^(1/3))*(b*coth(c + d*x)^4)^(1/3))/(d*coth(c + d*x)^(4/3)) - (3*b*coth(c + d*x)*(b*coth(c + d*x)^4)^(1/3))/(7*d) - (3*b*coth(c + d*x)^3*(b*coth(c + d*x)^4)^(1/3))/(13*d) - (b*(b*coth(c + d*x)^4)^(1/3)*log(1 - coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*coth(c + d*x)^(4/3)) + (b*(b*coth(c + d*x)^4)^(1/3)*log(1 + coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*coth(c + d*x)^(4/3)) - (3*b*(b*coth(c + d*x)^4)^(1/3)*tanh(c + d*x))/d]
    @test_int [(b*coth(c + d*x)^4)^(2/3),x, 14, (sqrt(3)*atan((1 - 2*coth(c + d*x)^(1/3))/sqrt(3))*(b*coth(c + d*x)^4)^(2/3))/(2*d*coth(c + d*x)^(8/3)) - (sqrt(3)*atan((1 + 2*coth(c + d*x)^(1/3))/sqrt(3))*(b*coth(c + d*x)^4)^(2/3))/(2*d*coth(c + d*x)^(8/3)) + (atanh(coth(c + d*x)^(1/3))*(b*coth(c + d*x)^4)^(2/3))/(d*coth(c + d*x)^(8/3)) - ((b*coth(c + d*x)^4)^(2/3)*log(1 - coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*coth(c + d*x)^(8/3)) + ((b*coth(c + d*x)^4)^(2/3)*log(1 + coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*coth(c + d*x)^(8/3)) - (3*(b*coth(c + d*x)^4)^(2/3)*tanh(c + d*x))/(5*d)]
    @test_int [(b*coth(c + d*x)^4)^(1/3), x, 14, -((sqrt(3)*atan((1 - 2*coth(c + d*x)^(1/3))/sqrt(3))*(b*coth(c + d*x)^4)^(1/3))/(2*d*coth(c + d*x)^(4/3))) + (sqrt(3)*atan((1 + 2*coth(c + d*x)^(1/3))/sqrt(3))*(b*coth(c + d*x)^4)^(1/3))/(2*d*coth(c + d*x)^(4/3)) + (atanh(coth(c + d*x)^(1/3))*(b*coth(c + d*x)^4)^(1/3))/(d*coth(c + d*x)^(4/3)) - ((b*coth(c + d*x)^4)^(1/3)*log(1 - coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*coth(c + d*x)^(4/3)) + ((b*coth(c + d*x)^4)^(1/3)*log(1 + coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*coth(c + d*x)^(4/3)) - (3*(b*coth(c + d*x)^4)^(1/3)*tanh(c + d*x))/d]
    @test_int [1/(b*coth(c + d*x)^4)^(1/3), x, 14, -((3*coth(c + d*x))/(d*(b*coth(c + d*x)^4)^(1/3))) + (sqrt(3)*atan((1 - 2*coth(c + d*x)^(1/3))/sqrt(3))*coth(c + d*x)^(4/3))/(2*d*(b*coth(c + d*x)^4)^(1/3)) - (sqrt(3)*atan((1 + 2*coth(c + d*x)^(1/3))/sqrt(3))*coth(c + d*x)^(4/3))/(2*d*(b*coth(c + d*x)^4)^(1/3)) + (atanh(coth(c + d*x)^(1/3))*coth(c + d*x)^(4/3))/(d*(b*coth(c + d*x)^4)^(1/3)) - (coth(c + d*x)^(4/3)*log(1 - coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*(b*coth(c + d*x)^4)^(1/3)) + (coth(c + d*x)^(4/3)*log(1 + coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*(b*coth(c + d*x)^4)^(1/3))]
    @test_int [1/(b*coth(c + d*x)^4)^(2/3), x, 14, -((3*coth(c + d*x))/(5*d*(b*coth(c + d*x)^4)^(2/3))) - (sqrt(3)*atan((1 - 2*coth(c + d*x)^(1/3))/sqrt(3))*coth(c + d*x)^(8/3))/(2*d*(b*coth(c + d*x)^4)^(2/3)) + (sqrt(3)*atan((1 + 2*coth(c + d*x)^(1/3))/sqrt(3))*coth(c + d*x)^(8/3))/(2*d*(b*coth(c + d*x)^4)^(2/3)) + (atanh(coth(c + d*x)^(1/3))*coth(c + d*x)^(8/3))/(d*(b*coth(c + d*x)^4)^(2/3)) - (coth(c + d*x)^(8/3)*log(1 - coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*(b*coth(c + d*x)^4)^(2/3)) + (coth(c + d*x)^(8/3)*log(1 + coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*d*(b*coth(c + d*x)^4)^(2/3))]
    @test_int [1/(b*coth(c + d*x)^4)^(4/3), x, 16, -((3*coth(c + d*x))/(b*d*(b*coth(c + d*x)^4)^(1/3))) + (sqrt(3)*atan((1 - 2*coth(c + d*x)^(1/3))/sqrt(3))*coth(c + d*x)^(4/3))/(2*b*d*(b*coth(c + d*x)^4)^(1/3)) - (sqrt(3)*atan((1 + 2*coth(c + d*x)^(1/3))/sqrt(3))*coth(c + d*x)^(4/3))/(2*b*d*(b*coth(c + d*x)^4)^(1/3)) + (atanh(coth(c + d*x)^(1/3))*coth(c + d*x)^(4/3))/(b*d*(b*coth(c + d*x)^4)^(1/3)) - (coth(c + d*x)^(4/3)*log(1 - coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*b*d*(b*coth(c + d*x)^4)^(1/3)) + (coth(c + d*x)^(4/3)*log(1 + coth(c + d*x)^(1/3) + coth(c + d*x)^(2/3)))/(4*b*d*(b*coth(c + d*x)^4)^(1/3)) - (3*tanh(c + d*x))/(7*b*d*(b*coth(c + d*x)^4)^(1/3)) - (3*tanh(c + d*x)^3)/(13*b*d*(b*coth(c + d*x)^4)^(1/3))]

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

    @test_int [(b*coth(c + d*x)^m)^n,x, 3, (coth(c + d*x)*(b*coth(c + d*x)^m)^n*HypergeometricFunctions._₂F₁(1, (1/2)*(1 + m*n), (1/2)*(3 + m*n), coth(c + d*x)^2))/(d*(1 + m*n))]

    @test_int [(b*coth(c + d*x)^m)^(3/2),x, 3, (2*b*coth(c + d*x)^(1 + m)*sqrt(b*coth(c + d*x)^m)*HypergeometricFunctions._₂F₁(1, (1/4)*(2 + 3*m), (3*(2 + m))/4, coth(c + d*x)^2))/(d*(2 + 3*m))]
    @test_int [(b*coth(c + d*x)^m)^(1/2), x, 3, (2*coth(c + d*x)*sqrt(b*coth(c + d*x)^m)*HypergeometricFunctions._₂F₁(1, (2 + m)/4, (6 + m)/4, coth(c + d*x)^2))/(d*(2 + m))]
    @test_int [1/(b*coth(c + d*x)^m)^(1/2), x, 3, (2*coth(c + d*x)*HypergeometricFunctions._₂F₁(1, (2 - m)/4, (6 - m)/4, coth(c + d*x)^2))/(d*(2 - m)*sqrt(b*coth(c + d*x)^m))]
    @test_int [1/(b*coth(c + d*x)^m)^(3/2), x, 3, (2*coth(c + d*x)^(1 - m)*HypergeometricFunctions._₂F₁(1, (1/4)*(2 - 3*m), (3*(2 - m))/4, coth(c + d*x)^2))/(b*d*(2 - 3*m)*sqrt(b*coth(c + d*x)^m))]

    @test_int [(b*coth(c + d*x)^m)^(4/3),x, 3, (3*b*coth(c + d*x)^(1 + m)*(b*coth(c + d*x)^m)^(1/3)*HypergeometricFunctions._₂F₁(1, (1/6)*(3 + 4*m), (1/6)*(9 + 4*m), coth(c + d*x)^2))/(d*(3 + 4*m))]
    @test_int [(b*coth(c + d*x)^m)^(2/3),x, 3, (3*coth(c + d*x)*(b*coth(c + d*x)^m)^(2/3)*HypergeometricFunctions._₂F₁(1, (1/6)*(3 + 2*m), (1/6)*(9 + 2*m), coth(c + d*x)^2))/(d*(3 + 2*m))]
    @test_int [(b*coth(c + d*x)^m)^(1/3), x, 3, (3*coth(c + d*x)*(b*coth(c + d*x)^m)^(1/3)*HypergeometricFunctions._₂F₁(1, (3 + m)/6, (9 + m)/6, coth(c + d*x)^2))/(d*(3 + m))]
    @test_int [1/(b*coth(c + d*x)^m)^(1/3), x, 3, (3*coth(c + d*x)*HypergeometricFunctions._₂F₁(1, (3 - m)/6, (9 - m)/6, coth(c + d*x)^2))/(d*(3 - m)*(b*coth(c + d*x)^m)^(1/3))]
    @test_int [1/(b*coth(c + d*x)^m)^(2/3), x, 3, (3*coth(c + d*x)*HypergeometricFunctions._₂F₁(1, (1/6)*(3 - 2*m), (1/6)*(9 - 2*m), coth(c + d*x)^2))/(d*(3 - 2*m)*(b*coth(c + d*x)^m)^(2/3))]
    @test_int [1/(b*coth(c + d*x)^m)^(4/3), x, 3, (3*coth(c + d*x)^(1 - m)*HypergeometricFunctions._₂F₁(1, (1/6)*(3 - 4*m), (1/6)*(9 - 4*m), coth(c + d*x)^2))/(b*d*(3 - 4*m)*(b*coth(c + d*x)^m)^(1/3))]

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

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

    #= ::Subsubsection::Closed:: =#
    #=a^2-b^2 = 0=#

    @test_int [(1 + coth(x))^5, x, 5, 16*x - 8*coth(x) - 2*(1 + coth(x))^2 - (2/3)*(1 + coth(x))^3 - (1/4)*(1 + coth(x))^4 + 16*log(sinh(x))]
    @test_int [(1 + coth(x))^4, x, 4, 8*x - 4*coth(x) - (1 + coth(x))^2 - (1/3)*(1 + coth(x))^3 + 8*log(sinh(x))]
    @test_int [(1 + coth(x))^3, x, 3, 4*x - 2*coth(x) - (1/2)*(1 + coth(x))^2 + 4*log(sinh(x))]
    @test_int [(1 + coth(x))^2, x, 2, 2*x - coth(x) + 2*log(sinh(x))]
    @test_int [1/(1 + coth(x)), x, 2, x/2 - 1/(2*(1 + coth(x)))]
    @test_int [1/(1 + coth(x))^2, x, 3, x/4 - 1/(4*(1 + coth(x))^2) - 1/(4*(1 + coth(x)))]
    @test_int [1/(1 + coth(x))^3, x, 4, x/8 - 1/(6*(1 + coth(x))^3) - 1/(8*(1 + coth(x))^2) - 1/(8*(1 + coth(x)))]
    @test_int [1/(1 + coth(x))^4, x, 5, x/16 - 1/(8*(1 + coth(x))^4) - 1/(12*(1 + coth(x))^3) - 1/(16*(1 + coth(x))^2) - 1/(16*(1 + coth(x)))]
    @test_int [1/(1 + coth(x))^5, x, 6, x/32 - 1/(10*(1 + coth(x))^5) - 1/(16*(1 + coth(x))^4) - 1/(24*(1 + coth(x))^3) - 1/(32*(1 + coth(x))^2) - 1/(32*(1 + coth(x)))]

    @test_int [(1 + coth(x))^(7/2), x, 5, 8*sqrt(2)*atanh(sqrt(1 + coth(x))/sqrt(2)) - 8*sqrt(1 + coth(x)) - (4/3)*(1 + coth(x))^(3/2) - (2/5)*(1 + coth(x))^(5/2)]
    @test_int [(1 + coth(x))^(5/2), x, 4, 4*sqrt(2)*atanh(sqrt(1 + coth(x))/sqrt(2)) - 4*sqrt(1 + coth(x)) - (2/3)*(1 + coth(x))^(3/2)]
    @test_int [(1 + coth(x))^(3/2), x, 3, 2*sqrt(2)*atanh(sqrt(1 + coth(x))/sqrt(2)) - 2*sqrt(1 + coth(x))]
    @test_int [(1 + coth(x))^(1/2), x, 2, sqrt(2)*atanh(sqrt(1 + coth(x))/sqrt(2))]
    @test_int [1/(1 + coth(x))^(1/2), x, 3, atanh(sqrt(1 + coth(x))/sqrt(2))/sqrt(2) - 1/sqrt(1 + coth(x))]
    @test_int [1/(1 + coth(x))^(3/2), x, 4, atanh(sqrt(1 + coth(x))/sqrt(2))/(2*sqrt(2)) - 1/(3*(1 + coth(x))^(3/2)) - 1/(2*sqrt(1 + coth(x)))]
    @test_int [1/(1 + coth(x))^(5/2), x, 5, atanh(sqrt(1 + coth(x))/sqrt(2))/(4*sqrt(2)) - 1/(5*(1 + coth(x))^(5/2)) - 1/(6*(1 + coth(x))^(3/2)) - 1/(4*sqrt(1 + coth(x)))]

    #= ::Subsubsection::Closed:: =#
    #=a^2-b^2 /= 0=#

    @test_int [(a + b*coth(c + d*x))^5, x, 5, a*(a^4 + 10*a^2*b^2 + 5*b^4)*x - (4*a*b^2*(a^2 + b^2)*coth(c + d*x))/d - (b*(3*a^2 + b^2)*(a + b*coth(c + d*x))^2)/(2*d) - (2*a*b*(a + b*coth(c + d*x))^3)/(3*d) - (b*(a + b*coth(c + d*x))^4)/(4*d) + (b*(5*a^4 + 10*a^2*b^2 + b^4)*log(sinh(c + d*x)))/d]
    @test_int [(a + b*coth(c + d*x))^4, x, 4, (a^4 + 6*a^2*b^2 + b^4)*x - (b^2*(3*a^2 + b^2)*coth(c + d*x))/d - (a*b*(a + b*coth(c + d*x))^2)/d - (b*(a + b*coth(c + d*x))^3)/(3*d) + (4*a*b*(a^2 + b^2)*log(sinh(c + d*x)))/d]
    @test_int [(a + b*coth(c + d*x))^3, x, 3, a*(a^2 + 3*b^2)*x - (2*a*b^2*coth(c + d*x))/d - (b*(a + b*coth(c + d*x))^2)/(2*d) + (b*(3*a^2 + b^2)*log(sinh(c + d*x)))/d]
    @test_int [(a + b*coth(c + d*x))^2, x, 2, (a^2 + b^2)*x - (b^2*coth(c + d*x))/d + (2*a*b*log(sinh(c + d*x)))/d]
    @test_int [1/(a + b*coth(c + d*x)),x, 2, (a*x)/(a^2 - b^2) - (b*log(b*cosh(c + d*x) + a*sinh(c + d*x)))/((a^2 - b^2)*d)]
    @test_int [1/(a + b*coth(c + d*x))^2,x, 3, ((a^2 + b^2)*x)/(a^2 - b^2)^2 + b/((a^2 - b^2)*d*(a + b*coth(c + d*x))) - (2*a*b*log(b*cosh(c + d*x) + a*sinh(c + d*x)))/((a^2 - b^2)^2*d)]
    @test_int [1/(a + b*coth(c + d*x))^3,x, 4, (a*(a^2 + 3*b^2)*x)/(a^2 - b^2)^3 + b/(2*(a^2 - b^2)*d*(a + b*coth(c + d*x))^2) + (2*a*b)/((a^2 - b^2)^2*d*(a + b*coth(c + d*x))) - (b*(3*a^2 + b^2)*log(b*cosh(c + d*x) + a*sinh(c + d*x)))/((a^2 - b^2)^3*d)]
    @test_int [1/(a + b*coth(c + d*x))^4,x, 5, ((a^4 + 6*a^2*b^2 + b^4)*x)/(a^2 - b^2)^4 + b/(3*(a^2 - b^2)*d*(a + b*coth(c + d*x))^3) + (a*b)/((a^2 - b^2)^2*d*(a + b*coth(c + d*x))^2) + (b*(3*a^2 + b^2))/((a^2 - b^2)^3*d*(a + b*coth(c + d*x))) - (4*a*b*(a^2 + b^2)*log(b*cosh(c + d*x) + a*sinh(c + d*x)))/((a^2 - b^2)^4*d)]

    @test_int [1/(4 + 6*coth(c + d*x)), x, 2, -(x/5) + (3*log(3*cosh(c + d*x) + 2*sinh(c + d*x)))/(10*d)]
    @test_int [1/(4 - 6*coth(c + d*x)), x, 2, -(x/5) - (3*log(3*cosh(c + d*x) - 2*sinh(c + d*x)))/(10*d)]

    @test_int [sqrt(a + b*coth(c + d*x)), x, 5, -((sqrt(a - b)*atanh(sqrt(a + b*coth(c + d*x))/sqrt(a - b)))/d) + (sqrt(a + b)*atanh(sqrt(a + b*coth(c + d*x))/sqrt(a + b)))/d]
    @test_int [1/sqrt(a + b*coth(c + d*x)), x, 5, -(atanh(sqrt(a + b*coth(c + d*x))/sqrt(a - b))/(sqrt(a - b)*d)) + atanh(sqrt(a + b*coth(c + d*x))/sqrt(a + b))/(sqrt(a + b)*d)]

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

    #= ::Subsubsection::Closed:: =#
    #=a^2-b^2 = 0=#

    @test_int [sinh(x)^4/(1 + coth(x)), x, 4, (5*x)/16 + 1/(32*(1 - coth(x))^2) + 1/(8*(1 - coth(x))) - 1/(24*(1 + coth(x))^3) - 3/(32*(1 + coth(x))^2) - 3/(16*(1 + coth(x)))]
    @test_int [sinh(x)^3/(1 + coth(x)), x, 3, -((4*cosh(x))/5) + (4*cosh(x)^3)/15 - sinh(x)^3/(5*(1 + coth(x)))]
    @test_int [sinh(x)^2/(1 + coth(x)), x, 4, -((3*x)/8) - 1/(8*(1 - coth(x))) + 1/(8*(1 + coth(x))^2) + 1/(4*(1 + coth(x)))]
    @test_int [sinh(x)^1/(1 + coth(x)), x, 2, (2*cosh(x))/3 - sinh(x)/(3*(1 + coth(x)))]
    @test_int [csch(x)^1/(1 + coth(x)), x, 1, -(csch(x)/(1 + coth(x)))]
    @test_int [csch(x)^2/(1 + coth(x)), x, 2, -log(1 + coth(x))]
    @test_int [csch(x)^3/(1 + coth(x)), x, 2, atanh(cosh(x)) - csch(x)]
    @test_int [csch(x)^4/(1 + coth(x)), x, 2, coth(x) - coth(x)^2/2]

    #= ::Subsubsection::Closed:: =#
    #=a^2-b^2 /= 0=#

    @test_int [sinh(x)^4/(a + b*coth(x)), x, 5, -(((3*a^2 + 9*a*b + 8*b^2)*log(1 - coth(x)))/(16*(a + b)^3)) + ((3*a^2 - 9*a*b + 8*b^2)*log(1 + coth(x)))/(16*(a - b)^3) - (b^5*log(a + b*coth(x)))/(a^2 - b^2)^3 - ((4*b^3 - a*(7 - (3*a^2)/b^2)*b^2*coth(x))*sinh(x)^2)/(8*(a^2 - b^2)^2) - ((b - a*coth(x))*sinh(x)^4)/(4*(a^2 - b^2))]
    @test_int [sinh(x)^3/(a + b*coth(x)), x, 9, -((b^4*atanh(((b + a*coth(x))*sinh(x))/sqrt(a^2 - b^2)))/(a^2 - b^2)^(5/2)) + (a*b^2*cosh(x))/(a^2 - b^2)^2 - (a*cosh(x))/(a^2 - b^2) + (a*cosh(x)^3)/(3*(a^2 - b^2)) - (b^3*sinh(x))/(a^2 - b^2)^2 - (b*sinh(x)^3)/(3*(a^2 - b^2))]
    @test_int [sinh(x)^2/(a + b*coth(x)), x, 4, ((a + 2*b)*log(1 - coth(x)))/(4*(a + b)^2) - ((a - 2*b)*log(1 + coth(x)))/(4*(a - b)^2) - (b^3*log(a + b*coth(x)))/(a^2 - b^2)^2 - ((b - a*coth(x))*sinh(x)^2)/(2*(a^2 - b^2))]
    @test_int [sinh(x)^1/(a + b*coth(x)), x, 5, -((b^2*atanh(((b + a*coth(x))*sinh(x))/sqrt(a^2 - b^2)))/(a^2 - b^2)^(3/2)) + (a*cosh(x))/(a^2 - b^2) - (b*sinh(x))/(a^2 - b^2)]
    @test_int [csch(x)^1/(a + b*coth(x)), x, 2, -(atanh(((b + a*coth(x))*sinh(x))/sqrt(a^2 - b^2))/sqrt(a^2 - b^2))]
    @test_int [csch(x)^2/(a + b*coth(x)), x, 2, -(log(a + b*coth(x))/b)]
    @test_int [csch(x)^3/(a + b*coth(x)), x, 5, (a*atanh(cosh(x)))/b^2 - (sqrt(a^2 - b^2)*atanh(((b + a*coth(x))*sinh(x))/sqrt(a^2 - b^2)))/b^2 - csch(x)/b]
    @test_int [csch(x)^4/(a + b*coth(x)), x, 3, (a*coth(x))/b^2 - coth(x)^2/(2*b) - ((a^2 - b^2)*log(a + b*coth(x)))/b^3]

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

    #= ::Subsubsection::Closed:: =#
    #=a^2-b^2 = 0=#

    @test_int [cosh(x)^4/(1 + coth(x)), x, 5, x/16 + 1/(32*(1 - coth(x))^2) - 1/(8*(1 - coth(x))) - 1/(24*(1 + coth(x))^3) + 5/(32*(1 + coth(x))^2) - 3/(16*(1 + coth(x)))]
    @test_int [cosh(x)^3/(1 + coth(x)), x, 9, cosh(x)^5/5 - sinh(x)^3/3 - sinh(x)^5/5]
    @test_int [cosh(x)^2/(1 + coth(x)), x, 5, x/8 - 1/(8*(1 - coth(x))) + 1/(8*(1 + coth(x))^2) - 1/(4*(1 + coth(x)))]
    @test_int [cosh(x)^1/(1 + coth(x)), x, 8, cosh(x)^3/3 - sinh(x)^3/3]
    @test_int [sech(x)^1/(1 + coth(x)), x, 8, atan(sinh(x)) + cosh(x) - sinh(x)]
    @test_int [sech(x)^2/(1 + coth(x)), x, 3, -log(1 + coth(x)) - log(tanh(x)) + tanh(x)]
    @test_int [sech(x)^3/(1 + coth(x)), x, 8, (-(1/2))*atan(sinh(x)) - sech(x) + (1/2)*sech(x)*tanh(x)]
    @test_int [sech(x)^4/(1 + coth(x)), x, 4, tanh(x)^2/2 - tanh(x)^3/3]

    @test_int [sech(x)^2*sqrt(1 + coth(x)), x, 4, atanh(sqrt(1 + coth(x))) + sqrt(1 + coth(x))*tanh(x)]

    #= ::Subsubsection::Closed:: =#
    #=a^2-b^2 /= 0=#

    @test_int [cosh(x)^4/(a + b*coth(x)), x, 5, -((a*(3*a + b)*log(1 - coth(x)))/(16*(a + b)^3)) + (a*(3*a - b)*log(1 + coth(x)))/(16*(a - b)^3) - (a^4*b*log(a + b*coth(x)))/(a^2 - b^2)^3 - ((4*b*(2*a^2 - b^2) - a*(5*a^2 - b^2)*coth(x))*sinh(x)^2)/(8*(a^2 - b^2)^2) - ((b - a*coth(x))*sinh(x)^4)/(4*(a^2 - b^2))]
    @test_int [cosh(x)^3/(a + b*coth(x)), x, 10, (a^3*b*atanh((a*cosh(x) + b*sinh(x))/sqrt(a^2 - b^2)))/(a^2 - b^2)^(5/2) - (a^2*b*cosh(x))/(a^2 - b^2)^2 - (b*cosh(x)^3)/(3*(a^2 - b^2)) + (a*b^2*sinh(x))/(a^2 - b^2)^2 + (a*sinh(x))/(a^2 - b^2) + (a*sinh(x)^3)/(3*(a^2 - b^2))]
    @test_int [cosh(x)^2/(a + b*coth(x)), x, 4, -((a*log(1 - coth(x)))/(4*(a + b)^2)) + (a*log(1 + coth(x)))/(4*(a - b)^2) - (a^2*b*log(a + b*coth(x)))/(a^2 - b^2)^2 - ((b - a*coth(x))*sinh(x)^2)/(2*(a^2 - b^2))]
    @test_int [cosh(x)^1/(a + b*coth(x)), x, 6, (a*b*atanh((a*cosh(x) + b*sinh(x))/sqrt(a^2 - b^2)))/(a^2 - b^2)^(3/2) - (b*cosh(x))/(a^2 - b^2) + (a*sinh(x))/(a^2 - b^2)]
    @test_int [sech(x)^1/(a + b*coth(x)), x, 6, atan(sinh(x))/a + (b*atanh((a*cosh(x) + b*sinh(x))/sqrt(a^2 - b^2)))/(a*sqrt(a^2 - b^2))]
    @test_int [sech(x)^2/(a + b*coth(x)), x, 3, -((b*log(a + b*coth(x)))/a^2) - (b*log(tanh(x)))/a^2 + tanh(x)/a]
    @test_int [sech(x)^3/(a + b*coth(x)), x, 9, atan(sinh(x))/(2*a) - (b^2*atan(sinh(x)))/a^3 + (b*sqrt(a^2 - b^2)*atanh((a*cosh(x) + b*sinh(x))/sqrt(a^2 - b^2)))/a^3 - (b*sech(x))/a^2 + (sech(x)*tanh(x))/(2*a)]
    @test_int [sech(x)^4/(a + b*coth(x)), x, 3, -((b*(a^2 - b^2)*log(a + b*coth(x)))/a^4) - (b*(a^2 - b^2)*log(tanh(x)))/a^4 + ((a^2 - b^2)*tanh(x))/a^3 + (b*tanh(x)^2)/(2*a^2) - tanh(x)^3/(3*a)]

    #= Following*hangs*Mathematica: =#
    @test_int [sech(x)/(im + 2*coth(x)), x, 6, (-im)*atan(sinh(x)) - (2*atanh((cosh(x) - 2*im*sinh(x))/sqrt(5)))/sqrt(5)]

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

    #= ::Subsubsection::Closed:: =#
    #=a^2-b^2 = 0=#

    @test_int [tanh(x)^4/(1 + coth(x)), x, 6, (5*x)/2 - 2*log(cosh(x)) - (5*tanh(x))/2 + tanh(x)^2 - (5*tanh(x)^3)/6 + tanh(x)^3/(2*(1 + coth(x)))]
    @test_int [tanh(x)^3/(1 + coth(x)), x, 5, -((3*x)/2) + 2*log(cosh(x)) + (3*tanh(x))/2 - tanh(x)^2 + tanh(x)^2/(2*(1 + coth(x)))]
    @test_int [tanh(x)^2/(1 + coth(x)), x, 4, (3*x)/2 - log(cosh(x)) - (3*tanh(x))/2 + tanh(x)/(2*(1 + coth(x)))]
    @test_int [tanh(x)^1/(1 + coth(x)), x, 4, -(x/2) + 1/(2*(1 + coth(x))) + log(cosh(x))]
    @test_int [tanh(x)^0/(1 + coth(x)), x, 2, x/2 - 1/(2*(1 + coth(x)))]
    @test_int [coth(x)^1/(1 + coth(x)), x, 2, x/2 + 1/(2*(1 + coth(x)))]
    @test_int [coth(x)^2/(1 + coth(x)), x, 3, -(x/2) - 1/(2*(1 + coth(x))) + log(sinh(x))]
    @test_int [coth(x)^3/(1 + coth(x)), x, 3, (3*x)/2 - (3*coth(x))/2 + coth(x)^2/(2*(1 + coth(x))) - log(sinh(x))]
    @test_int [coth(x)^4/(1 + coth(x)), x, 4, -((3*x)/2) + (3*coth(x))/2 - coth(x)^2 + coth(x)^3/(2*(1 + coth(x))) + 2*log(sinh(x))]

    @test_int [coth(x)*(1 + coth(x))^(3/2), x, 4, 2*sqrt(2)*atanh(sqrt(1 + coth(x))/sqrt(2)) - 2*sqrt(1 + coth(x)) - (2/3)*(1 + coth(x))^(3/2)]
    @test_int [coth(x)*sqrt(1 + coth(x)), x, 3, sqrt(2)*atanh(sqrt(1 + coth(x))/sqrt(2)) - 2*sqrt(1 + coth(x))]
    @test_int [coth(x)/sqrt(1 + coth(x)), x, 3, atanh(sqrt(1 + coth(x))/sqrt(2))/sqrt(2) + 1/sqrt(1 + coth(x))]
    @test_int [coth(x)/(1 + coth(x))^(3/2), x, 4, atanh(sqrt(1 + coth(x))/sqrt(2))/(2*sqrt(2)) + 1/(3*(1 + coth(x))^(3/2)) - 1/(2*sqrt(1 + coth(x)))]

    @test_int [coth(x)^2*(1 + coth(x))^(3/2), x, 4, 2*sqrt(2)*atanh(sqrt(1 + coth(x))/sqrt(2)) - 2*sqrt(1 + coth(x)) - (2/5)*(1 + coth(x))^(5/2)]
    @test_int [coth(x)^2*sqrt(1 + coth(x)), x, 3, sqrt(2)*atanh(sqrt(1 + coth(x))/sqrt(2)) - (2/3)*(1 + coth(x))^(3/2)]
    @test_int [coth(x)^2/sqrt(1 + coth(x)), x, 4, atanh(sqrt(1 + coth(x))/sqrt(2))/sqrt(2) - 1/sqrt(1 + coth(x)) - 2*sqrt(1 + coth(x))]
    @test_int [coth(x)^2/(1 + coth(x))^(3/2), x, 4, atanh(sqrt(1 + coth(x))/sqrt(2))/(2*sqrt(2)) - 1/(3*(1 + coth(x))^(3/2)) + 3/(2*sqrt(1 + coth(x)))]

    #= ::Subsubsection::Closed:: =#
    #=a^2-b^2 /= 0=#

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

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

    @test_int [x*csch(x)^2/(a + b*coth(x))^2, x, 3, -((a*x)/(b*(a^2 - b^2))) + x/(b*(a + b*coth(x))) + log(b*cosh(x) + a*sinh(x))/(a^2 - b^2)]

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*coth(a+2*log(x))^p=#

    @test_int [x^3*coth(a + 2*log(x)), x, 4, x^4/4 + log(1 - ℯ^(2*a)*x^4)/(2*ℯ^(2*a))]
    @test_int [x^2*coth(a + 2*log(x)), x, 5, x^3/3 + atan(ℯ^(a/2)*x)/ℯ^((3*a)/2) - atanh(ℯ^(a/2)*x)/ℯ^((3*a)/2)]
    @test_int [x^1*coth(a + 2*log(x)), x, 4, x^2/2 - atanh(ℯ^a*x^2)/ℯ^a]
    @test_int [x^0*coth(a + 2*log(x)), x, 5, x - atan(ℯ^(a/2)*x)/ℯ^(a/2) - atanh(ℯ^(a/2)*x)/ℯ^(a/2)]
    @test_int [coth(a + 2*log(x))/x^1, x, 2, log(sinh(a + 2*log(x)))/2]
    @test_int [coth(a + 2*log(x))/x^2, x, 5, x^(-1) + ℯ^(a/2)*atan(ℯ^(a/2)*x) - ℯ^(a/2)*atanh(ℯ^(a/2)*x)]
    @test_int [coth(a + 2*log(x))/x^3, x, 4, 1/(2*x^2) - ℯ^a*atanh(ℯ^a*x^2)]

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*coth(a+b*log(x))^p*with*m*symbolic=#

    @test_int [(e*x)^m*coth(a + 2*log(x))^1, x, 3, (e*x)^(1 + m)/(e*(1 + m)) - (2*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/4, (5 + m)/4, ℯ^(2*a)*x^4))/(e*(1 + m))]
    @test_int [(e*x)^m*coth(a + 2*log(x))^2, x, 4, (e*x)^(1 + m)/(e*(1 + m)) + (e*x)^(1 + m)/(e*(1 - ℯ^(2*a)*x^4)) - ((e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/4, (5 + m)/4, ℯ^(2*a)*x^4))/e]
    @test_int [(e*x)^m*coth(a + 2*log(x))^3, x, 5, ((3 + m)*(5 + m)*(e*x)^(1 + m))/(8*e*(1 + m)) - ((e*x)^(1 + m)*(1 + ℯ^(2*a)*x^4)^2)/(4*e*(1 - ℯ^(2*a)*x^4)^2) - ((e*x)^(1 + m)*(ℯ^(2*a)*(3 - m) - ℯ^(4*a)*(5 + m)*x^4))/(ℯ^(2*a)*(8*e*(1 - ℯ^(2*a)*x^4))) - ((9 + 2*m + m^2)*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/4, (5 + m)/4, ℯ^(2*a)*x^4))/(4*e*(1 + m))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*coth(a+b*log(x))^p*with*p*symbolic=#
    #==#

    @test_int [coth(a + b*log(x))^p, x, 3, (x*(-1 - ℯ^(2*a)*x^(2*b))^p*AppellF1(1/(2*b), p, -p, (1/2)*(2 + 1/b), ℯ^(2*a)*x^(2*b), (-ℯ^(2*a))*x^(2*b)))/(1 + ℯ^(2*a)*x^(2*b))^p]
    @test_int [(e*x)^m*coth(a + b*log(x))^p, x, 3, ((e*x)^(1 + m)*(-1 - ℯ^(2*a)*x^(2*b))^p*AppellF1((1 + m)/(2*b), p, -p, 1 + (1 + m)/(2*b), ℯ^(2*a)*x^(2*b), (-ℯ^(2*a))*x^(2*b)))/((1 + ℯ^(2*a)*x^(2*b))^p*(e*(1 + m)))]

    @test_int [coth(a + 1/2*log(x))^p, x, 2, -(((-1 - ℯ^(2*a)*x)^(1 + p)*HypergeometricFunctions._₂F₁(p, 1 + p, 2 + p, (1/2)*(1 + ℯ^(2*a)*x)))/(2^p*ℯ^(2*a)*(1 + p)))]
    @test_int [coth(a + 1/4*log(x))^p, x, 4, ((-1 - ℯ^(2*a)*sqrt(x))^(1 + p)*(1 - ℯ^(2*a)*sqrt(x))^(1 - p))/ℯ^(4*a) - (2^(1 - p)*p*(-1 - ℯ^(2*a)*sqrt(x))^(1 + p)*HypergeometricFunctions._₂F₁(p, 1 + p, 2 + p, (1/2)*(1 + ℯ^(2*a)*sqrt(x))))/(ℯ^(4*a)*(1 + p))]
    @test_int [coth(a + 1/6*log(x))^p, x, 5, (p*(-1 - ℯ^(2*a)*x^(1/3))^(1 + p)*(1 - ℯ^(2*a)*x^(1/3))^(1 - p))/ℯ^(6*a) + ((-1 - ℯ^(2*a)*x^(1/3))^(1 + p)*(1 - ℯ^(2*a)*x^(1/3))^(1 - p)*x^(1/3))/ℯ^(4*a) - ((1 + 2*p^2)*(-1 - ℯ^(2*a)*x^(1/3))^(1 + p)*HypergeometricFunctions._₂F₁(p, 1 + p, 2 + p, (1/2)*(1 + ℯ^(2*a)*x^(1/3))))/(2^p*ℯ^(6*a)*(1 + p))]
    @test_int [coth(a + 1/8*log(x))^p, x, 5, ((1/3)*(-1 - ℯ^(2*a)*x^(1/4))^(1 + p)*(1 - ℯ^(2*a)*x^(1/4))^(1 - p)*(ℯ^(4*a)*(3 + 2*p^2) + 2*ℯ^(6*a)*p*x^(1/4)))/ℯ^(12*a) + ((-1 - ℯ^(2*a)*x^(1/4))^(1 + p)*(1 - ℯ^(2*a)*x^(1/4))^(1 - p)*sqrt(x))/ℯ^(4*a) - (2^(2 - p)*p*(2 + p^2)*(-1 - ℯ^(2*a)*x^(1/4))^(1 + p)*HypergeometricFunctions._₂F₁(p, 1 + p, 2 + p, (1/2)*(1 + ℯ^(2*a)*x^(1/4))))/(ℯ^(8*a)*(3*(1 + p)))]

    @test_int [coth(a + 1*log(x))^p, x, 3, (x*(-1 - ℯ^(2*a)*x^2)^p*AppellF1(1/2, p, -p, 3/2, ℯ^(2*a)*x^2, (-ℯ^(2*a))*x^2))/(1 + ℯ^(2*a)*x^2)^p]
    @test_int [coth(a + 2*log(x))^p, x, 3, (x*(-1 - ℯ^(2*a)*x^4)^p*AppellF1(1/4, p, -p, 5/4, ℯ^(2*a)*x^4, (-ℯ^(2*a))*x^4))/(1 + ℯ^(2*a)*x^4)^p]
    @test_int [coth(a + 3*log(x))^p, x, 3, (x*(-1 - ℯ^(2*a)*x^6)^p*AppellF1(1/6, p, -p, 7/6, ℯ^(2*a)*x^6, (-ℯ^(2*a))*x^6))/(1 + ℯ^(2*a)*x^6)^p]

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

    @test_int [x^3*coth(d*(a + b*log(c*x^n))), x, 4, x^4/4 - (x^4*HypergeometricFunctions._₂F₁(1, 2/(b*d*n), 1 + 2/(b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/2]
    @test_int [x^2*coth(d*(a + b*log(c*x^n))), x, 4, x^3/3 - (2*x^3*HypergeometricFunctions._₂F₁(1, 3/(2*b*d*n), 1 + 3/(2*b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/3]
    @test_int [x^1*coth(d*(a + b*log(c*x^n))), x, 4, x^2/2 - x^2*HypergeometricFunctions._₂F₁(1, 1/(b*d*n), 1 + 1/(b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d))]
    @test_int [x^0*coth(d*(a + b*log(c*x^n))), x, 4, x - 2*x*HypergeometricFunctions._₂F₁(1, 1/(2*b*d*n), 1 + 1/(2*b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d))]
    @test_int [coth(d*(a + b*log(c*x^n)))/x^1, x, 2, log(sinh(a*d + b*d*log(c*x^n)))/(b*d*n)]
    @test_int [coth(d*(a + b*log(c*x^n)))/x^2, x, 4, -(1/x) + (2*HypergeometricFunctions._₂F₁(1, -(1/(2*b*d*n)), 1 - 1/(2*b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/x]
    @test_int [coth(d*(a + b*log(c*x^n)))/x^3, x, 4, -1/(2*x^2) + HypergeometricFunctions._₂F₁(1, -(1/(b*d*n)), 1 - 1/(b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d))/x^2]

    @test_int [x^3*coth(d*(a + b*log(c*x^n)))^2, x, 5, (1/4)*(1 + 4/(b*d*n))*x^4 + (x^4*(1 + ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/(b*d*n*(1 - ℯ^(2*a*d)*(c*x^n)^(2*b*d))) - (2*x^4*HypergeometricFunctions._₂F₁(1, 2/(b*d*n), 1 + 2/(b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/(b*d*n)]
    @test_int [x^2*coth(d*(a + b*log(c*x^n)))^2, x, 5, (1/3)*(1 + 3/(b*d*n))*x^3 + (x^3*(1 + ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/(b*d*n*(1 - ℯ^(2*a*d)*(c*x^n)^(2*b*d))) - (2*x^3*HypergeometricFunctions._₂F₁(1, 3/(2*b*d*n), 1 + 3/(2*b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/(b*d*n)]
    @test_int [x^1*coth(d*(a + b*log(c*x^n)))^2, x, 5, (1/2)*(1 + 2/(b*d*n))*x^2 + (x^2*(1 + ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/(b*d*n*(1 - ℯ^(2*a*d)*(c*x^n)^(2*b*d))) - (2*x^2*HypergeometricFunctions._₂F₁(1, 1/(b*d*n), 1 + 1/(b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/(b*d*n)]
    @test_int [x^0*coth(d*(a + b*log(c*x^n)))^2, x, 5, (1 + 1/(b*d*n))*x + (x*(1 + ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/(b*d*n*(1 - ℯ^(2*a*d)*(c*x^n)^(2*b*d))) - (2*x*HypergeometricFunctions._₂F₁(1, 1/(2*b*d*n), 1 + 1/(2*b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/(b*d*n)]
    @test_int [coth(d*(a + b*log(c*x^n)))^2/x^1, x, 3, -(coth(a*d + b*d*log(c*x^n))/(b*d*n)) + log(x)]
    @test_int [coth(d*(a + b*log(c*x^n)))^2/x^2, x, 5, -((1 - 1/(b*d*n))/x) + (1 + ℯ^(2*a*d)*(c*x^n)^(2*b*d))/(b*d*n*x*(1 - ℯ^(2*a*d)*(c*x^n)^(2*b*d))) - (2*HypergeometricFunctions._₂F₁(1, -(1/(2*b*d*n)), 1 - 1/(2*b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/(b*d*n*x)]
    @test_int [coth(d*(a + b*log(c*x^n)))^2/x^3, x, 5, (2 - b*d*n)/(2*b*d*n*x^2) + (1 + ℯ^(2*a*d)*(c*x^n)^(2*b*d))/(b*d*n*x^2*(1 - ℯ^(2*a*d)*(c*x^n)^(2*b*d))) - (2*HypergeometricFunctions._₂F₁(1, -(1/(b*d*n)), 1 - 1/(b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/(b*d*n*x^2)]

    @test_int [coth(a + b*log(c*x^n))^3/x, x, 3, -(coth(a + b*log(c*x^n))^2/(2*b*n)) + log(sinh(a + b*log(c*x^n)))/(b*n)]
    @test_int [coth(a + b*log(c*x^n))^4/x, x, 4, -(coth(a + b*log(c*x^n))/(b*n)) - coth(a + b*log(c*x^n))^3/(3*b*n) + log(x)]
    @test_int [coth(a + b*log(c*x^n))^5/x, x, 4, -(coth(a + b*log(c*x^n))^2/(2*b*n)) - coth(a + b*log(c*x^n))^4/(4*b*n) + log(sinh(a + b*log(c*x^n)))/(b*n)]

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

    @test_int [(e*x)^m*coth(d*(a + b*log(c*x^n)))^1, x, 4, (e*x)^(1 + m)/(e*(1 + m)) - (2*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/(2*b*d*n), 1 + (1 + m)/(2*b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/(e*(1 + m))]
    @test_int [(e*x)^m*coth(d*(a + b*log(c*x^n)))^2, x, 5, ((1 + m + b*d*n)*(e*x)^(1 + m))/(b*d*e*(1 + m)*n) + ((e*x)^(1 + m)*(1 + ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/(b*d*e*n*(1 - ℯ^(2*a*d)*(c*x^n)^(2*b*d))) - (2*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/(2*b*d*n), 1 + (1 + m)/(2*b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/(b*d*e*n)]
    @test_int [(e*x)^m*coth(d*(a + b*log(c*x^n)))^3, x, 6, ((1 + m + b*d*n)*(1 + m + 2*b*d*n)*(e*x)^(1 + m))/(2*b^2*d^2*e*(1 + m)*n^2) - ((e*x)^(1 + m)*(1 + ℯ^(2*a*d)*(c*x^n)^(2*b*d))^2)/(2*b*d*e*n*(1 - ℯ^(2*a*d)*(c*x^n)^(2*b*d))^2) + ((e*x)^(1 + m)*((ℯ^(2*a*d)*(1 + m - 2*b*d*n))/n + (ℯ^(4*a*d)*(1 + m + 2*b*d*n)*(c*x^n)^(2*b*d))/n))/(ℯ^(2*a*d)*(2*b^2*d^2*e*n*(1 - ℯ^(2*a*d)*(c*x^n)^(2*b*d)))) - ((1 + 2*m + m^2 + 2*b^2*d^2*n^2)*(e*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/(2*b*d*n), 1 + (1 + m)/(2*b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d)))/(b^2*d^2*e*(1 + m)*n^2)]

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

    @test_int [coth(d*(a + b*log(c*x^n)))^p, x, 4, (x*(-1 - ℯ^(2*a*d)*(c*x^n)^(2*b*d))^p*AppellF1(1/(2*b*d*n), p, -p, 1 + 1/(2*b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d), (-ℯ^(2*a*d))*(c*x^n)^(2*b*d)))/(1 + ℯ^(2*a*d)*(c*x^n)^(2*b*d))^p]
    @test_int [(e*x)^m*coth(d*(a + b*log(c*x^n)))^p, x, 4, ((e*x)^(1 + m)*(-1 - ℯ^(2*a*d)*(c*x^n)^(2*b*d))^p*AppellF1((1 + m)/(2*b*d*n), p, -p, 1 + (1 + m)/(2*b*d*n), ℯ^(2*a*d)*(c*x^n)^(2*b*d), (-ℯ^(2*a*d))*(c*x^n)^(2*b*d)))/((1 + ℯ^(2*a*d)*(c*x^n)^(2*b*d))^p*(e*(1 + m)))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*coth(a+b*log(c*x^n))^(p/2)=#

    @test_int [coth(a + b*log(c*x^n))^(5/2)/x, x, 7, -(atan(sqrt(coth(a + b*log(c*x^n))))/(b*n)) + atanh(sqrt(coth(a + b*log(c*x^n))))/(b*n) - (2*coth(a + b*log(c*x^n))^(3/2))/(3*b*n)]
    @test_int [coth(a + b*log(c*x^n))^(3/2)/x, x, 7, atan(sqrt(coth(a + b*log(c*x^n))))/(b*n) + atanh(sqrt(coth(a + b*log(c*x^n))))/(b*n) - (2*sqrt(coth(a + b*log(c*x^n))))/(b*n)]
    @test_int [sqrt(coth(a + b*log(c*x^n)))/x, x, 6, -(atan(sqrt(coth(a + b*log(c*x^n))))/(b*n)) + atanh(sqrt(coth(a + b*log(c*x^n))))/(b*n)]
    @test_int [1/(x*sqrt(coth(a + b*log(c*x^n)))), x, 6, atan(sqrt(coth(a + b*log(c*x^n))))/(b*n) + atanh(sqrt(coth(a + b*log(c*x^n))))/(b*n)]
    @test_int [1/(x*coth(a + b*log(c*x^n))^(3/2)), x, 7, -(atan(sqrt(coth(a + b*log(c*x^n))))/(b*n)) + atanh(sqrt(coth(a + b*log(c*x^n))))/(b*n) - 2/(b*n*sqrt(coth(a + b*log(c*x^n))))]
    @test_int [1/(x*coth(a + b*log(c*x^n))^(5/2)), x, 7, atan(sqrt(coth(a + b*log(c*x^n))))/(b*n) + atanh(sqrt(coth(a + b*log(c*x^n))))/(b*n) - 2/(3*b*n*coth(a + b*log(c*x^n))^(3/2))]

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

    @test_int [coth(x)^5/sqrt(a + b*coth(x)^2 + c*coth(x)^4), x, 8, ((b - 2*c)*atanh((b + 2*c*coth(x)^2)/(2*sqrt(c)*sqrt(a + b*coth(x)^2 + c*coth(x)^4))))/(4*c^(3/2)) + atanh((2*a + b + (b + 2*c)*coth(x)^2)/(2*sqrt(a + b + c)*sqrt(a + b*coth(x)^2 + c*coth(x)^4)))/(2*sqrt(a + b + c)) - sqrt(a + b*coth(x)^2 + c*coth(x)^4)/(2*c)]
    @test_int [coth(x)^3/sqrt(a + b*coth(x)^2 + c*coth(x)^4), x, 7, -(atanh((b + 2*c*coth(x)^2)/(2*sqrt(c)*sqrt(a + b*coth(x)^2 + c*coth(x)^4)))/(2*sqrt(c))) + atanh((2*a + b + (b + 2*c)*coth(x)^2)/(2*sqrt(a + b + c)*sqrt(a + b*coth(x)^2 + c*coth(x)^4)))/(2*sqrt(a + b + c))]
    @test_int [coth(x)/sqrt(a + b*coth(x)^2 + c*coth(x)^4), x, 4, atanh((2*a + b + (b + 2*c)*coth(x)^2)/(2*sqrt(a + b + c)*sqrt(a + b*coth(x)^2 + c*coth(x)^4)))/(2*sqrt(a + b + c))]
    @test_int [tanh(x)/sqrt(a + b*coth(x)^2 + c*coth(x)^4), x, 8, -(atanh((2*a + b*coth(x)^2)/(2*sqrt(a)*sqrt(a + b*coth(x)^2 + c*coth(x)^4)))/(2*sqrt(a))) + atanh((2*a + b + (b + 2*c)*coth(x)^2)/(2*sqrt(a + b + c)*sqrt(a + b*coth(x)^2 + c*coth(x)^4)))/(2*sqrt(a + b + c))]
    @test_int [tanh(x)^3/sqrt(a + b*coth(x)^2 + c*coth(x)^4), x, 11, -(atanh((2*a + b*coth(x)^2)/(2*sqrt(a)*sqrt(a + b*coth(x)^2 + c*coth(x)^4)))/(2*sqrt(a))) + (b*atanh((2*a + b*coth(x)^2)/(2*sqrt(a)*sqrt(a + b*coth(x)^2 + c*coth(x)^4))))/(4*a^(3/2)) + atanh((2*a + b + (b + 2*c)*coth(x)^2)/(2*sqrt(a + b + c)*sqrt(a + b*coth(x)^2 + c*coth(x)^4)))/(2*sqrt(a + b + c)) - (sqrt(a + b*coth(x)^2 + c*coth(x)^4)*tanh(x)^2)/(2*a)]

    #= [coth(x)^5*sqrt(a + b*coth(x)^2 + c*coth(x)^4), x, 0, 0] =#
    #= [coth(x)^3*sqrt(a + b*coth(x)^2 + c*coth(x)^4), x, 0, 0] =#
    @test_int [coth(x)*sqrt(a + b*coth(x)^2 + c*coth(x)^4), x, 8, -(((b + 2*c)*atanh((b + 2*c*coth(x)^2)/(2*sqrt(c)*sqrt(a + b*coth(x)^2 + c*coth(x)^4))))/(4*sqrt(c))) + (1/2)*sqrt(a + b + c)*atanh((2*a + b + (b + 2*c)*coth(x)^2)/(2*sqrt(a + b + c)*sqrt(a + b*coth(x)^2 + c*coth(x)^4))) - (1/2)*sqrt(a + b*coth(x)^2 + c*coth(x)^4)]
    #= [tanh(x)*sqrt(a + b*coth(x)^2 + c*coth(x)^4), x, 0, 0] =#
    #= [tanh(x)^3*sqrt(a + b*coth(x)^2 + c*coth(x)^4), x, 0, 0] =#

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*ℯ^(c*(a+b*x))*(coth(a*c+b*c*x)^2)^(m/2)=#

    @test_int [ℯ^(c*(a + b*x))*(coth(a*c + b*c*x)^2)^(5/2), x, 9, (ℯ^(c*(a + b*x))*sqrt(coth(a*c + b*c*x)^2)*tanh(a*c + b*c*x))/(b*c) - (4*ℯ^(c*(a + b*x))*sqrt(coth(a*c + b*c*x)^2)*tanh(a*c + b*c*x))/(b*c*(1 - ℯ^(2*c*(a + b*x)))^4) + (26*ℯ^(c*(a + b*x))*sqrt(coth(a*c + b*c*x)^2)*tanh(a*c + b*c*x))/(3*b*c*(1 - ℯ^(2*c*(a + b*x)))^3) - (55*ℯ^(c*(a + b*x))*sqrt(coth(a*c + b*c*x)^2)*tanh(a*c + b*c*x))/(6*b*c*(1 - ℯ^(2*c*(a + b*x)))^2) + (25*ℯ^(c*(a + b*x))*sqrt(coth(a*c + b*c*x)^2)*tanh(a*c + b*c*x))/(4*b*c*(1 - ℯ^(2*c*(a + b*x)))) - (15*atanh(ℯ^(c*(a + b*x)))*sqrt(coth(a*c + b*c*x)^2)*tanh(a*c + b*c*x))/(4*b*c)]
    @test_int [ℯ^(c*(a + b*x))*(coth(a*c + b*c*x)^2)^(3/2), x, 8, (ℯ^(c*(a + b*x))*sqrt(coth(a*c + b*c*x)^2)*tanh(a*c + b*c*x))/(b*c) - (2*ℯ^(c*(a + b*x))*sqrt(coth(a*c + b*c*x)^2)*tanh(a*c + b*c*x))/(b*c*(1 - ℯ^(2*c*(a + b*x)))^2) + (3*ℯ^(c*(a + b*x))*sqrt(coth(a*c + b*c*x)^2)*tanh(a*c + b*c*x))/(b*c*(1 - ℯ^(2*c*(a + b*x)))) - (3*atanh(ℯ^(c*(a + b*x)))*sqrt(coth(a*c + b*c*x)^2)*tanh(a*c + b*c*x))/(b*c)]
    @test_int [ℯ^(c*(a + b*x))*sqrt(coth(a*c + b*c*x)^2), x, 4, (ℯ^(c*(a + b*x))*sqrt(coth(a*c + b*c*x)^2)*tanh(a*c + b*c*x))/(b*c) - (2*atanh(ℯ^(c*(a + b*x)))*sqrt(coth(a*c + b*c*x)^2)*tanh(a*c + b*c*x))/(b*c)]
    @test_int [ℯ^(c*(a + b*x))/sqrt(coth(a*c + b*c*x)^2), x, 4, (ℯ^(c*(a + b*x))*coth(a*c + b*c*x))/(b*c*sqrt(coth(a*c + b*c*x)^2)) - (2*atan(ℯ^(c*(a + b*x)))*coth(a*c + b*c*x))/(b*c*sqrt(coth(a*c + b*c*x)^2))]
    @test_int [ℯ^(c*(a + b*x))/(coth(a*c + b*c*x)^2)^(3/2), x, 8, (ℯ^(c*(a + b*x))*coth(a*c + b*c*x))/(b*c*sqrt(coth(a*c + b*c*x)^2)) - (2*ℯ^(c*(a + b*x))*coth(a*c + b*c*x))/(b*c*(1 + ℯ^(2*c*(a + b*x)))^2*sqrt(coth(a*c + b*c*x)^2)) + (3*ℯ^(c*(a + b*x))*coth(a*c + b*c*x))/(b*c*(1 + ℯ^(2*c*(a + b*x)))*sqrt(coth(a*c + b*c*x)^2)) - (3*atan(ℯ^(c*(a + b*x)))*coth(a*c + b*c*x))/(b*c*sqrt(coth(a*c + b*c*x)^2))]
    @test_int [ℯ^(c*(a + b*x))/(coth(a*c + b*c*x)^2)^(5/2), x, 9, (ℯ^(c*(a + b*x))*coth(a*c + b*c*x))/(b*c*sqrt(coth(a*c + b*c*x)^2)) - (4*ℯ^(c*(a + b*x))*coth(a*c + b*c*x))/(b*c*(1 + ℯ^(2*c*(a + b*x)))^4*sqrt(coth(a*c + b*c*x)^2)) + (26*ℯ^(c*(a + b*x))*coth(a*c + b*c*x))/(3*b*c*(1 + ℯ^(2*c*(a + b*x)))^3*sqrt(coth(a*c + b*c*x)^2)) - (55*ℯ^(c*(a + b*x))*coth(a*c + b*c*x))/(6*b*c*(1 + ℯ^(2*c*(a + b*x)))^2*sqrt(coth(a*c + b*c*x)^2)) + (25*ℯ^(c*(a + b*x))*coth(a*c + b*c*x))/(4*b*c*(1 + ℯ^(2*c*(a + b*x)))*sqrt(coth(a*c + b*c*x)^2)) - (15*atan(ℯ^(c*(a + b*x)))*coth(a*c + b*c*x))/(4*b*c*sqrt(coth(a*c + b*c*x)^2))]

    #= ::Section::Closed:: =#
    #=Integrands*that*are*function*of*coth(c+d*x)=#

    @test_int [sin(coth(a + b*x))^3, x, 19, -((3*CosIntegral(1 - coth(a + b*x))*sin(1))/(8*b)) - (3*CosIntegral(1 + coth(a + b*x))*sin(1))/(8*b) + (CosIntegral(3 - 3*coth(a + b*x))*sin(3))/(8*b) + (CosIntegral(3 + 3*coth(a + b*x))*sin(3))/(8*b) - (cos(3)*SinIntegral(3 - 3*coth(a + b*x)))/(8*b) + (3*cos(1)*SinIntegral(1 - coth(a + b*x)))/(8*b) + (3*cos(1)*SinIntegral(1 + coth(a + b*x)))/(8*b) - (cos(3)*SinIntegral(3 + 3*coth(a + b*x)))/(8*b)]
    @test_int [sin(coth(a + b*x))^2, x, 13, (cos(2)*CosIntegral(2 - 2*coth(a + b*x)))/(4*b) - (cos(2)*CosIntegral(2 + 2*coth(a + b*x)))/(4*b) - log(1 - coth(a + b*x))/(4*b) + log(1 + coth(a + b*x))/(4*b) + (sin(2)*SinIntegral(2 - 2*coth(a + b*x)))/(4*b) - (sin(2)*SinIntegral(2 + 2*coth(a + b*x)))/(4*b)]
    @test_int [sin(coth(a + b*x))^1, x, 9, -((CosIntegral(1 - coth(a + b*x))*sin(1))/(2*b)) - (CosIntegral(1 + coth(a + b*x))*sin(1))/(2*b) + (cos(1)*SinIntegral(1 - coth(a + b*x)))/(2*b) + (cos(1)*SinIntegral(1 + coth(a + b*x)))/(2*b)]
    @test_int [csc(coth(a + b*x))^1, x, 3, (1/2)*Unintegrable((csc(coth(a + b*x))*csch(a + b*x)^2)/(-1 + coth(a + b*x)), x) - (1/2)*Unintegrable((csc(coth(a + b*x))*csch(a + b*x)^2)/(1 + coth(a + b*x)), x)]

    @test_int [cos(coth(a + b*x))^3, x, 19, -((cos(3)*CosIntegral(3 - 3*coth(a + b*x)))/(8*b)) - (3*cos(1)*CosIntegral(1 - coth(a + b*x)))/(8*b) + (3*cos(1)*CosIntegral(1 + coth(a + b*x)))/(8*b) + (cos(3)*CosIntegral(3 + 3*coth(a + b*x)))/(8*b) - (sin(3)*SinIntegral(3 - 3*coth(a + b*x)))/(8*b) - (3*sin(1)*SinIntegral(1 - coth(a + b*x)))/(8*b) + (3*sin(1)*SinIntegral(1 + coth(a + b*x)))/(8*b) + (sin(3)*SinIntegral(3 + 3*coth(a + b*x)))/(8*b)]
    @test_int [cos(coth(a + b*x))^2, x, 13, -((cos(2)*CosIntegral(2 - 2*coth(a + b*x)))/(4*b)) + (cos(2)*CosIntegral(2 + 2*coth(a + b*x)))/(4*b) - log(1 - coth(a + b*x))/(4*b) + log(1 + coth(a + b*x))/(4*b) - (sin(2)*SinIntegral(2 - 2*coth(a + b*x)))/(4*b) + (sin(2)*SinIntegral(2 + 2*coth(a + b*x)))/(4*b)]
    @test_int [cos(coth(a + b*x))^1, x, 9, -((cos(1)*CosIntegral(1 - coth(a + b*x)))/(2*b)) + (cos(1)*CosIntegral(1 + coth(a + b*x)))/(2*b) - (sin(1)*SinIntegral(1 - coth(a + b*x)))/(2*b) + (sin(1)*SinIntegral(1 + coth(a + b*x)))/(2*b)]
    @test_int [sec(coth(a + b*x))^1, x, 3, (1/2)*Unintegrable((csch(a + b*x)^2*sec(coth(a + b*x)))/(-1 + coth(a + b*x)), x) - (1/2)*Unintegrable((csch(a + b*x)^2*sec(coth(a + b*x)))/(1 + coth(a + b*x)), x)]
    end

max_stars_repo_path: test/6 Hyperbolic functions/6.4 Hyperbolic cotangent/6.4.2 Hyperbolic cotangent functions.jl


text:
const PRECHELT_REF = "[Prechelt, Lutz (1998): \"Early Stopping"*
    "- But When?\", in *Neural Networks: Tricks of the Trade*, "*
    "ed. G. Orr, Springer.](https://link.springer.com/chapter"*
    "/10.1007%2F3-540-49430-8_3)"

const STOPPING_DOC = "An early stopping criterion for loss-reporting "*
    "iterative algorithms. "

const CUSTOM_ALTERNATIVE_DOC = "For a customizable loss-based stopping "*
    "criterion, use [`WithLossDo`](@ref) or [`WithTrainingLossesDo`](@ref) "*
    "with the `stop_if_true=true` option. "

## NEVER

"""
    Never()

$STOPPING_DOC

Indicates early stopping is to be disabled.

See also [`NotANumber`](@ref), for stopping on encountering `NaN`.

"""
struct Never <: StoppingCriterion end

## OUT OF BOUNDS

"""
    InvalidValue()

$STOPPING_DOC

Stop if a loss (or training loss) is `NaN`, `Inf` or `-Inf` (or, more
precisely, if `isnan(loss)` or `isinf(loss)` is `true`).

$CUSTOM_ALTERNATIVE_DOC

"""
struct InvalidValue <: StoppingCriterion end

_isinf(x) = isinf(x)
_isinf(::Nothing) = false
_isnan(x) = isnan(x)
_isnan(::Nothing) = false

# state = `true` when `NaN`, `Inf` or `-Inf` has been encountered
update(::InvalidValue, loss, state=false) =
    state !== nothing && state || _isinf(loss) || _isnan(loss)
update_training(c::InvalidValue, loss, state) = update(c, loss, state)
done(::InvalidValue, state) = state !== nothing && state

message(::InvalidValue, state) = "Stopping early as `NaN`, "*
    "`Inf` or `-Inf` encountered. "

## TIME LIMIT

"""
    TimeLimit(; t=0.5)

$STOPPING_DOC

Stopping is triggered after `t` hours have elapsed since the stopping
criterion was initiated.

Any Julia built-in `Real` type can be used for `t`. Subtypes of
`Period` may also be used, as in `TimeLimit(t=Minute(30))`.

Internally, `t` is rounded to nearest millisecond.
``
"""
struct TimeLimit <: StoppingCriterion
    t::Millisecond
    function TimeLimit(t::Millisecond)
        t > Millisecond(0) ||
            throw(ArgumentError("Time limit `t` must be positive. "))
        return new(t)
    end
end
TimeLimit(t::T) where T <: Period = TimeLimit(convert(Millisecond, t))
# for t::T a "numeric" time in hours; assumes `round(Int, ::T)` implemented:
TimeLimit(t) = TimeLimit(round(Int, 3_600_000*t) |> Millisecond)
TimeLimit(; t =Minute(30)) = TimeLimit(t)

# state = time at initialization
update(::TimeLimit, loss, ::Nothing) = now()
update_training(::TimeLimit, loss, ::Nothing) = now()
done(criterion::TimeLimit, state) =
    state === nothing ? false : criterion.t < now() - state

## GL

# helper:

generalization_loss(E, E_opt) =  100*(E/abs(E_opt) - one(E_opt))

"""
    GL(; alpha=2.0)

$STOPPING_DOC

A stop is triggered when the (rescaled) generalization loss exceeds
the threshold `alpha`.

**Terminology.** Suppose ``E_1, E_2, ..., E_t`` are a sequence of
losses, for example, out-of-sample estimates of the loss associated
with some iterative machine learning algorithm. Then the
*generalization loss* at time `t`, is given by

`` GL_t = 100 (E_t - E_{opt}) \\over |E_{opt}|``

where ``E_{opt}`` is the minimum value of the sequence.

Reference: $PRECHELT_REF.

"""
struct GL <: StoppingCriterion
    alpha::Float64
    function GL(alpha)
        alpha > 0 ||
            throw(ArgumentError("Threshold `alpha` must be positive. "))
        return new(alpha)
    end
end
GL(; alpha=2.0) = GL(alpha)
update(::GL, loss, ::Nothing) = (loss=loss, min_loss=loss)
update(::GL, loss, state) = (loss=loss, min_loss=min(loss, state.min_loss))
function done(criterion::GL, state)
    if state === nothing
        return false
    else
        gl = generalization_loss(state.loss, state.min_loss)
        return  gl > criterion.alpha
    end
end

## PQ

# helpers:

function prepend(v, x, max_length)
    if length(v) < max_length
        vcat([x, ], v)
    else
        vcat([x, ], v[1:end-1])
    end
end

progress(losses::AbstractVector{T}) where T =
    1000*(mean(losses)/minimum(losses) - one(T))

_min(x, ::Nothing) = x
_min(::Nothing, x) = x
_min(x, y) = min(x, y)

"""
    PQ(; alpha=0.75, k=5, tol=eps(Float64))

A stopping criterion for training iterative supervised learners.

A stop is triggered when Prechelt's progress-modified generalization
loss exceeds the threshold ``PQ_T > alpha``, or if the training progress drops
below ``P_j ≤ tol``. Here `k` is the number of training (in-sample) losses used to
estimate the training progress.

## Context and explanation of terminology

The *training progress* at time ``j`` is defined by

`` P_j = 1000 |M - m|/|m| ``

where ``M`` is the mean of the last `k` training losses ``F_1, F_2, …, F_k``
and ``m`` is the minimum value of those losses.

The *progress-modified generalization loss* at time ``t`` is then given by

`` PQ_t = GL_t / P_t``

where ``GL_t`` is the generalization loss at time ``t``; see
[`GL`](@ref).

PQ will stop when the following are true:

1) At least `k` training samples have been collected via
   `done!(c::PQ, loss; training = true)` or `update_training(c::PQ, loss, state)`
2) The last update was an out-of-sample update.
   (`done!(::PQ, loss; training=true)` is always false)
3) The progress-modified generalization loss exceeds the threshold
   ``PQ_t > alpha`` **OR** the training progress stalls ``P_j ≤ tol``.

Reference: $PRECHELT_REF.

"""
struct PQ <: StoppingCriterion
    alpha::Float64
    k::Union{Int,Float64} # could be Inf
    tol::Float64
    function PQ(alpha, k, tol)
        alpha > 0 ||
            throw(ArgumentError("Threshold `alpha` must be positive. "))
        if k < 2 || !(k isa Int || isinf(k))
            throw(ArgumentError("`k` must `Int` or `Inf` and `k > 1`. "))
        end
        tol > 0 ||
            throw(ArgumentError("`tol` must be positive. "))

        return new(alpha, k, tol)
    end
end
PQ(; alpha=0.75, k=5, tol=eps(Float64)) = PQ(alpha, k, tol)

struct PQState{T}
    training_losses::Vector{T}
    waiting_for_out_of_sample::Bool
    loss::Union{Nothing,T}
    min_loss::Union{Nothing,T}
end

update_training(::PQ, loss, ::Nothing) = PQState([loss, ], true, nothing, nothing)
update(::PQ, loss::T, ::Nothing) where T = PQState(T[], false, loss, loss)

function update_training(criterion::PQ, loss, state)
    training_losses = prepend(state.training_losses, loss, criterion.k)
    return PQState(training_losses, true, state.loss, state.min_loss)
end

function update(::PQ, loss, state)
    min_loss = _min(loss, state.min_loss)
    return PQState(state.training_losses, false, loss, min_loss)
end

function done(criterion::PQ, state)
    state === nothing && return false
    state.loss === nothing && return false
    state.waiting_for_out_of_sample && return false
    length(state.training_losses) < criterion.k && return false
    GL = generalization_loss(state.loss, state.min_loss)
    P = progress(state.training_losses)
    P > criterion.tol || return true
    PQ = GL/P
    return  PQ > criterion.alpha
end

## PATIENCE

# This is UP_s in Prechelt 1998

"""
    Patience(; n=5)

$STOPPING_DOC

A stop is triggered by `n` consecutive increases in the loss.

Denoted "_UP_s" in $PRECHELT_REF.

$CUSTOM_ALTERNATIVE_DOC

"""
struct Patience <: StoppingCriterion
    n::Int
    function Patience(n::Int)
        n > 0 ||
            throw(ArgumentError("The patience level `n` must be positive. "))
        return new(n)
    end
end
Patience(; n=5) = Patience(n)

# Prechelt alias:
const UP = Patience

update(::Patience, loss, ::Nothing) = (loss=loss, n_increases=0)
@inline function update(::Patience, loss, state)
    old_loss, n = state
    if loss > old_loss
        n += 1
    else
        n = 0
    end
    return (loss=loss, n_increases=n)
end

done(criterion::Patience, state) =
    state === nothing ? false : state.n_increases == criterion.n

## NUMBER SINCE BEST

"""
    NumberSinceBest(; n=6)

$STOPPING_DOC

A stop is triggered when the number of calls to the control, since the
lowest value of the loss so far, is `n`.

$CUSTOM_ALTERNATIVE_DOC

"""
struct NumberSinceBest <: StoppingCriterion
    n::Int
    function NumberSinceBest(n::Int)
        n > 0 ||
            throw(ArgumentError("`n` must be positive. "))
        return new(n)
    end
end
NumberSinceBest(; n=6) = NumberSinceBest(n)

update(::NumberSinceBest, loss, ::Nothing) = (best=loss, number_since_best=0)
@inline function update(::NumberSinceBest, loss, state)
    best, number_since_best = state
    if loss < best
        best = loss
        number_since_best = 0
    else
        number_since_best += 1
    end
    return (best=best, number_since_best=number_since_best)
end

done(criterion::NumberSinceBest, state) =
    state === nothing ? false : state.number_since_best == criterion.n

# # NUMBER LIMIT

"""
    NumberLimit(; n=100)

$STOPPING_DOC

A stop is triggered by `n` consecutive loss updates, excluding
"training" loss updates.

If wrapped in a `stopper::EarlyStopper`, this is the number of calls
to `done!(stopper)`.

"""
struct NumberLimit <: StoppingCriterion
    n::Int
    function NumberLimit(n::Int)
        n > 0 ||
            throw(ArgumentError("`n` must be positive. "))
        return new(n)
    end
end
NumberLimit(; n=100) = NumberLimit(n)

update(criterion::NumberLimit, loss, ::Nothing) = 1
update(::NumberLimit, loss, state) = state+1
done(criterion::NumberLimit, state) =
    state === nothing ? false : state >= criterion.n

# ## THRESHOLD

"""
    Threshold(; value=0.0)

$STOPPING_DOC

A stop is triggered as soon as the loss drops below `value`.

$CUSTOM_ALTERNATIVE_DOC

"""
struct Threshold <: StoppingCriterion
    value::Float64
end
Threshold(; value=0.0) = Threshold(value)

update(::Threshold, loss, state) = loss
done(criterion::Threshold, state) =
    state === nothing ? false : state < criterion.value

"""
    Warmup(c::StoppingCriterion, n)

Wait for `n` updates before checking stopping criterion `c`
"""
struct Warmup{C} <: StoppingCriterion where {C <: StoppingCriterion}
    criterion::C
    n::Int
    function Warmup(criterion::C, n::N) where {C <: StoppingCriterion, N <: Integer}
        n > 0 || throw(ArgumentError("`n` must be positive. "))
        new{C}(criterion, Int(n))
    end
end

# Constructors for Warmup
Warmup() = Warmup(InvalidValue())   # Default for testing
Warmup(c; n = 1) = Warmup(c, n)     # Provide kwargs interface

# Initialize inner state for type-stability, and record first observation
update(c::Warmup, loss, ::Nothing) = update(c, loss)
update(criterion::Warmup, loss) = (1, update(criterion.criterion, loss))

# Catch uninitialized state
update_training(c::Warmup, loss, ::Nothing) = update_training(c, loss)
update_training(c::Warmup, loss) = (1, update_training(c.criterion, loss))

# Handle update vs update_training
update(c::Warmup, loss, state) = _update(update, c, loss, state)
update_training(c::Warmup, loss, state) =
    _update(update_training, c, loss, state)

# Dispatch update and update_training here
function _update(f::Function, criterion::Warmup, loss, state)
    n, inner = state
    n += 1
    if n <= criterion.n
        # Skip inner criterion
        return n, inner
    elseif n == criterion.n+1
        # First step of inner criterion
        return n, f(criterion.criterion, loss)
    else
        # Step inner criterion
        return n, f(criterion.criterion, loss, inner)
    end
end

function done(criterion::Warmup, state)
    # Only check if inner criterion is done after n updates
    state === nothing && return false
    return state[1] <= criterion.n ? false : done(criterion.criterion, state[2])
end

message(c::Warmup, state) = message(c.criterion, state[2])

## NOT A NUMBER (deprecated)

"""
    NotANumber()

$STOPPING_DOC

Stop if a loss of `NaN` is encountered.

**Now deprecated** in favour of [`InvalidValue`](@ref).

"""
struct NotANumber <: StoppingCriterion
    function NotANumber()
        Base.depwarn("`NotANumber()` is deprecated. Use `InvalidValue()` "*
                     "to trap `NaN`, `Inf` or `-Inf`. ", :NotANumber)
        return new()
    end
end

# state = `true` when NaN has been encountered
update(::NotANumber, loss, state) = state !== nothing && state || isnan(loss)
update_training(c::NotANumber, loss, state) = update(c, loss, state)

done(::NotANumber, state) = state !== nothing && state

message(::NotANumber, state) = "Stopping early as NaN encountered. "

ext: jl


text:
# test for updateBL function
# Claudia April 2015

include("../examples/tree_example.jl");
printEdges(net)

parts = edgesParts(net);
[isInternalEdge(e)?e.number:0 for e in net.edge]
i = 1
[n.number for n in parts[i].part1]
[n.number for n in parts[i].part2]
[n.number for n in parts[i].part3]
[n.number for n in parts[i].part4]

include("../src/types.jl")
include("../src/functions.jl")
net=readTopologyUpdate("1_astral.out");
printEdges(net)
net.names
df0 = readtable("HGT_truenet_expCF.csv")
d = readTableCF(df0); #expCF

parts = edgesParts(net);
df = makeTable(net,parts,d)
x = updateBL!(net,d)


text:
mutable struct SystemDescriptor <: PowerSystemCaseBuilderType
    name::AbstractString
    description::AbstractString
    category::Type{<:SystemCategory}
    raw_data::AbstractString
    build_function::Function
    download_function::Union{Nothing, Function}
    stats::Union{Nothing, SystemBuildStats}
end

function SystemDescriptor(;
    name,
    description,
    category,
    raw_data,
    build_function,
    download_function = nothing,
    stats = nothing,
)
    return SystemDescriptor(
        name,
        description,
        category,
        raw_data,
        build_function,
        download_function,
        stats,
    )
end

get_name(v::SystemDescriptor) = v.name
get_description(v::SystemDescriptor) = v.description
get_category(v::SystemDescriptor) = v.category
get_raw_data(v::SystemDescriptor) = v.raw_data
get_build_function(v::SystemDescriptor) = v.build_function
get_download_function(v::SystemDescriptor) = v.download_function
get_stats(v::SystemDescriptor) = v.stats

set_name!(v::SystemDescriptor, value::String) = v.name = value
set_description!(v::SystemDescriptor, value::String) = v.description = value
set_category!(v::SystemDescriptor, value::Type{<:SystemCategory}) = v.category = value
set_raw_data!(v::SystemDescriptor, value::String) = v.raw_data = value
set_build_function!(v::SystemDescriptor, value::Function) = v.build_function = value
set_download_function!(v::SystemDescriptor, value::Function) = v.download_function = value
set_stats!(v::SystemDescriptor, value::SystemBuildStats) = v.stats = value

update_stats!(v::SystemDescriptor, deserialize_time::Float64) =
    update_stats!(v.stats, deserialize_time)

max_stars_repo_name: NREL-SIIP/PowerSystemBuilder

max_issues_repo_name: NREL-SIIP/PowerSystemBuilder

max_forks_repo_name: NREL-SIIP/PowerSystemBuilder


text:
using Test

@testset "PolyVar and Monomial tests" begin
    @testset "PolyVar macro index set" begin
        n = 3
        @polyvar x[1:n] y z[1:n-1] u[1:n,1:n-1]
        @test x isa Vector{PolyVar{true}}
        @test y isa PolyVar{true}
        @test z isa Vector{PolyVar{true}}
        @test u isa Matrix{PolyVar{true}}
        @test length(x) == 3
        @test length(z) == 2
        @test size(u) == (3, 2)
        @test x[1] > x[2] > x[3] > y > z[1] > z[2]
        @test u[1, 1] > u[2, 1] > u[2, 2]

        @polyvar a[1:5, 1:3, 1:2]
        @test size(a) == (5, 3, 2)
    end
    @testset "PolyVar macro tuple return" begin
        vars = @polyvar x y z
        @test vars isa Tuple
        @test vars == (x, y, z)

        vars = @ncpolyvar x y z
        @test vars isa Tuple
        @test vars == (x, y, z)
    end

    @testset "variable_union_type" begin
        @polyvar x
        @test DynamicPolynomials.MP.variable_union_type(x) == PolyVar{true}
        @test DynamicPolynomials.MP.variable_union_type(x^2) == PolyVar{true}
        @test DynamicPolynomials.MP.variable_union_type(2x) == PolyVar{true}
        @test DynamicPolynomials.MP.variable_union_type(x + 1) == PolyVar{true}
        @ncpolyvar y
        @test DynamicPolynomials.MP.variable_union_type(y) == PolyVar{false}
        @test DynamicPolynomials.MP.variable_union_type(y^2) == PolyVar{false}
        @test DynamicPolynomials.MP.variable_union_type(2y) == PolyVar{false}
        @test DynamicPolynomials.MP.variable_union_type(y + 1) == PolyVar{false}
    end
    @testset "PolyVar" begin
        @test zeroterm(PolyVar{true}) == 0
        @test zero(PolyVar{true}) == 0
        @test one(PolyVar{false}) == 1
        @polyvar x
        @test zeroterm(x) isa Term{true, Int}
        @test zero(x) isa Polynomial{true, Int}
        @test one(x) isa Monomial{true}
    end
    @testset "Monomial" begin
        @test zeroterm(Monomial{false}) == 0
        @test zero(Monomial{false}) == 0
        @test one(Monomial{true}) == 1
        if VERSION ≥ v"0.7-"
            @test_throws ErrorException Monomial(2)
            @test (@inferred Monomial(1)) isa Monomial{true}
            @test Monomial(1) == 1
            @test_throws ErrorException convert(Monomial{true}, 2)
            @test (@inferred convert(Monomial{true}, 1)) isa Monomial{true}
            @test convert(Monomial{true}, 1) == 1
            @test_throws ErrorException convert(Monomial{false}, 2)
            @test (@inferred convert(Monomial{false}, 1)) isa Monomial{false}
            @test convert(Monomial{false}, 1) == 1
        end
        @polyvar x
        @test_throws ArgumentError Monomial{true}([x], [1,0])
        @test zeroterm(x^2) isa Term{true, Int}
        @test zero(x^2) isa Polynomial{true, Int}
        @test one(x^2) isa Monomial{true}

        @polyvar y
        @test Monomial([x, y], [1, 0]) == x
        @test x != Monomial([x, y], [0, 1])
    end
    @testset "MonomialVector" begin
        @polyvar x y
        @test_throws AssertionError MonomialVector{true}([x], [[1], [1,0]])
        @test_throws AssertionError monomials([y, x], 1:2) # should be [x, y]
        X = MonomialVector([x, 1, x*y])
        @test variables(X) == [x, y]
        @test X.Z == [[1, 1], [1, 0], [0, 0]]
        @test MonomialVector{true}([1]) isa MonomialVector{true}
        @test MonomialVector{false}([1]) isa MonomialVector{false}
        a = [1, 5, 3]
        b, Y = monovec(a, X)
        @test b === a
        @test Y === X
        σ, Y = sortmonovec(X)
        @test σ == 1:length(X)
        @test Y === X
        XX = @inferred mergemonovec([X, X])
        @test XX isa MonomialVector{true}
        @test X == XX

        @test 2 != MonomialVector([x, y], 1)
        @test x != MonomialVector([x, y], 1)
        @test MonomialVector([x, y], [[1, 0], [0, 0]]) == MonomialVector([x], [[1], [0]])
    end
    @testset "Non-commutative" begin
        @ncpolyvar x
        @test_throws ArgumentError Monomial{false}([x], [1,0])
        @test_throws AssertionError MonomialVector{false}([x], [[1], [1,0]])
    end
    @testset "NC PolyVar * Monomial" begin
        @ncpolyvar x y z
        m = y * Monomial([y, z, x, z], [0, 0, 2, 1])
        @test variables(m) == [y, z, x, z]
        @test m.z == [1, 0, 2, 1]
        m = x * Monomial([z, y, y, z], [0, 0, 2, 1])
        @test variables(m) == [z, x, y, y, z]
        @test m.z == [0, 1, 0, 2, 1]
        m = x * Monomial([y, z, y, z], [0, 0, 2, 1])
        @test variables(m) == [y, z, x, y, z]
        @test m.z == [0, 0, 1, 2, 1]
    end
    @testset "NC Monomial * PolyVar" begin
        @ncpolyvar x y z
        m = Monomial([x, z, x, y], [2, 1, 0, 0]) * y
        @test variables(m) == [x, z, x, y]
        @test m.z == [2, 1, 0, 1]
        m = Monomial([x, y, y, x], [2, 1, 0, 0]) * z
        @test variables(m) == [x, y, y, z, x]
        @test m.z == [2, 1, 0, 1, 0]
        m = Monomial([x, y, x, y], [2, 1, 0, 0]) * z
        @test variables(m) == [x, y, z, x, y]
        @test m.z == [2, 1, 1, 0, 0]
    end

    @testset "Evaluation" begin
        @polyvar x y
        @test (x^2*y)(3,2) == 18
        @test (x^2*y)((3,2)) == 18
        @test (x^2*y)([3,2]) == 18
        @test (x^2)(3) == 9
        @test (x)(3) == 3
    end
    @testset "TODO remove when added to MP" begin
        @polyvar x y
        @test x == DynamicPolynomials.MP.mapexponents!(div, x^1, x * y^2)
        for z in [x * y, x^2, y^2]
            @test y == DynamicPolynomials.MP.mapexponents_to!(z, -, x * y^2, x * y)
        end
    end
    # TODO add to MP
    @testset "Indexing with boolean" begin
        @polyvar x y
        X = monomials([x, y], 2)
        @test X[[true, false, true]] == monovec([x^2, y^2])
        X = monomials([x, y], 0:1)
        @test filter(mono -> degree(mono) == 1, X) == monovec([x, y])
        @test filter(mono -> degree(mono) == 0, X) == monovec([x^0])
    end
end

max_stars_repo_name: blegat/DynamicPolynomials

max_issues_repo_name: blegat/DynamicPolynomials

max_forks_repo_name: blegat/DynamicPolynomials


text:
using GLMakie
GLMakie.activate!()
using Makie
using JLD
##############################

wait_for_key(prompt) = (print(stdout, prompt); read(stdin, 1); nothing)

function replot_problem_results(problem, problem_result_dir)
    println("Loading data from $problem_result_dir")
    result_data = load(joinpath(problem_result_dir, "data.jld"))["result_data"]
    objs = Float64[]

    r1_x = result_data["init_solve"]["minimizer"]
    fig1 = TopOpt.TopOptProblems.Visualization.visualize(problem, topology=r1_x, undeformed_mesh_color=(:black, 1.0),
        display_supports=false)
    hidedecorations!(fig1.current_axis.x)
    save(joinpath(problem_result_dir, "r1_replot.pdf"), fig1, pt_per_unit = 1)
    push!(objs, result_data["init_solve"]["minimum"] )

    for k in keys(result_data)
        if !occursin("deflate", k)
            continue
        end
        df_result_minimizer = result_data[k]["minimizer"]
        df_result_minimum = result_data[k]["minimum"]
        push!(objs, df_result_minimum)

        fig2 = TopOpt.TopOptProblems.Visualization.visualize(problem, topology=df_result_minimizer[1:end-1],
            undeformed_mesh_color=(:black, 1.0), display_supports=false)
        hidedecorations!(fig2.current_axis.x)
        save(joinpath(problem_result_dir, "$(k)_replot.pdf"), fig2, pt_per_unit = 1)
    end

    objs ./= minimum(objs)
    obj_fig = Figure(resolution=(800,800), font="Arial")
    ax = Axis(obj_fig[1,1])
    # ax.xlabel = "deflation iterations"
    # ax.ylabel = "objectives"
    lines!(ax, 0:1:length(objs)-1, objs)
    save(joinpath(problem_result_dir, "_obj_history_replot.pdf"), obj_fig)
end

function plot_runtime(problem_result_dir)
    println("Loading data from $problem_result_dir")
    result_data = load(joinpath(problem_result_dir, "data.jld"))["result_data"]
    runtimes = Float64[]
    push!(runtimes, result_data["init_solve"]["runtime"] )
    for k in keys(result_data)
        if !occursin("deflate", k)
            continue
        end
        push!(runtimes, result_data[k]["runtime"])
    end
    println(runtimes)
    obj_fig = Figure(resolution=(800,800), font="Arial")
    ax = Axis(obj_fig[1,1])
    lines!(ax, 0:1:length(runtimes)-1, runtimes)
    save(joinpath(problem_result_dir, "_runtimes.png"), obj_fig)
    println("Initial runtime $(runtimes[1]), avr deflated runtime $(sum(runtimes[2:end])/(length(runtimes)-1))")
end


text:
"""
    sparsify(x, nnzrate)

Replace `AbstractMatrix` `x` with `SparseMatrixCSC` if at most `nnzrate` fraction of elements is non-zero.

```jldoctest
julia> n = ArrayNode([0 0; 0 0])
2×2 ArrayNode{Matrix{Int64}, Nothing}:
 0  0
 0  0

julia> Mill.mapdata(i -> sparsify(i, 0.05), n)
2×2 ArrayNode{SparseMatrixCSC{Int64, Int64}, Nothing}:
 ⋅  ⋅
 ⋅  ⋅
```

See also: [`Mill.mapdata`](@ref).
"""
sparsify(x, nnzrate) = x
sparsify(x::Matrix, nnzrate) = (mean(x .!= 0) < nnzrate) ? sparse(x) : x

# can be removed when https://github.com/FluxML/Flux.jl/pull/1357 is merged
function Base.:*(A::AbstractMatrix, B::Adjoint{Bool,<: Flux.OneHotArray})
    m = size(A,1)
    Y = similar(A, m, size(B,2))
    Y .= 0
    BT = B'
    for (j, ix) in enumerate(BT.indices)
        for i in 1:m
            @inbounds Y[i, ix] += A[i, j]
        end
    end
    Y
end

"""
    pred_lens(p, n)

Return a `Vector` of `Setfield.Lens`es for accessing all nodes/fields in `n` conforming to predicate `p`.

# Examples
```jldoctest
julia> n = ProductNode((BagNode(missing, bags([0:-1, 0:-1])), ArrayNode([1 2; 3 4])))
ProductNode with 2 obs
  ├── BagNode with 2 obs
  │     └── ∅
  └── ArrayNode(2×2 Array with Int64 elements)

julia> pred_lens(x -> x isa ArrayNode, n)
1-element Vector{Setfield.ComposedLens{Setfield.PropertyLens{:data}, Setfield.IndexLens{Tuple{Int64}}}}:
 (@lens _.data[2])
```

See also: [`list_lens`](@ref), [`find_lens`](@ref), [`findnonempty_lens`](@ref).
"""
pred_lens(p::Function, n) = _pred_lens(p, n)

"""
    list_lens(n)

Return a `Vector` of `Setfield.Lens`es for accessing all nodes/fields in `n`.

# Examples
```jldoctest
julia> n = ProductNode((BagNode(missing, bags([0:-1, 0:-1])), ArrayNode([1 2; 3 4])))
ProductNode with 2 obs
  ├── BagNode with 2 obs
  │     └── ∅
  └── ArrayNode(2×2 Array with Int64 elements)

julia> list_lens(n)
9-element Vector{Lens}:
 (@lens _)
 (@lens _.data[1])
 (@lens _.data[1].data)
 (@lens _.data[1].bags)
 (@lens _.data[1].metadata)
 (@lens _.data[2])
 (@lens _.data[2].data)
 (@lens _.data[2].metadata)
 (@lens _.metadata)
```

See also: [`pred_lens`](@ref), [`find_lens`](@ref), [`findnonempty_lens`](@ref).
"""
list_lens(n) = pred_lens(t -> true, n)

"""
    findnonempty_lens(n)

Return a `Vector` of `Setfield.Lens`es for accessing all nodes/fields in `n` that have at least one observation.

# Examples
```jldoctest
julia> n = ProductNode((BagNode(missing, bags([0:-1, 0:-1])), ArrayNode([1 2; 3 4])))
ProductNode with 2 obs
  ├── BagNode with 2 obs
  │     └── ∅
  └── ArrayNode(2×2 Array with Int64 elements)

julia> findnonempty_lens(n)
3-element Vector{Lens}:
 (@lens _)
 (@lens _.data[1])
 (@lens _.data[2])
```

See also: [`pred_lens`](@ref), [`list_lens`](@ref), [`find_lens`](@ref).
"""
findnonempty_lens(n) = pred_lens(t -> t isa AbstractMillNode && nobs(t) > 0, n)

"""
    find_lens(n, x)

Return a `Vector` of `Setfield.Lens`es for accessing all nodes/fields in `n` that return `true` when
compared to `x` using `Base.===`.

# Examples
```jldoctest
julia> n = ProductNode((BagNode(missing, bags([0:-1, 0:-1])), ArrayNode([1 2; 3 4])))
ProductNode with 2 obs
  ├── BagNode with 2 obs
  │     └── ∅
  └── ArrayNode(2×2 Array with Int64 elements)

julia> find_lens(n, n.data[1])
1-element Vector{Setfield.ComposedLens{Setfield.PropertyLens{:data}, Setfield.IndexLens{Tuple{Int64}}}}:
 (@lens _.data[1])
```

See also: [`pred_lens`](@ref), [`list_lens`](@ref), [`findnonempty_lens`](@ref).
"""
find_lens(n, x) = pred_lens(t -> t === x, n)

_pred_lens(p::Function, n) = p(n) ? [IdentityLens()] : Lens[]
function _pred_lens(p::Function, n::T) where T <: AbstractMillStruct
    res = [map(l -> PropertyLens{k}() ∘ l, _pred_lens(p, getproperty(n, k))) for k in fieldnames(T)]
    res = vcat(filter(!isempty, res)...)
    p(n) ? [IdentityLens(); res] : res
end
function _pred_lens(p::Function, n::Union{Tuple, NamedTuple})
    res = [map(l -> IndexLens(tuple(i)) ∘ l, _pred_lens(p, n[i])) for i in eachindex(n)]
    vcat(filter(!isempty, res)...)
end

"""
    code2lens(n, c)

Convert code `c` from [HierarchicalUtils.jl](@ref) traversal to `Setfield.Lens` such that they access
the same node in tree `n`.

# Examples
```jldoctest
julia> n = ProductNode((BagNode(missing, bags([0:-1, 0:-1])), ArrayNode([1 2; 3 4])));

julia> printtree(n; trav=true)
ProductNode with 2 obs [""]
  ├── BagNode with 2 obs ["E"]
  │     └── ∅ ["M"]
  └── ArrayNode(2×2 Array with Int64 elements) ["U"]

julia> code2lens(n, "U")
(@lens _.data[2])
```

See also: [`lens2code`](@ref).
"""
code2lens(n::AbstractMillStruct, c::AbstractString) = find_lens(n, n[c]) |> only

"""
    lens2code(n, l)

Convert `Setfield.Lens` l to code `c` from [HierarchicalUtils.jl](@ref) traversal such that they access
the same node in tree `n`.

# Examples
```jldoctest
julia> n = ProductNode((BagNode(missing, bags([0:-1, 0:-1])), ArrayNode([1 2; 3 4])));

julia> printtree(n; trav=true)
ProductNode with 2 obs [""]
  ├── BagNode with 2 obs ["E"]
  │     └── ∅ ["M"]
  └── ArrayNode(2×2 Array with Int64 elements) ["U"]

julia> lens2code(n, (@lens _.data[2]))
"U"

```

See also: [`code2lens`](@ref).
"""
lens2code(n::AbstractMillStruct, l::Lens) = HierarchicalUtils.find_traversal(n, get(n, l)) |> only

"""
    model_lens(m, l)

Convert `Setfield.Lens` `l` for a data node to a new lens for accessing the same location in model `m`.

# Examples
```jldoctest
julia> n = ProductNode((BagNode(ArrayNode(randn(2, 2)), bags([0:-1, 0:-1])), ArrayNode([1 2; 3 4])))
ProductNode with 2 obs
  ├── BagNode with 2 obs
  │     └── ArrayNode(2×2 Array with Float64 elements)
  └── ArrayNode(2×2 Array with Int64 elements)

julia> m = reflectinmodel(n)
ProductModel … ↦ ArrayModel(Dense(20, 10))
  ├── BagModel … ↦ BagCount([SegmentedMean(10); SegmentedMax(10)]) ↦ ArrayModel(Dense(21, 10))
  │     └── ArrayModel(Dense(2, 10))
  └── ArrayModel(Dense(2, 10))

julia> model_lens(m, (@lens _.data[2]))
(@lens _.ms[2])
```

See also: [`data_lens`](@ref).
"""
function model_lens(model, lens::ComposedLens)
    outerlens = model_lens(model, lens.outer)
    outerlens ∘ model_lens(get(model, outerlens), lens.inner)
end
model_lens(::ArrayModel, ::PropertyLens{:data}) = @lens _.m
model_lens(::BagModel, ::PropertyLens{:data}) = @lens _.im
model_lens(::ProductModel, ::PropertyLens{:data}) = @lens _.ms
model_lens(::Union{NamedTuple, Tuple}, lens::IndexLens) = lens
model_lens(::Union{AbstractMillModel, NamedTuple, Tuple}, lens::IdentityLens) = lens

"""
    data_lens(n, l)

Convert `Setfield.Lens` `l` for a model node to a new lens for accessing the same location in data node `n`.

# Examples
```jldoctest
julia> n = ProductNode((BagNode(ArrayNode(randn(2, 2)), bags([0:-1, 0:-1])), ArrayNode([1 2; 3 4])))
ProductNode with 2 obs
  ├── BagNode with 2 obs
  │     └── ArrayNode(2×2 Array with Float64 elements)
  └── ArrayNode(2×2 Array with Int64 elements)

julia> m = reflectinmodel(n)
ProductModel … ↦ ArrayModel(Dense(20, 10))
  ├── BagModel … ↦ BagCount([SegmentedMean(10); SegmentedMax(10)]) ↦ ArrayModel(Dense(21, 10))
  │     └── ArrayModel(Dense(2, 10))
  └── ArrayModel(Dense(2, 10))

julia> data_lens(n, (@lens _.ms[2]))
(@lens _.data[2])
```

See also: [`data_lens`](@ref).
"""
function data_lens(ds, lens::ComposedLens)
    outerlens = data_lens(ds, lens.outer)
    outerlens ∘ data_lens(get(ds, outerlens), lens.inner)
end
data_lens(::ArrayNode, ::PropertyLens{:m}) = @lens _.data
data_lens(::AbstractBagNode, ::PropertyLens{:im}) = @lens _.data
data_lens(::AbstractProductNode, ::PropertyLens{:ms}) = @lens _.data
data_lens(::Union{NamedTuple, Tuple}, lens::IndexLens) = lens
data_lens(::Union{AbstractMillNode, NamedTuple, Tuple}, lens::IdentityLens) = lens

"""
    replacein(n, old, new)

Replace in data node or model `n` each occurence of `old` by `new`.

Short description

# Examples
```jldoctest
julia> n = ProductNode((BagNode(ArrayNode(randn(2, 2)), bags([0:-1, 0:-1])), ArrayNode([1 2; 3 4])))
ProductNode with 2 obs
  ├── BagNode with 2 obs
  │     └── ArrayNode(2×2 Array with Float64 elements)
  └── ArrayNode(2×2 Array with Int64 elements)

julia> replacein(n, n.data[1], ArrayNode(maybehotbatch([1, 2], 1:2)))
ProductNode with 2 obs
  ├── ArrayNode(2×2 MaybeHotMatrix with Bool elements)
  └── ArrayNode(2×2 Array with Int64 elements)
```
"""
replacein(x, oldnode, newnode) = x
replacein(x::Tuple, oldnode, newnode) = tuple([replacein(m, oldnode, newnode) for m in x]...)
replacein(x::NamedTuple, oldnode, newnode) = (; [k => replacein(x[k], oldnode, newnode) for k in keys(x)]...)

function replacein(x::T, oldnode, newnode) where T <: AbstractMillStruct
    x === oldnode && return(newnode)
    fields = map(f -> replacein(getproperty(x, f), oldnode, newnode), fieldnames(T))
    n = nameof(T)
    p = parentmodule(T)
    eval(:($p.$n))(fields...)
end

function replacein(x::LazyNode{N}, oldnode, newnode) where {N}
    x === oldnode && return newnode
    LazyNode{N}(replacein(x.data, oldnode, newnode))
end

function replacein(x::LazyModel{N}, oldnode, newnode) where {N}
    x === oldnode && return newnode
    LazyModel{N}(replacein(x.m, oldnode, newnode))
end


text:
## Contour plot -- but this is too slow to be usable
## algorithm from MATLAB
## cf. http://www.mathworks.com/help/matlab/creating_plots/contour-plots.html
type Contourc
    contours
end

function contourc(f::Function, x, y; cs::Union(Nothing, Number)=nothing)

    fxy = [f(x,y) for x in x, y in y]

    ## we have edges 1,2,3,4
    ## with 1 on bottom, 2 on right, 3 top, 4 left. So
    ## exiting 1 goes to 3; 2->4, 3->1 and 4->2

    ## helper functions
    prop(c, z0, z1) = (c - z0)/(z1-z0)
    function interp_square(c, i,j)
        ## worry of i,j on boundary of 1:length(x), length(y)
        if (i < 1 || i >= length(x)) || (j < 1 || j >= length(y))
            return -1 * ones(4)
        else
            [prop(c, fxy[i,j], fxy[i+1,j]),
             prop(c, fxy[i+1,j], fxy[i+1, j+1]),
             prop(c, fxy[i,j+1], fxy[i+1, j+1]),
             prop(c, fxy[i,j], fxy[i, j+1])]
        end
    end
    insquare(x) =  any(0 .<= x .<= 1)

    function interp_point(edge, i, j, t)
        if edge == 1
            newx = x[i] + t*(x[i+1] - x[i])
            newy = y[j]
        elseif edge == 2
            newx = x[i+1]
            newy = y[j] + t*(y[j+1] - y[j])
        elseif edge == 3
            newx = x[i] + t*(x[i+1] - x[i])
            newy = y[j+1]
        else
            newx = x[i]
            newy = y[j] + t*(y[j+1] - y[j])
        end
        (newx, newy)
    end

    function which_next(edge, i, j)
        if edge == 1
            (i, j-1)
        elseif edge == 2
            (i+1, j)
        elseif edge == 3
            (i, j+1)
        else
            (i-1, j)
        end
    end

    function next_square(c, i, j, enter_edge, cx, cy, m)
##        println("Chasing $i $j")

        sq = interp_square(c, i, j)
        w = 0 .<= sq .<= 1
        ## check if 2 or more (saddle point) XXX

        ## what is next edge?
        if enter_edge == 1
            next_edge = setdiff((1:4)[w], 3)
        elseif enter_edge == 2
            next_edge = setdiff((1:4)[w], 4)
        elseif enter_edge == 3
            next_edge = setdiff((1:4)[w], 1)
        else
            next_edge = setdiff((1:4)[w], 2)
        end

        if length(next_edge) == 0
            return (cx, cy)
        end

        next_edge = next_edge[1]

        ##
        newx, newy = interp_point(next_edge, i, j, sq[next_edge])
        push!(cx, newx); push!(cy, newy)

        if m[i,j] == 1
            ## already visited
            return (cx, cy)
        end
        m[i,j] = 1
        next_square(c, which_next(next_edge, i,j)..., next_edge, cx, cy, m)
    end

    function chase_square(c, i, j, m)
##        println("chase $i $j")

        sq = interp_square(c, i, j)
        w = 0 .<= sq .<= 1

        ## chase both edges, then put together
        edges = [1:4][w]
        ## should be 2, might be more (saddle point)

        m[i,j] = 1                  # visited

        out = map(edges) do edge
            cx = Float64[]; cy = Float64[]
            t = sq[edge]
            newx, newy = interp_point(edge, i,j,t)
            push!(cx, newx); push!(cy, newy)
            next_square(c, which_next(edge, i,j)..., edge, cx, cy, m)
        end
        ## out is array of tuples
        if !isa(out[1], Nothing)
            ([reverse(out[1][1]), out[2][1]], [reverse(out[1][2]), out[2][2]])
        else
            nothing
        end

    end

## for each level to plot
    if isa(cs, Nothing)
        cs = linspace(minimum(fxy), maximum(fxy), 7+2)[2:8]
    else
        cs = [cs]
    end
    contours = Any[]

    for c in cs
        m = zeros(Int, size(fxy)...)
        c_contours = Any[]

        for i in 2:length(x)-1, j in 2:length(y)-1
            if m[i,j] == 1
                continue
            end
            sq = interp_square(c, i, j)
            if insquare(sq)
                path = chase_square(c, i,j, m)
##                println("Chased path:", path)
                if !isa(path, Nothing)
                    push!(c_contours, path)
                end
            else
                m[i,j] = 1              # visited
            end
        end
        push!(contours, (c, c_contours))
    end
    Contourc(contours)
end

## contour plot
function plot(f::Contourc; kwargs...)
    p = FramedPlot()
    ## out is array of tuples (c, c_countrous)
    for clevels in f.contours
        for contours in clevels[2]
            add(p, Curve(contours[1], contours[2]))
        end
    end
    for (k, v) in kwargs
        setattr(p, k, v)
    end

    p
end