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

    #= ::Package:: =#

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

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sin(e+f*x)^m*(a+a*im*tan(e+f*x))^n=#

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

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

    @test_int [sin(x)^4/(im + tan(x)), x, 5, -((im*x)/16) - 1/(32*(im - tan(x))^2) - im/(8*(im - tan(x))) + im/(24*(im + tan(x))^3) - 5/(32*(im + tan(x))^2) - (3*im)/(16*(im + tan(x)))]
    @test_int [sin(x)^3/(im + tan(x)), x, 9, (1/3)*im*cos(x)^3 - (1/5)*im*cos(x)^5 + sin(x)^5/5]
    @test_int [sin(x)^2/(im + tan(x)), x, 5, -((im*x)/8) - im/(8*(im - tan(x))) - 1/(8*(im + tan(x))^2) - im/(4*(im + tan(x)))]
    @test_int [sin(x)^1/(im + tan(x)), x, 8, (1/3)*im*cos(x)^3 + sin(x)^3/3]
    @test_int [csc(x)^1/(im + tan(x)), x, 8, im*atanh(cos(x)) - im*cos(x) + sin(x)]
    @test_int [csc(x)^2/(im + tan(x)), x, 3, im*x + im*cot(x) + log(cos(x)) + log(tan(x))]
    @test_int [csc(x)^3/(im + tan(x)), x, 8, (-(1/2))*im*atanh(cos(x)) - csc(x) + (1/2)*im*cot(x)*csc(x)]
    @test_int [csc(x)^4/(im + tan(x)), x, 4, (-(1/2))*cot(x)^2 + (1/3)*im*cot(x)^3]
    @test_int [csc(x)^5/(im + tan(x)), x, 9, (-(1/8))*im*atanh(cos(x)) - (1/8)*im*cot(x)*csc(x) - csc(x)^3/3 + (1/4)*im*cot(x)*csc(x)^3]
    @test_int [csc(x)^6/(im + tan(x)), x, 4, (-(1/2))*cot(x)^2 + (1/3)*im*cot(x)^3 - cot(x)^4/4 + (1/5)*im*cot(x)^5]

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

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

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

    @test_int [sin(c + d*x)^5*(a + b*tan(c + d*x)), x, 8, (b*atanh(sin(c + d*x)))/d - (a*cos(c + d*x))/d + (2*a*cos(c + d*x)^3)/(3*d) - (a*cos(c + d*x)^5)/(5*d) - (b*sin(c + d*x))/d - (b*sin(c + d*x)^3)/(3*d) - (b*sin(c + d*x)^5)/(5*d)]
    @test_int [sin(c + d*x)^4*(a + b*tan(c + d*x)), x, 6, (3*a*x)/8 - (b*log(cos(c + d*x)))/d - (cos(c + d*x)*sin(c + d*x)^3*(a + b*tan(c + d*x)))/(4*d) - (cos(c + d*x)*sin(c + d*x)*(3*a + 4*b*tan(c + d*x)))/(8*d)]
    @test_int [sin(c + d*x)^3*(a + b*tan(c + d*x)), x, 8, (b*atanh(sin(c + d*x)))/d - (a*cos(c + d*x))/d + (a*cos(c + d*x)^3)/(3*d) - (b*sin(c + d*x))/d - (b*sin(c + d*x)^3)/(3*d)]
    @test_int [sin(c + d*x)^2*(a + b*tan(c + d*x)), x, 5, (a*x)/2 - (b*log(cos(c + d*x)))/d - (cos(c + d*x)*sin(c + d*x)*(a + b*tan(c + d*x)))/(2*d)]
    @test_int [sin(c + d*x)^1*(a + b*tan(c + d*x)), x, 6, (b*atanh(sin(c + d*x)))/d - (a*cos(c + d*x))/d - (b*sin(c + d*x))/d]
    @test_int [csc(c + d*x)^1*(a + b*tan(c + d*x)), x, 4, -((a*atanh(cos(c + d*x)))/d) + (b*atanh(sin(c + d*x)))/d]
    @test_int [csc(c + d*x)^2*(a + b*tan(c + d*x)), x, 3, -((a*cot(c + d*x))/d) + (b*log(tan(c + d*x)))/d]
    @test_int [csc(c + d*x)^3*(a + b*tan(c + d*x)), x, 7, -((a*atanh(cos(c + d*x)))/(2*d)) + (b*atanh(sin(c + d*x)))/d - (b*csc(c + d*x))/d - (a*cot(c + d*x)*csc(c + d*x))/(2*d)]
    @test_int [csc(c + d*x)^4*(a + b*tan(c + d*x)), x, 3, -((a*cot(c + d*x))/d) - (b*cot(c + d*x)^2)/(2*d) - (a*cot(c + d*x)^3)/(3*d) + (b*log(tan(c + d*x)))/d]
    @test_int [csc(c + d*x)^5*(a + b*tan(c + d*x)), x, 9, -((3*a*atanh(cos(c + d*x)))/(8*d)) + (b*atanh(sin(c + d*x)))/d - (b*csc(c + d*x))/d - (3*a*cot(c + d*x)*csc(c + d*x))/(8*d) - (b*csc(c + d*x)^3)/(3*d) - (a*cot(c + d*x)*csc(c + d*x)^3)/(4*d)]
    @test_int [csc(c + d*x)^6*(a + b*tan(c + d*x)), x, 3, -((a*cot(c + d*x))/d) - (b*cot(c + d*x)^2)/d - (2*a*cot(c + d*x)^3)/(3*d) - (b*cot(c + d*x)^4)/(4*d) - (a*cot(c + d*x)^5)/(5*d) + (b*log(tan(c + d*x)))/d]

    @test_int [sin(c + d*x)^4*(a + b*tan(c + d*x))^2, x, 8, (3/8)*(a^2 - 5*b^2)*x - (2*a*b*log(cos(c + d*x)))/d + (b^2*tan(c + d*x))/d + (cos(c + d*x)^2*(7*b - 5*a*tan(c + d*x))*(a + b*tan(c + d*x)))/(8*d) + (cos(c + d*x)^3*sin(c + d*x)*(a + b*tan(c + d*x))^2)/(4*d)]
    @test_int [sin(c + d*x)^3*(a + b*tan(c + d*x))^2, x, 11, (2*a*b*atanh(sin(c + d*x)))/d - (a^2*cos(c + d*x))/d + (2*b^2*cos(c + d*x))/d + (a^2*cos(c + d*x)^3)/(3*d) - (b^2*cos(c + d*x)^3)/(3*d) + (b^2*sec(c + d*x))/d - (2*a*b*sin(c + d*x))/d - (2*a*b*sin(c + d*x)^3)/(3*d)]
    @test_int [sin(c + d*x)^2*(a + b*tan(c + d*x))^2, x, 6, (1/2)*(a^2 - 3*b^2)*x - (2*a*b*log(cos(c + d*x)))/d + (3*b^2*tan(c + d*x))/(2*d) - (cos(c + d*x)*sin(c + d*x)*(a + b*tan(c + d*x))^2)/(2*d)]
    @test_int [sin(c + d*x)^1*(a + b*tan(c + d*x))^2, x, 9, (2*a*b*atanh(sin(c + d*x)))/d - (a^2*cos(c + d*x))/d + (b^2*cos(c + d*x))/d + (b^2*sec(c + d*x))/d - (2*a*b*sin(c + d*x))/d]
    @test_int [csc(c + d*x)^1*(a + b*tan(c + d*x))^2, x, 6, -((a^2*atanh(cos(c + d*x)))/d) + (2*a*b*atanh(sin(c + d*x)))/d + (b^2*sec(c + d*x))/d]
    @test_int [csc(c + d*x)^2*(a + b*tan(c + d*x))^2, x, 3, -((a^2*cot(c + d*x))/d) + (2*a*b*log(tan(c + d*x)))/d + (b^2*tan(c + d*x))/d]
    @test_int [csc(c + d*x)^3*(a + b*tan(c + d*x))^2, x, 10, -((a^2*atanh(cos(c + d*x)))/(2*d)) - (b^2*atanh(cos(c + d*x)))/d + (2*a*b*atanh(sin(c + d*x)))/d - (2*a*b*csc(c + d*x))/d - (a^2*cot(c + d*x)*csc(c + d*x))/(2*d) + (b^2*sec(c + d*x))/d]
    @test_int [csc(c + d*x)^4*(a + b*tan(c + d*x))^2, x, 3, -(((a^2 + b^2)*cot(c + d*x))/d) - (a*b*cot(c + d*x)^2)/d - (a^2*cot(c + d*x)^3)/(3*d) + (2*a*b*log(tan(c + d*x)))/d + (b^2*tan(c + d*x))/d]
    @test_int [csc(c + d*x)^5*(a + b*tan(c + d*x))^2, x, 13, -((3*a^2*atanh(cos(c + d*x)))/(8*d)) - (3*b^2*atanh(cos(c + d*x)))/(2*d) + (2*a*b*atanh(sin(c + d*x)))/d - (2*a*b*csc(c + d*x))/d - (3*a^2*cot(c + d*x)*csc(c + d*x))/(8*d) - (2*a*b*csc(c + d*x)^3)/(3*d) - (a^2*cot(c + d*x)*csc(c + d*x)^3)/(4*d) + (3*b^2*sec(c + d*x))/(2*d) - (b^2*csc(c + d*x)^2*sec(c + d*x))/(2*d)]
    @test_int [csc(c + d*x)^6*(a + b*tan(c + d*x))^2, x, 3, -(((a^2 + 2*b^2)*cot(c + d*x))/d) - (2*a*b*cot(c + d*x)^2)/d - ((2*a^2 + b^2)*cot(c + d*x)^3)/(3*d) - (a*b*cot(c + d*x)^4)/(2*d) - (a^2*cot(c + d*x)^5)/(5*d) + (2*a*b*log(tan(c + d*x)))/d + (b^2*tan(c + d*x))/d]

    @test_int [sin(c + d*x)^3*(a + b*tan(c + d*x))^3, x, 16, (3*a^2*b*atanh(sin(c + d*x)))/d - (5*b^3*atanh(sin(c + d*x)))/(2*d) - (a^3*cos(c + d*x))/d + (6*a*b^2*cos(c + d*x))/d + (a^3*cos(c + d*x)^3)/(3*d) - (a*b^2*cos(c + d*x)^3)/d + (3*a*b^2*sec(c + d*x))/d - (3*a^2*b*sin(c + d*x))/d + (5*b^3*sin(c + d*x))/(2*d) - (a^2*b*sin(c + d*x)^3)/d + (5*b^3*sin(c + d*x)^3)/(6*d) + (b^3*sin(c + d*x)^3*tan(c + d*x)^2)/(2*d)]
    @test_int [sin(c + d*x)^2*(a + b*tan(c + d*x))^3, x, 7, (1/2)*a*(a^2 - 9*b^2)*x - (b*(3*a^2 - 2*b^2)*log(cos(c + d*x)))/d + (9*a*b^2*tan(c + d*x))/(2*d) + (b^3*tan(c + d*x)^2)/d - (cos(c + d*x)*sin(c + d*x)*(a + b*tan(c + d*x))^3)/(2*d)]
    @test_int [sin(c + d*x)^1*(a + b*tan(c + d*x))^3, x, 13, (3*a^2*b*atanh(sin(c + d*x)))/d - (3*b^3*atanh(sin(c + d*x)))/(2*d) - (a^3*cos(c + d*x))/d + (3*a*b^2*cos(c + d*x))/d + (3*a*b^2*sec(c + d*x))/d - (3*a^2*b*sin(c + d*x))/d + (3*b^3*sin(c + d*x))/(2*d) + (b^3*sin(c + d*x)*tan(c + d*x)^2)/(2*d)]
    @test_int [csc(c + d*x)^1*(a + b*tan(c + d*x))^3, x, 8, -((a^3*atanh(cos(c + d*x)))/d) + (3*a^2*b*atanh(sin(c + d*x)))/d - (b^3*atanh(sin(c + d*x)))/(2*d) + (3*a*b^2*sec(c + d*x))/d + (b^3*sec(c + d*x)*tan(c + d*x))/(2*d)]
    @test_int [csc(c + d*x)^2*(a + b*tan(c + d*x))^3, x, 3, -((a^3*cot(c + d*x))/d) + (3*a^2*b*log(tan(c + d*x)))/d + (3*a*b^2*tan(c + d*x))/d + (b^3*tan(c + d*x)^2)/(2*d)]
    @test_int [csc(c + d*x)^3*(a + b*tan(c + d*x))^3, x, 12, -((a^3*atanh(cos(c + d*x)))/(2*d)) - (3*a*b^2*atanh(cos(c + d*x)))/d + (3*a^2*b*atanh(sin(c + d*x)))/d + (b^3*atanh(sin(c + d*x)))/(2*d) - (3*a^2*b*csc(c + d*x))/d - (a^3*cot(c + d*x)*csc(c + d*x))/(2*d) + (3*a*b^2*sec(c + d*x))/d + (b^3*sec(c + d*x)*tan(c + d*x))/(2*d)]
    @test_int [csc(c + d*x)^4*(a + b*tan(c + d*x))^3, x, 3, -((a*(a^2 + 3*b^2)*cot(c + d*x))/d) - (3*a^2*b*cot(c + d*x)^2)/(2*d) - (a^3*cot(c + d*x)^3)/(3*d) + (b*(3*a^2 + b^2)*log(tan(c + d*x)))/d + (3*a*b^2*tan(c + d*x))/d + (b^3*tan(c + d*x)^2)/(2*d)]
    @test_int [csc(c + d*x)^5*(a + b*tan(c + d*x))^3, x, 17, -((3*a^3*atanh(cos(c + d*x)))/(8*d)) - (9*a*b^2*atanh(cos(c + d*x)))/(2*d) + (3*a^2*b*atanh(sin(c + d*x)))/d + (3*b^3*atanh(sin(c + d*x)))/(2*d) - (3*a^2*b*csc(c + d*x))/d - (3*b^3*csc(c + d*x))/(2*d) - (3*a^3*cot(c + d*x)*csc(c + d*x))/(8*d) - (a^2*b*csc(c + d*x)^3)/d - (a^3*cot(c + d*x)*csc(c + d*x)^3)/(4*d) + (9*a*b^2*sec(c + d*x))/(2*d) - (3*a*b^2*csc(c + d*x)^2*sec(c + d*x))/(2*d) + (b^3*csc(c + d*x)*sec(c + d*x)^2)/(2*d)]
    @test_int [csc(c + d*x)^6*(a + b*tan(c + d*x))^3, x, 3, -((a*(a^2 + 6*b^2)*cot(c + d*x))/d) - (b*(6*a^2 + b^2)*cot(c + d*x)^2)/(2*d) - (a*(2*a^2 + 3*b^2)*cot(c + d*x)^3)/(3*d) - (3*a^2*b*cot(c + d*x)^4)/(4*d) - (a^3*cot(c + d*x)^5)/(5*d) + (b*(3*a^2 + 2*b^2)*log(tan(c + d*x)))/d + (3*a*b^2*tan(c + d*x))/d + (b^3*tan(c + d*x)^2)/(2*d)]

    @test_int [sin(c + d*x)^3*(a + b*tan(c + d*x))^4, x, 19, (4*a^3*b*atanh(sin(c + d*x)))/d - (10*a*b^3*atanh(sin(c + d*x)))/d - (a^4*cos(c + d*x))/d + (12*a^2*b^2*cos(c + d*x))/d - (3*b^4*cos(c + d*x))/d + (a^4*cos(c + d*x)^3)/(3*d) - (2*a^2*b^2*cos(c + d*x)^3)/d + (b^4*cos(c + d*x)^3)/(3*d) + (6*a^2*b^2*sec(c + d*x))/d - (3*b^4*sec(c + d*x))/d + (b^4*sec(c + d*x)^3)/(3*d) - (4*a^3*b*sin(c + d*x))/d + (10*a*b^3*sin(c + d*x))/d - (4*a^3*b*sin(c + d*x)^3)/(3*d) + (10*a*b^3*sin(c + d*x)^3)/(3*d) + (2*a*b^3*sin(c + d*x)^3*tan(c + d*x)^2)/d]
    @test_int [sin(c + d*x)^2*(a + b*tan(c + d*x))^4, x, 7, (1/2)*(a^4 - 18*a^2*b^2 + 5*b^4)*x - (4*a*b*(a^2 - 2*b^2)*log(cos(c + d*x)))/d + (b^2*(18*a^2 - 5*b^2)*tan(c + d*x))/(2*d) + (4*a*b^3*tan(c + d*x)^2)/d + (5*b^4*tan(c + d*x)^3)/(6*d) - (cos(c + d*x)*sin(c + d*x)*(a + b*tan(c + d*x))^4)/(2*d)]
    @test_int [sin(c + d*x)^1*(a + b*tan(c + d*x))^4, x, 16, (4*a^3*b*atanh(sin(c + d*x)))/d - (6*a*b^3*atanh(sin(c + d*x)))/d - (a^4*cos(c + d*x))/d + (6*a^2*b^2*cos(c + d*x))/d - (b^4*cos(c + d*x))/d + (6*a^2*b^2*sec(c + d*x))/d - (2*b^4*sec(c + d*x))/d + (b^4*sec(c + d*x)^3)/(3*d) - (4*a^3*b*sin(c + d*x))/d + (6*a*b^3*sin(c + d*x))/d + (2*a*b^3*sin(c + d*x)*tan(c + d*x)^2)/d]
    @test_int [csc(c + d*x)^1*(a + b*tan(c + d*x))^4, x, 10, -((a^4*atanh(cos(c + d*x)))/d) + (4*a^3*b*atanh(sin(c + d*x)))/d - (2*a*b^3*atanh(sin(c + d*x)))/d + (6*a^2*b^2*sec(c + d*x))/d - (b^4*sec(c + d*x))/d + (b^4*sec(c + d*x)^3)/(3*d) + (2*a*b^3*sec(c + d*x)*tan(c + d*x))/d]
    @test_int [csc(c + d*x)^2*(a + b*tan(c + d*x))^4, x, 3, -((a^4*cot(c + d*x))/d) + (4*a^3*b*log(tan(c + d*x)))/d + (6*a^2*b^2*tan(c + d*x))/d + (2*a*b^3*tan(c + d*x)^2)/d + (b^4*tan(c + d*x)^3)/(3*d)]
    @test_int [csc(c + d*x)^3*(a + b*tan(c + d*x))^4, x, 14, -((a^4*atanh(cos(c + d*x)))/(2*d)) - (6*a^2*b^2*atanh(cos(c + d*x)))/d + (4*a^3*b*atanh(sin(c + d*x)))/d + (2*a*b^3*atanh(sin(c + d*x)))/d - (4*a^3*b*csc(c + d*x))/d - (a^4*cot(c + d*x)*csc(c + d*x))/(2*d) + (6*a^2*b^2*sec(c + d*x))/d + (b^4*sec(c + d*x)^3)/(3*d) + (2*a*b^3*sec(c + d*x)*tan(c + d*x))/d]
    @test_int [csc(c + d*x)^4*(a + b*tan(c + d*x))^4, x, 3, -((a^2*(a^2 + 6*b^2)*cot(c + d*x))/d) - (2*a^3*b*cot(c + d*x)^2)/d - (a^4*cot(c + d*x)^3)/(3*d) + (4*a*b*(a^2 + b^2)*log(tan(c + d*x)))/d + (b^2*(6*a^2 + b^2)*tan(c + d*x))/d + (2*a*b^3*tan(c + d*x)^2)/d + (b^4*tan(c + d*x)^3)/(3*d)]
    @test_int [csc(c + d*x)^5*(a + b*tan(c + d*x))^4, x, 21, -((3*a^4*atanh(cos(c + d*x)))/(8*d)) - (9*a^2*b^2*atanh(cos(c + d*x)))/d - (b^4*atanh(cos(c + d*x)))/d + (4*a^3*b*atanh(sin(c + d*x)))/d + (6*a*b^3*atanh(sin(c + d*x)))/d - (4*a^3*b*csc(c + d*x))/d - (6*a*b^3*csc(c + d*x))/d - (3*a^4*cot(c + d*x)*csc(c + d*x))/(8*d) - (4*a^3*b*csc(c + d*x)^3)/(3*d) - (a^4*cot(c + d*x)*csc(c + d*x)^3)/(4*d) + (9*a^2*b^2*sec(c + d*x))/d + (b^4*sec(c + d*x))/d - (3*a^2*b^2*csc(c + d*x)^2*sec(c + d*x))/d + (2*a*b^3*csc(c + d*x)*sec(c + d*x)^2)/d + (b^4*sec(c + d*x)^3)/(3*d)]
    @test_int [csc(c + d*x)^6*(a + b*tan(c + d*x))^4, x, 3, -(((a^4 + 12*a^2*b^2 + b^4)*cot(c + d*x))/d) - (2*a*b*(2*a^2 + b^2)*cot(c + d*x)^2)/d - (2*a^2*(a^2 + 3*b^2)*cot(c + d*x)^3)/(3*d) - (a^3*b*cot(c + d*x)^4)/d - (a^4*cot(c + d*x)^5)/(5*d) + (4*a*b*(a^2 + 2*b^2)*log(tan(c + d*x)))/d + (2*b^2*(3*a^2 + b^2)*tan(c + d*x))/d + (2*a*b^3*tan(c + d*x)^2)/d + (b^4*tan(c + d*x)^3)/(3*d)]
    @test_int [csc(c + d*x)^7*(a + b*tan(c + d*x))^4, x, 25, -((5*a^4*atanh(cos(c + d*x)))/(16*d)) - (45*a^2*b^2*atanh(cos(c + d*x)))/(4*d) - (5*b^4*atanh(cos(c + d*x)))/(2*d) + (4*a^3*b*atanh(sin(c + d*x)))/d + (10*a*b^3*atanh(sin(c + d*x)))/d - (4*a^3*b*csc(c + d*x))/d - (10*a*b^3*csc(c + d*x))/d - (5*a^4*cot(c + d*x)*csc(c + d*x))/(16*d) - (4*a^3*b*csc(c + d*x)^3)/(3*d) - (10*a*b^3*csc(c + d*x)^3)/(3*d) - (5*a^4*cot(c + d*x)*csc(c + d*x)^3)/(24*d) - (4*a^3*b*csc(c + d*x)^5)/(5*d) - (a^4*cot(c + d*x)*csc(c + d*x)^5)/(6*d) + (45*a^2*b^2*sec(c + d*x))/(4*d) + (5*b^4*sec(c + d*x))/(2*d) - (15*a^2*b^2*csc(c + d*x)^2*sec(c + d*x))/(4*d) - (3*a^2*b^2*csc(c + d*x)^4*sec(c + d*x))/(2*d) + (2*a*b^3*csc(c + d*x)^3*sec(c + d*x)^2)/d + (5*b^4*sec(c + d*x)^3)/(6*d) - (b^4*csc(c + d*x)^2*sec(c + d*x)^3)/(2*d)]

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

    @test_int [sin(c + d*x)^5/(a + b*tan(c + d*x)), x, 13, (a^5*b*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/((a^2 + b^2)^(7/2)*d) + (a^3*b^2*cos(c + d*x))/((a^2 + b^2)^3*d) + (a*b^2*cos(c + d*x))/((a^2 + b^2)^2*d) - (a*cos(c + d*x))/((a^2 + b^2)*d) - (a*b^2*cos(c + d*x)^3)/(3*(a^2 + b^2)^2*d) + (2*a*cos(c + d*x)^3)/(3*(a^2 + b^2)*d) - (a*cos(c + d*x)^5)/(5*(a^2 + b^2)*d) + (a^4*b*sin(c + d*x))/((a^2 + b^2)^3*d) + (a^2*b*sin(c + d*x)^3)/(3*(a^2 + b^2)^2*d) + (b*sin(c + d*x)^5)/(5*(a^2 + b^2)*d)]
    @test_int [sin(c + d*x)^4/(a + b*tan(c + d*x)), x, 8, (a*(3*a^4 - 6*a^2*b^2 - b^4)*x)/(8*(a^2 + b^2)^3) + (a^4*b*log(a*cos(c + d*x) + b*sin(c + d*x)))/((a^2 + b^2)^3*d) + (cos(c + d*x)^4*(b + a*tan(c + d*x)))/(4*(a^2 + b^2)*d) - (cos(c + d*x)^2*(4*b*(2*a^2 + b^2) + a*(5*a^2 + b^2)*tan(c + d*x)))/(8*(a^2 + b^2)^2*d)]
    @test_int [sin(c + d*x)^3/(a + b*tan(c + d*x)), x, 10, (a^3*b*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/((a^2 + b^2)^(5/2)*d) + (a*b^2*cos(c + d*x))/((a^2 + b^2)^2*d) - (a*cos(c + d*x))/((a^2 + b^2)*d) + (a*cos(c + d*x)^3)/(3*(a^2 + b^2)*d) + (a^2*b*sin(c + d*x))/((a^2 + b^2)^2*d) + (b*sin(c + d*x)^3)/(3*(a^2 + b^2)*d)]
    @test_int [sin(c + d*x)^2/(a + b*tan(c + d*x)), x, 7, (a*(a^2 - b^2)*x)/(2*(a^2 + b^2)^2) + (a^2*b*log(a*cos(c + d*x) + b*sin(c + d*x)))/((a^2 + b^2)^2*d) - (cos(c + d*x)^2*(b + a*tan(c + d*x)))/(2*(a^2 + b^2)*d)]
    @test_int [sin(c + d*x)^1/(a + b*tan(c + d*x)), x, 6, (a*b*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/((a^2 + b^2)^(3/2)*d) - (a*cos(c + d*x))/((a^2 + b^2)*d) + (b*sin(c + d*x))/((a^2 + b^2)*d)]
    @test_int [csc(c + d*x)^1/(a + b*tan(c + d*x)), x, 6, -(atanh(cos(c + d*x))/(a*d)) + (b*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/(a*sqrt(a^2 + b^2)*d)]
    @test_int [csc(c + d*x)^2/(a + b*tan(c + d*x)), x, 3, -(cot(c + d*x)/(a*d)) - (b*log(tan(c + d*x)))/(a^2*d) + (b*log(a + b*tan(c + d*x)))/(a^2*d)]
    @test_int [csc(c + d*x)^3/(a + b*tan(c + d*x)), x, 15, -(atanh(cos(c + d*x))/(2*a*d)) - (b^2*atanh(cos(c + d*x)))/(a^3*d) + (b*sqrt(a^2 + b^2)*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/(a^3*d) + (b*csc(c + d*x))/(a^2*d) - (cot(c + d*x)*csc(c + d*x))/(2*a*d)]
    @test_int [csc(c + d*x)^4/(a + b*tan(c + d*x)), x, 3, -(((a^2 + b^2)*cot(c + d*x))/(a^3*d)) + (b*cot(c + d*x)^2)/(2*a^2*d) - cot(c + d*x)^3/(3*a*d) - (b*(a^2 + b^2)*log(tan(c + d*x)))/(a^4*d) + (b*(a^2 + b^2)*log(a + b*tan(c + d*x)))/(a^4*d)]
    #= [csc(c + d*x)^5/(a + b*tan(c + d*x)), x, 25, -((3*atanh(cos(c + d*x)))/(8*a*d)) - (3*b^2*atanh(cos(c + d*x)))/(2*a^3*d) - (b^4*atanh(cos(c + d*x)))/(a^5*d) - (b*atanh(sin(c + d*x)))/(a^2*d) - (b^3*atanh(sin(c + d*x)))/(a^4*d) + (b*(a^2 + b^2)*atanh(sin(c + d*x)))/(a^4*d) + (b*(a^2 + b^2)^(3/2)*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/(a^5*d) + (b*csc(c + d*x))/(a^2*d) + (3*b^3*csc(c + d*x))/(2*a^4*d) - (3*cot(c + d*x)*csc(c + d*x))/(8*a*d) + (b*csc(c + d*x)^3)/(3*a^2*d) - (cot(c + d*x)*csc(c + d*x)^3)/(4*a*d) + (3*b^2*sec(c + d*x))/(2*a^3*d) + (b^4*sec(c + d*x))/(a^5*d) - (b^2*(a^2 + b^2)*sec(c + d*x))/(a^5*d) - (b^2*csc(c + d*x)^2*sec(c + d*x))/(2*a^3*d) - (b^3*csc(c + d*x)*sec(c + d*x)^2)/(2*a^4*d) + (b^3*sec(c + d*x)*tan(c + d*x))/(2*a^4*d)] =#
    @test_int [csc(c + d*x)^6/(a + b*tan(c + d*x)), x, 3, -(((a^2 + b^2)^2*cot(c + d*x))/(a^5*d)) + (b*(2*a^2 + b^2)*cot(c + d*x)^2)/(2*a^4*d) - ((2*a^2 + b^2)*cot(c + d*x)^3)/(3*a^3*d) + (b*cot(c + d*x)^4)/(4*a^2*d) - cot(c + d*x)^5/(5*a*d) - (b*(a^2 + b^2)^2*log(tan(c + d*x)))/(a^6*d) + (b*(a^2 + b^2)^2*log(a + b*tan(c + d*x)))/(a^6*d)]

    @test_int [sin(c + d*x)^6/(a + b*tan(c + d*x))^2, x, 9, ((5*a^8 - 80*a^6*b^2 + 50*a^4*b^4 + 8*a^2*b^6 + b^8)*x)/(16*(a^2 + b^2)^5) + (2*a^5*b*(a^2 - 3*b^2)*log(a*cos(c + d*x) + b*sin(c + d*x)))/((a^2 + b^2)^5*d) - (a^6*b)/((a^2 + b^2)^4*d*(a + b*tan(c + d*x))) - (cos(c + d*x)^6*(2*a*b + (a^2 - b^2)*tan(c + d*x)))/(6*(a^2 + b^2)^2*d) + (cos(c + d*x)^4*(12*a*b*(3*a^2 + b^2) + (13*a^4 - 18*a^2*b^2 - 7*b^4)*tan(c + d*x)))/(24*(a^2 + b^2)^3*d) - (cos(c + d*x)^2*(48*a^5*b + (11*a^6 - 43*a^4*b^2 - 7*a^2*b^4 - b^6)*tan(c + d*x)))/(16*(a^2 + b^2)^4*d)]
    @test_int [sin(c + d*x)^4/(a + b*tan(c + d*x))^2, x, 8, ((3*a^6 - 33*a^4*b^2 + 13*a^2*b^4 + b^6)*x)/(8*(a^2 + b^2)^4) + (2*a^3*b*(a^2 - 2*b^2)*log(a*cos(c + d*x) + b*sin(c + d*x)))/((a^2 + b^2)^4*d) - (a^4*b)/((a^2 + b^2)^3*d*(a + b*tan(c + d*x))) + (cos(c + d*x)^4*(2*a*b + (a^2 - b^2)*tan(c + d*x)))/(4*(a^2 + b^2)^2*d) - (cos(c + d*x)^2*(16*a^3*b + (5*a^4 - 12*a^2*b^2 - b^4)*tan(c + d*x)))/(8*(a^2 + b^2)^3*d)]
    @test_int [sin(c + d*x)^2/(a + b*tan(c + d*x))^2, x, 7, ((a^4 - 6*a^2*b^2 + b^4)*x)/(2*(a^2 + b^2)^3) + (2*a*b*(a^2 - b^2)*log(a*cos(c + d*x) + b*sin(c + d*x)))/((a^2 + b^2)^3*d) - (a^2*b)/((a^2 + b^2)^2*d*(a + b*tan(c + d*x))) - (cos(c + d*x)^2*(2*a*b + (a^2 - b^2)*tan(c + d*x)))/(2*(a^2 + b^2)^2*d)]
    @test_int [csc(c + d*x)^2/(a + b*tan(c + d*x))^2, x, 3, -(cot(c + d*x)/(a^2*d)) - (2*b*log(tan(c + d*x)))/(a^3*d) + (2*b*log(a + b*tan(c + d*x)))/(a^3*d) - b/(a^2*d*(a + b*tan(c + d*x)))]
    @test_int [csc(c + d*x)^4/(a + b*tan(c + d*x))^2, x, 3, -(((a^2 + 3*b^2)*cot(c + d*x))/(a^4*d)) + (b*cot(c + d*x)^2)/(a^3*d) - cot(c + d*x)^3/(3*a^2*d) - (2*b*(a^2 + 2*b^2)*log(tan(c + d*x)))/(a^5*d) + (2*b*(a^2 + 2*b^2)*log(a + b*tan(c + d*x)))/(a^5*d) - (b*(a^2 + b^2))/(a^4*d*(a + b*tan(c + d*x)))]
    @test_int [csc(c + d*x)^6/(a + b*tan(c + d*x))^2, x, 3, -(((a^2 + b^2)*(a^2 + 5*b^2)*cot(c + d*x))/(a^6*d)) + (2*b*(a^2 + b^2)*cot(c + d*x)^2)/(a^5*d) - ((2*a^2 + 3*b^2)*cot(c + d*x)^3)/(3*a^4*d) + (b*cot(c + d*x)^4)/(2*a^3*d) - cot(c + d*x)^5/(5*a^2*d) - (2*b*(a^2 + b^2)*(a^2 + 3*b^2)*log(tan(c + d*x)))/(a^7*d) + (2*b*(a^2 + b^2)*(a^2 + 3*b^2)*log(a + b*tan(c + d*x)))/(a^7*d) - (b*(a^2 + b^2)^2)/(a^6*d*(a + b*tan(c + d*x)))]
    #=
    @test_int [sin(c + d*x)^5/(a + b*tan(c + d*x))^2, x, 0, (5*sec(c + d*x)^2*(a*cos(c + d*x) + b*sin(c + d*x))*(a*(a^2 + b^2) + 2*b*sqrt(a^2 + b^2)*atanh((-b + a*tan((1/2)*(c + d*x)))/sqrt(a^2 + b^2))*(a*cos(c + d*x) + b*sin(c + d*x))))/(64*(a^2 + b^2)^2*d*(a + b*tan(c + d*x))^2) + (1/(64*(a^2 + b^2)^3*d*(a + b*tan(c + d*x))^2))*(sec(c + d*x)^2*(a*cos(c + d*x) + b*sin(c + d*x))*(6*b*(-3*a^2 + b^2)*sqrt(a^2 + b^2)*atanh((-b + a*tan((1/2)*(c + d*x)))/sqrt(a^2 + b^2))*(a*cos(c + d*x) + b*sin(c + d*x)) - (a^2 + b^2)*(3*a*(a^2 - 3*b^2) + 2*a*(a^2 + b^2)*cos(2*(c + d*x)) - 2*b*(a^2 + b^2)*sin(2*(c + d*x))))) - (3*sec(c + d*x)^2*(a*cos(c + d*x) + b*sin(c + d*x))^2*((10*b*(5*a^4 - 10*a^2*b^2 + b^4)*atanh((-b + a*tan((1/2)*(c + d*x)))/sqrt(a^2 + b^2)))/(a^2 + b^2)^(7/2) + (8*(a^4 - 6*a^2*b^2 + b^4)*cos(c + d*x))/(a^2 + b^2)^3 - (4*(a^2 - b^2)*cos(3*(c + d*x)))/(3*(a^2 + b^2)^2) - (32*a*b*(a^2 - b^2)*sin(c + d*x))/(a^2 + b^2)^3 + (a*(a^4 - 10*a^2*b^2 + 5*b^4))/((a^2 + b^2)^3*(a*cos(c + d*x) + b*sin(c + d*x))) + (8*a*b*sin(3*(c + d*x)))/(3*(a^2 + b^2)^2)))/(64*d*(a + b*tan(c + d*x))^2) + (1/(64*d*(a + b*tan(c + d*x))^2))*(sec(c + d*x)^2*(a*cos(c + d*x) + b*sin(c + d*x))^2*((14*b*(-7*a^6 + 35*a^4*b^2 - 21*a^2*b^4 + b^6)*atanh((-b + a*tan((1/2)*(c + d*x)))/sqrt(a^2 + b^2)))/(a^2 + b^2)^(9/2) - (12*(a^6 - 15*a^4*b^2 + 15*a^2*b^4 - b^6)*cos(c + d*x))/(a^2 + b^2)^4 + (8*(a^4 - 6*a^2*b^2 + b^4)*cos(3*(c + d*x)))/(3*(a^2 + b^2)^3) - (4*(a^2 - b^2)*cos(5*(c + d*x)))/(5*(a^2 + b^2)^2) + (24*a*b*(3*a^4 - 10*a^2*b^2 + 3*b^4)*sin(c + d*x))/(a^2 + b^2)^4 - (a*(a^6 - 21*a^4*b^2 + 35*a^2*b^4 - 7*b^6))/((a^2 + b^2)^4*(a*cos(c + d*x) + b*sin(c + d*x))) - (32*a*b*(a^2 - b^2)*sin(3*(c + d*x)))/(3*(a^2 + b^2)^3) + (8*a*b*sin(5*(c + d*x)))/(5*(a^2 + b^2)^2)))]
    @test_int [sin(c + d*x)^3/(a + b*tan(c + d*x))^2, x, 31, (2*a^4*b*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/((a^2 + b^2)^(7/2)*d) - (3*a^2*b^3*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/((a^2 + b^2)^(7/2)*d) + (4*a^2*b^2*cos(c + d*x))/((a^2 + b^2)^3*d) - (a^2*cos(c + d*x))/((a^2 + b^2)^2*d) + (a^2*cos(c + d*x)^3)/(3*(a^2 + b^2)^2*d) - (b^2*cos(c + d*x)^3)/(3*(a^2 + b^2)^2*d) + (2*a^3*b*sin(c + d*x))/((a^2 + b^2)^3*d) - (2*a*b^3*sin(c + d*x))/((a^2 + b^2)^3*d) + (2*a*b*sin(c + d*x)^3)/(3*(a^2 + b^2)^2*d) + (a^3*b^2)/((a^2 + b^2)^3*d*(a*cos(c + d*x) + b*sin(c + d*x)))]
    @test_int [sin(c + d*x)^1/(a + b*tan(c + d*x))^2, x, 14, (2*a^2*b*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/((a^2 + b^2)^(5/2)*d) - (b^3*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/((a^2 + b^2)^(5/2)*d) - (a^2*cos(c + d*x))/((a^2 + b^2)^2*d) + (b^2*cos(c + d*x))/((a^2 + b^2)^2*d) + (2*a*b*sin(c + d*x))/((a^2 + b^2)^2*d) + (a*b^2)/((a^2 + b^2)^2*d*(a*cos(c + d*x) + b*sin(c + d*x)))]
    @test_int [csc(c + d*x)^1/(a + b*tan(c + d*x))^2, x, 9, -(atanh(cos(c + d*x))/(a^2*d)) + (b*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/((a^2 + b^2)^(3/2)*d) + (b*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/(a^2*sqrt(a^2 + b^2)*d) + b^2/(a*(a^2 + b^2)*d*(a*cos(c + d*x) + b*sin(c + d*x)))]
    @test_int [csc(c + d*x)^3/(a + b*tan(c + d*x))^2, x, 17, -(atanh(cos(c + d*x))/(2*a^2*d)) - (3*b^2*atanh(cos(c + d*x)))/(a^4*d) - (b*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/(a^2*sqrt(a^2 + b^2)*d) + (3*b*sqrt(a^2 + b^2)*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/(a^4*d) + (2*b*csc(c + d*x))/(a^3*d) - (cot(c + d*x)*csc(c + d*x))/(2*a^2*d) + b^2/(a^3*d*(a*cos(c + d*x) + b*sin(c + d*x)))]
    =#

    @test_int [sin(c + d*x)^6/(a + b*tan(c + d*x))^3, x, 9, (a*(5*a^8 - 180*a^6*b^2 + 390*a^4*b^4 - 68*a^2*b^6 - 3*b^8)*x)/(16*(a^2 + b^2)^6) + (a^4*b*(3*a^4 - 22*a^2*b^2 + 15*b^4)*log(a*cos(c + d*x) + b*sin(c + d*x)))/((a^2 + b^2)^6*d) - (a^6*b)/(2*(a^2 + b^2)^4*d*(a + b*tan(c + d*x))^2) - (2*a^5*b*(a^2 - 3*b^2))/((a^2 + b^2)^5*d*(a + b*tan(c + d*x))) - (cos(c + d*x)^6*(b*(3*a^2 - b^2) + a*(a^2 - 3*b^2)*tan(c + d*x)))/(6*(a^2 + b^2)^3*d) + (cos(c + d*x)^4*(6*b*(9*a^4 - 4*a^2*b^2 - b^4) + a*(13*a^4 - 62*a^2*b^2 - 3*b^4)*tan(c + d*x)))/(24*(a^2 + b^2)^4*d) - (a*cos(c + d*x)^2*(24*a^3*b*(3*a^2 - 5*b^2) + (11*a^6 - 119*a^4*b^2 + 65*a^2*b^4 + 3*b^6)*tan(c + d*x)))/(16*(a^2 + b^2)^5*d)]
    @test_int [sin(c + d*x)^4/(a + b*tan(c + d*x))^3, x, 8, (3*a*(a^6 - 25*a^4*b^2 + 35*a^2*b^4 - 3*b^6)*x)/(8*(a^2 + b^2)^5) + (3*a^2*b*(a^4 - 5*a^2*b^2 + 2*b^4)*log(a*cos(c + d*x) + b*sin(c + d*x)))/((a^2 + b^2)^5*d) - (a^4*b)/(2*(a^2 + b^2)^3*d*(a + b*tan(c + d*x))^2) - (2*a^3*b*(a^2 - 2*b^2))/((a^2 + b^2)^4*d*(a + b*tan(c + d*x))) + (cos(c + d*x)^4*(b*(3*a^2 - b^2) + a*(a^2 - 3*b^2)*tan(c + d*x)))/(4*(a^2 + b^2)^3*d) - (a*cos(c + d*x)^2*(24*a*b*(a^2 - b^2) + (5*a^4 - 34*a^2*b^2 + 9*b^4)*tan(c + d*x)))/(8*(a^2 + b^2)^4*d)]
    @test_int [sin(c + d*x)^2/(a + b*tan(c + d*x))^3, x, 7, (a*(a^4 - 14*a^2*b^2 + 9*b^4)*x)/(2*(a^2 + b^2)^4) + (b*(3*a^4 - 8*a^2*b^2 + b^4)*log(a*cos(c + d*x) + b*sin(c + d*x)))/((a^2 + b^2)^4*d) - (a^2*b)/(2*(a^2 + b^2)^2*d*(a + b*tan(c + d*x))^2) - (2*a*b*(a^2 - b^2))/((a^2 + b^2)^3*d*(a + b*tan(c + d*x))) - (cos(c + d*x)^2*(b*(3*a^2 - b^2) + a*(a^2 - 3*b^2)*tan(c + d*x)))/(2*(a^2 + b^2)^3*d)]
    @test_int [csc(c + d*x)^2/(a + b*tan(c + d*x))^3, x, 3, -(cot(c + d*x)/(a^3*d)) - (3*b*log(tan(c + d*x)))/(a^4*d) + (3*b*log(a + b*tan(c + d*x)))/(a^4*d) - b/(2*a^2*d*(a + b*tan(c + d*x))^2) - (2*b)/(a^3*d*(a + b*tan(c + d*x)))]
    @test_int [csc(c + d*x)^4/(a + b*tan(c + d*x))^3, x, 3, -(((a^2 + 6*b^2)*cot(c + d*x))/(a^5*d)) + (3*b*cot(c + d*x)^2)/(2*a^4*d) - cot(c + d*x)^3/(3*a^3*d) - (b*(3*a^2 + 10*b^2)*log(tan(c + d*x)))/(a^6*d) + (b*(3*a^2 + 10*b^2)*log(a + b*tan(c + d*x)))/(a^6*d) - (b*(a^2 + b^2))/(2*a^4*d*(a + b*tan(c + d*x))^2) - (2*b*(a^2 + 2*b^2))/(a^5*d*(a + b*tan(c + d*x)))]
    @test_int [csc(c + d*x)^6/(a + b*tan(c + d*x))^3, x, 3, -(((a^4 + 12*a^2*b^2 + 15*b^4)*cot(c + d*x))/(a^7*d)) + (b*(3*a^2 + 5*b^2)*cot(c + d*x)^2)/(a^6*d) - (2*(a^2 + 3*b^2)*cot(c + d*x)^3)/(3*a^5*d) + (3*b*cot(c + d*x)^4)/(4*a^4*d) - cot(c + d*x)^5/(5*a^3*d) - (b*(3*a^4 + 20*a^2*b^2 + 21*b^4)*log(tan(c + d*x)))/(a^8*d) + (b*(3*a^4 + 20*a^2*b^2 + 21*b^4)*log(a + b*tan(c + d*x)))/(a^8*d) - (b*(a^2 + b^2)^2)/(2*a^6*d*(a + b*tan(c + d*x))^2) - (2*b*(a^2 + b^2)*(a^2 + 3*b^2))/(a^7*d*(a + b*tan(c + d*x)))]
    #=
    @test_int [sin(c + d*x)^5/(a + b*tan(c + d*x))^3, x, 0, (3*sec(c + d*x)^3*(a*cos(c + d*x) + b*sin(c + d*x))*(6*a*b*sqrt(a^2 + b^2)*atanh((-b + a*tan((1/2)*(c + d*x)))/sqrt(a^2 + b^2))*(a*cos(c + d*x) + b*sin(c + d*x))^2 + (a^2 + b^2)*(a*(2*a^2 - b^2)*cos(c + d*x) + b*(a^2 - 2*b^2)*sin(c + d*x))))/(64*(a^2 + b^2)^3*d*(a + b*tan(c + d*x))^3) + (1/(64*d*(a + b*tan(c + d*x))^3))*(sec(c + d*x)^3*(a*cos(c + d*x) + b*sin(c + d*x))^3*(-((70*a*b*(3*a^4 - 10*a^2*b^2 + 3*b^4)*atanh((-b + a*tan((1/2)*(c + d*x)))/sqrt(a^2 + b^2)))/(a^2 + b^2)^(9/2)) - (24*a*(a^4 - 10*a^2*b^2 + 5*b^4)*cos(c + d*x))/(a^2 + b^2)^4 + (8*a*(a^2 - 3*b^2)*cos(3*(c + d*x)))/(3*(a^2 + b^2)^3) + (24*b*(5*a^4 - 10*a^2*b^2 + b^4)*sin(c + d*x))/(a^2 + b^2)^4 - (b*(-3*a^4 + 10*a^2*b^2 - 3*b^4)*sin(c + d*x))/((a^2 + b^2)^3*(a*cos(c + d*x) + b*sin(c + d*x))^2) - (6*a^6 - 87*a^4*b^2 + 80*a^2*b^4 - 3*b^6)/((a^2 + b^2)^4*(a*cos(c + d*x) + b*sin(c + d*x))) + (8*b*(-3*a^2 + b^2)*sin(3*(c + d*x)))/(3*(a^2 + b^2)^3))) + (sec(c + d*x)^3*(a*cos(c + d*x) + b*sin(c + d*x))^3*((60*a*b*(a^2 - b^2)*atanh((-b + a*tan((1/2)*(c + d*x)))/sqrt(a^2 + b^2)))/(a^2 + b^2)^(7/2) + (2*(5*a^3*(a^2 - 5*b^2)*cos(c + d*x) + a*(a^2 + b^2)^2*cos(3*(c + d*x)) - b*(-5*b^2*(-5*a^2 + b^2)*sin(c + d*x) + (a^2 + b^2)^2*sin(3*(c + d*x)))))/((a^2 + b^2)^3*(a*cos(c + d*x) + b*sin(c + d*x))^2)))/(64*d*(a + b*tan(c + d*x))^3) - (1/(320*(a^2 + b^2)^6*d*(a + b*tan(c + d*x))^3))*(sec(c + d*x)^3*(a*cos(c + d*x) + b*sin(c + d*x))*(1260*a*b*sqrt(a^2 + b^2)*(a^6 - 7*a^4*b^2 + 7*a^2*b^4 - b^6)*atanh((-b + a*tan((1/2)*(c + d*x)))/sqrt(a^2 + b^2))*(a*cos(c + d*x) + b*sin(c + d*x))^2 + (a^2 + b^2)*(105*a*(a^8 - 26*a^6*b^2 + 56*a^4*b^4 - 14*a^2*b^6 - b^8)*cos(c + d*x) + 21*a*(a^2 + b^2)^2*(a^4 - 10*a^2*b^2 + 5*b^4)*cos(3*(c + d*x)) - 3*a^9*cos(5*(c + d*x)) + 18*a^5*b^4*cos(5*(c + d*x)) + 24*a^3*b^6*cos(5*(c + d*x)) + 9*a*b^8*cos(5*(c + d*x)) + a^9*cos(7*(c + d*x)) + 4*a^7*b^2*cos(7*(c + d*x)) + 6*a^5*b^4*cos(7*(c + d*x)) + 4*a^3*b^6*cos(7*(c + d*x)) + a*b^8*cos(7*(c + d*x)) - 105*a^8*b*sin(c + d*x) - 1470*a^6*b^3*sin(c + d*x) + 5880*a^4*b^5*sin(c + d*x) - 2730*a^2*b^7*sin(c + d*x) + 105*b^9*sin(c + d*x) - 105*a^8*b*sin(3*(c + d*x)) + 294*a^4*b^5*sin(3*(c + d*x)) + 168*a^2*b^7*sin(3*(c + d*x)) - 21*b^9*sin(3*(c + d*x)) + 9*a^8*b*sin(5*(c + d*x)) + 24*a^6*b^3*sin(5*(c + d*x)) + 18*a^4*b^5*sin(5*(c + d*x)) - 3*b^9*sin(5*(c + d*x)) - a^8*b*sin(7*(c + d*x)) - 4*a^6*b^3*sin(7*(c + d*x)) - 6*a^4*b^5*sin(7*(c + d*x)) - 4*a^2*b^7*sin(7*(c + d*x)) - b^9*sin(7*(c + d*x)))))]
    @test_int [sin(c + d*x)^3/(a + b*tan(c + d*x))^3, x, 99, (3*a^5*b*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/((a^2 + b^2)^(9/2)*d) - (23*a^3*b^3*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/(2*(a^2 + b^2)^(9/2)*d) + (3*a*b^5*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/((a^2 + b^2)^(9/2)*d) + (9*a^3*b^2*cos(c + d*x))/((a^2 + b^2)^4*d) - (3*a*b^4*cos(c + d*x))/((a^2 + b^2)^4*d) - (a^3*cos(c + d*x))/((a^2 + b^2)^3*d) + (a^3*cos(c + d*x)^3)/(3*(a^2 + b^2)^3*d) - (a*b^2*cos(c + d*x)^3)/((a^2 + b^2)^3*d) + (3*a^4*b*sin(c + d*x))/((a^2 + b^2)^4*d) - (9*a^2*b^3*sin(c + d*x))/((a^2 + b^2)^4*d) + (b^3*sin(c + d*x))/((a^2 + b^2)^3*d) + (a^2*b*sin(c + d*x)^3)/((a^2 + b^2)^3*d) - (b^3*sin(c + d*x)^3)/(3*(a^2 + b^2)^3*d) + (a^3*b^3*(b*cos(c + d*x) - a*sin(c + d*x)))/(2*(a^2 + b^2)^4*d*(a*cos(c + d*x) + b*sin(c + d*x))^2) + (3*a^4*b^2)/((a^2 + b^2)^4*d*(a*cos(c + d*x) + b*sin(c + d*x))) - (3*a^2*b^4)/((a^2 + b^2)^4*d*(a*cos(c + d*x) + b*sin(c + d*x)))]
    @test_int [sin(c + d*x)^1/(a + b*tan(c + d*x))^3, x, 0, -((1/(4*(a^2 + b^2)^(7/2)*d*(a + b*tan(c + d*x))^3))*(sec(c + d*x)^3*(a*cos(c + d*x) + b*sin(c + d*x))*(12*a*b*(2*a^2 - 3*b^2)*atanh((-b + a*tan((1/2)*(c + d*x)))/sqrt(a^2 + b^2))*(a*cos(c + d*x) + b*sin(c + d*x))^2 + sqrt(a^2 + b^2)*(a*(3*a^4 - 26*a^2*b^2 + b^4)*cos(c + d*x) + a*(a^2 + b^2)^2*cos(3*(c + d*x)) - 2*b*(a^4 + 13*a^2*b^2 - 3*b^4 + (a^2 + b^2)^2*cos(2*(c + d*x)))*sin(c + d*x)))))]
    @test_int [csc(c + d*x)^1/(a + b*tan(c + d*x))^3, x, 0, (1/(2*a^3*d*(a + b*tan(c + d*x))^3))*(sec(c + d*x)^2*(a*cos(c + d*x) + b*sin(c + d*x))*((3*a*b^2*(2*a^2 + b^2)*sec(c + d*x)*(a*cos(c + d*x) + b*sin(c + d*x)))/(a^2 + b^2)^2 - (2*b*(6*a^4 + 5*a^2*b^2 + 2*b^4)*atanh((-b + a*tan((1/2)*(c + d*x)))/sqrt(a^2 + b^2))*sec(c + d*x)*(a*cos(c + d*x) + b*sin(c + d*x))^2)/(a^2 + b^2)^(5/2) - 2*log(cos((1/2)*(c + d*x)))*sec(c + d*x)*(a*cos(c + d*x) + b*sin(c + d*x))^2 + 2*log(sin((1/2)*(c + d*x)))*sec(c + d*x)*(a*cos(c + d*x) + b*sin(c + d*x))^2 - (a*b^3*tan(c + d*x))/(a^2 + b^2)))]
    @test_int [csc(c + d*x)^3/(a + b*tan(c + d*x))^3, x, 20, -(atanh(cos(c + d*x))/(2*a^3*d)) - (6*b^2*atanh(cos(c + d*x)))/(a^5*d) + (b^3*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/(2*a^3*(a^2 + b^2)^(3/2)*d) - (3*b*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/(a^3*sqrt(a^2 + b^2)*d) + (6*b*sqrt(a^2 + b^2)*atanh((b*cos(c + d*x) - a*sin(c + d*x))/sqrt(a^2 + b^2)))/(a^5*d) + (3*b*csc(c + d*x))/(a^4*d) - (cot(c + d*x)*csc(c + d*x))/(2*a^3*d) + (b^3*(b*cos(c + d*x) - a*sin(c + d*x)))/(2*a^3*(a^2 + b^2)*d*(a*cos(c + d*x) + b*sin(c + d*x))^2) + (3*b^2)/(a^4*d*(a*cos(c + d*x) + b*sin(c + d*x)))]
    =#

    @test_int [sin(c + d*x)^4/(a + b*tan(c + d*x))^4, x, 8, ((3*a^8 - 132*a^6*b^2 + 370*a^4*b^4 - 132*a^2*b^6 + 3*b^8)*x)/(8*(a^2 + b^2)^6) + (4*a*b*(a^2 - b^2)*(a^4 - 8*a^2*b^2 + b^4)*log(a*cos(c + d*x) + b*sin(c + d*x)))/((a^2 + b^2)^6*d) - (a^4*b)/(3*(a^2 + b^2)^3*d*(a + b*tan(c + d*x))^3) - (a^3*b*(a^2 - 2*b^2))/((a^2 + b^2)^4*d*(a + b*tan(c + d*x))^2) - (3*a^2*b*(a^4 - 5*a^2*b^2 + 2*b^4))/((a^2 + b^2)^5*d*(a + b*tan(c + d*x))) + (cos(c + d*x)^4*(4*a*b*(a^2 - b^2) + (a^4 - 6*a^2*b^2 + b^4)*tan(c + d*x)))/(4*(a^2 + b^2)^4*d) - (cos(c + d*x)^2*(16*a*b*(2*a^4 - 5*a^2*b^2 + b^4) + (5*a^6 - 65*a^4*b^2 + 55*a^2*b^4 - 3*b^6)*tan(c + d*x)))/(8*(a^2 + b^2)^5*d)]
    @test_int [sin(c + d*x)^2/(a + b*tan(c + d*x))^4, x, 7, ((a^6 - 25*a^4*b^2 + 35*a^2*b^4 - 3*b^6)*x)/(2*(a^2 + b^2)^5) + (4*a*b*(a^4 - 5*a^2*b^2 + 2*b^4)*log(a*cos(c + d*x) + b*sin(c + d*x)))/((a^2 + b^2)^5*d) - (a^2*b)/(3*(a^2 + b^2)^2*d*(a + b*tan(c + d*x))^3) - (a*b*(a^2 - b^2))/((a^2 + b^2)^3*d*(a + b*tan(c + d*x))^2) - (b*(3*a^4 - 8*a^2*b^2 + b^4))/((a^2 + b^2)^4*d*(a + b*tan(c + d*x))) - (cos(c + d*x)^2*(4*a*b*(a^2 - b^2) + (a^4 - 6*a^2*b^2 + b^4)*tan(c + d*x)))/(2*(a^2 + b^2)^4*d)]
    @test_int [csc(c + d*x)^2/(a + b*tan(c + d*x))^4, x, 3, -(cot(c + d*x)/(a^4*d)) - (4*b*log(tan(c + d*x)))/(a^5*d) + (4*b*log(a + b*tan(c + d*x)))/(a^5*d) - b/(3*a^2*d*(a + b*tan(c + d*x))^3) - b/(a^3*d*(a + b*tan(c + d*x))^2) - (3*b)/(a^4*d*(a + b*tan(c + d*x)))]
    @test_int [csc(c + d*x)^4/(a + b*tan(c + d*x))^4, x, 3, -(((a^2 + 10*b^2)*cot(c + d*x))/(a^6*d)) + (2*b*cot(c + d*x)^2)/(a^5*d) - cot(c + d*x)^3/(3*a^4*d) - (4*b*(a^2 + 5*b^2)*log(tan(c + d*x)))/(a^7*d) + (4*b*(a^2 + 5*b^2)*log(a + b*tan(c + d*x)))/(a^7*d) - (b*(a^2 + b^2))/(3*a^4*d*(a + b*tan(c + d*x))^3) - (b*(a^2 + 2*b^2))/(a^5*d*(a + b*tan(c + d*x))^2) - (b*(3*a^2 + 10*b^2))/(a^6*d*(a + b*tan(c + d*x)))]
    @test_int [csc(c + d*x)^6/(a + b*tan(c + d*x))^4, x, 3, -(((a^4 + 20*a^2*b^2 + 35*b^4)*cot(c + d*x))/(a^8*d)) + (2*b*(2*a^2 + 5*b^2)*cot(c + d*x)^2)/(a^7*d) - (2*(a^2 + 5*b^2)*cot(c + d*x)^3)/(3*a^6*d) + (b*cot(c + d*x)^4)/(a^5*d) - cot(c + d*x)^5/(5*a^4*d) - (4*b*(a^4 + 10*a^2*b^2 + 14*b^4)*log(tan(c + d*x)))/(a^9*d) + (4*b*(a^4 + 10*a^2*b^2 + 14*b^4)*log(a + b*tan(c + d*x)))/(a^9*d) - (b*(a^2 + b^2)^2)/(3*a^6*d*(a + b*tan(c + d*x))^3) - (b*(a^2 + b^2)*(a^2 + 3*b^2))/(a^7*d*(a + b*tan(c + d*x))^2) - (b*(3*a^4 + 20*a^2*b^2 + 21*b^4))/(a^8*d*(a + b*tan(c + d*x)))]

    #= Hangs*Mathematica*6 & 7: =#
    @test_int [csc(x)/(1 + tan(x)), x, 6, -atanh(cos(x)) + atanh((cos(x) - sin(x))/sqrt(2))/sqrt(2)]

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

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(d*sin(e+f*x))^(m/2)*(a+b*tan(e+f*x))^n=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(d*sin(e+f*x))^(m/2)*(a+b*tan(e+f*x))^(n/2)=#

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

    @test_int [sin(c + d*x)^m*(a + b*tan(c + d*x))^3, x, 8, (a^3*cos(c + d*x)*HypergeometricFunctions._₂F₁(1/2, (1 + m)/2, (3 + m)/2, sin(c + d*x)^2)*sin(c + d*x)^(1 + m))/(d*(1 + m)*sqrt(cos(c + d*x)^2)) + (3*a^2*b*HypergeometricFunctions._₂F₁(1, (2 + m)/2, (4 + m)/2, sin(c + d*x)^2)*sin(c + d*x)^(2 + m))/(d*(2 + m)) + (3*a*b^2*sqrt(cos(c + d*x)^2)*HypergeometricFunctions._₂F₁(3/2, (3 + m)/2, (5 + m)/2, sin(c + d*x)^2)*sec(c + d*x)*sin(c + d*x)^(3 + m))/(d*(3 + m)) + (b^3*HypergeometricFunctions._₂F₁(2, (4 + m)/2, (6 + m)/2, sin(c + d*x)^2)*sin(c + d*x)^(4 + m))/(d*(4 + m))]
    @test_int [sin(c + d*x)^m*(a + b*tan(c + d*x))^2, x, 6, (a^2*cos(c + d*x)*HypergeometricFunctions._₂F₁(1/2, (1 + m)/2, (3 + m)/2, sin(c + d*x)^2)*sin(c + d*x)^(1 + m))/(d*(1 + m)*sqrt(cos(c + d*x)^2)) + (2*a*b*HypergeometricFunctions._₂F₁(1, (2 + m)/2, (4 + m)/2, sin(c + d*x)^2)*sin(c + d*x)^(2 + m))/(d*(2 + m)) + (b^2*sqrt(cos(c + d*x)^2)*HypergeometricFunctions._₂F₁(3/2, (3 + m)/2, (5 + m)/2, sin(c + d*x)^2)*sec(c + d*x)*sin(c + d*x)^(3 + m))/(d*(3 + m))]
    @test_int [sin(c + d*x)^m*(a + b*tan(c + d*x))^1, x, 5, (a*cos(c + d*x)*HypergeometricFunctions._₂F₁(1/2, (1 + m)/2, (3 + m)/2, sin(c + d*x)^2)*sin(c + d*x)^(1 + m))/(d*(1 + m)*sqrt(cos(c + d*x)^2)) + (b*HypergeometricFunctions._₂F₁(1, (2 + m)/2, (4 + m)/2, sin(c + d*x)^2)*sin(c + d*x)^(2 + m))/(d*(2 + m))]
    @test_int [sin(c + d*x)^m/(a + b*tan(c + d*x))^1, x, 14, (2^(1 + m)*HypergeometricFunctions._₂F₁((1 + m)/2, 1 + m, (3 + m)/2, -tan((1/2)*(c + d*x))^2)*tan((1/2)*(c + d*x))*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/(a*d*(1 + m)) + (2^(1 + m)*b*AppellF1((2 + m)/2, 1 + m, 1, (4 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b - sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))^2*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/(sqrt(a^2 + b^2)*(b - sqrt(a^2 + b^2))*d*(2 + m)) - (2^(1 + m)*b*AppellF1((2 + m)/2, 1 + m, 1, (4 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b + sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))^2*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/(sqrt(a^2 + b^2)*(b + sqrt(a^2 + b^2))*d*(2 + m)) + (2^(1 + m)*a*b*AppellF1((3 + m)/2, 1 + m, 1, (5 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b - sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))^3*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/(sqrt(a^2 + b^2)*(b - sqrt(a^2 + b^2))^2*d*(3 + m)) - (2^(1 + m)*a*b*AppellF1((3 + m)/2, 1 + m, 1, (5 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b + sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))^3*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/(sqrt(a^2 + b^2)*(b + sqrt(a^2 + b^2))^2*d*(3 + m))]
    #= [sin(c + d*x)^m/(a + b*tan(c + d*x))^2, x, 52, -((2^(1 + m)*b^2*AppellF1((1 + m)/2, 1 + m, 1, (3 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b - sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)^(3/2)*(b - sqrt(a^2 + b^2))*d*(1 + m))) - (2^(2 + m)*b*(1 - (2*b)/sqrt(a^2 + b^2))*AppellF1((1 + m)/2, 1 + m, 1, (3 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b - sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/(a^2*(b - sqrt(a^2 + b^2))*d*(1 + m)) + (2^(1 + m)*b^2*AppellF1((1 + m)/2, 1 + m, 1, (3 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b + sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)^(3/2)*(b + sqrt(a^2 + b^2))*d*(1 + m)) - (2^(2 + m)*b*(1 + (2*b)/sqrt(a^2 + b^2))*AppellF1((1 + m)/2, 1 + m, 1, (3 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b + sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/(a^2*(b + sqrt(a^2 + b^2))*d*(1 + m)) + (2^(1 + m)*b^2*(b - sqrt(a^2 + b^2))^2*AppellF1((1 + m)/2, 1 + m, 2, (3 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(a^2 + 2*b*(b - sqrt(a^2 + b^2))))*tan((1/2)*(c + d*x))*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)*(a^2 + 2*b*(b - sqrt(a^2 + b^2)))^2*d*(1 + m)) + (2^(1 + m)*b^2*(b + sqrt(a^2 + b^2))^2*AppellF1((1 + m)/2, 1 + m, 2, (3 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(a^2 + 2*b*(b + sqrt(a^2 + b^2))))*tan((1/2)*(c + d*x))*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)*(a^2 + 2*b*(b + sqrt(a^2 + b^2)))^2*d*(1 + m)) + (2^(1 + m)*HypergeometricFunctions._₂F₁((1 + m)/2, 1 + m, (3 + m)/2, -tan((1/2)*(c + d*x))^2)*tan((1/2)*(c + d*x))*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/(a^2*d*(1 + m)) - (2^(1 + m)*a*b^2*AppellF1((2 + m)/2, 1 + m, 1, (4 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b - sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))^2*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)^(3/2)*(b - sqrt(a^2 + b^2))^2*d*(2 + m)) - (2^(2 + m)*b*(1 - (2*b)/sqrt(a^2 + b^2))*AppellF1((2 + m)/2, 1 + m, 1, (4 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b - sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))^2*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/(a*(b - sqrt(a^2 + b^2))^2*d*(2 + m)) - (2^(2 + m)*b^3*AppellF1((2 + m)/2, 1 + m, 1, (4 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b - sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))^2*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/(a*(a^2 + b^2)^(3/2)*(b - sqrt(a^2 + b^2))*d*(2 + m)) + (2^(1 + m)*a*b^2*AppellF1((2 + m)/2, 1 + m, 1, (4 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b + sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))^2*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)^(3/2)*(b + sqrt(a^2 + b^2))^2*d*(2 + m)) - (2^(2 + m)*b*(1 + (2*b)/sqrt(a^2 + b^2))*AppellF1((2 + m)/2, 1 + m, 1, (4 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b + sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))^2*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/(a*(b + sqrt(a^2 + b^2))^2*d*(2 + m)) + (2^(2 + m)*b^3*AppellF1((2 + m)/2, 1 + m, 1, (4 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b + sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))^2*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/(a*(a^2 + b^2)^(3/2)*(b + sqrt(a^2 + b^2))*d*(2 + m)) + (2^(2 + m)*a*b^2*(b - sqrt(a^2 + b^2))*AppellF1((2 + m)/2, 1 + m, 2, (4 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(a^2 + 2*b*(b - sqrt(a^2 + b^2))))*tan((1/2)*(c + d*x))^2*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)*(a^2 + 2*b*(b - sqrt(a^2 + b^2)))^2*d*(2 + m)) + (2^(2 + m)*b^3*(b - sqrt(a^2 + b^2))^2*AppellF1((2 + m)/2, 1 + m, 2, (4 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(a^2 + 2*b*(b - sqrt(a^2 + b^2))))*tan((1/2)*(c + d*x))^2*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/(a*(a^2 + b^2)*(a^2 + 2*b*(b - sqrt(a^2 + b^2)))^2*d*(2 + m)) + (2^(2 + m)*a*b^2*(b + sqrt(a^2 + b^2))*AppellF1((2 + m)/2, 1 + m, 2, (4 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(a^2 + 2*b*(b + sqrt(a^2 + b^2))))*tan((1/2)*(c + d*x))^2*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)*(a^2 + 2*b*(b + sqrt(a^2 + b^2)))^2*d*(2 + m)) + (2^(2 + m)*b^3*(b + sqrt(a^2 + b^2))^2*AppellF1((2 + m)/2, 1 + m, 2, (4 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(a^2 + 2*b*(b + sqrt(a^2 + b^2))))*tan((1/2)*(c + d*x))^2*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/(a*(a^2 + b^2)*(a^2 + 2*b*(b + sqrt(a^2 + b^2)))^2*d*(2 + m)) - (2^(2 + m)*b^3*AppellF1((3 + m)/2, 1 + m, 1, (5 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b - sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))^3*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)^(3/2)*(b - sqrt(a^2 + b^2))^2*d*(3 + m)) + (2^(2 + m)*b^3*AppellF1((3 + m)/2, 1 + m, 1, (5 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(b + sqrt(a^2 + b^2))^2)*tan((1/2)*(c + d*x))^3*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)^(3/2)*(b + sqrt(a^2 + b^2))^2*d*(3 + m)) + (2^(1 + m)*a^2*b^2*AppellF1((3 + m)/2, 1 + m, 2, (5 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(a^2 + 2*b*(b - sqrt(a^2 + b^2))))*tan((1/2)*(c + d*x))^3*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)*(a^2 + 2*b*(b - sqrt(a^2 + b^2)))^2*d*(3 + m)) + (2^(3 + m)*b^3*(b - sqrt(a^2 + b^2))*AppellF1((3 + m)/2, 1 + m, 2, (5 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(a^2 + 2*b*(b - sqrt(a^2 + b^2))))*tan((1/2)*(c + d*x))^3*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)*(a^2 + 2*b*(b - sqrt(a^2 + b^2)))^2*d*(3 + m)) + (2^(1 + m)*a^2*b^2*AppellF1((3 + m)/2, 1 + m, 2, (5 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(a^2 + 2*b*(b + sqrt(a^2 + b^2))))*tan((1/2)*(c + d*x))^3*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)*(a^2 + 2*b*(b + sqrt(a^2 + b^2)))^2*d*(3 + m)) + (2^(3 + m)*b^3*(b + sqrt(a^2 + b^2))*AppellF1((3 + m)/2, 1 + m, 2, (5 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(a^2 + 2*b*(b + sqrt(a^2 + b^2))))*tan((1/2)*(c + d*x))^3*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)*(a^2 + 2*b*(b + sqrt(a^2 + b^2)))^2*d*(3 + m)) + (2^(2 + m)*a*b^3*AppellF1((4 + m)/2, 1 + m, 2, (6 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(a^2 + 2*b*(b - sqrt(a^2 + b^2))))*tan((1/2)*(c + d*x))^4*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)*(a^2 + 2*b*(b - sqrt(a^2 + b^2)))^2*d*(4 + m)) + (2^(2 + m)*a*b^3*AppellF1((4 + m)/2, 1 + m, 2, (6 + m)/2, -tan((1/2)*(c + d*x))^2, (a^2*tan((1/2)*(c + d*x))^2)/(a^2 + 2*b*(b + sqrt(a^2 + b^2))))*tan((1/2)*(c + d*x))^4*(tan((1/2)*(c + d*x))/(1 + tan((1/2)*(c + d*x))^2))^m*(1 + tan((1/2)*(c + d*x))^2)^m)/((a^2 + b^2)*(a^2 + 2*b*(b + sqrt(a^2 + b^2)))^2*d*(4 + m))] =#

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

    @test_int [sin(c + d*x)^m*(a + b*tan(c + d*x))^n, x, 0, CannotIntegrate(sin(c + d*x)^m*(a + b*tan(c + d*x))^n, x)]

    @test_int [sin(c + d*x)^4*(a + b*tan(c + d*x))^n, x, 7, -(((a*b^2*n*(5*a^2 + b^2*(3 + 2*n)) + sqrt(-b^2)*(3*a^4 + a^2*b^2*(6 + 6*n - n^2) + b^4*(3 + 4*n + n^2)))*HypergeometricFunctions._₂F₁(1, 1 + n, 2 + n, (a + b*tan(c + d*x))/(a - sqrt(-b^2)))*(a + b*tan(c + d*x))^(1 + n))/(16*b*(a^2 + b^2)^2*(a - sqrt(-b^2))*d*(1 + n))) - ((a*b^2*n*(5*a^2 + b^2*(3 + 2*n)) - sqrt(-b^2)*(3*a^4 + a^2*b^2*(6 + 6*n - n^2) + b^4*(3 + 4*n + n^2)))*HypergeometricFunctions._₂F₁(1, 1 + n, 2 + n, (a + b*tan(c + d*x))/(a + sqrt(-b^2)))*(a + b*tan(c + d*x))^(1 + n))/(16*b*(a^2 + b^2)^2*(a + sqrt(-b^2))*d*(1 + n)) + (cos(c + d*x)^4*(b + a*tan(c + d*x))*(a + b*tan(c + d*x))^(1 + n))/(4*(a^2 + b^2)*d) - (cos(c + d*x)^2*(a + b*tan(c + d*x))^(1 + n)*(b*(a^2*(7 - n) + b^2*(5 + n)) + a*(5*a^2 + b^2*(3 + 2*n))*tan(c + d*x)))/(8*(a^2 + b^2)^2*d)]
    @test_int [sin(c + d*x)^2*(a + b*tan(c + d*x))^n, x, 6, -(((a*b^2*n + sqrt(-b^2)*(a^2 + b^2*(1 + n)))*HypergeometricFunctions._₂F₁(1, 1 + n, 2 + n, (a + b*tan(c + d*x))/(a - sqrt(-b^2)))*(a + b*tan(c + d*x))^(1 + n))/(4*b*(a^2 + b^2)*(a - sqrt(-b^2))*d*(1 + n))) - ((a*b^2*n - sqrt(-b^2)*(a^2 + b^2*(1 + n)))*HypergeometricFunctions._₂F₁(1, 1 + n, 2 + n, (a + b*tan(c + d*x))/(a + sqrt(-b^2)))*(a + b*tan(c + d*x))^(1 + n))/(4*b*(a^2 + b^2)*(a + sqrt(-b^2))*d*(1 + n)) - (cos(c + d*x)^2*(b + a*tan(c + d*x))*(a + b*tan(c + d*x))^(1 + n))/(2*(a^2 + b^2)*d)]
    @test_int [csc(c + d*x)^2*(a + b*tan(c + d*x))^n, x, 2, (b*HypergeometricFunctions._₂F₁(2, 1 + n, 2 + n, 1 + (b*tan(c + d*x))/a)*(a + b*tan(c + d*x))^(1 + n))/(a^2*d*(1 + n))]
    @test_int [csc(c + d*x)^4*(a + b*tan(c + d*x))^n, x, 4, (b*(2 - n)*cot(c + d*x)^2*(a + b*tan(c + d*x))^(1 + n))/(6*a^2*d) - (cot(c + d*x)^3*(a + b*tan(c + d*x))^(1 + n))/(3*a*d) + (b*(6*a^2 + b^2*(2 - 3*n + n^2))*HypergeometricFunctions._₂F₁(2, 1 + n, 2 + n, 1 + (b*tan(c + d*x))/a)*(a + b*tan(c + d*x))^(1 + n))/(6*a^4*d*(1 + n))]

    @test_int [sin(c + d*x)^3*(a + b*tan(c + d*x))^n, x, 0, CannotIntegrate(sin(c + d*x)^3*(a + b*tan(c + d*x))^n, x)]
    @test_int [sin(c + d*x)^1*(a + b*tan(c + d*x))^n, x, 0, CannotIntegrate(sin(c + d*x)*(a + b*tan(c + d*x))^n, x)]
    @test_int [csc(c + d*x)^1*(a + b*tan(c + d*x))^n, x, 0, CannotIntegrate(csc(c + d*x)*(a + b*tan(c + d*x))^n, x)]
    @test_int [csc(c + d*x)^3*(a + b*tan(c + d*x))^n, x, 0, CannotIntegrate(csc(c + d*x)^3*(a + b*tan(c + d*x))^n, x)]

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

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

    #= ::Subsection:: =#
    #=Integrands*of*the*form*csc(e+f*x)^(m/2)*(a+b*tan(e+f*x))^(n/2)=#
    end

max_stars_repo_path: test/4 Trig functions/4.3 Tangent/4.3.1.3 (d sin)^m (a+b tan)^n.jl


text:
module Day22

export get_inputs, get_solution1, get_solution2

## Input getting
function get_inputs()
    test_input1 = test_input2 = read_input(joinpath(@__DIR__, "test_input1.txt"))
    test_output1 = 306
    test_output2 = 291
    data = read_input(joinpath(@__DIR__, "input.txt"))
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

function read_input(io)
    str = read(io, String)
    decks = split(str, r"(\n\n|\r\n\r\n)*Player [1-2]:", keepempty=false)
    return [parse.(Int, split(deck, r"(\n|\r\n)", keepempty=false)) for deck in decks]
end

## Solution functions
function combat(deck1, deck2, recursive=false)
    previous_rounds = Vector{Int}[]
    while true
        if recursive
            deck1 in previous_rounds && return (; deck=deck1, id=1)
            push!(previous_rounds, copy(deck1))
        end

        card1 = popfirst!(deck1)
        card2 = popfirst!(deck2)

        if recursive && card1<=length(deck1) && card2<=length(deck2)
            out = combat(deck1[1:card1], deck2[1:card2], recursive)
            out.id==1 ? push!(deck1, card1, card2) : push!(deck2, card2, card1)
        elseif card1 > card2
            push!(deck1, card1, card2)
        else
            push!(deck2, card2, card1)
        end

        isempty(deck1) && return (deck=deck2, id=2)
        isempty(deck2) && return (deck=deck1, id=1)
    end
end

get_score(winning_deck) = winning_deck' * reverse(only(axes(winning_deck)))

get_solution(decks, recursive=false) = get_score(combat(deepcopy(decks)..., recursive).deck)

# Part 1
get_solution1(data) = get_solution(data)

# Part 2
get_solution2(data) = get_solution(data, true)

end

ext: jl


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Descriptions
# ==============================================================================
#
#   Function calls to create the precompilation statements using SnoopCompiler.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

function precompilation_input()
    include("../test/runtests.jl")
end


text:
using ChaosTools, DynamicalSystemsBase
using Test
using ChaosTools.Distances: Cityblock, Euclidean

test_value = (val, vmin, vmax) -> @test vmin <= val <= vmax

println("\nTesting nonlinear timeseries analysis...")
@testset "Lyapunov from data" begin
    ds = Systems.henon()
    data = trajectory(ds, 100000)
    x = data[:, 1] # some "recorded" timeseries
    ks = 1:20
    @testset "meth = $meth" for meth in
        [NeighborNumber(1), NeighborNumber(4), WithinRange(0.01)]
        @testset "distance = $di" for di in [Euclidean(), Cityblock()]
            for D in [2, 4]
                R = embed(x, D, 1)
                E = lyapunov_from_data(R, ks,
                refstates = 1:1000, distance=di, ntype=meth)
                λ = linear_region(ks, E)[2]
                test_value(λ, 0.3, 0.5)
            end
        end
    end
end

@testset "Broomhead-King" begin
    using Random
    Random.seed!(1234);
    ds = Systems.gissinger(ones(3)) # standard initial condition

    data = trajectory(ds, 1000.0, dt = 0.05)
    x = data[1:end-1, 1]
    s = x .+ 0.01rand(length(x))

    Ux, Σx = broomhead_king(x, 40)
    Us, Σs = broomhead_king(s, 40)
    R = embed(x, 2, 30)

    for j in 10:40
        @test Σx[j] < Σs[j]
    end

    newcoords = Dataset(Us[:, 1], Us[:, 2])
    newcoordsclean = Dataset(Ux[:, 1], Ux[:, 2])

    Dnew = generalized_dim(newcoords)
    DR = generalized_dim(R)
    DC = generalized_dim(newcoordsclean)
    @test abs(Dnew - 0.1 - DR) < 0.2 # subtract 0.1 for "added dimensionality"
    @test abs(DC - DR) < 0.2
end


text:
function dirclustercoeffs_test()
    (A,xy,labels) = load_matrix_network_metadata("celegans")
    (cc, cccyc, ccmid, ccin, ccout, nf) = dirclustercoeffs(A, true, true)
    if length(cc) != 202
        error("dirclustercoeffs failed")
    end
    (maxval, maxind) = findmax(cc)
    if maxind != 113
        error("dirclustercoeffs failed")
    end

    return true
end


text:
"""
    Sparse(container, values, indices)

Create a sparse container to make gradient calculations efficient. `s::Sparse` represents
the value `a` as defined below:

    a = zero(s.container)
    for (idx, val) in zip(s.indices, s.values)
        a[idx] .+= val
    end

except when there are repeated indices in `idx`, the corresponding values get added rather
than being overwritten. See https://github.com/JuliaLang/julia/issues/31392.
"""
struct Sparse{T,N} <: AbstractArray{T,N}
    container
    values::Vector{Any}
    indices::Vector{Any}
end

Sparse(a::AbstractArray{T,N},v::Vector{Any},i::Vector{Any}) where {T,N} = Sparse{T,N}(a,v,i)

# To add a Sparse to an Array without allocating extra space, we need to use:
# a .+= s  OR  a .= a .+ s
# Both of which translate to:
# materialize!(a, broadcasted(+, a, s))

import Base: size, copyto!
using Base.Broadcast: Broadcasted

# This is used in broadcasted
size(b::Sparse,d...)=size(b.container,d...)

# This is used in materialize!
function copyto!(a::AbstractArray, bc::Broadcasted{S,A,F,X}) where
    {S, A, F <: Union{typeof(+),typeof(-)}, X <: Tuple{Any,Sparse}}
    (b,c) = bc.args
    if !(size(a) == size(b) == size(c.container))
        a .= bc.f.(b, full(c))
        return a
    end
    a === b || copyto!(a, b)
    F <: typeof(-) && (c = -c)
    addto!(a, c)
    return a
end

# These are used in Knet/src/update.jl:
import LinearAlgebra: axpy!, norm, lmul!
axpy!(a::Number, x::Sparse, y::AbstractArray) = addto!(y, a*x)
lmul!(a::Number, x::Sparse{T,N}) where {T,N} = Sparse{T,N}(x.container, Any[ a*v for v in x.values ], x.indices)

# This does not give the correct result when there are repeated indices, but should be good enough for gclip
norm(x::Sparse) = sqrt(sum(abs2, norm(v) for v in x.values))

# Convert to regular array -- use this as last resort
full(b::Sparse)=addto!(zeroslike(b.container), b)
zeroslike(a::AbstractArray{T}) where T = (isbitstype(T) ? zero(a) : Array{Any}(nothing,size(a)))
full(x)=x

# Arithmetic with numbers
import Base: *, +, -, /
import Base.Broadcast: broadcasted
*(s::Sparse, n::Number) = Sparse(s.container, Any[ v*n for v in s.values ], s.indices)
*(n::Number, s::Sparse) = Sparse(s.container, Any[ v*n for v in s.values ], s.indices)
/(s::Sparse, n::Number) = Sparse(s.container, Any[ v/n for v in s.values ], s.indices)
broadcasted(::typeof(*), s::Sparse, n::Number) = Sparse(s.container, Any[ v.*n for v in s.values ], s.indices)
broadcasted(::typeof(*), n::Number, s::Sparse) = Sparse(s.container, Any[ v.*n for v in s.values ], s.indices)
broadcasted(::typeof(/), s::Sparse, n::Number) = Sparse(s.container, Any[ v./n for v in s.values ], s.indices)

# Arithmetic with arrays (can use addto! which overwrites its first argument)
+(a::AbstractArray, s::Sparse) = addto!(copy(a), s)
+(s::Sparse, a::AbstractArray) = addto!(copy(a), s)
-(a::AbstractArray, s::Sparse) = addto!(copy(a), -s)
-(s::Sparse, a::AbstractArray) = addto!(-a, s)
-(s::Sparse) = -1*s

# Issue #114: we may need to add multiple gradients
function +(a::Sparse, b::Sparse)
    @assert matches(a.container, b.container) "$(summary.((a.container, b.container)))"
    Sparse(a.container, Any[ a.values; b.values ], Any[ a.indices; b.indices ])
end

# Do we need these?
# sum(b::Sparse)=sum(sum(v) for v in b.values)
# zero(b::Sparse)=Sparse(b.container,Any[],Any[])
# ones(b::Sparse)=ones(b.container)
# length(b::Sparse)=length(b.container)

# We do not create Sparse for these types any more:
# zeroslike(a::AbstractDict)=empty(a)
# zeroslike(a::Tuple)=ntuple(i->nothing, length(a))
# zeroslike(a::Sparse)=zeroslike(a.container)
# zeroslike(a::T) where {T<:Number} = T(0)   # This comes up if people use getindex on a single number


text:
#Testbed for synthetic spike recovery from example background

using nnlearn, BGHMM, HMMBase, Distributed, Distributions, Random, Serialization, Distributions

Random.seed!(786)

#CONSTANTS
no_obs=1000
obsl=140

hmm=HMM{Univariate,Float64}([0.4016518533961019, 0.2724399569450827, 0.3138638675018568, 0.012044322156962559], [3.016523036789942e-9 2.860631288328858e-6 0.2299906524188302 0.7700064839333549; 3.0102278323431375e-11 0.7494895424906354 0.23378615437778671 0.016724303101477486; 8.665894321573098e-17 0.2789950381410553 0.7141355461949104 0.006869415664033568; 0.006872526597796038 0.016052425322133648 0.017255179541768192 0.9598198685383041], [Categorical([0.1582723599684065, 0.031949729618356536, 0.653113286526948, 0.15666462388628763]), Categorical([0.4610499748798372, 0.2613013005680122, 0.15161801872560146, 0.12603070582654768]), Categorical([0.08601960130086236, 0.13869192872427524, 0.26945182329686523, 0.5058366466779973]), Categorical([0.3787366527613563, 0.11034604356978756, 0.1119586976058099, 0.3989586060630392])])

struc_sig=[.1 .7 .1 .1
           .1 .1 .1 .7
           .1 .7 .1 .1]
periodicity=8
struc_frac_obs=.35

tata_box=[.05 .05 .05 .85
          .85 .05 .05 .05
          .05 .05 .05 .85
          .85 .05 .05 .05
          .425 .075 .075 .425
          .85 .05 .05 .05
          .425 .075 .075 .425]
tata_frac_obs=.7
tata_recur_range=1:4

combined_ensemble = "/bench/PhD/NGS_binaries/nnlearn/combined_ensemble"

#JOB CONSTANTS
const position_size = 141
const ensemble_size = 100
const no_sources = 3
const source_min_bases = 3
const source_max_bases = 12
@assert source_min_bases < source_max_bases
const source_length_range= source_min_bases:source_max_bases
const mixing_prior = .3
@assert mixing_prior >= 0 && mixing_prior <= 1
const models_to_permute = ensemble_size * 3

job_sets=[
([
    ("PS", (no_sources)),
    ("PM", (no_sources)),
    ("PSFM", (no_sources)),
    ("PSFM", (no_sources*10, .8, 1.)),
    ("FM", ()),
    # ("DM", (no_sources)),
    # ("SM", (no_sources)),
    ("RD", (no_sources)),
    ("RI", (no_sources)),
    ("EM", (no_sources))
],[.025, .025, .025, .025, .775, .075, .025, .025]),
([
    ("PS", (no_sources)),
    ("PM", (no_sources)),
    ("PSFM", (no_sources)),
    ("PSFM", (no_sources, 0., 1.)),
    ("PSFM", (no_sources, 0.8, .0)),
    ("FM", ()),
    ("DM", (no_sources)),
    ("SM", (no_sources)),
    ("RD", (no_sources)),
    ("RI", (no_sources*10)),
    ("EM", (no_sources))
],[.05, .15, .15, .10, .05, .20, 0.0, 0.00, 0.20, 0.10, .0]),
]
job_limit=4

const prior_wt=1.1

#FUNCTIONS
function setup_obs(hmm, no_obs, obsl)
    obs=vec(UInt8.(rand(hmm,obsl)[2]))
    for o in 2:no_obs
        obs=hcat(obs, vec(UInt8.(rand(hmm,obsl)[2])))
    end
    obs=vcat(obs,zeros(UInt8,1,no_obs))
    return obs
end

function spike_irreg!(obs, source, frac_obs, recur)
    truth=Vector{UInt8}()
    for o in 1:size(obs,2)
        if rand()<frac_obs
            push!(truth, o)
            for r in rand(recur)
                rand()<.5 && (source=nnlearn.revcomp_pwm(source))
                pos=rand(1:size(obs,1)-1)
                pwm_ctr=1
                while pos<=size(obs,1)-1&& pwm_ctr<=size(source,1)
                    obs[pos,o]=rand(Categorical(source[pwm_ctr,:]))
                    pos+=1
                    pwm_ctr+=1
                end
            end
        end
    end
    return truth
end

function spike_struc!(obs, source, frac_obs, periodicity)
    truth=Vector{UInt8}()
    truthpos=Vector{Vector{UInt8}}()
    for o in 1:size(obs,2)
        if rand()<frac_obs
            push!(truth, o)
            rand()<.5 && (source=nnlearn.revcomp_pwm(source))
            pos=rand(1:periodicity)
            posvec=Vector{UInt8}()
            while pos<=size(obs,1)
                pos_ctr=pos
                pwm_ctr=1
                while pos_ctr<=size(obs,1)-1&&pwm_ctr<=size(source,1)
                    push!(posvec,pos_ctr)
                    base=rand(Categorical(source[pwm_ctr,:]))
                    obs[pos_ctr,o]=base
                    pos_ctr+=1
                    pwm_ctr+=1
                end
                pos+=periodicity+pwm_ctr
            end
            push!(truthpos,posvec)
        end
    end
    return truth,truthpos
end

function get_BGHMM_lhs(obs,hmm)
    lh_mat=zeros(size(obs,1)-1,size(obs,2))
    for o in 1:size(obs,2)
        obso=zeros(UInt8,size(obs,1),1)
        obso[1:end,1] = obs[:,o]
        lh_mat[:,o]=BGHMM.get_BGHMM_symbol_lh(Matrix(transpose(obso)),hmm)
    end
    return lh_mat
end

@info "Setting up synthetic observation set..."
obs=setup_obs(hmm, no_obs, obsl)
struc_truth,struc_postruth=spike_struc!(obs, struc_sig, struc_frac_obs, periodicity)
irreg_truth=spike_irreg!(obs, tata_box, tata_frac_obs, tata_recur_range)

@info "Calculating background likelihood matrix..."
bg_lhs=get_BGHMM_lhs(obs,hmm)

@info "Assembling source priors..."
#prior_array= [struc_sig, tata_box]
prior_array= Vector{Matrix{Float64}}()
source_priors = nnlearn.assemble_source_priors(no_sources, prior_array, prior_wt, source_length_range)

@info "Assembling ensemble..."
path=randstring()
isfile(string(path,'/',"ens")) ? (ens = deserialize(string(path,'/',"ens"))) :
    (ens = nnlearn.Bayes_IPM_ensemble(path, ensemble_size, source_priors, (falses(0,0), mixing_prior), bg_lhs, obs, source_length_range))

job_set_thresh=[-Inf,ens.naive_lh]
param_set=(job_sets,job_set_thresh,job_limit)

@info "Spawning worker pool..."
remote_machine = "18.222.227.156"
remote_pool=addprocs([(remote_machine, 1)], tunnel=true, topology=:master_worker)
librarians=[1]

@everywhere using nnlearn

@info "Converging ensemble..."
final_logZ = nnlearn.ns_converge!(ens, [param_set], models_to_permute, librarians, remote_pool, 0.01, backup=(true,5), model_display=no_sources)

rm(path,recursive=true)


text:
# NOTE: this mapping is valid for Level 3/Version 2, it *may* not be valid for
# other versions.  See
# https://github.com/sbmlteam/libsbml/blob/d4bc12abc4e72e451a0a0f2be4b0b6101ac94160/src/sbml/UnitKind.c#L46-L85
const UNITFUL_KIND_STRING = Dict(
    "ampere" => 1.0 * u"A", # UNIT_KIND_AMPERE
    "avogadro" => ustrip(u"mol^-1", Unitful.Na), # UNIT_KIND_AVOGADRO
    "becquerel" => 1.0 * u"Bq", # UNIT_KIND_BECQUEREL
    "candela" => 1.0 * u"cd", # UNIT_KIND_CANDELA
    "Celsius" => 1.0 * u"°C", # UNIT_KIND_CELSIUS
    "coulomb" => 1.0 * u"C", # UNIT_KIND_COULOMB
    "dimensionless" => 1, # UNIT_KIND_DIMENSIONLESS
    "farad" => 1.0 * u"F", # UNIT_KIND_FARAD
    "gram" => 1.0 * u"g", # UNIT_KIND_GRAM
    "gray" => 1.0 * u"Gy", # UNIT_KIND_GRAY
    "henry" => 1.0 * u"H", # UNIT_KIND_HENRY
    "hertz" => 1.0 * u"Hz", # UNIT_KIND_HERTZ
    "item" => 1, # UNIT_KIND_ITEM
    "joule" => 1.0 * u"J", # UNIT_KIND_JOULE
    "katal" => 1.0 * u"kat", # UNIT_KIND_KATAL
    "kelvin" => 1.0 * u"K", # UNIT_KIND_KELVIN
    "kilogram" => 1.0 * u"kg", # UNIT_KIND_KILOGRAM
    "liter" => 1.0 * u"L", # UNIT_KIND_LITER
    "litre" => 1.0 * u"L", # UNIT_KIND_LITRE
    "lumen" => 1.0 * u"lm", # UNIT_KIND_LUMEN
    "lux" => 1.0 * u"lx", # UNIT_KIND_LUX
    "meter" => 1.0 * u"m", # UNIT_KIND_METER
    "metre" => 1.0 * u"m", # UNIT_KIND_METRE
    "mole" => 1.0 * u"mol", # UNIT_KIND_MOLE
    "newton" => 1.0 * u"N", # UNIT_KIND_NEWTON
    "ohm" => 1.0 * u"Ω", # UNIT_KIND_OHM
    "pascal" => 1.0 * u"Pa", # UNIT_KIND_PASCAL
    "radian" => 1.0 * u"rad", # UNIT_KIND_RADIAN
    "second" => 1.0 * u"s", # UNIT_KIND_SECOND
    "siemens" => 1.0 * u"S", # UNIT_KIND_SIEMENS
    "sievert" => 1.0 * u"Sv", # UNIT_KIND_SIEVERT
    "steradian" => 1.0 * u"sr", # UNIT_KIND_STERADIAN
    "tesla" => 1.0 * u"T", # UNIT_KIND_TESLA
    "volt" => 1.0 * u"V", # UNIT_KIND_VOLT
    "watt" => 1.0 * u"W", # UNIT_KIND_WATT
    "weber" => 1.0 * u"W", # UNIT_KIND_WEBER
    "(Invalid UnitKind)" => 1, # UNIT_KIND_INVALID (let's treat is as a dimensionless quantity)
)

"""
    unitful(u::UnitPart)

Converts an [`UnitPart`](@ref) to a corresponding Unitful unit.

The conversion is done according to the formula from
[SBML L3v2 core manual release 2](http://sbml.org/Special/specifications/sbml-level-3/version-2/core/release-2/sbml-level-3-version-2-release-2-core.pdf)(section 4.4.2).
"""
unitful(u::UnitPart) =
    (u.multiplier * UNITFUL_KIND_STRING[u.kind] * exp10(u.scale))^u.exponent

"""
    unitful(units::Vector{UnitPart})

Converts a SBML unit (i.e., a vector of [`UnitPart`](@ref)s) to a corresponding
Unitful unit.
"""
unitful(u::Vector{UnitPart}) = prod(unitful.(u))

"""
    unitful(m::Model, val::Tuple{Float64,String})

Computes a properly unitful value from a value-unit pair stored in the model
`m`.
"""
unitful(m::Model, val::Tuple{Float64,String}) = unitful(m.units[val[2]]) * val[1]

"""
    unitful(m::Model, val::Tuple{Float64, String}, default_unit::Number)

Overload of [`unitful`](@ref) that uses the `default_unit` if the unit is not
found in the model.

# Example
```
julia> SBML.unitful(mdl, (10.0,"firkin"), 90 * u"lb")
990.0 lb
```
"""
unitful(m::Model, val::Tuple{Float64,String}, default_unit::Number) =
    mayfirst(maylift(unitful, get(m.units, val[2], nothing)), default_unit) * val[1]

"""
    unitful(m::Model, val::Tuple{Float64, String}, default_unit::String)

Overload of [`unitful`](@ref) that allows specification of the `default_unit` by
string ID.
"""
unitful(m::Model, val::Tuple{Float64,String}, default_unit::String) =
    unitful(m, val, unitful(m.units[default_unit]))


text:
import GR
using Triplot

include("common.jl")

function plot_pcolor(fname, n, f)
    x,y,t = cartmesh(n)
    z = f.(x,y)

    GR.beginprint(fname)
    GR.setcolormap(1)
    tripcolor(x,y,z,t)
    GR.endprint()
end

plot_pcolor(joinpath(@__DIR__, "diag_pcolor.png"), 5, (x,y) -> x + y)
plot_pcolor(joinpath(@__DIR__, "peaks_pcolor.png"), 100, peaks_unit)
plot_pcolor(joinpath(@__DIR__, "sine_pcolor.png"), 100, (x,y) -> sin(2*pi*x)*cos(2*pi*y))


text:
function tests_invert_bkwds()
    model = MetNets.toy_model()

    bkwd_idxs = MetNets.bkwds_bounded(model)
    rename = (rxn) -> string(rxn, "_REV")
    new_model = MetNets.invert_bkwds(model; rename)

    for rxn_i in bkwd_idxs
        @test all(new_model.S[:,rxn_i] .== -model.S[:,rxn_i])
        @test new_model.lb[rxn_i] == -model.ub[rxn_i]
        @test new_model.ub[rxn_i] == -model.lb[rxn_i]
        @test new_model.rxns[rxn_i] == rename(model.rxns[rxn_i])
        @test new_model.rxnNames[rxn_i] == model.rxnNames[rxn_i]
    end
end
tests_invert_bkwds()


text:
export Shape
export Line, DashLine, Rectangle

abstract type Shape end

struct Line <: Shape
    len::Int
    background::Char
end

Line(len::Int; background=' ') = Line(len, background)

function DisplayRow(line::Line)
    fill_char = Char(0x2500)

    return DisplayRow(
        line.len,
        line.background,
        fill(fill_char, line.len)
    )
end

struct DashLine <: Shape
    len::Int
    background::Char
end

DashLine(len::Int; background=' ') = DashLine(len, background)

function DisplayRow(line::DashLine)
    fill_char = '-'

    return DisplayRow(
        line.len,
        line.background,
        fill(fill_char, line.len)
    )
end

struct Rectangle <: Shape
    h::Int
    w::Int
    background::Char
    border::Bool
end

Rectangle(h::Int, w::Int; background=' ', border=true) = Rectangle(h, w, background, border)
function Rectangle(size::Tuple{Int, Int}; background=' ', border=true)
    return Rectangle(size[1], size[2], background, border)
end

function DisplayArray(rectangle::Rectangle)
    h, w, background = rectangle.h, rectangle.w, rectangle.background
    array = DisplayArray(h, w, background=background)
    if rectangle.border

        # top, bottom ─
        line = Line(w, background)
        array.content[1] = DisplayRow(line)
        array.content[end] = DisplayRow(line)

        # side │
        foreach(i->array[i, 1]=Char(0x2502), 2:(h-1))
        foreach(i->array[i, end]=Char(0x2502), 2:(h-1))

        # corner
        array[1, 1] = Char(0x250C) # ┌
        array[1, end] = Char(0x2510) # ┐
        array[end, 1] = Char(0x2514) # └
        array[end, end] = Char(0x2518) # ┘
    end

    return array
end


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

# This implements an efficient sampler from a discrete probability distribution
# which changes in time.
# See Wong and Easton 1980: http://dx.doi.org/10.1137/0209009

module DynamicSamplers

using Random
using ExtractMacro
using ..Common

import Base: getindex, setindex!, @propagate_inbounds, copy!, copy
import Random: rand

export DynamicSampler, refresh!, getel

mutable struct DynamicSampler
    v::Vec
    ps::Vec
    z::Float64
    N::Int
    levs::Int
    tinds::Vector{Int}
    tpos::Vector{Int}
    function DynamicSampler(N::Int)
        levs = ceil(Int, log2(N))
        N2 = 2^levs
        v2 = zeros(N2)
        ps = zeros(N2 - 1)
        tinds, tpos = buildtable(levs)
        return new(v2, ps, 0.0, N, levs, tinds, tpos)
    end
    function DynamicSampler(v)
        isempty(v) && throw(ArgumentError("input must be non-empty"))
        N = length(v)
        levs = ceil(Int, log2(N))
        N2 = 2^levs
        v2 = zeros(N2)
        for (i,x) in enumerate(v)
            v2[i] = x
        end
        n = findfirst(v2 .< 0)
        n ≡ nothing || throw(ArgumentError("negative element given: v[$n] = $(v2[n])"))
        ps = zeros(N2 - 1)
        tinds, tpos = buildtable(levs)

        dynsmp = new(v2, ps, 0.0, N, levs, tinds, tpos)
        refresh!(dynsmp)

        return dynsmp
    end
end

function buildtable(levs::Int)
    N2 = 2^levs

    tinds = Array{Int}(undef, N2+1)
    tpos = empty!(Array{Int}(undef, N2 * levs))

    j0 = 1
    for i = 1:N2
        tinds[i] = j0
        j1 = j0
        k = 0
        u = 1 << (levs-1)
        off = 1
        @inbounds for lev = 1:levs
            if (i-1) & u == 0
                push!(tpos, off + k)
                j1 += 1
                k *= 2
            else
                k = 2k + 1
            end
            u >>>= 1
            off *= 2
        end
        j0 = j1
    end
    tinds[N2+1] = j0
    return tinds, tpos
end

function refresh!(dynsmp::DynamicSampler)
    @extract dynsmp : v ps N levs tinds tpos

    dynsmp.z = sum(v)
    fill!(ps, 0.0)
    @inbounds for i = 1:N
        x = v[i]
        @simd for j = tinds[i]:(tinds[i+1]-1)
            k = tpos[j]
            ps[k] += v[i]
        end
    end
    return dynsmp
end

function copy!(dest::DynamicSampler, src::DynamicSampler)
    @extract src : v ps z N levs
    dest.N == N || throw(ArgumentError("dest N=$(dest.N), src N=$N"))
    dest.levs == levs || throw(ArgumentError("dest levs=$(dest.levs), src levs=$levs"))

    copy!(dest.v, v)
    copy!(dest.ps, ps)
    dest.z = src.z
    return dest
end

copy(dynsmp::DynamicSampler) = copy!(DynamicSampler(dynsmp.N), dynsmp)

# for checking purposes
function getel_naive(v::Vec, z::Float64, x::Float64)
    @assert 0 ≤ x < 1
    #@extract dynsmp : v z N
    N = length(v)
    x *= z
    i = 0
    c = 0.0
    @inbounds for i = 1:N
        c += v[i]
        c ≥ x && break
    end
    return i
end

getel_naive(dynsmp::DynamicSampler, x::Float64) = getel_naive(dynsmp.v, dynsmp.z, x)

function getel(dynsmp::DynamicSampler, x::Float64)
    #@assert 0 ≤ x < 1
    @extract dynsmp : v ps z N levs
    x *= z

    k = 0
    off = 1
    @inbounds for lev = 1:levs
        p = ps[off+k]
        k *= 2
        if x > p
            x -= p
            k += 1
        end
        off *= 2
    end
    #@assert k < N
    #@assert v[k+1] > 0
    return k+1
end

rand(rng::AbstractRNG, dynsmp::DynamicSampler) = getel(dynsmp, rand(rng))
rand(dynsmp::DynamicSampler) = rand(Random.GLOBAL_RNG, dynsmp)

@inline getindex(dynsmp::DynamicSampler, i::Int) = dynsmp.v[i]

@propagate_inbounds function setindex!(dynsmp::DynamicSampler, x::Float64, i::Int)
    @extract dynsmp : v ps N levs tinds tpos
    @boundscheck 1 ≤ i ≤ N || throw(BoundsError(dynsmp, i))

    @inbounds d = x - v[i]
    @inbounds v[i] = x
    dynsmp.z += d

    @inbounds @simd for j = tinds[i]:(tinds[i+1]-1)
        k = tpos[j]
        ps[k] += d
    end
    return dynsmp

    ## following is the version without using
    ## the tinds/tpos tables, for reference
    ## (may be useful to save memory)

    # k = 0
    # u = 1 << (levs-1)
    # i -= 1
    # off = 1
    # @inbounds for lev = 1:levs
    #     if i & u == 0
    #         ps[off+k] += d
    #         k *= 2
    #     else
    #         k = 2k + 1
    #     end
    #     u >>>= 1
    #     off *= 2
    # end
    # return dynsmp
end

end # module


text:
# Blatantly lifted from BinDeps.jl because for some reason the build
# script sometimes cannot find BinDeps.
# @windows_only unpack_cmd has been modified to use zip in addition to 7z.

function splittarpath(path)
    path,extension = splitext(path)
    base_filename,secondary_extension = splitext(path)
    if extension == ".tgz" || extension == ".tbz" || extension == ".zip" && !isempty(secondary_extension)
        base_filename *= secondary_extension
        secondary_extension = ""
    end
    (base_filename,extension,secondary_extension)
end

downloadcmd_candidates = @static if is_windows()
  (:powershell, :curl, :wget, :fetch)
else
  (:curl, :wget, :fetch)
end

downloadcmd = nothing
function download_cmd(url::AbstractString, filename::AbstractString)
    global downloadcmd
    whichcmd = @static if is_windows() "where" else "which" end
    if downloadcmd === nothing
        for checkcmd in downloadcmd_candidates
            try
                if success(`$whichcmd $checkcmd`)
                    downloadcmd = checkcmd
                    break
                end
            catch
                continue # don't bail if one of these fails
            end
        end
    end
    if downloadcmd == :wget
        return `wget -O $filename $url`
    elseif downloadcmd == :curl
        return `curl -o $filename -L $url`
    elseif downloadcmd == :fetch
        return `fetch -f $filename $url`
    elseif downloadcmd == :powershell
        return `powershell -Command "(new-object net.webclient).DownloadFile(\"$url\", \"$filename\")"`
    else
        error("No download agent available; install curl, wget, or fetch.")
    end
end

@static if is_unix()
    function unpack_cmd(file,directory,extension,secondary_extension)
        if (extension == ".gz" && secondary_extension == ".tar") || extension == ".tgz"
            return (`tar xzf $file --directory=$directory`)
        elseif (extension == ".bz2" && secondary_extension == ".tar") || extension == ".tbz"
            return (`tar xjf $file --directory=$directory`)
        elseif extension == ".xz" && secondary_extension == ".tar"
            return (`unxz -c $file `|>`tar xv --directory=$directory`)
        elseif extension == ".tar"
            return (`tar xf $file --directory=$directory`)
        elseif extension == ".zip"
            return (`unzip -x $file -d $directory`)
        elseif extension == ".gz"
            return (`mkdir $directory` |> `cp $file $directory` |> `gzip -d $directory/$file`)
        end
        error("I don't know how to unpack $file")
    end
end

@static if is_windows()
    has_7z  = nothing
    has_zip = nothing
    function unpack_cmd(file,directory,extension,secondary_extension)
        global has_7z
        global has_zip

        if has_7z === nothing
            has_7z  = success(`where 7z`)
            has_zip = success(`where unzip`)
        end

        if has_7z
            if((extension == ".gz" || extension == ".xz" || extension == ".bz2") && secondary_extension == ".tar") ||
                extension == ".tgz" || extension == ".tbz"
                return (`7z x $file -y -so`|>`7z x -si -y -ttar -o$directory`)
            elseif extension == ".zip" || extension == ".7z"
                return (`7z x $file -y -o$directory`)
            else
                error("I don't know how to unpack $file")
            end
        elseif has_zip
            if extension == ".zip"
                return (`unzip $file -o -d$directory`)
            else
                error("I don't know how to unpack $file")
            end
        elseif extension == ".zip"
            rm(directory,recursive=true)
            return (`powershell -Command "Add-Type -AssemblyName System.IO.Compression.FileSystem; [System.IO.Compression.ZipFile]::ExtractToDirectory(\"$file\",\"$directory\")"`)
        else
            error("I don't know how to unpack $file")
        end
    end
end


text:
function simulate!(sim::Simulation)

    T, p = size(sim.y)
    a = size(sim.m, 1)
    d = convert(Int, (a - 1)/2)
    y = zeros(T+d, p)
    F = ones(a)

    for t in 1:T
        if t > d
            # Lagged values
            F[2:end] = vec(sim.y[t-1:-1:t-d, :])
            # Simulate Ω
            sim.Ω[t, :, :] = rand(MatrixNormal(sim.m, sim.P, sim.S))
            # Simulate y
            sim.y[t, :] = rand(MvNormal(vec(F' * sim.Ω[t, :, :]), sim.S))
        else
            sim.Ω[t, :, :] = sim.m
        end
    end
    #return sim
end


text:
# Functions for finding the Fermi level and occupation numbers for bands

import Roots

"""
Find the Fermi level.
"""
function find_fermi_level(basis, energies)
    find_occupation(basis, energies).εF
end

"""
Find the occupation and Fermi level.
"""
function find_occupation(basis::PlaneWaveBasis{T}, energies;
                         temperature=basis.model.temperature) where {T}
    @assert basis.model.spin_polarization in (:none, :spinless)
    model = basis.model
    n_electrons = model.n_electrons
    smearing = model.smearing

    # Maximum occupation per state
    filled_occ = filled_occupation(model)

    # The goal is to find εF so that
    # n_i = filled_occ * f((εi-εF)/T)
    # sum_i n_i = n_electrons
    # If temperature is zero, (εi-εF)/T = ±∞.
    # The occupation function is required to give 1 and 0 respectively in these cases.
    compute_occupation(εF) = [filled_occ * Smearing.occupation.(smearing, (ε .- εF) ./ temperature) for ε in energies]
    compute_n_elec(εF) = sum(basis.kweights .* sum.(compute_occupation(εF)))

    # assume that we can get the required number of electrons by filling every state
    @assert filled_occ*sum(basis.kweights .* length.(energies)) ≥ n_electrons - sqrt(eps(T))

    # Get rough bounds to bracket εF
    min_ε = minimum(minimum.(energies)) - 1
    max_ε = maximum(maximum.(energies)) + 1
    if temperature != 0
        @assert compute_n_elec(min_ε) < n_electrons < compute_n_elec(max_ε)
    end

    # Compute εF by bisection
    if compute_n_elec(max_ε) ≈ n_electrons
        # This branch takes care of the case of insulators at zero
        # temperature with as many bands as electrons; there it is
        # possible that compute_n_elec(max_ε) ≈ n_electrons but
        # compute_n_elec(max_ε) < n_electrons, so that bisection does
        # not work
        εF = max_ε
    else
        # Just use bisection here; note that with MP smearing there might
        # be multiple possible Fermi levels. This could be sped up with more
        # advanced methods (e.g. false position), but more care has to be
        # taken with convergence criteria and the like
        εF = Roots.find_zero(εF -> compute_n_elec(εF) - n_electrons, (min_ε, max_ε),
                             Roots.Bisection(), atol=eps(T))
    end

    if !isapprox(compute_n_elec(εF), n_electrons)
        if temperature == 0
            error("Unable to find non-fractional occupations that have the " *
                  "correct number of electrons. You should add a temperature.")
        else
            error("This should not happen, debug me!")
        end
    end

    minocc = maximum(minimum.(compute_occupation(εF)))
    if temperature > 0 && minocc > 1e-5
        @warn "One kpoint has a high minimum occupation $minocc. You should probably increase the number of bands."
    end

    (occupation=compute_occupation(εF), εF=εF)
end

"""
Find Fermi level and occupation for the given parameters, assuming a band gap
and zero temperature. This function is for DEBUG purposes only, and the
finite-temperature version with 0 temperature should be preferred.
"""
function find_occupation_bandgap(basis, energies)
    @assert basis.model.spin_polarization in (:none, :spinless)
    n_bands = length(energies[1])
    @assert all(e -> length(e) == n_bands, energies)
    n_electrons = basis.model.n_electrons
    T = eltype(basis)
    @assert basis.model.temperature == 0

    filled_occ = filled_occupation(basis.model)
    n_fill = div(n_electrons, filled_occ)
    @assert filled_occ * n_fill == n_electrons
    @assert n_bands ≥ n_fill

    # We need to fill n_fill states with occupation filled_occ
    # Find HOMO and LUMO
    HOMO = -Inf # highest occupied energy state
    LUMO = Inf  # lowest unoccupied energy state
    occupation = similar(basis.kpoints, Vector{T})
    for ik in 1:length(occupation)
        occupation[ik] = zeros(T, n_bands)
        occupation[ik][1:n_fill] .= filled_occ
        HOMO = max(HOMO, energies[ik][n_fill])
        if n_fill < n_bands
            LUMO = min(LUMO, energies[ik][n_fill + 1])
        end
    end
    @assert sum(basis.kweights .* sum.(occupation)) ≈ n_electrons

    # Put Fermi level slightly above HOMO energy, to ensure that HOMO < εF
    εF = nextfloat(HOMO)
    if εF > LUMO
        @warn("`find_occupation_zero_temperature` assumes an insulator, but the " *
              "system seems metallic. Try specifying a temperature and a smearing function.",
              HOMO, LUMO)
    end

    (occupation=occupation, εF=εF)
end


text:
module FieldArraysTests

using Gridap.Fields
using Gridap.Arrays
using Gridap.TensorValues

using Test

# Testing the default interface for field arrays

function result(f,x)
  T = return_type(testitem(f),testitem(x))
  r = zeros(T,size(x)...,size(f)...)
  for j in CartesianIndices(f)
    for i in CartesianIndices(x)
      r[i,j] = evaluate(f[j],x[i])
    end
  end
  r
end

p = Point(1.0,2.0)
np = 4
x = fill(p,np)
z = fill(p,0)

v = VectorValue{2,Float64}[(1,1),(4,2),(3,5)]
f = MockField.(v)

fp = v
∇fp = fill(zero(TensorValue{2,2,Float64}),length(v))
∇∇fp = fill(zero(ThirdOrderTensorValue{2,2,2,Float64,6}),length(v))
test_field_array(f,p,fp)
test_field_array(f,p,fp,grad=∇fp)
test_field_array(f,p,fp,grad=∇fp,gradgrad=∇∇fp)

test_field_array(f,x,result(f,x))
test_field_array(f,x,result(f,x),grad=result(∇.(f),x))
test_field_array(f,x,result(f,x),grad=result(∇.(f),x),gradgrad=result(∇∇.(f),x))

test_field_array(f,z,result(f,z))
test_field_array(f,z,result(f,z),grad=result(∇.(f),z))
test_field_array(f,z,result(f,z),grad=result(∇.(f),z),gradgrad=result(∇∇.(f),z))

#using BenchmarkTools
#
#c = return_cache(f,p)
#@btime evaluate!($c,$f,$p)
#c = return_cache(f,x)
#@btime evaluate!($c,$f,$x)
#
#∇f = Broadcasting(∇)(f)
#c = return_cache(∇f,p)
#@btime evaluate!($c,$∇f,$p)
#c = return_cache(∇f,x)
#@btime evaluate!($c,$∇f,$x)
#
#∇∇f = Broadcasting(∇∇)(f)
#c = return_cache(∇∇f,p)
#@btime evaluate!($c,$∇∇f,$p)
#c = return_cache(∇∇f,x)
#@btime evaluate!($c,$∇∇f,$x)

# Integration

v = VectorValue{2,Float64}[(1,1),(4,2),(3,5)]
f = MockField.(v)
ϕfun(x) = 2*x
ϕ = GenericField(ϕfun)
w = ones(size(x))

i = integrate(f,x,w)
@test i == reshape(sum(evaluate(f,x).*w,dims=1),size(f))

i = integrate(f,x,w,∇(ϕ))
@test i == reshape(sum(evaluate(f,x).*w.*meas.(∇(ϕ)(x)),dims=1),size(f))

#using BenchmarkTools
#c = return_cache(integrate,f,x,w)
#@btime evaluate!($c,$integrate,$f,$x,$w)
#
#J = ∇(ϕ)
#c = return_cache(integrate,f,x,w,J)
#@btime evaluate!($c,$integrate,$f,$x,$w,$J)

# testing 0-length array of fields

v = VectorValue{2,Float64}[]
f = MockField.(v)

fp = v
∇fp = fill(zero(TensorValue{2,2,Float64}),length(v))
∇∇fp = fill(zero(ThirdOrderTensorValue{2,2,2,Float64,6}),length(v))
test_field_array(f,p,fp)
test_field_array(f,p,fp,grad=∇fp)
test_field_array(f,p,fp,grad=∇fp,gradgrad=∇∇fp)

test_field_array(f,x,result(f,x))
test_field_array(f,x,result(f,x),grad=result(∇.(f),x))
test_field_array(f,x,result(f,x),grad=result(∇.(f),x),gradgrad=result(∇∇.(f),x))

test_field_array(f,z,result(f,z))
test_field_array(f,z,result(f,z),grad=result(∇.(f),z))
test_field_array(f,z,result(f,z),grad=result(∇.(f),z),gradgrad=result(∇∇.(f),z))

# MockFieldArray (this mimics how polynomial bases will be implemented)

v = VectorValue{2,Float64}[(1,1),(4,2),(3,5)]
f = MockFieldArray(v)

fp = v
∇fp = fill(zero(TensorValue{2,2,Float64}),length(v))
∇∇fp = fill(zero(ThirdOrderTensorValue{2,2,2,Float64,6}),length(v))
test_field_array(f,p,fp)
test_field_array(f,p,fp,grad=∇fp)

test_field_array(f,x,repeat(transpose(v),np))
test_field_array(f,x,repeat(transpose(v),np),grad=zeros(TensorValue{2,2,Float64},np,length(v)))

#using BenchmarkTools
#
#c = return_cache(f,p)
#@btime evaluate!($c,$f,$p)
#c = return_cache(f,x)
#@btime evaluate!($c,$f,$x)
#
#∇f = Broadcasting(∇)(f)
#c = return_cache(∇f,p)
#@btime evaluate!($c,$∇f,$p)
#c = return_cache(∇f,x)
#@btime evaluate!($c,$∇f,$x)
#
#∇∇f = Broadcasting(∇∇)(f)
#c = return_cache(∇∇f,p)
#@btime evaluate!($c,$∇∇f,$p)
#c = return_cache(∇∇f,x)
#@btime evaluate!($c,$∇∇f,$x)

# Transpose

v = VectorValue{2,Float64}[(1,1),(4,2),(3,5)]
b = MockField.(v)

f = transpose(b)

fp = transpose(v)
∇fp = zeros(TensorValue{2,2,Float64,4},1,length(v))

test_field_array(f,p,fp)
test_field_array(f,p,fp,grad=∇fp)
test_field_array(f,x,result(f,x))
test_field_array(f,x,result(f,x),grad=result(∇.(f),x))
test_field_array(f,z,result(f,z))
test_field_array(f,z,result(f,z),grad=result(∇.(f),z))

#using BenchmarkTools
#
#@btime transpose($b)
#
#c = return_cache(f,p)
#@btime evaluate!($c,$f,$p)
#c = return_cache(f,x)
#@btime evaluate!($c,$f,$x)
#
#∇f = Broadcasting(∇)(f)
#c = return_cache(∇f,p)
#@btime evaluate!($c,$∇f,$p)
#c = return_cache(∇f,x)
#@btime evaluate!($c,$∇f,$x)

# Composition

afun(x) = 2*x
a = GenericField(afun)

nf = 4
bfun(x) = x[1] + 3
bi = GenericField(bfun)
b = fill(bi,nf)

f = Broadcasting(∘)(b,a)

fp = map(i->i(a(p)),b)
∇fp = map(i->∇(a)(p)⋅∇(i)(a(p)),b)

test_field_array(f,p,fp)
test_field_array(f,p,fp,grad=∇fp)
test_field_array(f,x,result(f,x))
test_field_array(f,x,result(f,x),grad=result(∇.(f),x))
test_field_array(f,z,result(f,z))
test_field_array(f,z,result(f,z),grad=result(∇.(f),z))

#using BenchmarkTools
#
##@btime Broadcasting(∘)($b,$a)
#
#c = return_cache(f,p)
#@btime evaluate!($c,$f,$p)
#c = return_cache(f,x)
#@btime evaluate!($c,$f,$x)
#
#∇f = Broadcasting(∇)(f)
#c = return_cache(∇f,p)
#@btime evaluate!($c,$∇f,$p)
#c = return_cache(∇f,x)
#@btime evaluate!($c,$∇f,$x)

# Broadcasting operations

afun(x) = x[1]+2
a = GenericField(afun)

v = VectorValue{2,Float64}[(1,1),(4,2),(3,5)]
b = MockField.(v)

f = Broadcasting(Operation(*))(a,b)

fp = Broadcasting(*)(a(p),v)
∇fp = Broadcasting(⊗)(∇(a)(p),v)

test_field_array(f,p,fp)
test_field_array(f,p,fp,grad=∇fp)
test_field_array(f,x,result(f,x))
test_field_array(f,x,result(f,x),grad=result(∇.(f),x))
test_field_array(f,z,result(f,z))
test_field_array(f,z,result(f,z),grad=result(∇.(f),z))

#using BenchmarkTools
#
#c = return_cache(f,p)
#@btime evaluate!($c,$f,$p)
#c = return_cache(f,x)
#@btime evaluate!($c,$f,$x)
#
#∇f = Broadcasting(∇)(f)
#c = return_cache(∇f,p)
#@btime evaluate!($c,$∇f,$p)
#c = return_cache(∇f,x)
#@btime evaluate!($c,$∇f,$x)

# MockFieldArray (this mimics how polynomial bases will be implemented)

b = MockFieldArray(v)

f1 = Broadcasting(Operation(*))(a,b)
test_field_array(f1,p,fp)
test_field_array(f1,p,fp,grad=∇fp)
test_field_array(f1,x,result(f,x))
test_field_array(f1,x,result(f,x),grad=result(∇.(f),x))
test_field_array(f1,z,result(f,z))
test_field_array(f1,z,result(f,z),grad=result(∇.(f),z))

# transpose of an array of fields

avals = [1.0,0.5,1.2]
a = MockField.(avals)

f = transpose(a)

fp = transpose(evaluate(a,p))
∇fp = transpose(evaluate(Broadcasting(∇)(a),p))
∇∇fp = transpose(evaluate(Broadcasting(∇∇)(a),p))

test_field_array(f,p,fp)
test_field_array(f,p,fp,grad=∇fp,gradgrad=∇∇fp)
test_field_array(f,x,result(f,x))
test_field_array(f,x,result(f,x),grad=result(∇.(f),x),gradgrad=result(∇∇.(f),x))
test_field_array(f,z,result(f,z))
test_field_array(f,z,result(f,z),grad=result(∇.(f),z),gradgrad=result(∇∇.(f),z))

test_map(evaluate(Broadcasting(∇),a),Broadcasting(∇),a)
test_map(evaluate(Broadcasting(∇∇),a),Broadcasting(∇∇),a)

#using BenchmarkTools
#
#c = return_cache(f,p)
#@btime evaluate!($c,$f,$p)
#c = return_cache(f,x)
#@btime evaluate!($c,$f,$x)
#
#∇f = Broadcasting(∇)(f)
#c = return_cache(∇f,p)
#@btime evaluate!($c,$∇f,$p)
#c = return_cache(∇f,x)
#@btime evaluate!($c,$∇f,$x)

# transpose(i_to_f)*i_to_vals

avals = [1.0,0.5,1.2,3.9]
a = MockField.(avals)
b = VectorValue{2,Float64}[(1,1),(4,2),(3,5),(1,2)]

f = linear_combination(b,a)
@test isa(f,Fields.LinearCombinationField)

fp = transpose(b)*avals
∇fp = zero(TensorValue{2,2,Float64,4})

test_field(f,p,fp)
test_field(f,p,fp,grad=∇fp)
test_field(f,x,f.(x))
test_field(f,x,f.(x),grad=∇(f).(x))
test_field(f,z,f.(z))
test_field(f,z,f.(z),grad=∇(f).(z))

#using BenchmarkTools
#
#@btime linear_combination($b,$a)
#c = return_cache(f,p)
#@btime evaluate!($c,$f,$p)
#c = return_cache(f,x)
#@btime evaluate!($c,$f,$x)
#
#∇f = Broadcasting(∇)(f)
#c = return_cache(∇f,p)
#@btime evaluate!($c,$∇f,$p)
#c = return_cache(∇f,x)
#@btime evaluate!($c,$∇f,$x)

# ij_to_vals*j_to_f

avals =rand(4)
a = MockField.(avals)
b = zeros(VectorValue{2,Float64},4,3)

f = linear_combination(b,a)
#@show typeof(f)
@test isa(f,Fields.LinearCombinationFieldVector)

fp = transpose(b)*avals
∇fp = zeros(TensorValue{2,2,Float64,4},3)

test_field_array(f,p,fp)
test_field_array(f,p,fp,grad=∇fp)
test_field_array(f,x,result(f,x))
test_field_array(f,x,result(f,x),grad=result(∇.(f),x))
test_field_array(f,z,result(f,z))
test_field_array(f,z,result(f,z),grad=result(∇.(f),z))

#using BenchmarkTools
#@btime linear_combination($b,$a)
#
#c = return_cache(f,p)
#@btime evaluate!($c,$f,$p)
#c = return_cache(f,x)
#@btime evaluate!($c,$f,$x)
#
#∇f = Broadcasting(∇)(f)
#c = return_cache(∇f,p)
#@btime evaluate!($c,$∇f,$p)
#c = return_cache(∇f,x)
#@btime evaluate!($c,$∇f,$x)

## Test MockField
#
#np = 4
#p = Point(1.0,2.0)
#x = fill(p,np)
#npp = (np,np)
#
## v = 3.0
#d = 2
#v = VectorValue{d}(1.0,1.0)
#f = MockField{d}(v)
#fx = fill(v,np)
#
#nf = 3
#nff = (nf,nf)
#
#fa, p = test_field_array(f,p,nf,grad=true)
#
#c = return_cache(fa,p)
## @btime evaluate!($c,$fa,$p)
#c = return_cache(f,p)
## @btime evaluate!($c,$f,$p)
#
#bfa, x = test_broadcast_field_array(f,p,nf,np,grad=true)
#
#c = return_cache(bfa,x)
## @btime evaluate!($c,$bfa,$x);
#
#∇f = gradient(f)
#
#∇fa, p = test_field_array(∇f,p,nf)
#bfa, x = test_broadcast_field_array(∇f,p,nf,np)
#
#@test gradient.(fa) == ∇fa
#
#c = return_cache(bfa,x)
## @btime evaluate!($c,$bfa,$x)
#
#bf = f
#
#fa, x = test_field_array(bf,p,nf,grad=true)
#bfa, x = test_broadcast_field_array(bf,p,nf,np,grad=true)
#
#c = return_cache(bfa,x)
## @btime evaluate!($c,$bfa,$x)
#c = return_cache(fa,p)
## @btime evaluate!($c,$fa,$p)
#
##
#
#q(x) = 2*x
#∇q = gradient(q)
#
#qf = GenericField(q)
#
#af, x = test_field_array(qf,p,nf,grad=true)
#bfa, x = test_broadcast_field_array(qf,p,nf,np,grad=true)
#
#c = return_cache(bfa,x)
## @btime evaluate!($c,$bfa,$x)
#c = return_cache(af,p)
## @btime evaluate!($c,$af,$p)
#
#bqf = qf
#
#af, x = test_field_array(bqf,p,nf,grad=true)
#bfa, x = test_broadcast_field_array(bqf,p,nf,np,grad=true)
#
#c = return_cache(bfa,x)
## @btime evaluate!($c,$bfa,$x)
#c = return_cache(af,p)
## @btime evaluate!($c,$af,$p)
#
#
## GenericField (constant)
#
#v = 1.0
## v = VectorValue(4.0,3.0)
#vf = GenericField(v)
#
#af, x = test_field_array(vf,p,nf,hessian=true)
#bfa, x = test_broadcast_field_array(vf,p,nf,np,hessian=true)
#
#c = return_cache(bfa,x)
## @btime evaluate!($c,$bfa,$x)
#c = return_cache(af,p)
## @btime evaluate!($c,$af,$p)
#c = return_cache(vf,p)
## @btime evaluate!($c,$vf,$p)
#
#v = VectorValue(4.0,3.0)
#vf = GenericField(v)
#
#af, x = test_field_array(vf,p,nf,grad=true)
#bfa, x = test_broadcast_field_array(vf,p,nf,np,grad=true)
#
#c = return_cache(bfa,x)
## @btime evaluate!($c,$bfa,$x)
#c = return_cache(af,p)
## @btime evaluate!($c,$af,$p)
#
## ZeroField
#
#zvf = ZeroField(vf)
#
#af, x = test_field_array(zvf,p,nf,grad=true)
#bfa, x = test_broadcast_field_array(zvf,p,nf,np,grad=true)
#
#c = return_cache(bfa,x)
## @btime evaluate!($c,$bfa,$x)
#c = return_cache(af,p)
## @btime evaluate!($c,$af,$p)
#
#bzvf = zvf
#
#af, x = test_field_array(bzvf,p,nf,grad=true)
#bfa, x = test_broadcast_field_array(bzvf,p,nf,np,grad=true)
#
#c = return_cache(bfa,x)
## @btime evaluate!($c,$bfa,$x)
#c = return_cache(af,p)
## @btime evaluate!($c,$af,$p)
#
## Type unstable array
## @santiagobadia : Anything to do here?
## We would need to define promotion rules, convert, etc.
## It seems not feasible but I think this case is not of
## practical interest
## On the other hand, we can always optimize these methods
## for our particular problem at hand
#
#bfa = [zvf, vf, qf, zvf, vf, qf, zvf, vf, qf, zvf, vf, qf]
#
#c = return_cache(bfa,p)
## @btime evaluate!($c,$bfa,$p)
#
#
#q(x) = 2*x
#h(x) = 4*x
#qf = GenericField(q)
#hf = GenericField(h)
#
#bqf = qf
#bhf = hf
#
#af = [qf, hf, hf, qf]
#
#c = return_cache(af,p)
## @btime evaluate!($c,$af,$p)
#
#bfa = [bqf, bhf, bhf, bqf]
#c = return_cache(bfa,x)
## @btime evaluate!($c,$bfa,$x)

end # module


text:
using ControlSystems;
using DifferentialEquations;
using PyPlot;
using Polynomials;

M=2;m=0.1;l=0.5;g=9.81;
A=[0 1 0 0;(m+M)/M/l*g 0 0 0;0 0 0 1;-m*g/M 0 0 0];
B=[0;-1/M/l;0;1/M];
C=[1.0 0 0 0];
D=0.0;
sys=ss(A,B,C,D,0);

eigvals(A)

rank(ctrb(sys))

rank(obsv(sys))

u(x,t)=[0.0 0.0 0.0 0.0]*x;
x0=[0.1;0.0;0.0;0.0];
tspan=0:0.01:10;
y, t, x, uout = lsim(sys,u,tspan,x0=x0);

plot(t,x[:,1],label="theta")
plot(t,x[:,2],label="theta_dot")
plot(t,x[:,3],label="pos")
plot(t,x[:,4],label="vel")
xlabel("time (t)")
legend()
grid(true)

mutable struct Param
    A::Array{Float64,2}
    B::Array{Float64,1}
    C::Array{Float64,2}
    D::Float64
    K::Array{Float64,2}
    L::Array{Float64,2}
end

K=[0.0 0.0 0.0 0.0];
L=[0.0 0.0;0.0 0.0;0.0 0.0;0.0 0.0];
param= Param(A,B,C,D,K,L);

function f!(dx,x,param,t)
    A=param.A;
    B=param.B;
    u=0.0;
    temp=A*x+B*u;
#     make sure that you assign to each element
    dx[1]=temp[1];
    dx[2]=temp[2];
    dx[3]=temp[3];
    dx[4]=temp[4];

end;

prob=ODEProblem(f!,x0,(0.0,10.0),param);
sol=solve(prob,Tsit5(),dtmax=0.01);

n=length(sol.t);
theta=zeros(n); [theta[i]=sol.u[i][1] for i=1:n];
theta_dot=zeros(n); [theta_dot[i]=sol.u[i][2] for i=1:n];
pos=zeros(n); [pos[i]=sol.u[i][3] for i=1:n];
vel=zeros(n); [vel[i]=sol.u[i][4] for i=1:n];

plot(sol.t,theta,label="theta")
plot(sol.t,theta_dot,label="theta_dot")
plot(sol.t,pos,label="pos")
plot(sol.t,vel,label="vel")
xlabel("time (t)")
legend()
grid(true)

v1=nullspace([-2*eye(4)-sys.A sys.B]);
v2=nullspace([-3*eye(4)-sys.A sys.B]);
v3=nullspace([-4*eye(4)-sys.A sys.B]);
v4=nullspace([-5*eye(4)-sys.A sys.B]);

K=[v1[5,1] v2[5,1] v3[5,1] v4[5,1]]*inv([v1[1:4,1] v2[1:4,1] v3[1:4,1] v4[1:4,1]])

sys_new=sys.A-sys.B*K

eigvals(sys_new)

K1=place(sys,[-2,-3,-4,-5])

u(x,t)=-K*x;
y, t, x, uout = lsim(sys,u,tspan,x0=x0);

plot(t,x[:,1],label="theta")
plot(t,x[:,2],label="theta_dot")
plot(t,x[:,3],label="pos")
plot(t,x[:,4],label="vel")
xlabel("time (t)")
legend()
grid(true)

C=[1.0 0 0 0;0.0 0 1 0]
rank(obsv(sys.A,C))

function f!(dx,x,param,t)
    A=param.A;
    B=param.B;
    K=param.K;
    L=param.L;

    xtrue=x[1:4];
    xhat=x[5:8];
    u=-K*xhat;

    # observer
    temp_hat=A*xhat+B*u[1]+L*C*(xhat-xtrue);

    dx[5]=temp_hat[1];
    dx[6]=temp_hat[2];
    dx[7]=temp_hat[3];
    dx[8]=temp_hat[4];

    # true physical system
    temp_true=A*xtrue+B*u[1];

    dx[1]=temp_true[1];
    dx[2]=temp_true[2];
    dx[3]=temp_true[3];
    dx[4]=temp_true[4];

end;

v1=nullspace([-2*eye(4)-A' C'])
v2=nullspace([-3*eye(4)-A' C'])
v3=nullspace([-4*eye(4)-A' C'])
v4=nullspace([-5*eye(4)-A' C'])

L=[v1[5:6,1] v2[5:6,1] v3[5:6,1] v4[5:6,1]]*inv( [v1[1:4,1] v2[1:4,1] v3[1:4,1] v4[1:4,1]]);
L=L';
eigvals(A-L*C)

# L=place(A',C',[-10,-11,-12,-13])
# L=[22.9338 1.0388;151.3808 11.9770;0.9570 23.0662;10.4925 132.2187];
param.C=C;
param.K=K;
param.L=-L;

x0=[0.1,0.0,0.0,0.0,0.0,0.0,0.0,0.0];
prob=ODEProblem(f!,x0,(0.0,10.0),param);
sol=solve(prob,Tsit5(),dtmax=0.01);

n=length(sol.t);
theta=zeros(n); [theta[i]=sol.u[i][1] for i=1:n];
theta_dot=zeros(n); [theta_dot[i]=sol.u[i][2] for i=1:n];
pos=zeros(n); [pos[i]=sol.u[i][3] for i=1:n];
vel=zeros(n); [vel[i]=sol.u[i][4] for i=1:n];
theta_hat=zeros(n); [theta_hat[i]=sol.u[i][5] for i=1:n];
theta_dot_hat=zeros(n); [theta_dot_hat[i]=sol.u[i][6] for i=1:n];
pos_hat=zeros(n); [pos_hat[i]=sol.u[i][7] for i=1:n];
vel_hat=zeros(n); [vel_hat[i]=sol.u[i][8] for i=1:n];

figure(figsize=(10,10))
subplot(221)
plot(sol.t,theta,label="theta")
plot(sol.t,theta_hat,label="theta_hat")
xlabel("time (t)")
legend()
grid(true)
subplot(222)
plot(sol.t,theta_dot,label="theta_dot")
plot(sol.t,theta_dot_hat,label="theta_dot_hat")
xlabel("time (t)")
legend()
grid(true)
subplot(223)
plot(sol.t,pos,label="pos")
plot(sol.t,pos_hat,label="pos_hat")
xlabel("time (t)")
legend()
grid(true)
subplot(224)
plot(sol.t,vel,label="vel")
plot(sol.t,vel_hat,label="vel_hat")
xlabel("time (t)")
legend()
grid(true)

# place(A',C',[-10,-15,-20,-25])
L=[7.7024 -0.7207;34.5649 -2.5235;-0.7029 6.2976;-2.9390 9.0361];
param.L=-L;

prob=ODEProblem(f!,x0,(0.0,10.0),param);
sol=solve(prob,Tsit5(),dtmax=0.01);

n=length(sol.t);
theta=zeros(n); [theta[i]=sol.u[i][1] for i=1:n];
theta_dot=zeros(n); [theta_dot[i]=sol.u[i][2] for i=1:n];
pos=zeros(n); [pos[i]=sol.u[i][3] for i=1:n];
vel=zeros(n); [vel[i]=sol.u[i][4] for i=1:n];
theta_hat=zeros(n); [theta_hat[i]=sol.u[i][5] for i=1:n];
theta_dot_hat=zeros(n); [theta_dot_hat[i]=sol.u[i][6] for i=1:n];
pos_hat=zeros(n); [pos_hat[i]=sol.u[i][7] for i=1:n];
vel_hat=zeros(n); [vel_hat[i]=sol.u[i][8] for i=1:n];

figure(figsize=(10,10))
subplot(221)
plot(sol.t,theta,label="theta")
plot(sol.t,theta_hat,label="theta_hat")
xlabel("time (t)")
legend()
grid(true)
subplot(222)
plot(sol.t,theta_dot,label="theta_dot")
plot(sol.t,theta_dot_hat,label="theta_dot_hat")
xlabel("time (t)")
legend()
grid(true)
subplot(223)
plot(sol.t,pos,label="pos")
plot(sol.t,pos_hat,label="pos_hat")
xlabel("time (t)")
legend()
grid(true)
subplot(224)
plot(sol.t,vel,label="vel")
plot(sol.t,vel_hat,label="vel_hat")
xlabel("time (t)")
legend()
grid(true)

K=lqr(sys,diagm([1,1,100,100]),1)

param.K=K;

prob=ODEProblem(f!,x0,(0.0,10.0),param);
sol=solve(prob,Tsit5(),dtmax=0.01);

n=length(sol.t);
theta=zeros(n); [theta[i]=sol.u[i][1] for i=1:n];
theta_dot=zeros(n); [theta_dot[i]=sol.u[i][2] for i=1:n];
pos=zeros(n); [pos[i]=sol.u[i][3] for i=1:n];
vel=zeros(n); [vel[i]=sol.u[i][4] for i=1:n];
theta_hat=zeros(n); [theta_hat[i]=sol.u[i][5] for i=1:n];
theta_dot_hat=zeros(n); [theta_dot_hat[i]=sol.u[i][6] for i=1:n];
pos_hat=zeros(n); [pos_hat[i]=sol.u[i][7] for i=1:n];
vel_hat=zeros(n); [vel_hat[i]=sol.u[i][8] for i=1:n];

figure(figsize=(10,10))
subplot(221)
plot(sol.t,theta,label="theta")
plot(sol.t,theta_hat,label="theta_hat")
xlabel("time (t)")
legend()
grid(true)
subplot(222)
plot(sol.t,theta_dot,label="theta_dot")
plot(sol.t,theta_dot_hat,label="theta_dot_hat")
xlabel("time (t)")
legend()
grid(true)
subplot(223)
plot(sol.t,pos,label="pos")
plot(sol.t,pos_hat,label="pos_hat")
xlabel("time (t)")
legend()
grid(true)
subplot(224)
plot(sol.t,vel,label="vel")
plot(sol.t,vel_hat,label="vel_hat")
xlabel("time (t)")
legend()
grid(true)

max_stars_repo_name: akshayhiregoudar/Control-System-Design

max_issues_repo_name: akshayhiregoudar/Control-System-Design

max_forks_repo_name: akshayhiregoudar/Control-System-Design