text:
@testset "4.3.11 (e x)^m (a+b tan(c+d x^n))^p" begin
    (a, b, c, d, x, ) = @variables a b c d x

    #= ::Package:: =#

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

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

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [x^3*(a + b*tan(c + d*x^2)), x, 7, (a*x^4)/4 + (im/4)*b*x^4 - (b*x^2*log(1 + ℯ^((2*im)*(c + d*x^2))))/(2*d) + ((im/4)*b*PolyLog(2, -ℯ^((2*im)*(c + d*x^2))))/d^2]
    @test_int [x^2*(a + b*tan(c + d*x^2)), x, 2, (a*x^3)/3 + b*Unintegrable(x^2*tan(c + d*x^2), x)]
    @test_int [x^1*(a + b*tan(c + d*x^2)), x, 4, (a*x^2)/2 - (b*log(cos(c + d*x^2)))/(2*d)]
    @test_int [a + b*tan(c + d*x^2), x, 1, a*x + b*Unintegrable(tan(c + d*x^2), x)]
    @test_int [(a + b*tan(c + d*x^2))/x^1, x, 2, a*log(x) + b*Unintegrable(tan(c + d*x^2)/x, x)]
    @test_int [(a + b*tan(c + d*x^2))/x^2, x, 2, -(a/x) + b*Unintegrable(tan(c + d*x^2)/x^2, x)]

    @test_int [x^3*(a + b*tan(c + d*x^2))^2, x, 10, (a^2*x^4)/4 + (im/2)*a*b*x^4 - (b^2*x^4)/4 - (a*b*x^2*log(1 + ℯ^((2*im)*(c + d*x^2))))/d + (b^2*log(cos(c + d*x^2)))/(2*d^2) + ((im/2)*a*b*PolyLog(2, -ℯ^((2*im)*(c + d*x^2))))/d^2 + (b^2*x^2*tan(c + d*x^2))/(2*d)]
    @test_int [x^2*(a + b*tan(c + d*x^2))^2, x, 0, Unintegrable(x^2*(a + b*tan(c + d*x^2))^2, x)]
    @test_int [x^1*(a + b*tan(c + d*x^2))^2, x, 3, ((a^2 - b^2)*x^2)/2 - (a*b*log(cos(c + d*x^2)))/d + (b^2*tan(c + d*x^2))/(2*d)]
    @test_int [x^0*(a + b*tan(c + d*x^2))^2, x, 0, Unintegrable((a + b*tan(c + d*x^2))^2, x)]
    @test_int [(a + b*tan(c + d*x^2))^2/x^1, x, 0, Unintegrable((a + b*tan(c + d*x^2))^2/x, x)]
    @test_int [(a + b*tan(c + d*x^2))^2/x^2, x, 0, Unintegrable((a + b*tan(c + d*x^2))^2/x^2, x)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [x^3/(a + b*tan(c + d*x^2)), x, 5, x^4/(4*(a + im*b)) + (b*x^2*log(1 + ((a^2 + b^2)*ℯ^(2*im*(c + d*x^2)))/(a + im*b)^2))/(2*(a^2 + b^2)*d) - (im*b*PolyLog(2, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^2)))/(a + im*b)^2)))/(4*(a^2 + b^2)*d^2)]
    @test_int [x^2/(a + b*tan(c + d*x^2)), x, 0, Unintegrable(x^2/(a + b*tan(c + d*x^2)), x)]
    @test_int [x^1/(a + b*tan(c + d*x^2)), x, 3, (a*x^2)/(2*(a^2 + b^2)) + (b*log(a*cos(c + d*x^2) + b*sin(c + d*x^2)))/(2*(a^2 + b^2)*d)]
    @test_int [x^0/(a + b*tan(c + d*x^2)), x, 0, Unintegrable((a + b*tan(c + d*x^2))^(-1), x)]
    @test_int [1/(x*(a + b*tan(c + d*x^2))), x, 0, Unintegrable(1/(x*(a + b*tan(c + d*x^2))), x)]
    @test_int [1/(x^2*(a + b*tan(c + d*x^2))), x, 0, Unintegrable(1/(x^2*(a + b*tan(c + d*x^2))), x)]

    @test_int [x^3/(a + b*tan(c + d*x^2))^2, x, 6, -(x^4/(4*(a^2 + b^2))) + (b + 2*a*d*x^2)^2/(8*a*(a + im*b)*(a^2 + b^2)*d^2) + (b*(b + 2*a*d*x^2)*log(1 + ((a^2 + b^2)*ℯ^(2*im*(c + d*x^2)))/(a + im*b)^2))/(2*(a^2 + b^2)^2*d^2) - (im*a*b*PolyLog(2, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^2)))/(a + im*b)^2)))/(2*(a^2 + b^2)^2*d^2) - (b*x^2)/(2*(a^2 + b^2)*d*(a + b*tan(c + d*x^2)))]
    @test_int [x^2/(a + b*tan(c + d*x^2))^2, x, 0, Unintegrable(x^2/(a + b*tan(c + d*x^2))^2, x)]
    @test_int [x^1/(a + b*tan(c + d*x^2))^2, x, 4, ((a^2 - b^2)*x^2)/(2*(a^2 + b^2)^2) + (a*b*log(a*cos(c + d*x^2) + b*sin(c + d*x^2)))/((a^2 + b^2)^2*d) - b/(2*(a^2 + b^2)*d*(a + b*tan(c + d*x^2)))]
    @test_int [x^0/(a + b*tan(c + d*x^2))^2, x, 0, Unintegrable((a + b*tan(c + d*x^2))^(-2), x)]
    @test_int [1/(x*(a + b*tan(c + d*x^2))^2), x, 0, Unintegrable(1/(x*(a + b*tan(c + d*x^2))^2), x)]
    @test_int [1/(x^2*(a + b*tan(c + d*x^2))^2), x, 0, Unintegrable(1/(x^2*(a + b*tan(c + d*x^2))^2), x)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*tan(c+d*x^(1/2)))^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [x^3*(a + b*tan(c + d*sqrt(x))), x, 13, (a*x^4)/4 + (im/4)*b*x^4 - (2*b*x^(7/2)*log(1 + ℯ^((2*im)*(c + d*sqrt(x)))))/d + ((7*im)*b*x^3*PolyLog(2, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^2 - (21*b*x^(5/2)*PolyLog(3, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^3 - (((105*im)/2)*b*x^2*PolyLog(4, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^4 + (105*b*x^(3/2)*PolyLog(5, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^5 + (((315*im)/2)*b*x*PolyLog(6, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^6 - (315*b*sqrt(x)*PolyLog(7, -ℯ^((2*im)*(c + d*sqrt(x)))))/(2*d^7) - (((315*im)/4)*b*PolyLog(8, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^8]
    @test_int [x^2*(a + b*tan(c + d*sqrt(x))), x, 11, (a*x^3)/3 + (im/3)*b*x^3 - (2*b*x^(5/2)*log(1 + ℯ^((2*im)*(c + d*sqrt(x)))))/d + ((5*im)*b*x^2*PolyLog(2, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^2 - (10*b*x^(3/2)*PolyLog(3, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^3 - ((15*im)*b*x*PolyLog(4, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^4 + (15*b*sqrt(x)*PolyLog(5, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^5 + (((15*im)/2)*b*PolyLog(6, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^6]
    @test_int [x^1*(a + b*tan(c + d*sqrt(x))), x, 9, (a*x^2)/2 + (im/2)*b*x^2 - (2*b*x^(3/2)*log(1 + ℯ^((2*im)*(c + d*sqrt(x)))))/d + ((3*im)*b*x*PolyLog(2, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^2 - (3*b*sqrt(x)*PolyLog(3, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^3 - (((3*im)/2)*b*PolyLog(4, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^4]
    @test_int [a + b*tan(c + d*sqrt(x)), x, 6, a*x + im*b*x - (2*b*sqrt(x)*log(1 + ℯ^((2*im)*(c + d*sqrt(x)))))/d + (im*b*PolyLog(2, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^2]
    @test_int [(a + b*tan(c + d*sqrt(x)))/x^1, x, 2, a*log(x) + b*Unintegrable(tan(c + d*sqrt(x))/x, x)]
    @test_int [(a + b*tan(c + d*sqrt(x)))/x^2, x, 2, -(a/x) + b*Unintegrable(tan(c + d*sqrt(x))/x^2, x)]

    @test_int [x^2*(a + b*tan(c + d*sqrt(x)))^2, x, 20, ((-2*im)*b^2*x^(5/2))/d + (a^2*x^3)/3 + ((2*im)/3)*a*b*x^3 - (b^2*x^3)/3 + (10*b^2*x^2*log(1 + ℯ^((2*im)*(c + d*sqrt(x)))))/d^2 - (4*a*b*x^(5/2)*log(1 + ℯ^((2*im)*(c + d*sqrt(x)))))/d - ((20*im)*b^2*x^(3/2)*PolyLog(2, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^3 + ((10*im)*a*b*x^2*PolyLog(2, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^2 + (30*b^2*x*PolyLog(3, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^4 - (20*a*b*x^(3/2)*PolyLog(3, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^3 + ((30*im)*b^2*sqrt(x)*PolyLog(4, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^5 - ((30*im)*a*b*x*PolyLog(4, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^4 - (15*b^2*PolyLog(5, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^6 + (30*a*b*sqrt(x)*PolyLog(5, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^5 + ((15*im)*a*b*PolyLog(6, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^6 + (2*b^2*x^(5/2)*tan(c + d*sqrt(x)))/d]
    @test_int [x^1*(a + b*tan(c + d*sqrt(x)))^2, x, 16, ((-2*im)*b^2*x^(3/2))/d + (a^2*x^2)/2 + im*a*b*x^2 - (b^2*x^2)/2 + (6*b^2*x*log(1 + ℯ^((2*im)*(c + d*sqrt(x)))))/d^2 - (4*a*b*x^(3/2)*log(1 + ℯ^((2*im)*(c + d*sqrt(x)))))/d - ((6*im)*b^2*sqrt(x)*PolyLog(2, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^3 + ((6*im)*a*b*x*PolyLog(2, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^2 + (3*b^2*PolyLog(3, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^4 - (6*a*b*sqrt(x)*PolyLog(3, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^3 - ((3*im)*a*b*PolyLog(4, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^4 + (2*b^2*x^(3/2)*tan(c + d*sqrt(x)))/d]
    @test_int [x^0*(a + b*tan(c + d*sqrt(x)))^2, x, 10, a^2*x + (2*im)*a*b*x - b^2*x - (4*a*b*sqrt(x)*log(1 + ℯ^((2*im)*(c + d*sqrt(x)))))/d + (2*b^2*log(cos(c + d*sqrt(x))))/d^2 + ((2*im)*a*b*PolyLog(2, -ℯ^((2*im)*(c + d*sqrt(x)))))/d^2 + (2*b^2*sqrt(x)*tan(c + d*sqrt(x)))/d]
    @test_int [(a + b*tan(c + d*sqrt(x)))^2/x^1, x, 0, Unintegrable((a + b*tan(c + d*sqrt(x)))^2/x, x)]
    @test_int [(a + b*tan(c + d*sqrt(x)))^2/x^2, x, 0, Unintegrable((a + b*tan(c + d*sqrt(x)))^2/x^2, x)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [x^3/(a + b*tan(c + d*sqrt(x))), x, 11, x^4/(4*(a + im*b)) + (2*b*x^(7/2)*log(1 + ((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2))/((a^2 + b^2)*d) - (7*im*b*x^3*PolyLog(2, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/((a^2 + b^2)*d^2) + (21*b*x^(5/2)*PolyLog(3, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/((a^2 + b^2)*d^3) + (105*im*b*x^2*PolyLog(4, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/(2*(a^2 + b^2)*d^4) - (105*b*x^(3/2)*PolyLog(5, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/((a^2 + b^2)*d^5) - (315*im*b*x*PolyLog(6, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/(2*(a^2 + b^2)*d^6) + (315*b*sqrt(x)*PolyLog(7, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/(2*(a^2 + b^2)*d^7) + (315*im*b*PolyLog(8, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/(4*(a^2 + b^2)*d^8)]
    @test_int [x^2/(a + b*tan(c + d*sqrt(x))), x, 9, x^3/(3*(a + im*b)) + (2*b*x^(5/2)*log(1 + ((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2))/((a^2 + b^2)*d) - (5*im*b*x^2*PolyLog(2, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/((a^2 + b^2)*d^2) + (10*b*x^(3/2)*PolyLog(3, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/((a^2 + b^2)*d^3) + (15*im*b*x*PolyLog(4, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/((a^2 + b^2)*d^4) - (15*b*sqrt(x)*PolyLog(5, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/((a^2 + b^2)*d^5) - (15*im*b*PolyLog(6, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/(2*(a^2 + b^2)*d^6)]
    @test_int [x^1/(a + b*tan(c + d*sqrt(x))), x, 7, x^2/(2*(a + im*b)) + (2*b*x^(3/2)*log(1 + ((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2))/((a^2 + b^2)*d) - (3*im*b*x*PolyLog(2, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/((a^2 + b^2)*d^2) + (3*b*sqrt(x)*PolyLog(3, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/((a^2 + b^2)*d^3) + (3*im*b*PolyLog(4, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/(2*(a^2 + b^2)*d^4)]
    @test_int [x^0/(a + b*tan(c + d*sqrt(x))), x, 5, x/(a + im*b) + (2*b*sqrt(x)*log(1 + ((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2))/((a^2 + b^2)*d) - (im*b*PolyLog(2, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/((a^2 + b^2)*d^2)]
    @test_int [1/(x*(a + b*tan(c + d*sqrt(x)))), x, 0, Unintegrable(1/(x*(a + b*tan(c + d*sqrt(x)))), x)]
    @test_int [1/(x^2*(a + b*tan(c + d*sqrt(x)))), x, 0, Unintegrable(1/(x^2*(a + b*tan(c + d*sqrt(x)))), x)]

    @test_int [x^2/(a + b*tan(c + d*sqrt(x)))^2, x, 28, -((4*im*b^2*x^(5/2))/((a^2 + b^2)^2*d)) + (4*b^2*x^(5/2))/((a + im*b)*(im*a + b)^2*d*(im*a - b + (im*a + b)*ℯ^(2*im*(c + d*sqrt(x))))) + x^3/(3*(a - im*b)^2) + (4*b*x^3)/(3*(im*a - b)*(a - im*b)^2) - (4*b^2*x^3)/(3*(a^2 + b^2)^2) + (10*b^2*x^2*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)))/((a^2 + b^2)^2*d^2) + (4*b*x^(5/2)*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)))/((a - im*b)^2*(a + im*b)*d) - (4*im*b^2*x^(5/2)*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)))/((a^2 + b^2)^2*d) - (20*im*b^2*x^(3/2)*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a^2 + b^2)^2*d^3) + (10*b*x^2*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((im*a - b)*(a - im*b)^2*d^2) - (10*b^2*x^2*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a^2 + b^2)^2*d^2) + (30*b^2*x*PolyLog(3, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a^2 + b^2)^2*d^4) + (20*b*x^(3/2)*PolyLog(3, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a - im*b)^2*(a + im*b)*d^3) - (20*im*b^2*x^(3/2)*PolyLog(3, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a^2 + b^2)^2*d^3) + (30*im*b^2*sqrt(x)*PolyLog(4, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a^2 + b^2)^2*d^5) - (30*b*x*PolyLog(4, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((im*a - b)*(a - im*b)^2*d^4) + (30*b^2*x*PolyLog(4, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a^2 + b^2)^2*d^4) - (15*b^2*PolyLog(5, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a^2 + b^2)^2*d^6) - (30*b*sqrt(x)*PolyLog(5, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a - im*b)^2*(a + im*b)*d^5) + (30*im*b^2*sqrt(x)*PolyLog(5, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a^2 + b^2)^2*d^5) + (15*b*PolyLog(6, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((im*a - b)*(a - im*b)^2*d^6) - (15*b^2*PolyLog(6, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a^2 + b^2)^2*d^6)]
    @test_int [x^1/(a + b*tan(c + d*sqrt(x)))^2, x, 22, -((4*im*b^2*x^(3/2))/((a^2 + b^2)^2*d)) + (4*b^2*x^(3/2))/((a + im*b)*(im*a + b)^2*d*(im*a - b + (im*a + b)*ℯ^(2*im*(c + d*sqrt(x))))) + x^2/(2*(a - im*b)^2) + (2*b*x^2)/((im*a - b)*(a - im*b)^2) - (2*b^2*x^2)/(a^2 + b^2)^2 + (6*b^2*x*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)))/((a^2 + b^2)^2*d^2) + (4*b*x^(3/2)*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)))/((a - im*b)^2*(a + im*b)*d) - (4*im*b^2*x^(3/2)*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)))/((a^2 + b^2)^2*d) - (6*im*b^2*sqrt(x)*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a^2 + b^2)^2*d^3) + (6*b*x*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((im*a - b)*(a - im*b)^2*d^2) - (6*b^2*x*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a^2 + b^2)^2*d^2) + (3*b^2*PolyLog(3, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a^2 + b^2)^2*d^4) + (6*b*sqrt(x)*PolyLog(3, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a - im*b)^2*(a + im*b)*d^3) - (6*im*b^2*sqrt(x)*PolyLog(3, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a^2 + b^2)^2*d^3) - (3*b*PolyLog(4, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((im*a - b)*(a - im*b)^2*d^4) + (3*b^2*PolyLog(4, -(((a - im*b)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b))))/((a^2 + b^2)^2*d^4)]
    @test_int [x^0/(a + b*tan(c + d*sqrt(x)))^2, x, 6, (b + 2*a*d*sqrt(x))^2/(2*a*(a + im*b)*(a^2 + b^2)*d^2) - x/(a^2 + b^2) + (2*b*(b + 2*a*d*sqrt(x))*log(1 + ((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2))/((a^2 + b^2)^2*d^2) - (2*im*a*b*PolyLog(2, -(((a^2 + b^2)*ℯ^(2*im*(c + d*sqrt(x))))/(a + im*b)^2)))/((a^2 + b^2)^2*d^2) - (2*b*sqrt(x))/((a^2 + b^2)*d*(a + b*tan(c + d*sqrt(x))))]
    @test_int [1/(x*(a + b*tan(c + d*sqrt(x)))^2), x, 0, Unintegrable(1/(x*(a + b*tan(c + d*sqrt(x)))^2), x)]
    @test_int [1/(x^2*(a + b*tan(c + d*sqrt(x)))^2), x, 0, Unintegrable(1/(x^2*(a + b*tan(c + d*sqrt(x)))^2), x)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*tan(c+d*x^(1/3)))^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [x^2*(a + b*tan(c + d*x^(1/3))), x, 14, (a*x^3)/3 + (im/3)*b*x^3 - (3*b*x^(8/3)*log(1 + ℯ^((2*im)*(c + d*x^(1/3)))))/d + ((12*im)*b*x^(7/3)*PolyLog(2, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^2 - (42*b*x^2*PolyLog(3, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^3 - ((126*im)*b*x^(5/3)*PolyLog(4, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^4 + (315*b*x^(4/3)*PolyLog(5, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^5 + ((630*im)*b*x*PolyLog(6, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^6 - (945*b*x^(2/3)*PolyLog(7, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^7 - ((945*im)*b*x^(1/3)*PolyLog(8, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^8 + (945*b*PolyLog(9, -ℯ^((2*im)*(c + d*x^(1/3)))))/(2*d^9)]
    @test_int [x^1*(a + b*tan(c + d*x^(1/3))), x, 11, (a*x^2)/2 + (im/2)*b*x^2 - (3*b*x^(5/3)*log(1 + ℯ^((2*im)*(c + d*x^(1/3)))))/d + (((15*im)/2)*b*x^(4/3)*PolyLog(2, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^2 - (15*b*x*PolyLog(3, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^3 - (((45*im)/2)*b*x^(2/3)*PolyLog(4, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^4 + (45*b*x^(1/3)*PolyLog(5, -ℯ^((2*im)*(c + d*x^(1/3)))))/(2*d^5) + (((45*im)/4)*b*PolyLog(6, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^6]
    @test_int [a + b*tan(c + d*x^(1/3)), x, 7, a*x + im*b*x - (3*b*x^(2/3)*log(1 + ℯ^((2*im)*(c + d*x^(1/3)))))/d + ((3*im)*b*x^(1/3)*PolyLog(2, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^2 - (3*b*PolyLog(3, -ℯ^((2*im)*(c + d*x^(1/3)))))/(2*d^3)]
    @test_int [(a + b*tan(c + d*x^(1/3)))/x^1, x, 2, a*log(x) + b*Unintegrable(tan(c + d*x^(1/3))/x, x)]
    @test_int [(a + b*tan(c + d*x^(1/3)))/x^2, x, 2, -(a/x) + b*Unintegrable(tan(c + d*x^(1/3))/x^2, x)]

    @test_int [x^2*(a + b*tan(c + d*x^(1/3)))^2, x, 26, ((-3*im)*b^2*x^(8/3))/d + (a^2*x^3)/3 + ((2*im)/3)*a*b*x^3 - (b^2*x^3)/3 + (24*b^2*x^(7/3)*log(1 + ℯ^((2*im)*(c + d*x^(1/3)))))/d^2 - (6*a*b*x^(8/3)*log(1 + ℯ^((2*im)*(c + d*x^(1/3)))))/d - ((84*im)*b^2*x^2*PolyLog(2, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^3 + ((24*im)*a*b*x^(7/3)*PolyLog(2, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^2 + (252*b^2*x^(5/3)*PolyLog(3, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^4 - (84*a*b*x^2*PolyLog(3, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^3 + ((630*im)*b^2*x^(4/3)*PolyLog(4, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^5 - ((252*im)*a*b*x^(5/3)*PolyLog(4, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^4 - (1260*b^2*x*PolyLog(5, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^6 + (630*a*b*x^(4/3)*PolyLog(5, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^5 - ((1890*im)*b^2*x^(2/3)*PolyLog(6, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^7 + ((1260*im)*a*b*x*PolyLog(6, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^6 + (1890*b^2*x^(1/3)*PolyLog(7, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^8 - (1890*a*b*x^(2/3)*PolyLog(7, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^7 + ((945*im)*b^2*PolyLog(8, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^9 - ((1890*im)*a*b*x^(1/3)*PolyLog(8, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^8 + (945*a*b*PolyLog(9, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^9 + (3*b^2*x^(8/3)*tan(c + d*x^(1/3)))/d]
    @test_int [x^1*(a + b*tan(c + d*x^(1/3)))^2, x, 20, ((-3*im)*b^2*x^(5/3))/d + (a^2*x^2)/2 + im*a*b*x^2 - (b^2*x^2)/2 + (15*b^2*x^(4/3)*log(1 + ℯ^((2*im)*(c + d*x^(1/3)))))/d^2 - (6*a*b*x^(5/3)*log(1 + ℯ^((2*im)*(c + d*x^(1/3)))))/d - ((30*im)*b^2*x*PolyLog(2, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^3 + ((15*im)*a*b*x^(4/3)*PolyLog(2, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^2 + (45*b^2*x^(2/3)*PolyLog(3, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^4 - (30*a*b*x*PolyLog(3, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^3 + ((45*im)*b^2*x^(1/3)*PolyLog(4, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^5 - ((45*im)*a*b*x^(2/3)*PolyLog(4, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^4 - (45*b^2*PolyLog(5, -ℯ^((2*im)*(c + d*x^(1/3)))))/(2*d^6) + (45*a*b*x^(1/3)*PolyLog(5, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^5 + (((45*im)/2)*a*b*PolyLog(6, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^6 + (3*b^2*x^(5/3)*tan(c + d*x^(1/3)))/d]
    @test_int [x^0*(a + b*tan(c + d*x^(1/3)))^2, x, 14, ((-3*im)*b^2*x^(2/3))/d + a^2*x + (2*im)*a*b*x - b^2*x + (6*b^2*x^(1/3)*log(1 + ℯ^((2*im)*(c + d*x^(1/3)))))/d^2 - (6*a*b*x^(2/3)*log(1 + ℯ^((2*im)*(c + d*x^(1/3)))))/d - ((3*im)*b^2*PolyLog(2, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^3 + ((6*im)*a*b*x^(1/3)*PolyLog(2, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^2 - (3*a*b*PolyLog(3, -ℯ^((2*im)*(c + d*x^(1/3)))))/d^3 + (3*b^2*x^(2/3)*tan(c + d*x^(1/3)))/d]
    @test_int [(a + b*tan(c + d*x^(1/3)))^2/x^1, x, 0, Unintegrable((a + b*tan(c + d*x^(1/3)))^2/x, x)]
    @test_int [(a + b*tan(c + d*x^(1/3)))^2/x^2, x, 0, Unintegrable((a + b*tan(c + d*x^(1/3)))^2/x^2, x)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [x^2/(a + b*tan(c + d*x^(1/3))), x, 12, x^3/(3*(a + im*b)) + (3*b*x^(8/3)*log(1 + ((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2))/((a^2 + b^2)*d) - (12*im*b*x^(7/3)*PolyLog(2, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/((a^2 + b^2)*d^2) + (42*b*x^2*PolyLog(3, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/((a^2 + b^2)*d^3) + (126*im*b*x^(5/3)*PolyLog(4, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/((a^2 + b^2)*d^4) - (315*b*x^(4/3)*PolyLog(5, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/((a^2 + b^2)*d^5) - (630*im*b*x*PolyLog(6, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/((a^2 + b^2)*d^6) + (945*b*x^(2/3)*PolyLog(7, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/((a^2 + b^2)*d^7) + (945*im*b*x^(1/3)*PolyLog(8, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/((a^2 + b^2)*d^8) - (945*b*PolyLog(9, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/(2*(a^2 + b^2)*d^9)]
    @test_int [x^1/(a + b*tan(c + d*x^(1/3))), x, 9, x^2/(2*(a + im*b)) + (3*b*x^(5/3)*log(1 + ((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2))/((a^2 + b^2)*d) - (15*im*b*x^(4/3)*PolyLog(2, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/(2*(a^2 + b^2)*d^2) + (15*b*x*PolyLog(3, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/((a^2 + b^2)*d^3) + (45*im*b*x^(2/3)*PolyLog(4, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/(2*(a^2 + b^2)*d^4) - (45*b*x^(1/3)*PolyLog(5, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/(2*(a^2 + b^2)*d^5) - (45*im*b*PolyLog(6, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/(4*(a^2 + b^2)*d^6)]
    @test_int [x^0/(a + b*tan(c + d*x^(1/3))), x, 6, x/(a + im*b) + (3*b*x^(2/3)*log(1 + ((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2))/((a^2 + b^2)*d) - (3*im*b*x^(1/3)*PolyLog(2, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/((a^2 + b^2)*d^2) + (3*b*PolyLog(3, -(((a^2 + b^2)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)^2)))/(2*(a^2 + b^2)*d^3)]
    @test_int [1/(x*(a + b*tan(c + d*x^(1/3)))), x, 0, Unintegrable(1/(x*(a + b*tan(c + d*x^(1/3)))), x)]
    @test_int [1/(x^2*(a + b*tan(c + d*x^(1/3)))), x, 0, Unintegrable(1/(x^2*(a + b*tan(c + d*x^(1/3)))), x)]

    @test_int [x^2/(a + b*tan(c + d*x^(1/3)))^2, x, 37, -((6*im*b^2*x^(8/3))/((a^2 + b^2)^2*d)) + (6*b^2*x^(8/3))/((a + im*b)*(im*a + b)^2*d*(im*a - b + (im*a + b)*ℯ^(2*im*(c + d*x^(1/3))))) + x^3/(3*(a - im*b)^2) + (4*b*x^3)/(3*(im*a - b)*(a - im*b)^2) - (4*b^2*x^3)/(3*(a^2 + b^2)^2) + (24*b^2*x^(7/3)*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)))/((a^2 + b^2)^2*d^2) + (6*b*x^(8/3)*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)))/((a - im*b)^2*(a + im*b)*d) - (6*im*b^2*x^(8/3)*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)))/((a^2 + b^2)^2*d) - (84*im*b^2*x^2*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^3) + (24*b*x^(7/3)*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((im*a - b)*(a - im*b)^2*d^2) - (24*b^2*x^(7/3)*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^2) + (252*b^2*x^(5/3)*PolyLog(3, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^4) + (84*b*x^2*PolyLog(3, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a - im*b)^2*(a + im*b)*d^3) - (84*im*b^2*x^2*PolyLog(3, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^3) + (630*im*b^2*x^(4/3)*PolyLog(4, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^5) - (252*b*x^(5/3)*PolyLog(4, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((im*a - b)*(a - im*b)^2*d^4) + (252*b^2*x^(5/3)*PolyLog(4, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^4) - (1260*b^2*x*PolyLog(5, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^6) - (630*b*x^(4/3)*PolyLog(5, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a - im*b)^2*(a + im*b)*d^5) + (630*im*b^2*x^(4/3)*PolyLog(5, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^5) - (1890*im*b^2*x^(2/3)*PolyLog(6, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^7) + (1260*b*x*PolyLog(6, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((im*a - b)*(a - im*b)^2*d^6) - (1260*b^2*x*PolyLog(6, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^6) + (1890*b^2*x^(1/3)*PolyLog(7, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^8) + (1890*b*x^(2/3)*PolyLog(7, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a - im*b)^2*(a + im*b)*d^7) - (1890*im*b^2*x^(2/3)*PolyLog(7, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^7) + (945*im*b^2*PolyLog(8, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^9) - (1890*b*x^(1/3)*PolyLog(8, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((im*a - b)*(a - im*b)^2*d^8) + (1890*b^2*x^(1/3)*PolyLog(8, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^8) - (945*b*PolyLog(9, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a - im*b)^2*(a + im*b)*d^9) + (945*im*b^2*PolyLog(9, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^9)]
    @test_int [x^1/(a + b*tan(c + d*x^(1/3)))^2, x, 28, -((6*im*b^2*x^(5/3))/((a^2 + b^2)^2*d)) + (6*b^2*x^(5/3))/((a + im*b)*(im*a + b)^2*d*(im*a - b + (im*a + b)*ℯ^(2*im*(c + d*x^(1/3))))) + x^2/(2*(a - im*b)^2) + (2*b*x^2)/((im*a - b)*(a - im*b)^2) - (2*b^2*x^2)/(a^2 + b^2)^2 + (15*b^2*x^(4/3)*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)))/((a^2 + b^2)^2*d^2) + (6*b*x^(5/3)*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)))/((a - im*b)^2*(a + im*b)*d) - (6*im*b^2*x^(5/3)*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)))/((a^2 + b^2)^2*d) - (30*im*b^2*x*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^3) + (15*b*x^(4/3)*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((im*a - b)*(a - im*b)^2*d^2) - (15*b^2*x^(4/3)*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^2) + (45*b^2*x^(2/3)*PolyLog(3, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^4) + (30*b*x*PolyLog(3, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a - im*b)^2*(a + im*b)*d^3) - (30*im*b^2*x*PolyLog(3, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^3) + (45*im*b^2*x^(1/3)*PolyLog(4, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^5) - (45*b*x^(2/3)*PolyLog(4, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((im*a - b)*(a - im*b)^2*d^4) + (45*b^2*x^(2/3)*PolyLog(4, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^4) - (45*b^2*PolyLog(5, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/(2*(a^2 + b^2)^2*d^6) - (45*b*x^(1/3)*PolyLog(5, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a - im*b)^2*(a + im*b)*d^5) + (45*im*b^2*x^(1/3)*PolyLog(5, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^5) + (45*b*PolyLog(6, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/(2*(im*a - b)*(a - im*b)^2*d^6) - (45*b^2*PolyLog(6, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/(2*(a^2 + b^2)^2*d^6)]
    @test_int [x^0/(a + b*tan(c + d*x^(1/3)))^2, x, 19, -((6*im*b^2*x^(2/3))/((a^2 + b^2)^2*d)) + (6*b^2*x^(2/3))/((a + im*b)*(im*a + b)^2*d*(im*a - b + (im*a + b)*ℯ^(2*im*(c + d*x^(1/3))))) + x/(a - im*b)^2 + (4*b*x)/((im*a - b)*(a - im*b)^2) - (4*b^2*x)/(a^2 + b^2)^2 + (6*b^2*x^(1/3)*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)))/((a^2 + b^2)^2*d^2) + (6*b*x^(2/3)*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)))/((a - im*b)^2*(a + im*b)*d) - (6*im*b^2*x^(2/3)*log(1 + ((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b)))/((a^2 + b^2)^2*d) - (3*im*b^2*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^3) + (6*b*x^(1/3)*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((im*a - b)*(a - im*b)^2*d^2) - (6*b^2*x^(1/3)*PolyLog(2, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^2) + (3*b*PolyLog(3, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a - im*b)^2*(a + im*b)*d^3) - (3*im*b^2*PolyLog(3, -(((a - im*b)*ℯ^(2*im*(c + d*x^(1/3))))/(a + im*b))))/((a^2 + b^2)^2*d^3)]
    @test_int [1/(x*(a + b*tan(c + d*x^(1/3)))^2), x, 0, Unintegrable(1/(x*(a + b*tan(c + d*x^(1/3)))^2), x)]
    @test_int [1/(x^2*(a + b*tan(c + d*x^(1/3)))^2), x, 0, Unintegrable(1/(x^2*(a + b*tan(c + d*x^(1/3)))^2), x)]
    end

max_stars_repo_path: test/4 Trig functions/4.3 Tangent/4.3.11 (e x)^m (a+b tan(c+d x^n))^p.jl


text:
#**************************************************************************************
# test_Output_Interface.jl
# ======================== part of the GeoEfficiency.jl package.
#
#
#
#**************************************************************************************

@testset "Output Interface" begin
	@test G.checkResultsDirs() == nothing

	pnt = Point(1)
	cylDet = Detector(5, 10); wellDet = Detector(5, 4, 3.2, 2)
	@testset "function `calc` on CylDetector, WellDetector" for
	SrcRadius = Real[1, 1//2, e, pi, 1.0],
	SrcLength = Real[1, 1//2, e, pi, 1.0]

		@test calc(cylDet,  (pnt, SrcRadius , SrcLength))    == nothing
		@test calc(wellDet, (pnt, SrcRadius , SrcLength))    == nothing
     end #testset

	@testset "function `calcN`" for
	cylDet = [Detector(5,10), Detector(eps(),0)], #, wellDet= Detector(5, 4, 3, 2)
	consol_input = ["4 0 1 2 ", "4 0 1 2 " * "d " * "4 0 1 2 ", "4 0 1 2 " * "n " * "10 5 0 " * "4 0 1 2 "]
		@test exec_consol_unattended(calcN, consol_input, Fn_ARGs =[cylDet])  == nothing
		@test exec_consol_unattended(calcN, "10 5 0 " * consol_input)      == nothing
 	end #testset

info("test `_batch` & `batch`...")
	@testset "function `batch`" begin
		@test G._batch(Val{true},  CylDetector(eps(0.1)), [0.0], [0.0], [0.0], [0.0])[2][end] ≈ 0.5
		@test G._batch(Val{false}, CylDetector(eps(0.2)), [0.0], [0.0], [0.0], [0.0])[2][end] ≈ 0.5
		@test isnan(G._batch(Val{true}, CylDetector(eps(0.3)), [0.0], [1.0], [0.0],[0.0])[2][end])
		@test isnan(G._batch(Val{false}, CylDetector(eps(0.4)), [0.0], [1.0], [0.0],[0.0])[2][end])

		acylDetector = CylDetector(eps(0.5)); path = batch(acylDetector, [0.0])
		@test contains(path, G.id(acylDetector))
		@test readcsv(path)[2,end] ≈ 0.5

		path = batch(acylDetector, [0.0], [0.0], [0.0],[0.0],false)
		@test contains(path, G.id(acylDetector))
		@test readcsv(path)[2,end] ≈ 0.5

		path = batch([acylDetector], [0.0])
		@test contains(contains ,path, G.id(acylDetector))
		every_path::Vector{String} = path

		path = batch([acylDetector], [0.0], [0.0], [0.0],[0.0],false)
		@test contains(contains ,path, G.id(acylDetector))
		append!(every_path, path)

		aBDetector = BoreDetector(eps(0.5), eps(0.4), eps(0.2)); path = batch([aBDetector], [0.0])
		@test contains(contains ,path, G.id(aBDetector))
		append!(every_path, path)

		aWDetector = WellDetector(eps(0.5), eps(0.4), eps(0.2), eps(0.1)); path = batch([aWDetector], [0.0])
		@test contains(contains, path, G.id(aWDetector))
		append!(every_path, path)

		path = batch([aWDetector], [0.0], [0.0], [0.0],[0.0],false)
		@test contains(contains ,path, G.id(aWDetector))
		append!(every_path, path)

		path = batch([acylDetector, aWDetector], [0.0])
		@test contains(contains ,path, G.id(aWDetector))
		append!(every_path, path)

		path = batch([acylDetector, aWDetector], [0.0], [0.0], [0.0],[0.0],false)
		@test contains(contains ,path, G.id(acylDetector))
		@test contains(contains ,path, G.id(aWDetector))
		append!(every_path, path)

		path = batch([acylDetector, aBDetector], [0.0])
		@test contains(contains ,path, G.id(acylDetector))
		append!(every_path, path)

		path = batch([acylDetector, aBDetector], [0.0], [0.0], [0.0],[0.0],false)
		@test contains(contains ,path, G.id(acylDetector))
		@test contains(contains ,path, G.id(aBDetector))
		append!(every_path, path)

		path = batch([aBDetector, aWDetector], [0.0])
		@test contains(contains ,path, G.id(aBDetector))
		@test contains(contains ,path, G.id(aWDetector))
		append!(every_path, path)

		path = batch([aBDetector, aWDetector], [0.0], [0.0], [0.0],[0.0],false)
		@test contains(contains ,path, G.id(aBDetector))
		@test contains(contains ,path, G.id(aWDetector))
		append!(every_path, path)

		acylDetector = CylDetector(eps(0.6)); path = batch([acylDetector, aBDetector, aWDetector], [0.0])
		@test contains(contains ,path, G.id(acylDetector))
		@test contains(contains ,path, G.id(aBDetector))
		@test contains(contains ,path, G.id(aWDetector))
	chmod(path[1], 0o100444); chmod(path[2], 0o100444); chmod(path[3], 0o100444);
	append!(every_path, path)
		path = batch([acylDetector, aBDetector, aWDetector], [0.0])
		@test contains(contains ,path, "_" * G.id(acylDetector))
		@test contains(contains ,path, "_" * G.id(aBDetector))
		@test contains(contains ,path, "_" * G.id(aWDetector))
	append!(every_path, path)
	chmod(path[1],0o777); chmod(path[2], 0o777); chmod(path[3], 0o777); # 0o100666

		acylDetector = CylDetector(eps(0.7)); path = batch([acylDetector, aBDetector, aWDetector], [0.0], [0.0], [0.0],[0.0],false)
		@test contains(contains ,path, G.id(acylDetector))
		@test contains(contains ,path, G.id(aBDetector))
		@test contains(contains ,path, G.id(aWDetector))
	append!(every_path, path)
	chmod(path[1], 0o100444); chmod(path[2], 0o100444); chmod(path[3], 0o100444);
		path = batch([acylDetector, aBDetector, aWDetector], [0.0], [0.0], [0.0],[0.0],false)
		@test contains(contains ,path, "_" * G.id(acylDetector))
		@test contains(contains ,path, "_" * G.id(aBDetector))
		@test contains(contains ,path, "_" * G.id(aWDetector))
	append!(every_path, path)
	chmod(path[1],0o777); chmod(path[2], 0o777); chmod(path[3], 0o777);

		@test append!(every_path, batch([eps() 0], [0.0]))|> eltype === String
		@test append!(every_path, batch([eps() 0], [0.0], [0.0], [0.0], [0.0],false)) |> eltype === String
		@test append!(every_path, batch([1.0 0], [0.0]))|> eltype === String
		@test append!(every_path, batch([1.0 0], [0.0], [0.0], [0.0], [0.0],false)) |> eltype === String
		@test append!(every_path, batch([1//2 0.0], [0.0]))|> eltype === String
		@test append!(every_path, batch([1//2 0.0], [0.0], [0.0], [0.0], [0.0],false)) |> eltype === String
		@test append!(every_path, batch([1//2 0.0], [0.0]))|> eltype === String
		@test append!(every_path, batch([1//2 0.0], [0.0], [0.0], [0.0], [0.0],false)) |> eltype === String
		@test append!(every_path, batch([e pi], [0.0]))|> eltype === String
		@test append!(every_path, batch([e pi], [0.0], [0.0], [0.0], [0.0],false)) |> eltype === String

		@test append!(every_path, batch([5.0 4 3], [0.0]))|> eltype === String
		@test append!(every_path, batch([5.0 4 3], [0.0], [0.0],[0.0],[0.0],false)) |> eltype === String
		@test append!(every_path, batch([5.0 4 3//1], [0.0]))|> eltype === String
		@test append!(every_path, batch([5.0 4 3//1], [0.0], [0.0],[0.0],[0.0],false)) |> eltype === String
		@test append!(every_path, batch([5.0 4 pi], [0.0]))|> eltype === String
		@test append!(every_path, batch([5.0 4 pi], [0.0], [0.0],[0.0], [0.0], false)) |> eltype === String

		@test append!(every_path, batch([5.0 4 3 2], [0.0]))|> eltype === String
		@test append!(every_path, batch([5.0 4 3 2], [0.0], [0.0], [0.0],[0.0],false))|> eltype === String

		@test append!(every_path, batch([acylDetector, aWDetector], [0.0]))|> eltype === String
		@test append!(every_path, batch([acylDetector, aWDetector], [0.0], [0.0], [0.0],[0.0],false))|> eltype === String

		@test append!(every_path, batch([acylDetector, aBDetector], [0.0]))|> eltype === String
		@test append!(every_path, batch([acylDetector, aBDetector], [0.0], [0.0], [0.0],[0.0],false))|> eltype === String

		@test append!(every_path, batch([aBDetector, aWDetector], [0.0]))|> eltype === String
		@test append!(every_path, batch([aBDetector, aWDetector], [0.0], [0.0], [0.0],[0.0],false))|> eltype === String

		@test append!(every_path, batch([acylDetector, aBDetector, aWDetector], [0.0]))|> eltype === String
		@test append!(every_path, batch([acylDetector, aBDetector, aWDetector], [0.0], [0.0], [0.0],[0.0],false))|> eltype === String

		try
		G.detector_info_from_csvFile()
		if  [0.0] != G.read_from_csvFile(G.srcHeights, G.datadir)
			setSrcToPoint(true);
			@test append!(every_path, batch())|> eltype === String

			if [0.0] != G.read_from_csvFile(G.srcRadii, G.datadir)
				setSrcToPoint(false);
				@test append!(every_path, batch())|> eltype === String
			end
		end
	end

	rm.(every_path,  force=true)
	#rm.(batch([aWDetector], [0.0]))
	#rm.(batch([aWDetector], [0.0], [0.0], [0.0],[0.0],false))
	for cr = 0.0:0.1:0.7
		rm.(batch([Detector(eps(cr))], [0.0], [0.0], [0.0],[0.0],false) ; force=true)
		rm.(batch([Detector(eps(cr))], [0.0]) ; force=true)
	end
	end  #begin_testset
println()
end  #begin_testset

ext: jl


text:
# This file shows the proposed expansion of Protocol macro

# Protocol Iterable
#   start::Int64
#
#   function next(t :~ Iterable) # no default implementation
#
#   function secondnext(t :~ Iterable) # default implementation
#     next(t)
#     next(t)
#   end
# end

struct Iterable end
struct NotIterable end

is_Iterable(::Any) = NotIterable()

function next(::Iterable, t)
  error("Not implemented")
end

function secondnext(::Iterable, t) # default implementation
  next(Iterable, t)
  next(Iterable, t)
end

next(t) = next(is_iterable(t), t)
secondnext(t) = secondnext(is_iterable(t), t)

next(::NotIterable, t) = println("Argument does not implement or declare Protocol Iterable")
secondnext(::NotIterable, t) = println("Argument does not implement or declare Protocol Iterable")

# We now define a type...
struct T
  start::Int64
  c
end

# impl Iterable for T
#
#   function next(t::T)
#     t.c
#   end
#
# end

if ! hasfield(T, :start)
  throw(ProtocolError("Type T has no field start"))
end

if ! hasfieldtype(T, :start, Int64)  # Not sure how to implement this...
  throw(ProtocolError("Field T.start needs to be of type Int64"))
end

# Taken directly from the body of the implementation
function next(::Iterable, t::T)
  t.c
end

function secondnext(::Iterable, t::T) # default implementation taken from the Protocol
  next(t)
  next(t)
end

# Check that the body of the implementation does what it should:
if ! hasmethod(next, Tuple{Iterable, T})
  throw(ProtocolError("Protocol function next not implemented for Type T"))
end

is_Iterable(::T) = Iterable()

# function get_third(b :~ Iterable)
#   secondnext(b)
#   next(b)
# end

function get_third(::Iterable, b)
  secondnext(b)
  next(b)
end
get_third(b) = get_third(is_Iterable(b), b)
get_third(::NotIterable, b) = throw(ProtocolError("First argument doesn't implement Iterable."))

# we get an object of type T
t = T(1, 5)

# and call get_third on it:
get_third(t)

# and calling it this way:
get_third(1)
# raises "Type of argument arr does not implement protocol Iterable"


text:
@reexport module _Combinatorics

export combinations, multinomial, nthperm, permutations

import Combinatorics

using Brainstorm: flatten
using Nemo: QQ, ZZ, PowerSeriesRing, coeff, divexact, fac
using Combinatorics: combinations, factorial, nthperm, permutations
using Multicombinations: multicombinations

# Multinomial coefficient where n = sum(k)
# https://github.com/jiahao/Combinatorics.jl/blob/master/src/Combinatorics.jl
multinomial(k) = begin
    s = big(0)
    result = big(1)
    for i in k
        s += i
        result *= binomial(s, i)
    end
    result
end

include("permutation.jl")
include("combination.jl")

include("multipermutation.jl")
include("multicombination.jl")

include("subset.jl")

include("integerpartition.jl")
include("integersum.jl")

end


text:
# Min Trust LAMB (https://github.com/Smerity/pytorch-lamb/blob/master/pytorch_lamb/lamb.py)
#
# It has been proposed in Large Batch Optimization for Deep Learning: Training BERT in 76 minutes.
# https://arxiv.org/abs/1904.00962
#
# Adapted from BertAdam https://github.com/OsmanMutlu/BERT.jl/blob/master/src/optimizer.jl

import Knet: update!
using AutoGrad: full
using LinearAlgebra
using CUDA

warmup_cosine(x, warmup=0.002) = x < warmup ? x/warmup : 0.5 * (1.0 + cos(π * x))
warmup_constant(x, warmup=0.002) = x < warmup ? x/warmup : 1.0
warmup_linear(x, warmup=0.002) = x < warmup ? x/warmup : 1.0 - x

mutable struct MinTrustLAMB
    lr::AbstractFloat
    beta1::AbstractFloat
    beta2::AbstractFloat
    eps::AbstractFloat
    w_decay::AbstractFloat
    min_trust::AbstractFloat
    gclip::AbstractFloat
    t::Int
    schedule
    warmup
    t_total
    fstm
    scndm
end

MinTrustLAMB(; min_trust=0.25, lr=1e-3, gclip=0.25, beta1=0.9, beta2=0.999, eps=1e-6, w_decay=0.0, schedule="warmup_constant", warmup=-1, t_total=-1) = MinTrustLAMB(lr, beta1, beta2, eps, w_decay, min_trust, gclip, 0, schedule, warmup, t_total, nothing, nothing)

function initlamb!(model, t_total; min_trust=0.25, lr=0.001, warmup=0.1, schedule="warmup_constant")
    for par in params(model)
        if length(size(value(par))) === 1
            par.opt = MinTrustLAMB(;min_trust=min_trust, lr=lr, warmup=warmup, schedule=schedule, t_total=t_total, w_decay=1.2e-6)
        else
            par.opt = MinTrustLAMB(;min_trust=min_trust, lr=lr, warmup=warmup, schedule=schedule, t_total=t_total)
        end
    end
end

for T in (Array{Float32},Array{Float64},KnetArray{Float32},KnetArray{Float64},CuArray{Float32},CuArray{Float64})
    @eval begin
        function update!(w::$T, g, p::MinTrustLAMB)
            Knet.Train20.gclip!(g, p.gclip)
            g = full(g)
            if p.fstm===nothing; p.fstm=zero(w); p.scndm=zero(w); end
            p.t += 1

            # Decay the first and second moment running average coefficient
            # scndm : v_t, fstm : m_t
            lmul!(p.beta1, p.fstm)
            axpy!(1-p.beta1, g, p.fstm)
            lmul!(p.beta2, p.scndm)
            axpy!(1-p.beta2, g .* g, p.scndm)

            # Get Learning rate from scheduler
            if p.t_total !== -1
                schedule_func = eval(Meta.parse(p.schedule))
                lr_scheduled = p.lr * schedule_func(p.t/p.t_total, p.warmup)
            else
                lr_scheduled = p.lr
            end

            # Calculate Opt. Step
            if p.w_decay > 0.0
                adam_step = (p.fstm ./ (sqrt.(p.scndm) .+ p.eps)) .+ (p.w_decay * w)
            else
                adam_step = (p.fstm ./ (sqrt.(p.scndm) .+ p.eps))
            end

            # Calculate Trust Ratio
            weight_norm = clamp(norm(w), 0, 10)
            adam_norm = norm(adam_step)
            if weight_norm == 0 || adam_norm == 0
                trust_ratio = 1
            else
                trust_ratio = weight_norm / adam_norm
            end
            trust_ratio = max(trust_ratio, p.min_trust)

            # Update weight matrix
            axpy!(-lr_scheduled * trust_ratio, adam_step, w)
        end
    end
end;

nothing


text:
module VOC

import ...BBOX, ...TruthLabel, ...LabelledImageDataset

import LightXML: content,
                 get_elements_by_tagname,
                 parse_file,
                 find_element,
                 root

include(joinpath(@__DIR__, "..", "..", "..", "src", "common.jl"))

const namesdic = Dict("aeroplane"=>1,"bicycle"=>2,"bird"=>3, "boat"=>4,
            "bottle"=>5,"bus"=>6,"car"=>7,"cat"=>8,"chair"=>9,
            "cow"=>10,"diningtable"=>11,"dog"=>12,"horse"=>13,"motorbike"=>14,
            "person"=>15,"pottedplant"=>16,"sheep"=>17,"sofa"=>18,"train"=>19,"tvmonitor"=>20)

labels_dir = joinpath(
    datasets_dir,
    "voc2007",
    "VOCdevkit",
    "VOC2007",
    "Annotations",
)
images_dir = joinpath(datasets_dir, "voc2007", "VOCdevkit", "VOC2007", "JPEGImages")

"""
    populate()

Populate VOC dataset as a `LabelledImageDataset`.
"""
function populate()
    if !isdir(joinpath(@__DIR__, "VOCdevkit", "VOC2007"))
        @error """VOC dataset is not downloaded. Download with `YOLO.download_dataset("voc2007")` and try again."""
        return nothing
    end
    @info "Populating VOC dataset"
    image_paths, label_paths = collectTruthLabelImagePairs(labels_dir, images_dir)
    labels = loadTruthLabel.(label_paths)
    objectcounts = countobjects(labels)

    return LabelledImageDataset(
        name = "VOC",
        objects = flipdict(namesdic),
        objectcounts = objectcounts,
        image_size_lims = (500, 500),
        images_dir = images_dir,
        labels_dir = labels_dir,
        image_paths = image_paths,
        label_paths = label_paths,
        labels = labels,
    )
end

"""
    countobjects(labels::Vector{Vector{TruthLabel}})

Collect objects and count object occurances across entire dataset
"""
function countobjects(labels::Vector{Vector{TruthLabel}})
    occurances = Int[]
    for label in labels
        append!(occurances, map(x -> x.class, label))
    end
    objects = sort(unique(occurances))
    objectcounts = map(x -> sum(occurances .== x), objects)
    return objectcounts
end

"""
    collectTruthLabelImagePairs(labels_dir::String, images_dir::String)

Collects all labels and images' directories, checks for matched label and image files.
"""
function collectTruthLabelImagePairs(labels_dir::String, images_dir::String)
    labels, images = String[], String[]
    excluded = 0
    for (root, dirs, files) in walkdir(labels_dir)
        filter!(x -> endswith(x,".xml") && !startswith(x,"."), files)
        label_filepaths = joinpath.(labels_dir, files)
        image_filepaths = joinpath.(
            images_dir,
            map(x -> string(first(split(x, ".xml")), ".jpg"), files),
        )
        image_exists = isfile.(image_filepaths)
        append!(labels, label_filepaths[image_exists])
        append!(images, image_filepaths[image_exists])
        excluded += sum(.!image_exists) #Count missing images.
    end
    (excluded > 0) && @info "$(excluded) label files excluded due to missing image files"
    return images, labels
end

"""
    loadTruthLabel(label_file::String)

Loads label file into `TruthLabel` struct.
"""
function loadTruthLabel(label_file::String)
    xdoc = parse_file(label_file)
    xroot = root(xdoc)
    ces = get_elements_by_tagname(xroot, "size")
    image_width = parse(Int32, content(find_element(ces[1], "width")))
    image_height = parse(Int32, content(find_element(ces[1], "height")))
    ces = get_elements_by_tagname(xroot, "object")
    YLs = TruthLabel[]
    for i = 1:length(ces)
        name = content(find_element(ces[i], "name"))
        difficult = parse(Int32, content(find_element(ces[i], "difficult")))
        if difficult == 0 #TODO allow different difficulties
            bbox_loc = find_element(ces[i], "bndbox")
            xmin = parse(Int32, content(find_element(bbox_loc, "xmin"))) /
                   image_width
            xmax = parse(Int32, content(find_element(bbox_loc, "xmax"))) /
                   image_width
            ymin = parse(Int32, content(find_element(bbox_loc, "ymin"))) /
                   image_height
            ymax = parse(Int32, content(find_element(bbox_loc, "ymax"))) /
                   image_height
            class = get(namesdic,name,-1)
            class == -1 && @error "Class name mismatch found : $(name) isn't in class list"
            push!(
                YLs,
                TruthLabel(
                    bbox = BBOX(
                        x = xmin,
                        y = ymin,
                        w = xmax - xmin,
                        h = ymax - ymin,
                    ),
                    class = class,
                    difficult = difficult,
                ),
            )
        end
    end
    return YLs
end
end #module


text:
import PredictMD

PredictMD.probability_calibration_scores_and_fractions(
    [0, 0, 0], [0., 0., 0.];
    window = -1,
    )


text:
# ---
# title: Custom scales
# cover: assets/custom_scales.png
# description: Custom palettes and custom attributes
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie
using Colors
set_aog_theme!() #src

# Sometimes, there is no default palettes for a specific attribute. In that
# case, the user can pass their own. TODO: allow legend to use custom attribute
# of plot, such as the arrowhead or the arrowcolor and pass correct legend symbol.

x=repeat(1:20, inner=20)
y=repeat(1:20, outer=20)
u=cos.(x)
v=sin.(y)
c=rand(Bool, length(x))
d=rand(Bool, length(x))
df = (; x, y, u, v, c, d)
colors = [colorant"#E24A33", colorant"#348ABD"]
heads = ['▲', '●']
plt = data(df) *
    mapping(:x, :y, :u, :v) *
    mapping(arrowhead=:c => nonnumeric) *
    mapping(arrowcolor=:d => nonnumeric) *
    visual(Arrows, arrowsize=10, lengthscale=0.3)
draw(plt; palettes=(arrowcolor=colors, arrowhead=heads))

# To associate specific attribute values to specific data values, use pairs.
# Missing keys will cycle over values that are not pairs.

x = rand(100)
y = rand(100)
z = rand(["a", "b", "c", "d"], 100)
df = (; x, y, z)
plt = data(df) * mapping(:x, :y, color=:z)
colors = ["a" => colorant"#E24A33", "c" => colorant"#348ABD", colorant"#988ED5", colorant"#777777"]
fg = draw(plt; palettes=(color=colors,))

# save cover image #src
mkpath("assets") #src
save("assets/custom_scales.png", fg) #src


text:
# Command: Naviagte to
"""
	navigate!(session::Session, url::AbstractString)

The command causes the user agent to navigate the current top-level browsing context to a new location.
"""
function navigate!(session::Session, url::AbstractString)
    try
        isvalid(URI(url))
        @unpack addr, id = session
        response = HTTP.post(
            "$addr/session/$id/url",
            [("Content-Type" => "application/json")],
            """{"url": "$url"}""",
        )
        nothing
    catch
        throw(ArgumentError("url does not seem valid: $url"))
    end
end


text:
@doc raw"""
    Stiefel{n,k,𝔽} <: AbstractEmbeddedManifold{𝔽,DefaultIsometricEmbeddingType}

The Stiefel manifold consists of all $n × k$, $n ≥ k$ unitary matrices, i.e.

````math
\operatorname{St}(n,k) = \bigl\{ p ∈ 𝔽^{n × k}\ \big|\ p^{\mathrm{H}}p = I_k \bigr\},
````

where $𝔽 ∈ \{ℝ, ℂ\}$,
$\cdot^{\mathrm{H}}$ denotes the complex conjugate transpose or Hermitian, and
$I_k ∈ ℝ^{k × k}$ denotes the $k × k$ identity matrix.

The tangent space at a point $p ∈ \mathcal M$ is given by

````math
T_p \mathcal M = \{ X ∈ 𝔽^{n × k} : p^{\mathrm{H}}X + \overline{X^{\mathrm{H}}p} = 0_k\},
````

where $0_k$ is the $k × k$ zero matrix and $\overline{\cdot}$ the (elementwise) complex conjugate.

This manifold is modeled as an embedded manifold to the [`Euclidean`](@ref), i.e.
several functions like the [`inner`](@ref inner(::Euclidean, ::Any...)) product and the
[`zero_vector`](@ref zero_vector(::Euclidean, ::Any...)) are inherited from the embedding.

The manifold is named after
[Eduard L. Stiefel](https://en.wikipedia.org/wiki/Eduard_Stiefel) (1909–1978).

# Constructor
    Stiefel(n, k, field = ℝ)

Generate the (real-valued) Stiefel manifold of $n × k$ dimensional orthonormal matrices.
"""
struct Stiefel{n,k,𝔽} <: AbstractEmbeddedManifold{𝔽,DefaultIsometricEmbeddingType} end

Stiefel(n::Int, k::Int, field::AbstractNumbers=ℝ) = Stiefel{n,k,field}()

@doc raw"""
    PadeRetraction{m} <: AbstractRetractionMethod

A retraction based on the Padé approximation of order $m$
"""
struct PadeRetraction{m} <: AbstractRetractionMethod end

function PadeRetraction(m::Int)
    (m < 1) && error(
        "The Padé based retraction is only available for positive orders, not for order $m.",
    )
    return PadeRetraction{m}()
end
@doc raw"""
    CayleyRetraction <: AbstractRetractionMethod

A retraction based on the Cayley transform, which is realized by using the
[`PadeRetraction`](@ref)`{1}`.
"""
const CayleyRetraction = PadeRetraction{1}

function allocation_promotion_function(::Stiefel{n,k,ℂ}, ::Any, ::Tuple) where {n,k}
    return complex
end

@doc raw"""
    check_point(M::Stiefel, p; kwargs...)

Check whether `p` is a valid point on the [`Stiefel`](@ref) `M`=$\operatorname{St}(n,k)$, i.e. that it has the right
[`AbstractNumbers`](@ref) type and $p^{\mathrm{H}}p$ is (approximately) the identity, where $\cdot^{\mathrm{H}}$ is the
complex conjugate transpose. The settings for approximately can be set with `kwargs...`.
"""
function check_point(M::Stiefel{n,k,𝔽}, p; kwargs...) where {n,k,𝔽}
    mpv = invoke(check_point, Tuple{supertype(typeof(M)),typeof(p)}, M, p; kwargs...)
    mpv === nothing || return mpv
    c = p' * p
    if !isapprox(c, one(c); kwargs...)
        return DomainError(
            norm(c - one(c)),
            "The point $(p) does not lie on $(M), because p'p is not the unit matrix.",
        )
    end
    return nothing
end

@doc raw"""
    check_vector(M::Stiefel, p, X; kwargs...)

Checks whether `X` is a valid tangent vector at `p` on the [`Stiefel`](@ref)
`M`=$\operatorname{St}(n,k)$, i.e. the [`AbstractNumbers`](@ref) fits and
it (approximately) holds that $p^{\mathrm{H}}X + \overline{X^{\mathrm{H}}p} = 0$,
where $\cdot^{\mathrm{H}}$ denotes the Hermitian and $\overline{\cdot}$ the (elementwise) complex conjugate.
The settings for approximately can be set with `kwargs...`.
"""
function check_vector(M::Stiefel{n,k,𝔽}, p, X; kwargs...) where {n,k,𝔽}
    mpv = invoke(
        check_vector,
        Tuple{supertype(typeof(M)),typeof(p),typeof(X)},
        M,
        p,
        X;
        kwargs...,
    )
    mpv === nothing || return mpv
    if !isapprox(p' * X, -conj(X' * p); kwargs...)
        return DomainError(
            norm(p' * X + conj(X' * p)),
            "The matrix $(X) is does not lie in the tangent space of $(p) on the Stiefel manifold of dimension ($(n),$(k)), since p'X + X'p is not the zero matrix.",
        )
    end
    return nothing
end

decorated_manifold(::Stiefel{N,K,𝔽}) where {N,K,𝔽} = Euclidean(N, K; field=𝔽)

@doc raw"""
    inverse_retract(M::Stiefel, p, q, ::PolarInverseRetraction)

Compute the inverse retraction based on a singular value decomposition
for two points `p`, `q` on the [`Stiefel`](@ref) manifold `M`.
This follows the folloing approach: From the Polar retraction we know that

````math
\operatorname{retr}_p^{-1}q = qs - t
````

if such a symmetric positive definite $k × k$ matrix exists. Since $qs - t$
is also a tangent vector at $p$ we obtain

````math
p^{\mathrm{H}}qs + s(p^{\mathrm{H}}q)^{\mathrm{H}} + 2I_k = 0,
````
which can either be solved by a Lyapunov approach or a continuous-time
algebraic Riccati equation.

This implementation follows the Lyapunov approach.
"""
inverse_retract(::Stiefel, ::Any, ::Any, ::PolarInverseRetraction)

@doc raw"""
    inverse_retract(M::Stiefel, p, q, ::QRInverseRetraction)

Compute the inverse retraction based on a qr decomposition
for two points `p`, `q` on the [`Stiefel`](@ref) manifold `M` and return
the resulting tangent vector in `X`. The computation follows Algorithm 1
in [^KanekoFioriTanaka2013].

[^KanekoFioriTanaka2013]:
    > T. Kaneko, S. Fiori, T. Tanaka: "Empirical Arithmetic Averaging over the
    > Compact Stiefel AbstractManifold", IEEE Transactions on Signal Processing, 2013,
    > doi: [10.1109/TSP.2012.2226167](https://doi.org/10.1109/TSP.2012.2226167).
"""
inverse_retract(::Stiefel, ::Any, ::Any, ::QRInverseRetraction)

function _stiefel_inv_retr_qr_mul_by_r_generic!(::Stiefel{n,k}, X, q, R, A) where {n,k}
    @inbounds for i in 1:k
        b = zeros(eltype(R), i)
        b[i] = 1
        b[1:(end - 1)] = -transpose(R[1:(i - 1), 1:(i - 1)]) * A[i, 1:(i - 1)]
        R[1:i, i] = A[1:i, 1:i] \ b
    end
    #TODO: replace with this once it's supported by StaticArrays
    #return mul!(X, q, UpperTriangular(R))
    return mul!(X, q, R)
end

function _stiefel_inv_retr_qr_mul_by_r!(::Stiefel{n,1}, X, q, A, ::Type) where {n}
    @inbounds R = SMatrix{1,1}(inv(A[1, 1]))
    return mul!(X, q, R)
end
function _stiefel_inv_retr_qr_mul_by_r!(
    M::Stiefel{n,1},
    X,
    q,
    A::StaticArray,
    ::Type{ElT},
) where {n,ElT}
    return invoke(
        _stiefel_inv_retr_qr_mul_by_r!,
        Tuple{Stiefel{n,1},typeof(X),typeof(q),AbstractArray,typeof(ElT)},
        M,
        X,
        q,
        A,
        ElT,
    )
end
function _stiefel_inv_retr_qr_mul_by_r!(::Stiefel{n,2}, X, q, A, ::Type{ElT}) where {n,ElT}
    R11 = inv(A[1, 1])
    @inbounds R =
        hcat(SA[R11, zero(ElT)], A[SOneTo(2), SOneTo(2)] \ SA[-R11 * A[2, 1], one(ElT)])

    #TODO: replace with this once it's supported by StaticArrays
    #return mul!(X, q, UpperTriangular(R))
    return mul!(X, q, R)
end
function _stiefel_inv_retr_qr_mul_by_r!(
    M::Stiefel{n,2},
    X,
    q,
    A::StaticArray,
    ::Type{ElT},
) where {n,ElT}
    return invoke(
        _stiefel_inv_retr_qr_mul_by_r!,
        Tuple{Stiefel{n,2},typeof(X),typeof(q),AbstractArray,typeof(ElT)},
        M,
        X,
        q,
        A,
        ElT,
    )
end
function _stiefel_inv_retr_qr_mul_by_r!(
    M::Stiefel{n,k},
    X,
    q,
    A::StaticArray,
    ::Type{ElT},
) where {n,k,ElT}
    R = zeros(MMatrix{k,k,ElT})
    return _stiefel_inv_retr_qr_mul_by_r_generic!(M, X, q, R, A)
end
function _stiefel_inv_retr_qr_mul_by_r!(
    M::Stiefel{n,k},
    X,
    q,
    A,
    ::Type{ElT},
) where {n,k,ElT}
    R = zeros(ElT, k, k)
    return _stiefel_inv_retr_qr_mul_by_r_generic!(M, X, q, R, A)
end

function inverse_retract!(::Stiefel, X, p, q, ::PolarInverseRetraction)
    A = p' * q
    H = -2 * one(p' * p)
    B = lyap(A, H)
    mul!(X, q, B)
    X .-= p
    return X
end
function inverse_retract!(M::Stiefel{n,k}, X, p, q, ::QRInverseRetraction) where {n,k}
    A = p' * q
    @boundscheck size(A) === (k, k)
    ElT = typeof(one(eltype(p)) * one(eltype(q)))
    _stiefel_inv_retr_qr_mul_by_r!(M, X, q, A, ElT)
    X .-= p
    return X
end

function Base.isapprox(M::Stiefel, p, X, Y; kwargs...)
    return isapprox(sqrt(inner(M, p, zero_vector(M, p), X - Y)), 0; kwargs...)
end
Base.isapprox(::Stiefel, p, q; kwargs...) = isapprox(norm(p - q), 0; kwargs...)

@doc raw"""
    manifold_dimension(M::Stiefel)

Return the dimension of the [`Stiefel`](@ref) manifold `M`=$\operatorname{St}(n,k,𝔽)$.
The dimension is given by

````math
\begin{aligned}
\dim \mathrm{St}(n, k, ℝ) &= nk - \frac{1}{2}k(k+1)\\
\dim \mathrm{St}(n, k, ℂ) &= 2nk - k^2\\
\dim \mathrm{St}(n, k, ℍ) &= 4nk - k(2k-1)
\end{aligned}
````
"""
manifold_dimension(::Stiefel{n,k,ℝ}) where {n,k} = n * k - div(k * (k + 1), 2)
manifold_dimension(::Stiefel{n,k,ℂ}) where {n,k} = 2 * n * k - k * k
manifold_dimension(::Stiefel{n,k,ℍ}) where {n,k} = 4 * n * k - k * (2k - 1)

@doc raw"""
    retract(::Stiefel, p, X, ::CayleyRetraction)

Compute the retraction on the [`Stiefel`](@ref) that is based on the Cayley transform[^Zhu2017].
Using
````math
  W_{p,X} = \operatorname{P}_pXp^{\mathrm{H}} - pX^{\mathrm{H}}\operatorname{P_p}
  \quad\text{where}
  \operatorname{P}_p = I - \frac{1}{2}pp^{\mathrm{H}}
````
the formula reads
````math
    \operatorname{retr}_pX = \Bigl(I - \frac{1}{2}W_{p,X}\Bigr)^{-1}\Bigl(I + \frac{1}{2}W_{p,X}\Bigr)p.
````

It is implemented as the case $m=1$ of the [`PadeRetraction`](@ref).

[^Zhu2017]:
    > X. Zhu:
    > A Riemannian conjugate gradient method for optimizazion on the Stiefel manifold,
    > Computational Optimization and Applications 67(1), pp. 73–110, 2017.
    > doi [10.1007/s10589-016-9883-4](https://doi.org/10.1007/s10589-016-9883-4).
"""
retract(::Stiefel, ::Any, ::Any, ::CayleyRetraction)

@doc raw"""
    retract(M::Stiefel, p, X, ::PadeRetraction{m})

Compute the retraction on the [`Stiefel`](@ref) manifold `M` based on the Padé approximation of order $m$[^ZhuDuan2018].
Let $p_m$ and $q_m$ be defined for any matrix $A ∈ ℝ^{n×x}$ as
````math
  p_m(A) = \sum_{k=0}^m \frac{(2m-k)!m!}{(2m)!(m-k)!}\frac{A^k}{k!}
````
and
````math
  q_m(A) = \sum_{k=0}^m \frac{(2m-k)!m!}{(2m)!(m-k)!}\frac{(-A)^k}{k!}
````
respectively. Then the Padé approximation (of the matrix exponential $\exp(A)$) reads
````math
  r_m(A) = q_m(A)^{-1}p_m(A)
````
Defining further
````math
  W_{p,X} = \operatorname{P}_pXp^{\mathrm{H}} - pX^{\mathrm{H}}\operatorname{P_p}
  \quad\text{where}
  \operatorname{P}_p = I - \frac{1}{2}pp^{\mathrm{H}}
````
the retraction reads
````math
  \operatorname{retr}_pX = r_m(W_{p,X})p
````
[^ZhuDuan2018]:
    > X. Zhu, C. Duan:
    > On matrix exponentials and their approximations related to optimization on the Stiefel manifold,
    > Optimizazion Letters 13(5), pp. 1069–1083, 2018.
    > doi [10.1007/s11590-018-1341-z](https://doi.org/10.1007/s11590-018-1341-z).
"""
retract(::Stiefel, ::Any, ::Any, ::PadeRetraction)

@doc raw"""
    retract(M::Stiefel, p, X, ::PolarRetraction)

Compute the SVD-based retraction [`PolarRetraction`](@ref) on the
[`Stiefel`](@ref) manifold `M`. With $USV = p + X$ the retraction reads

````math
\operatorname{retr}_p X = U\bar{V}^\mathrm{H}.
````
"""
retract(::Stiefel, ::Any, ::Any, ::PolarRetraction)

@doc raw"""
    retract(M::Stiefel, p, X, ::QRRetraction)

Compute the QR-based retraction [`QRRetraction`](@ref) on the
[`Stiefel`](@ref) manifold `M`. With $QR = p + X$ the retraction reads

````math
\operatorname{retr}_p X = QD,
````

where $D$ is a $n × k$ matrix with

````math
D = \operatorname{diag}\bigl(\operatorname{sgn}(R_{ii}+0,5)_{i=1}^k \bigr),
````

where $\operatorname{sgn}(p) = \begin{cases}
1 & \text{ for } p > 0,\\
0 & \text{ for } p = 0,\\
-1& \text{ for } p < 0.
\end{cases}$
"""
retract(::Stiefel, ::Any, ::Any, ::QRRetraction)

_qrfac_to_q(qrfac) = Matrix(qrfac.Q)
_qrfac_to_q(qrfac::StaticArrays.QR) = qrfac.Q

function retract!(::Stiefel, q, p, X, ::PadeRetraction{m}) where {m}
    Pp = I - 1 // 2 * p * p'
    WpX = Pp * X * p' - p * X' * Pp
    pm = zeros(eltype(WpX), size(WpX))
    qm = zeros(eltype(WpX), size(WpX))
    WpXk = similar(WpX)
    copyto!(WpXk, factorial(m) / factorial(2 * m) * I) # factorial factor independent of k
    for k in 0:m
        # incrementally build (2m-k)!/(m-k)!(k)! for k > 0, i.e.
        # remove factor (2m-k+1) in the nominator, (m-k+1) in the denominator and multiply by 1/k
        WpXk .*= (k == 0 ? 2 : (m - k + 1) / ((2 * m - k + 1) * k))
        pm .+= WpXk
        if k % 2 == 0
            qm .+= WpXk
        else
            qm .-= WpXk
        end
        WpXk *= WpX
    end
    return copyto!(q, (qm \ pm) * p)
end
function retract!(::Stiefel, q, p, X, ::PolarRetraction)
    s = svd(p + X)
    return mul!(q, s.U, s.Vt)
end
function retract!(::Stiefel, q, p, X, ::QRRetraction)
    qrfac = qr(p + X)
    d = diag(qrfac.R)
    D = Diagonal(sign.(sign.(d .+ 0.5)))
    return mul!(q, _qrfac_to_q(qrfac), D)
end

@doc raw"""
    representation_size(M::Stiefel)

Returns the representation size of the [`Stiefel`](@ref) `M`=$\operatorname{St}(n,k)$,
i.e. `(n,k)`, which is the matrix dimensions.
"""
@generated representation_size(::Stiefel{n,k}) where {n,k} = (n, k)

Base.show(io::IO, ::CayleyRetraction) = print(io, "CayleyRetraction()")
Base.show(io::IO, ::PadeRetraction{m}) where {m} = print(io, "PadeRetraction($(m))")
Base.show(io::IO, ::Stiefel{n,k,F}) where {n,k,F} = print(io, "Stiefel($(n), $(k), $(F))")

"""
    uniform_distribution(M::Stiefel{n,k,ℝ}, p)

Uniform distribution on given (real-valued) [`Stiefel`](@ref) `M`.
Specifically, this is the normalized Haar and Hausdorff measure on `M`.
Generated points will be of similar type as `p`.

The implementation is based on Section 2.5.1 in [^Chikuse2003];
see also Theorem 2.2.1(iii) in [^Chikuse2003].

[^Chikuse2003]:
    > Y. Chikuse: "Statistics on Special Manifolds", Springer New York, 2003,
    > doi: [10.1007/978-0-387-21540-2](https://doi.org/10.1007/978-0-387-21540-2).
"""
function uniform_distribution(M::Stiefel{n,k,ℝ}, p) where {n,k}
    μ = Distributions.Zeros(n, k)
    σ = one(eltype(p))
    Σ1 = Distributions.PDMats.ScalMat(n, σ)
    Σ2 = Distributions.PDMats.ScalMat(k, σ)
    d = MatrixNormal(μ, Σ1, Σ2)

    return ProjectedPointDistribution(M, d, project!, p)
end

@doc raw"""
    vector_transport_direction(::Stiefel, p, X, d, ::DifferentiatedRetractionVectorTransport{CayleyRetraction})

Compute the vector transport given by the differentiated retraction of the [`CayleyRetraction`](@ref), cf. [^Zhu2017] Equation (17).

The formula reads
````math
\operatorname{T}_{p,d}(X) =
\Bigl(I - \frac{1}{2}W_{p,d}\Bigr)^{-1}W_{p,X}\Bigl(I - \frac{1}{2}W_{p,d}\Bigr)^{-1}p,
````
with
````math
  W_{p,X} = \operatorname{P}_pXp^{\mathrm{H}} - pX^{\mathrm{H}}\operatorname{P_p}
  \quad\text{where}
  \operatorname{P}_p = I - \frac{1}{2}pp^{\mathrm{H}}
````

Since this is the differentiated retraction as a vector transport, the result will be in the
tangent space at $q=\operatorname{retr}_p(d)$ using the [`CayleyRetraction`](@ref).
"""
vector_transport_direction(
    M::Stiefel,
    p,
    X,
    d,
    ::DifferentiatedRetractionVectorTransport{CayleyRetraction},
)

@doc raw"""
    vector_transport_direction(M::Stiefel, p, X, d, DifferentiatedRetractionVectorTransport{PolarRetraction})

Compute the vector transport by computing the push forward of
[`retract(::Stiefel, ::Any, ::Any, ::PolarRetraction)`](@ref) Section 3.5 of [^Zhu2017]:

```math
T_{p,d}^{\text{Pol}}(X) = q*Λ + (I-qq^{\mathrm{T}})X(1+d^\mathrm{T}d)^{-\frac{1}{2}},
```

where $q = \operatorname{retr}^{\mathrm{Pol}}_p(d)$, and $Λ$ is the unique solution of the Sylvester equation

```math
    Λ(I+d^\mathrm{T}d)^{\frac{1}{2}} + (I + d^\mathrm{T}d)^{\frac{1}{2}} = q^\mathrm{T}X - X^\mathrm{T}q
```
"""
vector_transport_direction(
    ::Stiefel,
    ::Any,
    ::Any,
    ::Any,
    ::DifferentiatedRetractionVectorTransport{PolarRetraction},
)
@doc raw"""
    vector_transport_direction(M::Stiefel, p, X, d, DifferentiatedRetractionVectorTransport{QRRetraction})

Compute the vector transport by computing the push forward of the
[`retract(::Stiefel, ::Any, ::Any, ::QRRetraction)`](@ref),
See [^AbsilMahonySepulchre2008], p. 173, or Section 3.5 of [^Zhu2017].
```math
T_{p,d}^{\text{QR}}(X) = q*\rho_{\mathrm{s}}(q^\mathrm{T}XR^{-1}) + (I-qq^{\mathrm{T}})XR^{-1},
```
where $q = \operatorname{retr}^{\mathrm{QR}}_p(d)$, $R$ is the $R$ factor of the QR
decomposition of $p + d$, and
```math
\bigl( \rho_{\mathrm{s}}(A) \bigr)_{ij}
= \begin{cases}
A_{ij}&\text{ if } i > j\\
0 \text{ if } i = j\\
-A_{ji} \text{ if } i < j.\\
\end{cases}
```
[^AbsilMahonySepulchre2008]:
    >Absil, P.-A., Mahony, R. and Sepulchre R.,
    > _Optimization Algorithms on Matrix Manifolds_
    > Princeton University Press, 2008,
    > doi: [10.1515/9781400830244](https://doi.org/10.1515/9781400830244)
    > [open access](http://press.princeton.edu/chapters/absil/)
"""
vector_transport_direction(
    ::Stiefel,
    ::Any,
    ::Any,
    ::Any,
    ::DifferentiatedRetractionVectorTransport{QRRetraction},
)

function vector_transport_direction!(
    ::Stiefel,
    Y,
    p,
    X,
    d,
    ::DifferentiatedRetractionVectorTransport{CayleyRetraction},
)
    Pp = I - 1 // 2 * p * p'
    Wpd = Pp * d * p' - p * d' * Pp
    WpX = Pp * X * p' - p * X' * Pp
    q1 = I - 1 // 2 * Wpd
    return copyto!(Y, (q1 \ WpX) * (q1 \ p))
end

function vector_transport_direction!(
    M::Stiefel,
    Y,
    p,
    X,
    d,
    ::DifferentiatedRetractionVectorTransport{PolarRetraction},
)
    q = retract(M, p, d, PolarRetraction())
    Iddsqrt = sqrt(I + d' * d)
    Λ = sylvester(Iddsqrt, Iddsqrt, -q' * X + X' * q)
    return copyto!(Y, q * Λ + (X - q * (q' * X)) / Iddsqrt)
end
function vector_transport_direction!(
    M::Stiefel,
    Y,
    p,
    X,
    d,
    ::DifferentiatedRetractionVectorTransport{QRRetraction},
)
    q = retract(M, p, d, QRRetraction())
    rf = UpperTriangular(qr(p + d).R)
    Xrf = X / rf
    qtXrf = q' * Xrf
    return copyto!(
        Y,
        q * (UpperTriangular(qtXrf) - UpperTriangular(qtXrf)') + Xrf - q * qtXrf,
    )
end

@doc raw"""
    vector_transport_to(M::Stiefel, p, X, q, DifferentiatedRetractionVectorTransport{PolarRetraction})

Compute the vector transport by computing the push forward of the
[`retract(M::Stiefel, ::Any, ::Any, ::PolarRetraction)`](@ref), see
Section 4 of [^HuangGallivanAbsil2015] or  Section 3.5 of [^Zhu2017]:

```math
T_{q\gets p}^{\text{Pol}}(X) = q*Λ + (I-qq^{\mathrm{T}})X(1+d^\mathrm{T}d)^{-\frac{1}{2}},
```

where $d = \bigl( \operatorname{retr}^{\mathrm{Pol}}_p\bigr)^{-1}(q)$,
and $Λ$ is the unique solution of the Sylvester equation

```math
    Λ(I+d^\mathrm{T}d)^{\frac{1}{2}} + (I + d^\mathrm{T}d)^{\frac{1}{2}} = q^\mathrm{T}X - X^\mathrm{T}q
```
[^HuangGallivanAbsil2015]:
    > Huang, W., Gallivan, K. A., and Absil, P.-A.:
    > _A Broyden class of quasi-Newton methods for Riemannian optimization_
    > SIAM Journal of Optimization, 2015, Vol. 25, No. 3, pp. 1660–1685
    > doi: [10.1137/140955483](https://doi.org/10.1137/140955483)
    > pdf: [tech. report](https://www.math.fsu.edu/~whuang2/pdf/RBroydenBasic_techrep.pdf)
"""
vector_transport_to(
    ::Stiefel,
    ::Any,
    ::Any,
    ::Any,
    ::DifferentiatedRetractionVectorTransport{PolarRetraction},
)

@doc raw"""
    vector_transport_to(M::Stiefel, p, X, q, DifferentiatedRetractionVectorTransport{QRRetraction})

Compute the vector transport by computing the push forward of the
[`retract(M::Stiefel, ::Any, ::Any, ::QRRetraction)`](@ref),
see [^AbsilMahonySepulchre2008], p. 173, or Section 3.5 of [^Zhu2017].

```math
T_{q \gets p}^{\text{QR}}(X) = q*\rho_{\mathrm{s}}(q^\mathrm{T}XR^{-1}) + (I-qq^{\mathrm{T}})XR^{-1},
```
where $d = \bigl(\operatorname{retr}^{\mathrm{QR}}\bigr)^{-1}_p(q)$, $R$ is the $R$ factor of the QR
decomposition of $p+X$, and
```math
\bigl( \rho_{\mathrm{s}}(A) \bigr)_{ij}
= \begin{cases}
A_{ij}&\text{ if } i > j\\
0 \text{ if } i = j\\
-A_{ji} \text{ if } i < j.\\
\end{cases}
```
"""
vector_transport_to(
    ::Stiefel,
    ::Any,
    ::Any,
    ::Any,
    ::DifferentiatedRetractionVectorTransport{QRRetraction},
)

@doc raw"""
    vector_transport_to(M::Stiefel, p, X, q, ::ProjectionTransport)

Compute a vector transport by projection, i.e. project `X` from the tangent space at `x` by
projection it onto the tangent space at `q`.
"""
vector_transport_to(::Stiefel, ::Any, ::Any, ::Any, ::ProjectionTransport)

function vector_transport_to!(
    M::Stiefel,
    Y,
    p,
    X,
    q,
    ::DifferentiatedRetractionVectorTransport{PolarRetraction},
)
    d = inverse_retract(M, p, q, PolarInverseRetraction())
    Iddsqrt = sqrt(I + d' * d)
    Λ = sylvester(Iddsqrt, Iddsqrt, -q' * X + X' * q)
    return copyto!(Y, q * Λ + (X - q * (q' * X)) / Iddsqrt)
end
function vector_transport_to!(
    M::Stiefel,
    Y,
    p,
    X,
    q,
    ::DifferentiatedRetractionVectorTransport{QRRetraction},
)
    d = inverse_retract(M, p, q, QRInverseRetraction())
    rf = UpperTriangular(qr(p + d).R)
    Xrf = X / rf
    qtXrf = q' * Xrf
    return copyto!(
        Y,
        q * (UpperTriangular(qtXrf) - UpperTriangular(qtXrf)') + Xrf - q * qtXrf,
    )
end


text:
function parse_script(document_body)
    lines = split(document_body, '\n')

    doc_line = r"(^#'.*)|(^#%%.*)|(^# %%.*)"
    doc_start = r"(^#')|(^#%%)|(^# %%)"
    opt_line = r"(^#\+.*$)|(^#%%\+.*$)|(^# %%\+.*$)"
    opt_start = r"(^#\+)|(^#%%\+)|(^# %%\+)"

    content = ""
    state = :code
    doc_no = 0
    code_no = 0
    start_line = 1

    options = OptionDict()
    option_string = ""

    chunks = WeaveChunk[]
    for (line_no, line) in enumerate(lines)
        if (m = match(doc_line, line)) !== nothing && (m = match(opt_line, line)) === nothing
            line = replace(line, doc_start => "", count = 1)
            startswith(line, ' ') && (line = replace(line, ' ' => "", count = 1))
            if state === :code && !isempty(strip(content))
                push!(chunks, CodeChunk(string('\n', strip(content)), code_no += 1, start_line, option_string, options))
                content = ""
                start_line = line_no
            end
            state = :doc
        elseif (m = match(opt_line, line)) !== nothing
            start_line = line_no
            if state === :code && !isempty(strip(content))
                push!(chunks, CodeChunk(string('\n', strip(content)), code_no += 1, start_line, option_string, options))
                content = ""
            end
            if state === :doc && !isempty(strip(content))
                iszero(doc_no) || (content = string('\n', content)) # Add whitespace to doc chunk. Needed for markdown output
                push!(chunks, DocChunk(content, doc_no += 1, start_line))
                content = ""
            end

            option_string = replace(line, opt_start => "", count = 1)
            options = parse_options(option_string)
            haskey(options, :label) && (options[:name] = options[:label])
            haskey(options, :name) || (options[:name] = nothing)

            state = :code
            continue
        elseif state === :doc # && strip(line) != "" && strip(content) != ""
            state = :code
            iszero(doc_no) || (content = string('\n', content)) # Add whitespace to doc chunk. Needed for markdown output
            push!(chunks, DocChunk(content, doc_no += 1, start_line))
            content = ""

            options = Dict{Symbol,Any}()
            start_line = line_no
        end
        content *= string(line, '\n')
    end

    # Handle the last chunk
    chunk = state === :code ?
        CodeChunk(string('\n', strip(content)), code_no, start_line, option_string, options) :
        DocChunk(content, doc_no, start_line)
    push!(chunks, chunk)

    return Dict(), chunks
end


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 1a5c4360-c5ec-11eb-073d-8dc280b34894
begin
	using CSV
	using DataFrames
	using GLM
	using Plots
	using PlutoUI
end

# ╔═╡ abf556e0-810d-443c-87b1-d0f95f9c080d
data = CSV.read("D:\\Workspace\\Machine Learning\\Canada GDP Analysis\\Data\\canada_gdp.csv",DataFrame)

# ╔═╡ 02178664-f416-4505-932a-019ca57f17cb
scatter(data.year,data.GDP)

# ╔═╡ e6a9282e-b580-4fd4-bb16-e83a8e5e0105
reg = fit(LinearModel,@formula(GDP~year),data)

# ╔═╡ b55b5c14-5761-4297-8f56-c1a6fb86e67c
linear_regression = predict(reg)

# ╔═╡ 5c0c577a-3571-4d5f-90da-c5c94942673a
begin
	scatter(data.year,data.GDP)
	plot!(data.year,linear_regression)
end

# ╔═╡ e8c03cd5-12a6-45d3-80b4-20feae154e32
md"Enter a year to predict: $@bind y TextField()"

# ╔═╡ d383695f-5bb5-4b8b-8876-d856405ec8aa
pyear = parse(Int64,y)

# ╔═╡ fbb51392-e820-4317-9334-e0c722d9a423
pGDP = coef(reg)[2]*pyear+coef(reg)[1]

# ╔═╡ Cell order:
# ╠═1a5c4360-c5ec-11eb-073d-8dc280b34894
# ╠═abf556e0-810d-443c-87b1-d0f95f9c080d
# ╠═02178664-f416-4505-932a-019ca57f17cb
# ╠═e6a9282e-b580-4fd4-bb16-e83a8e5e0105
# ╠═b55b5c14-5761-4297-8f56-c1a6fb86e67c
# ╠═5c0c577a-3571-4d5f-90da-c5c94942673a
# ╠═e8c03cd5-12a6-45d3-80b4-20feae154e32
# ╠═d383695f-5bb5-4b8b-8876-d856405ec8aa
# ╠═fbb51392-e820-4317-9334-e0c722d9a423

max_stars_repo_path: Canada GDP Analysis/main.jl

max_stars_repo_name: chidhvilash/GDP-PerCapita-Analysis

max_issues_repo_path: Canada GDP Analysis/main.jl

max_issues_repo_name: chidhvilash/GDP-PerCapita-Analysis

max_forks_repo_path: Canada GDP Analysis/main.jl

max_forks_repo_name: chidhvilash/GDP-PerCapita-Analysis


text:
export UpdateEloF

"""
    UpdateEloF <: UpdateRule

 Update using standard Elo where
* K is constant, but
* performance model distribution can be set (with constant parameters)
* extra factors are added to ratings, with a scale parameter

## Parameters
* `r0::Float64`: default rating
* `K::Float64`: the value of the gain parameter (use EloF2 for variable K)
* `dist::ContinuousUnivariateDistribution`: the model to use for performance as a function of rating, e.g, Normal or Logistic
* `factor_scale::Float64`
* the distribution can have its own parameters, e.g., `β::Float64`, the spread parameter for the Logistic

```
"""
struct UpdateEloF <: UpdateRule
    r0::Float64
    K::Float64
    dist::ContinuousUnivariateDistribution
    factor_scale::Float64
    function UpdateEloF(r0, K, dist, factor_scale)
        @check_args(UpdateEloF, K >= zero(K))
        @check_args(UpdateEloF, r0 >= zero(r0))
        @check_args(UpdateEloF, factor_scale >= zero(factor_scale))
        new(r0, K, dist, factor_scale)
    end
end
# defaults imply the distribution is  Logistic(μ=0.0, θ=100.0)
#   https://en.wikipedia.org/wiki/Logistic_distribution
#   https://juliastats.github.io/Distributions.jl/latest/univariate.html#Continuous-Distributions-1
#      mean=mu=0.0, sdev=π θ / √3
default_θ = 400 / log(10) # this results in the standard values used for Elo in Chess and many other cases
UpdateEloF(;r0=1500.0, K=32.0, θ=default_θ, factor_scale=default_θ) = UpdateEloF( r0, K, Logistic(0.0, θ), factor_scale )
# default factor_scale is one "standard deviation"

function update_info( rule::UpdateEloF )
    info = Dict(
                :name => "EloF",
                :reference => "\"Whos's #1\", Langville and Meyer, p.?? and ??",
                :computation => "simultaneous",
                :state_model => "recursive",
                :input => "outcome,factors",
                :output => "probabilistic",
                :model => "single",
                :ties => true,
                :factors => true,
                :parameters => ["r0, default rating", "K, gain", "dist(), performance model", "factor_scale"],
                :record => false
                )
    return info
end
update_info( ::Type{UpdateEloF} ) =  update_info( UpdateEloF() )

function update_ratings( rule::UpdateEloF,
                         input_ratings::RatingsList,
                         input_competitions::DataFrame)
    n = size(input_competitions,1)
    m = length( input_ratings.players )
    I = player_indexes( input_ratings.players )

    # set initial ratings if they aren't already
    old_r = copy( input_ratings.ratings )
    for p in input_ratings.players
        if !haskey( old_r, p )
            old_r[p] = rule.r0
        end
    end
    new_r = copy(old_r) # start with the old ratings, in case one isn't modified at all

    # perform update for each outcome
    d = input_competitions # just an abbreviation
    for i=1:n
        pA = d[i,PlayerA]
        pB = d[i,PlayerB]
        old_rA = old_r[ pA ]
        old_rB = old_r[ pB ]
        (ΔrA, ΔrB) = update( rule, old_rA, old_rB, d[i,Outcome], rule.factor_scale*d[i,FactorA], rule.factor_scale*d[i,FactorB])
        new_r[ pA ] += ΔrA
        new_r[ pB ] += ΔrB
    end

    # output ratings list
    output_ratings = RatingsList(input_ratings.players, new_r )
    return output_ratings
end

# see Elo.jl for definition of `update`


text:
#=
Example 3:
  Fit the following analytical model with 3 parameters:
  (p1  +  p2 * x  +  p3 * x^2)  *  cos(x)

A `FuncWrap` component is used to wrap a common Julia function
=#
using Random, DataFitting

# Define the analytic model function
f(x, p1, p2, p3) = @. (p1  +  p2 * x  +  p3 * x^2)  *  cos(x)

# "True" parameter values:
params = [1, 2, 3];

# "True" physical quantities
x = 0:0.1:15
y = f(x, params...);

# Add some noise to simulate a measurement process
noise = 0.1 .* y
y .+= noise .* randn(length(x));

# Wrap empirical measures and uncertainties in a `Measure` object
data = Measures(y, noise)

# Prepare the model with a `FuncWrap` component.  Provide the guess
# values with the `params` vector.
model = Model(Domain(x), :comp1 => DataFitting.FuncWrap(f, params...))

# Fit the model to empirical data
result = fit!(model, data)

# Print best fit parameters
for i in 1:length(result[:comp1].p)
    println("p[$i] = ", result[:comp1].p[i].val, " ± ", result[:comp1].p[i].unc)
end

# Plot data and best fit model
using Gnuplot
@gp    x data.val data.unc "w yerr" :-
@gp :- x model[1] "w line"


text:
const Coefficients{T1,T2} = NamedTuple{(:all,:dim), Tuple{T1,T2}}
const Residuals{T1,T2,T3} = NamedTuple{(:main,:dist,:weighted), Tuple{T1,T2,T3}}

"""
```
__allocate_coeff__(T, nfeatures, nclasses)
```

Allocate a `NamedTuple` storing data similar in shape to model coefficients.
"""
function __allocate_coeff__(T, p, c)
    arr = similar(Matrix{T}, p, c-1)
    return (
        all=arr,                                # full coefficients vector
        dim=[view(arr, :, j) for j in 1:c-1]  # views into parameters for each class
    )
end

"""
```
__allocate_res__(T, nsamples, nfeatures, nclasses)
```

Allocate a `NamedTuple` storing data similar in shape to model residuals.
"""
function __allocate_res__(T, n, p, c)
    arr_main = similar(Matrix{T}, n, c-1)
    arr_dist = similar(Matrix{T}, p, c-1)
    arr_weighted = similar(Matrix{T}, n, c-1)

    # views into main residuals, R = Y - BX
    main = (
        all=arr_main,
        dim=[view(arr_main, :, j) for j in 1:c-1],
        sample=[view(arr_main, i, :) for i in 1:n],
    )

    # views into distance residuals, P - B
    dist = (
        all=arr_dist,
        dim=[view(arr_dist, :, j) for j in 1:c-1]
    )

    # views into weighted residuals, Z - BX
    weighted = (
        all=arr_weighted,
        dim=[view(arr_weighted, :, j) for j in 1:c-1],
        sample=[view(arr_weighted, i, :) for i in 1:n],
    )

    return (main=main, dist=dist, weighted=weighted,)
end

function __allocate_problem_arrays__(::Nothing, ::Type{T}, n, p, c, intercept::Bool) where T<:Real
    coeff = __allocate_coeff__(T, p+intercept, c)
    coeff_prev = __allocate_coeff__(T, p+intercept, c)
    proj = __allocate_coeff__(T, p+intercept, c)
    res = __allocate_res__(T, n, p+intercept, c)
    grad = __allocate_coeff__(T, p+intercept, c)
    return coeff, coeff_prev, proj, res, grad
end

function __allocate_problem_arrays__(::Kernel, ::Type{T}, n, p, c, intercept::Bool) where T<:Real
    coeff = __allocate_coeff__(T, n+intercept, c)
    coeff_prev = __allocate_coeff__(T, n+intercept, c)
    proj = __allocate_coeff__(T, n+intercept, c)
    res = __allocate_res__(T, n, n+intercept, c)
    grad = __allocate_coeff__(T, n+intercept, c)
    return coeff, coeff_prev, proj, res, grad
end

function vertex_encoding(::Type{T}, labels, class) where T<:Real
    # dimensions
    c = length(class)

    # enumerate vertices
    a = ( 1 + sqrt(c) ) / ( (c-1)^(3/2) )
    b = sqrt( c / (c-1) )
    vertex = Vector{Vector{T}}(undef, c)
    vertex[1] = 1 / sqrt(c-1) * ones(c-1)
    for j in 2:c
        v = -a * ones(c-1)
        v[j-1] += b
        vertex[j] = v
    end

    # encoding: map labels to vertices and vice-versa
    label2vertex = Dict(class_j => vertex[j] for (j, class_j) in enumerate(class))
    vertex2label = Dict(vertex[j] => class_j for (j, class_j) in enumerate(class))

    # create response matrix based on encoding
    Y = create_Y(T, labels, label2vertex)

    return Y, vertex, label2vertex, vertex2label
end

function create_Y(::Type{T}, labels, label2vertex) where T<:Real
    n, c = length(labels), length(label2vertex)
    Y = Matrix{T}(undef, n, c-1)
    @views for (i, label_i) in enumerate(labels)
        Y[i,:] .= label2vertex[label_i]
    end
    return Y
end

function create_X_and_K(kernel::Kernel, ::Type{T}, data, intercept) where T<:Real
    K = kernelmatrix(kernel, data, obsdim=1)
    intercept && (K = [K ones(size(K, 1))])
    X = copy(data)
    return X, K
end

function create_X_and_K(kernel::Nothing, ::Type{T}, data, intercept) where T<:Real
    X = intercept ? [data ones(size(data, 1))] : copy(data)
    return X, nothing
end

"""
    MVDAProblem{T} where T<:AbstractFloat

Representation of a vertex discriminant analysis problem using floating-point type `T`.
Each unique class label is mapped to a unique vertex of a regular simplex. A `(c-1)`-simplex is used to encode `c` classes.

**Note**: The matrix `X` is used as the *design matrix* for linear problems (i.e. `kernel isa Nothing`); otherwise the matrix `K` is used for the design (i.e. `kernel isa Kernel`).
In either case, the *design matrix* will contain a column of `1`s when `intecept=true`.

# Fields

- `n`: Number of samples/instances in problem.
- `p`: Number of features in each sample/instance.
- `c`: Total number of classes represented in vertex space.

- `Y`: Response matrix. Each row `Y[i,:]` corresponds to the vertex assigned to sample `i`.
- `X`: Data matrix. Each row `X[i,:]` corresponds to a sample `i` with `p` features.
- `K`: Kernel matrix. Each row `K[i,:]` corresponds to a sample `i`.

- `vertex`: List of vertices used to encode classes.
- `label2vertex`: An associative map translating a value from original label space to vertex space.
- `vertex2label`: An associative map trasnlating a value from vertex space to the original label space.
- `intercept`: Indicates whether the design matrix/model include an intercept term (`intercept=true`).

- `kernel`: A `Kernel` object from KernelFunctions.jl. See also: [`Kernel`](@ref).
- `coeff`: A `NamedTuple` containing the current estimates for the full model (`coeff.all`) and views along each dimension in vertex space (`coeff.dim`).
- `coeff_prev`: Similar to `coeff`, but contains previous estimates.
- `proj`: Similar to `coeff`, but contains projection of `coeff`.
- `res`: A `NamedTuple` containing various residuals.
- `grad`: Similar to `coeff`, but contains gradient with respect to `coeff`.
"""
struct MVDAProblem{T<:AbstractFloat,kernT,matT<:AbstractMatrix{T},labelT,viewT,R1,R2,R3}
    ##### dimensions #####
    n::Int
    p::Int
    c::Int

    ##### data #####
    Y::matT
    X::matT
    K::Union{Nothing,matT}
    vertex::Vector{Vector{T}}
    label2vertex::Dict{labelT,Vector{T}}
    vertex2label::Dict{Vector{T},labelT}
    intercept::Bool

    ##### model #####
    kernel::kernT
    coeff::Coefficients{matT,viewT}
    coeff_prev::Coefficients{matT,viewT}

    ##### quadratic surrogate #####
    proj::Coefficients{matT,viewT}
    res::Residuals{R1,R2,R3}
    grad::Coefficients{matT,viewT}

    ##### Inner Constructor #####
    function MVDAProblem{T}(n, p, c, Y, X, K,
        vertex, label2vertex, vertex2label, intercept,
        kernel, coeff, coeff_prev,
        proj, res, grad) where T <: Real
        # get type parameters
        kernT = typeof(kernel)
        matT = typeof(Y)
        labelT = keytype(label2vertex)
        viewT = typeof(coeff.dim)
        R1 = typeof(res.main)
        R2 = typeof(res.dist)
        R3 = typeof(res.weighted)

        new{T,kernT,matT,labelT,viewT,R1,R2,R3}(
            n, p, c,
            Y, X, K, vertex, label2vertex, vertex2label, intercept,
            kernel, coeff, coeff_prev,
            proj, res, grad,
        )
    end
end

"""
    MVDAProblem(labels, data; [intercept=true], [kernel=nothing])

Create a `MVDAProblem` instance from the labeled dataset `(label, data)`.

The `label` information should enter as an iterable object, and `data` should be a `n × p` matrix with samples/instances aligned along rows (e.g. `data[i,:]` is sample `i`).

!!! note

    Defaults to linear classifier, `kernel=nothing`.
    Specifying a `Kernel` requires `using KernelFunctions` first.

# Keyword Arguments

- `intercept`: Should the model include an intercept term?
- `kernel`: How should the data be transformed?.
"""
function MVDAProblem(labels, data; intercept::Bool=true, kernel::Union{Nothing,Kernel}=nothing)
    # get problem info
    class = sort!(unique(labels))
    (n, p), c = size(data), length(class)
    T = Float64 # TODO

    # create design matrices
    X, K = create_X_and_K(kernel, T, data, intercept)

    # assign classes to vertices and create response matrix
    Y, vertex, label2vertex, vertex2label = vertex_encoding(T, labels, class)

    # allocate data structures for coefficients, projections, residuals, and gradient
    coeff, coeff_prev, proj, res, grad = __allocate_problem_arrays__(kernel, T, n, p, c, intercept)

    return MVDAProblem{T}(
        n, p, c,
        Y, X, K, vertex, label2vertex, vertex2label, intercept,
        kernel, coeff, coeff_prev,
        proj, res, grad,
    )
end

"""
    change_data(problem::MVDAProblem, labels::AbstractVector, data)

Create a new `MVDAProblem` instance from the labeled dataset `(label, data)` using the vertex encoding from the reference `problem`.
"""
function change_data(problem::MVDAProblem, labels::AbstractVector, data)
    @unpack label2vertex = problem          # extract encoding-dependent fields + problem info
    T = floattype(problem)
    Y = create_Y(T, labels, label2vertex)   # create response matrix
    change_data(problem, Y, data)           # dispatch
end

"""
    change_data(problem::MVDAProblem, Y::AbstractMatrix, data)

Create a new `MVDAProblem` instance from the labeled dataset `(Y, data)` using the vertex encoding from the reference `problem`.

Assumes `Y` already contains vertex assignments consistent with the encoding in `problem`.
"""
function change_data(problem::MVDAProblem, Y::AbstractMatrix, data)
    # extract encoding-dependent fields + problem info
    @unpack vertex, label2vertex, vertex2label, intercept, kernel = problem
    has_intercept = all(isequal(1), data[:,end])
    n, p, c = size(Y, 1), has_intercept ? size(data, 2)-1 : size(data, 2), problem.c
    T = floattype(problem)

    # create new design matrices
    X, K = create_X_and_K(kernel, T, data, intercept)

    # allocate data structures for coefficients, projections, residuals, and gradient
    coeff, coeff_prev, proj, res, grad = __allocate_problem_arrays__(kernel, T, n, p, c, intercept)

    return MVDAProblem{T}(
        n, p, c,
        Y, X, K, vertex, label2vertex, vertex2label, intercept,
        kernel, coeff, coeff_prev,
        proj, res, grad,
    )
end

"""
Return the floating-point type used for model coefficients.
"""
floattype(::MVDAProblem{T}) where T = T

"""
Return the design matrix used for fitting a classifier.

Uses `problem.X` when `problem.kernel isa Nothing` and `problem.K` when `problem.kernel isa Kernel` from KernelFunctions.jl.
"""
get_design_matrix(problem::MVDAProblem) = __get_design_matrix__(problem.kernel, problem) # dispatch
__get_design_matrix__(::Nothing, problem::MVDAProblem) = problem.X  # linear case
__get_design_matrix__(::Kernel, problem::MVDAProblem) = problem.K   # nonlinear case

"""
Returns the number of samples, number of features, and number of categories, respectively.
"""
probdims(problem::MVDAProblem) = (problem.n, problem.p, problem.c)

"""
    predict(problem::MVDAProblem, x)

When `x` is a vector, predict the vertex value of a sample/instance `x` based on the fitted model in `problem`.
Otherwise if `x` is a matrix then each sample is assumed to be aligned along rows (e.g. `x[i,:]` is sample `i`).

See also: [`classify`](@ref)
"""
predict(problem::MVDAProblem, x) = __predict__(problem.kernel, problem, x)

function __predict__(::Nothing, problem::MVDAProblem, x::AbstractVector)
    @unpack p, proj, intercept = problem
    B = view(proj.all, 1:p, :)
    B0 = view(proj.all, p+intercept, :)
    y = B' * x
    intercept && (y += B0)
    return y
end

function __predict__(::Nothing, problem::MVDAProblem, X::AbstractMatrix)
    @unpack p, proj, intercept = problem
    B = view(proj.all, 1:p, :)
    B0 = view(proj.all, p+intercept, :)
    Y = X * B
    intercept && (Y .+= B0')
    return Y
end

function __predict__(::Kernel, problem::MVDAProblem, x::AbstractVector)
    c = problem.c
    κ = problem.kernel
    Γ = problem.proj.all
    X = problem.X

    ϕ = zeros(floattype(problem), c-1)
    for (i, xᵢ) in enumerate(eachrow(X))
        k = κ(xᵢ, x)
        γ = Γ[i,:]
        @inbounds for j in eachindex(ϕ)
            ϕ[j] = muladd(k, γ[j], ϕ[j])
        end
    end

    # ...and accumulate intercept term, ϕⱼ += Γ₀ⱼ.
    if problem.intercept
        @inbounds for j in eachindex(ϕ)
            ϕ[j] += Γ[end,j]
        end
    end

    return ϕ
end

function __predict__(::Kernel, problem::MVDAProblem, X::AbstractMatrix)
    @unpack c = problem
    n = size(X, 1)
    Y = Matrix{floattype(problem)}(undef, n, c-1)
    nthreads = BLAS.get_num_threads()
    BLAS.set_num_threads(1)
    @batch per=core for i in 1:n
        Y[i,:] .= predict(problem, X[i,:])
    end
    BLAS.set_num_threads(nthreads)
    return Y
end

"""
    classify(problem::MVDAProblem, x)

Classify the samples/instances in `x` based on the model in `problem`.

If `x` is a vector then it is treated as an instance.
Otherwise if `x` is a matrix then each sample is assumed to be aligned along rows (e.g. `x[i,:]` is sample `i`).
See also: [`predict`](@ref)
"""
classify(problem::MVDAProblem, x) = __classify__(problem, predict(problem, x))

function __classify__(problem::MVDAProblem, y::AbstractVector)
    @unpack vertex, vertex2label = problem
    distances = [norm(y - vertex[j]) for j in eachindex(vertex)]
    j = argmin(distances)
    return problem.vertex2label[vertex[j]]
end

function __classify__(problem::MVDAProblem, Y::AbstractMatrix)
    n = size(Y, 1)
    L = Vector{valtype(problem.vertex2label)}(undef, n)
    nthreads = BLAS.get_num_threads()
    BLAS.set_num_threads(1)
    @batch per=core for i in eachindex(L)
        L[i] = __classify__(problem, Y[i,:])
    end
    BLAS.set_num_threads(nthreads)
    return L
end

function Base.show(io::IO, problem::MVDAProblem)
    n, p, c = probdims(problem)
    T = floattype(problem)
    kernT = typeof(problem.kernel)
    respT = eltype(problem.Y)
    matT = eltype(problem.X)
    labelT = keytype(problem.label2vertex)
    kernel_info = kernT <: Nothing ? "linear classifier" : "nonlinear classifier ($(kernT))"
    print(io, "MVDAProblem{$(T)}")
    print(io, "\n  ∘ $(kernel_info)")
    print(io, "\n  ∘ $(n) sample(s) ($(respT))")
    print(io, "\n  ∘ $(p) feature(s) ($(matT))")
    print(io, "\n  ∘ $(c) categories ($(labelT))")
    print(io, "\n  ∘ intercept? $(problem.intercept)")
end

function maximal_deadzone(problem::MVDAProblem)
    c = problem.c
    ifelse(c == 2, 0.5, 1//2 * sqrt(2*c/(c-1)))
end

function set_initial_coefficients!(::Nothing, train_coeff, coeff, idx)
    copyto!(train_coeff, coeff)
end

function set_initial_coefficients!(::Kernel, train_coeff, coeff, idx)
    @views for (i, idx_i) in enumerate(idx)
        train_coeff[i, :] .= coeff[idx_i, :]
    end
end

function set_initial_coefficients!(train_problem::MVDAProblem, problem::MVDAProblem, idx)
    set_initial_coefficients!(train_problem.kernel, train_problem.coeff.all, problem.coeff.all, idx)
    set_initial_coefficients!(train_problem.kernel, train_problem.coeff_prev.all, problem.coeff_prev.all, idx)
end


text:
# This file was generated, do not modify it. # hide
@from_network MyTwoModelStack(regressor1=model1,
                              regressor2=model2,
                              judge=judge) <= yhat

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
module CodecBGZF

export
    BGZFCompressor,
    BGZFCompressorStream,
    BGZFDecompressor,
    BGZFDecompressorStream

import CodecZlib
import TranscodingStreams:
    TranscodingStreams,
    TranscodingStream,
    Memory,
    Error,
    initialize,
    finalize

# BGZF blocks are no larger than 64 KiB before and after compression.
const BGZF_MAX_BLOCK_SIZE = UInt(64 * 1024)

# BGZF_MAX_BLOCK_SIZE minus "margin for safety"
# NOTE: Data block will become slightly larger after deflation when bytes are
# randomly distributed.
const BGZF_SAFE_BLOCK_SIZE = UInt(BGZF_MAX_BLOCK_SIZE - 256)

include("compressor.jl")
include("decompressor.jl")

end # module


text:
"""
Julia package for implementing bundle methods
The current version has implemented a proximal bundle method.

minimize
  sum_{j=1}^N f_j(x)
subject to
  x in a feasible set of dimension n

f_j may be evaluated by calling evaluate_f.
"""

module BundleMethod

using JuMP

include("Parameters.jl")
include("Abstract.jl")
include("Model.jl")
include("Proximal.jl")
include("TrustRegion.jl")

end


text:
eye(n) = 1.0 * Matrix(I, n ,n);

export FilterParam, KRLSFilter

struct FilterParam
    order::Int
    fALD::Float64
    iDtMax::Int
end

mutable struct KRLSFilter
    mKinv::Array
    vAlpha::Array
    mP::Array
    mDict::Array
    m::Int
    mode::Char
    params::FilterParam
    function KRLSFilter(vx::Array, vy::Array, mode::Char, params::FilterParam)
        # Frst Check vx and vy shape, should be a vertical vector
        ktt = vx * vx'
        mKinv = 1 ./ ktt
        vAlpha = vx * vy'
        mP = eye(length(vx))
        mDict = vx
        m = length(vx)
        new(mKinv, vAlpha, mP, mDict, m, mode, params)
    end
end

@doc raw"""
	update(vx::Array, vy::Array, filter::KRLSFilter)

Update Filter intrinsic parameters for new observations.

# Examples
```julia-repl
julia> update(vx, vy, filter)
```
"""
function update(vx::Array, vy::Array, filter::KRLSFilter)
    at = filter.mKinv * kernel(vx, filter.m, filter.mode, filter.mDict);
    ktt = [filter.mDict vx]' * vx;
    deltat = abs(ktt[end] - dot(kernel(vx, filter.m, filter.mode, filter.mDict), at))
    # Dot is required here to ensure it produces a scalar

    rowC = size(filter.mDict, 2);

    if deltat > filter.params.fALD && rowC < filter.params.iDtMax
        filter.mDict = [filter.mDict vx];
        filter.mKinv = 1 / deltat * [deltat * filter.mKinv + at * at' -at; -at' 1];
        (r_mp, c_mp) = size(filter.mP);
        filter.mP = [filter.mP zeros(r_mp, 1); zeros(1, c_mp) 1];

        filter.vAlpha = [filter.vAlpha + at / deltat * (vy' - kernel(vx, filter.m, filter.mode,
            filter.mDict)' * filter.vAlpha); 1 / deltat * (vy' - kernel(vx, filter.m, filter.mode,
            filter.mDict)' * filter.vAlpha)];
        filter.m = filter.m + 1;
    else
        qt = filter.mP * at / (1 + (at' * filter.mP * at)[1]);
        filter.mP = filter.mP - qt * at' * filter.mP;
        filter.vAlpha = filter.vAlpha + filter.mKinv * qt * (vy' - kernel(vx, filter.m, filter.mode,
            filter.mDict)' * filter.vAlpha);
    end

    return filter;
end

function predict(vx::Array, filter::KRLSFilter)
    (kernel(vx, filter.m, filter.mode, filter.mDict)' * filter.vAlpha)';
end


text:
using Alazar
using Base.Test

@test Alazar.alazar_exception(512) == "ApiNoError"
@test Alazar.alazar_exception(579) == "ApiWaitTimeout"

p = PageAlignedVector{Int}(512)
@test eltype(p) == Int
@test Integer(pointer(p)) % Base.Mmap.PAGESIZE == 0
@test (p[:] = 1) == 1
@test @inferred(p[1]) == 1
@test all(p .== 1)
@test length(p) == 512
@test size(p) == (512,)
@test Base.IndexStyle(p) == Base.IndexLinear()

# Should assert n_buf >= 0
@test_throws AssertionError DMABufferVector(PageAlignedVector{Int}, -1, 10)

# Should assert bytes_buf * n_buf % sizeof(eltype(T)) == 0
@test_throws AssertionError DMABufferVector(PageAlignedVector{Int}, 3, 3)

# Should throw an error if bytes_buf % Base.Mmap.PAGESIZE != 0 && n_buf > 1
@test_throws ErrorException DMABufferVector(PageAlignedVector{Int}, 64, 10)

# Should throw an error if a non-page-aligned AbstractVector type is used.
# However, sometimes a Vector may be page-aligned, by chance. Not sure how to ensure
# this test will pass in general.
# @test_throws ErrorException DMABufferVector(Vector{Int}, 512, 1)

# Should throw an error if anything but an AbstractVector type is used.
@test_throws MethodError DMABufferVector(Array{Int,2}, 512, 1)
@test_throws MethodError DMABufferVector(Array{Int}, 512, 1)

dma = DMABufferVector(PageAlignedVector{Int}, Base.Mmap.PAGESIZE, 10)
@test @inferred(DMABufferVector(PageAlignedVector{Int}, Base.Mmap.PAGESIZE, 10)) isa
    DMABufferVector{Ptr{Int}, PageAlignedVector{Int}}
@test eltype(dma) == Ptr{Int}
@test @inferred(dma[1]) isa Ptr{Int}
@test length(dma) == 10
@test size(dma) == (10,)
@test Base.IndexStyle(dma) == Base.IndexLinear()


text: touch(joinpath(@__DIR__, "deps.jl"))

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
using Plots
using DifferentialEquations
using DiffEqUncertainty
#----------------------------------------------
# Specify and Solve model:

function fitzhugh_nagumo!(du, u, p, t)

    v, w = u
    a, b, I, R, τ = p

    du[1] = v - (v^3)/3 - w  + R * I

    du[2] = (v + a - (b * w)) / τ
end

u0 = [-2.0, 1.0]
p = [0.7, 0.8, 0.34, 1.0, 10]
tspan = (0.0, 100.0)
prob = ODEProblem(fitzhugh_nagumo!,u0,tspan,p)
sol = solve(prob)

plot(sol, vars =(1))

cb = ProbIntsUncertainty(0.1,1)
cb2 = ProbIntsUncertainty(0.05,1)
ensemble_prob = EnsembleProblem(prob)
sim = solve(ensemble_prob,Euler(),trajectories=100,callback=cb,dt=1/10)
sim2 = solve(ensemble_prob,Euler(),trajectories=50,callback=cb2,dt=1/100)

plot(sim, vars=(0,1), color="grey")
plot!(sim2, vars=(0,1), color="lightsalmon2")
plot!(sol, vars =(1), color="red3", linewidth =3)


text:
file_dir = Pkg.dir("D4M")*"/examples/2Apps/1EntityAnalysis/Entity.jld";

using JLD

E = load(file_dir)["E"]

E = logical(E)

x = "LOCATION/new york,";
p = "PERSON/michael chang,";
F = ( nocol(E[:,x]) & nocol(E[:,p]))' * E;
print(F' > 1 )

Fn = F./ sum(E,1)
print((Fn' > 0.02))