text:
@testset "4.2.1.2 (g sin)^p (a+b cos)^m" begin
    (a, b, c, d, e, x, ) = @variables a b c d e x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(g*sin(e+f*x))^p*(a+b*cos(e+f*x))^m=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(g*sin(e+f*x))^p*(a+b*cos(e+f*x))^m*when*a^2-b^2=0=#

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

    @test_int [sin(x)^4/(a + a*cos(x)), x, 3, x/(2*a) - (cos(x)*sin(x))/(2*a) - sin(x)^3/(3*a)]
    @test_int [sin(x)^3/(a + a*cos(x)), x, 2, -(cos(x)/a) + cos(x)^2/(2*a)]
    @test_int [sin(x)^2/(a + a*cos(x)), x, 2, x/a - sin(x)/a]
    @test_int [sin(x)^1/(a + a*cos(x)), x, 2, -(log(1 + cos(x))/a)]
    @test_int [sin(x)^0/(a + a*cos(x)), x, 1, sin(x)/(a + a*cos(x))]
    @test_int [csc(x)^1/(a + a*cos(x)), x, 4, -(atanh(cos(x))/(2*a)) + 1/(2*(a + a*cos(x)))]
    @test_int [csc(x)^2/(a + a*cos(x)), x, 3, -((2*cot(x))/(3*a)) + csc(x)/(3*(a + a*cos(x)))]
    @test_int [csc(x)^3/(a + a*cos(x)), x, 4, -((3*atanh(cos(x)))/(8*a)) - 1/(8*(a - a*cos(x))) + a/(8*(a + a*cos(x))^2) + 1/(4*(a + a*cos(x)))]
    @test_int [csc(x)^4/(a + a*cos(x)), x, 3, -((4*cot(x))/(5*a)) - (4*cot(x)^3)/(15*a) + csc(x)^3/(5*(a + a*cos(x)))]

    @test_int [sin(2*x)/(1 + cos(2*x)), x, 2, -log(cos(x)), (-(1/2))*log(1 + cos(2*x))]
    @test_int [sin(2*x)/(1 - cos(2*x)), x, 2, log(sin(x)), (1/2)*log(1 - cos(2*x))]

    @test_int [sin(x)/(1 + cos(x))^2, x, 2, 1/(1 + cos(x))]
    @test_int [sin(x)/(1 - cos(x))^2, x, 2, -(1/(1 - cos(x)))]
    @test_int [sin(x)^2/(1 + cos(x))^2, x, 2, -x + (2*sin(x))/(1 + cos(x))]
    @test_int [sin(x)^2/(1 - cos(x))^2, x, 2, -x - (2*sin(x))/(1 - cos(x))]
    @test_int [sin(x)^3/(1 + cos(x))^2, x, 3, cos(x) - 2*log(1 + cos(x))]
    @test_int [sin(x)^3/(1 - cos(x))^2, x, 3, cos(x) + 2*log(1 - cos(x))]

    @test_int [sin(x)/(1 + cos(x))^3, x, 2, 1/(2*(1 + cos(x))^2)]
    @test_int [sin(x)/(1 - cos(x))^3, x, 2, -(1/(2*(1 - cos(x))^2))]
    @test_int [sin(x)^2/(1 + cos(x))^3, x, 1, sin(x)^3/(3*(1 + cos(x))^3)]
    @test_int [sin(x)^2/(1 - cos(x))^3, x, 1, -(sin(x)^3/(3*(1 - cos(x))^3))]
    @test_int [sin(x)^3/(1 + cos(x))^3, x, 3, 2/(1 + cos(x)) + log(1 + cos(x))]
    @test_int [sin(x)^3/(1 - cos(x))^3, x, 3, -(2/(1 - cos(x))) - log(1 - cos(x))]

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

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(g*sin(e+f*x))^(p/2)*(a+a*cos(e+f*x))^m=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(g*sin(e+f*x))^(p/2)*(a+a*cos(e+f*x))^(m/2)=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(g*sin(e+f*x))^p*(a+b*cos(e+f*x))^m=#

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

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

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

    @test_int [sin(x)^4/(a + b*cos(x)), x, 5, -((a*(2*a^2 - 3*b^2)*x)/(2*b^4)) + (2*(a - b)^(3/2)*(a + b)^(3/2)*atan((sqrt(a - b)*tan(x/2))/sqrt(a + b)))/b^4 + ((2*(a^2 - b^2) - a*b*cos(x))*sin(x))/(2*b^3) - sin(x)^3/(3*b)]
    @test_int [sin(x)^3/(a + b*cos(x)), x, 3, -((a*cos(x))/b^2) + cos(x)^2/(2*b) + ((a^2 - b^2)*log(a + b*cos(x)))/b^3]
    @test_int [sin(x)^2/(a + b*cos(x)), x, 4, (a*x)/b^2 - (2*sqrt(a - b)*sqrt(a + b)*atan((sqrt(a - b)*tan(x/2))/sqrt(a + b)))/b^2 - sin(x)/b]
    @test_int [sin(x)^1/(a + b*cos(x)), x, 2, -(log(a + b*cos(x))/b)]
    @test_int [sin(x)^0/(a + b*cos(x)), x, 2, (2*atan((sqrt(a - b)*tan(x/2))/sqrt(a + b)))/(sqrt(a - b)*sqrt(a + b))]
    @test_int [csc(x)^1/(a + b*cos(x)), x, 6, log(1 - cos(x))/(2*(a + b)) - log(1 + cos(x))/(2*(a - b)) + (b*log(a + b*cos(x)))/(a^2 - b^2)]
    @test_int [csc(x)^2/(a + b*cos(x)), x, 4, -((2*b^2*atan((sqrt(a - b)*tan(x/2))/sqrt(a + b)))/((a - b)^(3/2)*(a + b)^(3/2))) + ((b - a*cos(x))*csc(x))/(a^2 - b^2)]
    @test_int [csc(x)^3/(a + b*cos(x)), x, 4, ((b - a*cos(x))*csc(x)^2)/(2*(a^2 - b^2)) + ((a + 2*b)*log(1 - cos(x)))/(4*(a + b)^2) - ((a - 2*b)*log(1 + cos(x)))/(4*(a - b)^2) - (b^3*log(a + b*cos(x)))/(a^2 - b^2)^2]
    @test_int [csc(x)^4/(a + b*cos(x)), x, 5, (2*b^4*atan((sqrt(a - b)*tan(x/2))/sqrt(a + b)))/((a - b)^(5/2)*(a + b)^(5/2)) - ((3*b^3 + a*(2*a^2 - 5*b^2)*cos(x))*csc(x))/(3*(a^2 - b^2)^2) + ((b - a*cos(x))*csc(x)^3)/(3*(a^2 - b^2))]

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(g*sin(e+f*x))^(p/2)*(a+b*cos(e+f*x))^m=#

    #= ::Subsubsection::Closed:: =#
    #=m>0=#

    @test_int [(a + b*cos(c + d*x))*(e*sin(c + d*x))^(7/2), x, 5, (10*a*e^4*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(21*d*sqrt(e*sin(c + d*x))) - (10*a*e^3*cos(c + d*x)*sqrt(e*sin(c + d*x)))/(21*d) - (2*a*e*cos(c + d*x)*(e*sin(c + d*x))^(5/2))/(7*d) + (2*b*(e*sin(c + d*x))^(9/2))/(9*d*e)]
    @test_int [(a + b*cos(c + d*x))*(e*sin(c + d*x))^(5/2), x, 4, (6*a*e^2*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(5*d*sqrt(sin(c + d*x))) - (2*a*e*cos(c + d*x)*(e*sin(c + d*x))^(3/2))/(5*d) + (2*b*(e*sin(c + d*x))^(7/2))/(7*d*e)]
    @test_int [(a + b*cos(c + d*x))*(e*sin(c + d*x))^(3/2), x, 4, (2*a*e^2*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(3*d*sqrt(e*sin(c + d*x))) - (2*a*e*cos(c + d*x)*sqrt(e*sin(c + d*x)))/(3*d) + (2*b*(e*sin(c + d*x))^(5/2))/(5*d*e)]
    @test_int [(a + b*cos(c + d*x))*(e*sin(c + d*x))^(1/2), x, 3, (2*a*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(d*sqrt(sin(c + d*x))) + (2*b*(e*sin(c + d*x))^(3/2))/(3*d*e)]
    @test_int [(a + b*cos(c + d*x))/(e*sin(c + d*x))^(1/2), x, 3, (2*a*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(d*sqrt(e*sin(c + d*x))) + (2*b*sqrt(e*sin(c + d*x)))/(d*e)]
    @test_int [(a + b*cos(c + d*x))/(e*sin(c + d*x))^(3/2), x, 4, (-2*b)/(d*e*sqrt(e*sin(c + d*x))) - (2*a*cos(c + d*x))/(d*e*sqrt(e*sin(c + d*x))) - (2*a*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(d*e^2*sqrt(sin(c + d*x)))]
    @test_int [(a + b*cos(c + d*x))/(e*sin(c + d*x))^(5/2), x, 4, (-2*b)/(3*d*e*(e*sin(c + d*x))^(3/2)) - (2*a*cos(c + d*x))/(3*d*e*(e*sin(c + d*x))^(3/2)) + (2*a*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(3*d*e^2*sqrt(e*sin(c + d*x)))]
    @test_int [(a + b*cos(c + d*x))/(e*sin(c + d*x))^(7/2), x, 5, (-2*b)/(5*d*e*(e*sin(c + d*x))^(5/2)) - (2*a*cos(c + d*x))/(5*d*e*(e*sin(c + d*x))^(5/2)) - (6*a*cos(c + d*x))/(5*d*e^3*sqrt(e*sin(c + d*x))) - (6*a*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(5*d*e^4*sqrt(sin(c + d*x)))]

    @test_int [(a + b*cos(c + d*x))^2*(e*sin(c + d*x))^(7/2), x, 6, (10*(11*a^2 + 2*b^2)*e^4*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(231*d*sqrt(e*sin(c + d*x))) - (10*(11*a^2 + 2*b^2)*e^3*cos(c + d*x)*sqrt(e*sin(c + d*x)))/(231*d) - (2*(11*a^2 + 2*b^2)*e*cos(c + d*x)*(e*sin(c + d*x))^(5/2))/(77*d) + (26*a*b*(e*sin(c + d*x))^(9/2))/(99*d*e) + (2*b*(a + b*cos(c + d*x))*(e*sin(c + d*x))^(9/2))/(11*d*e)]
    @test_int [(a + b*cos(c + d*x))^2*(e*sin(c + d*x))^(5/2), x, 5, (2*(9*a^2 + 2*b^2)*e^2*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(15*d*sqrt(sin(c + d*x))) - (2*(9*a^2 + 2*b^2)*e*cos(c + d*x)*(e*sin(c + d*x))^(3/2))/(45*d) + (22*a*b*(e*sin(c + d*x))^(7/2))/(63*d*e) + (2*b*(a + b*cos(c + d*x))*(e*sin(c + d*x))^(7/2))/(9*d*e)]
    @test_int [(a + b*cos(c + d*x))^2*(e*sin(c + d*x))^(3/2), x, 5, (2*(7*a^2 + 2*b^2)*e^2*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(21*d*sqrt(e*sin(c + d*x))) - (2*(7*a^2 + 2*b^2)*e*cos(c + d*x)*sqrt(e*sin(c + d*x)))/(21*d) + (18*a*b*(e*sin(c + d*x))^(5/2))/(35*d*e) + (2*b*(a + b*cos(c + d*x))*(e*sin(c + d*x))^(5/2))/(7*d*e)]
    @test_int [(a + b*cos(c + d*x))^2*(e*sin(c + d*x))^(1/2), x, 4, (2*(5*a^2 + 2*b^2)*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(5*d*sqrt(sin(c + d*x))) + (14*a*b*(e*sin(c + d*x))^(3/2))/(15*d*e) + (2*b*(a + b*cos(c + d*x))*(e*sin(c + d*x))^(3/2))/(5*d*e)]
    @test_int [(a + b*cos(c + d*x))^2/(e*sin(c + d*x))^(1/2), x, 4, (2*(3*a^2 + 2*b^2)*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(3*d*sqrt(e*sin(c + d*x))) + (10*a*b*sqrt(e*sin(c + d*x)))/(3*d*e) + (2*b*(a + b*cos(c + d*x))*sqrt(e*sin(c + d*x)))/(3*d*e)]
    @test_int [(a + b*cos(c + d*x))^2/(e*sin(c + d*x))^(3/2), x, 4, (-2*(b + a*cos(c + d*x))*(a + b*cos(c + d*x)))/(d*e*sqrt(e*sin(c + d*x))) - (2*(a^2 + 2*b^2)*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(d*e^2*sqrt(sin(c + d*x))) - (2*a*b*(e*sin(c + d*x))^(3/2))/(d*e^3)]
    @test_int [(a + b*cos(c + d*x))^2/(e*sin(c + d*x))^(5/2), x, 4, (-2*(b + a*cos(c + d*x))*(a + b*cos(c + d*x)))/(3*d*e*(e*sin(c + d*x))^(3/2)) + (2*(a^2 - 2*b^2)*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(3*d*e^2*sqrt(e*sin(c + d*x))) - (2*a*b*sqrt(e*sin(c + d*x)))/(3*d*e^3)]
    @test_int [(a + b*cos(c + d*x))^2/(e*sin(c + d*x))^(7/2), x, 5, (-2*(b + a*cos(c + d*x))*(a + b*cos(c + d*x)))/(5*d*e*(e*sin(c + d*x))^(5/2)) - (2*a*b)/(5*d*e^3*sqrt(e*sin(c + d*x))) - (2*(3*a^2 - 2*b^2)*cos(c + d*x))/(5*d*e^3*sqrt(e*sin(c + d*x))) - (2*(3*a^2 - 2*b^2)*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(5*d*e^4*sqrt(sin(c + d*x)))]

    @test_int [(a + b*cos(c + d*x))^3*(e*sin(c + d*x))^(7/2), x, 7, (10*a*(11*a^2 + 6*b^2)*e^4*Elliptic.F((1/2)*(c - pi/2 + d*x), 2)*sqrt(sin(c + d*x)))/(231*d*sqrt(e*sin(c + d*x))) - (10*a*(11*a^2 + 6*b^2)*e^3*cos(c + d*x)*sqrt(e*sin(c + d*x)))/(231*d) - (2*a*(11*a^2 + 6*b^2)*e*cos(c + d*x)*(e*sin(c + d*x))^(5/2))/(77*d) + (2*b*(177*a^2 + 44*b^2)*(e*sin(c + d*x))^(9/2))/(1287*d*e) + (34*a*b*(a + b*cos(c + d*x))*(e*sin(c + d*x))^(9/2))/(143*d*e) + (2*b*(a + b*cos(c + d*x))^2*(e*sin(c + d*x))^(9/2))/(13*d*e)]
    @test_int [(a + b*cos(c + d*x))^3*(e*sin(c + d*x))^(5/2), x, 6, (2*a*(3*a^2 + 2*b^2)*e^2*Elliptic.E((1/2)*(c - pi/2 + d*x), 2)*sqrt(e*sin(c + d*x)))/(5*d*sqrt(sin(c + d*x))) - (2*a*(3*a^2 + 2*b^2)*e*cos(c + d*x)*(e*sin(c + d*x))^(3/2))/(15*d) + (2*b*(43*a^2 + 12*b^2)*(e*sin(c + d*x))^(7/2))/(231*d*e) + (10*a*b*(a + b*cos(c + d*x))*(e*sin(c + d*x))^(7/2))/(33*d*e) + (2*b*(a + b*cos(c + d*x))^2*(e*sin(c + d*x))^(7/2))/(11*d*e)]
    @test_int [(a + b*cos(c + d*x))^3*(e*sin(c + d*x))^(3/2), x, 6, (2*a*(7*a^2 + 6*b^2)*e^2*Elliptic.F((1/2)*(c - pi/2 + d*x), 2)*sqrt(sin(c + d*x)))/(21*d*sqrt(e*sin(c + d*x))) - (2*a*(7*a^2 + 6*b^2)*e*cos(c + d*x)*sqrt(e*sin(c + d*x)))/(21*d) + (2*b*(89*a^2 + 28*b^2)*(e*sin(c + d*x))^(5/2))/(315*d*e) + (26*a*b*(a + b*cos(c + d*x))*(e*sin(c + d*x))^(5/2))/(63*d*e) + (2*b*(a + b*cos(c + d*x))^2*(e*sin(c + d*x))^(5/2))/(9*d*e)]
    @test_int [(a + b*cos(c + d*x))^3*(e*sin(c + d*x))^(1/2), x, 5, (2*a*(5*a^2 + 6*b^2)*Elliptic.E((1/2)*(c - pi/2 + d*x), 2)*sqrt(e*sin(c + d*x)))/(5*d*sqrt(sin(c + d*x))) + (2*b*(57*a^2 + 20*b^2)*(e*sin(c + d*x))^(3/2))/(105*d*e) + (22*a*b*(a + b*cos(c + d*x))*(e*sin(c + d*x))^(3/2))/(35*d*e) + (2*b*(a + b*cos(c + d*x))^2*(e*sin(c + d*x))^(3/2))/(7*d*e)]
    @test_int [(a + b*cos(c + d*x))^3/(e*sin(c + d*x))^(1/2), x, 5, (2*a*(a^2 + 2*b^2)*Elliptic.F((1/2)*(c - pi/2 + d*x), 2)*sqrt(sin(c + d*x)))/(d*sqrt(e*sin(c + d*x))) + (2*b*(11*a^2 + 4*b^2)*sqrt(e*sin(c + d*x)))/(5*d*e) + (6*a*b*(a + b*cos(c + d*x))*sqrt(e*sin(c + d*x)))/(5*d*e) + (2*b*(a + b*cos(c + d*x))^2*sqrt(e*sin(c + d*x)))/(5*d*e)]
    @test_int [(a + b*cos(c + d*x))^3/(e*sin(c + d*x))^(3/2), x, 5, -((2*(b + a*cos(c + d*x))*(a + b*cos(c + d*x))^2)/(d*e*sqrt(e*sin(c + d*x)))) - (2*a*(a^2 + 6*b^2)*Elliptic.E((1/2)*(c - pi/2 + d*x), 2)*sqrt(e*sin(c + d*x)))/(d*e^2*sqrt(sin(c + d*x))) - (2*b*(3*a^2 + 4*b^2)*(e*sin(c + d*x))^(3/2))/(3*d*e^3) - (2*a*b*(a + b*cos(c + d*x))*(e*sin(c + d*x))^(3/2))/(d*e^3)]
    @test_int [(a + b*cos(c + d*x))^3/(e*sin(c + d*x))^(5/2), x, 5, -((2*(b + a*cos(c + d*x))*(a + b*cos(c + d*x))^2)/(3*d*e*(e*sin(c + d*x))^(3/2))) + (2*a*(a^2 - 6*b^2)*Elliptic.F((1/2)*(c - pi/2 + d*x), 2)*sqrt(sin(c + d*x)))/(3*d*e^2*sqrt(e*sin(c + d*x))) - (2*b*(a^2 + 4*b^2)*sqrt(e*sin(c + d*x)))/(3*d*e^3) - (2*a*b*(a + b*cos(c + d*x))*sqrt(e*sin(c + d*x)))/(3*d*e^3)]
    @test_int [(a + b*cos(c + d*x))^3/(e*sin(c + d*x))^(7/2), x, 5, -((2*(b + a*cos(c + d*x))*(a + b*cos(c + d*x))^2)/(5*d*e*(e*sin(c + d*x))^(5/2))) + (2*(a + b*cos(c + d*x))*(a*b - (3*a^2 - 4*b^2)*cos(c + d*x)))/(5*d*e^3*sqrt(e*sin(c + d*x))) - (6*a*(a^2 - 2*b^2)*Elliptic.E((1/2)*(c - pi/2 + d*x), 2)*sqrt(e*sin(c + d*x)))/(5*d*e^4*sqrt(sin(c + d*x))) - (2*b*(3*a^2 - 4*b^2)*(e*sin(c + d*x))^(3/2))/(5*d*e^5)]
    @test_int [(a + b*cos(c + d*x))^3/(e*sin(c + d*x))^(9/2), x, 5, -((2*(b + a*cos(c + d*x))*(a + b*cos(c + d*x))^2)/(7*d*e*(e*sin(c + d*x))^(7/2))) - (2*(a + b*cos(c + d*x))*(a*b + (5*a^2 - 4*b^2)*cos(c + d*x)))/(21*d*e^3*(e*sin(c + d*x))^(3/2)) + (2*a*(5*a^2 - 6*b^2)*Elliptic.F((1/2)*(c - pi/2 + d*x), 2)*sqrt(sin(c + d*x)))/(21*d*e^4*sqrt(e*sin(c + d*x))) - (2*b*(5*a^2 - 4*b^2)*sqrt(e*sin(c + d*x)))/(21*d*e^5)]

    #= ::Subsubsection::Closed:: =#
    #=m<0=#

    @test_int [(e*sin(c + d*x))^(11/2)/(a + b*cos(c + d*x)), x, 15, ((-a^2 + b^2)^(9/4)*e^(11/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(b^(11/2)*d) + ((-a^2 + b^2)^(9/4)*e^(11/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(b^(11/2)*d) + (2*a*(21*a^4 - 49*a^2*b^2 + 33*b^4)*e^6*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(21*b^6*d*sqrt(e*sin(c + d*x))) - (a*(a^2 - b^2)^3*e^6*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(b^6*(a^2 - b*(b - sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) - (a*(a^2 - b^2)^3*e^6*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(b^6*(a^2 - b*(b + sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) - (2*e^5*(21*(a^2 - b^2)^2 - a*b*(7*a^2 - 12*b^2)*cos(c + d*x))*sqrt(e*sin(c + d*x)))/(21*b^5*d) + (2*e^3*(7*(a^2 - b^2) - 5*a*b*cos(c + d*x))*(e*sin(c + d*x))^(5/2))/(35*b^3*d) - (2*e*(e*sin(c + d*x))^(9/2))/(9*b*d)]
    @test_int [(e*sin(c + d*x))^(9/2)/(a + b*cos(c + d*x)), x, 14, -(((-a^2 + b^2)^(7/4)*e^(9/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(b^(9/2)*d)) + ((-a^2 + b^2)^(7/4)*e^(9/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(b^(9/2)*d) + (a*(a^2 - b^2)^2*e^5*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(b^5*(b - sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) + (a*(a^2 - b^2)^2*e^5*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(b^5*(b + sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) - (2*a*(5*a^2 - 8*b^2)*e^4*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(5*b^4*d*sqrt(sin(c + d*x))) + (2*e^3*(5*(a^2 - b^2) - 3*a*b*cos(c + d*x))*(e*sin(c + d*x))^(3/2))/(15*b^3*d) - (2*e*(e*sin(c + d*x))^(7/2))/(7*b*d)]
    @test_int [(e*sin(c + d*x))^(7/2)/(a + b*cos(c + d*x)), x, 14, ((-a^2 + b^2)^(5/4)*e^(7/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(b^(7/2)*d) + ((-a^2 + b^2)^(5/4)*e^(7/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(b^(7/2)*d) - (2*a*(3*a^2 - 4*b^2)*e^4*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(3*b^4*d*sqrt(e*sin(c + d*x))) + (a*(a^2 - b^2)^2*e^4*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(b^4*(a^2 - b*(b - sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) + (a*(a^2 - b^2)^2*e^4*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(b^4*(a^2 - b*(b + sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) + (2*e^3*(3*(a^2 - b^2) - a*b*cos(c + d*x))*sqrt(e*sin(c + d*x)))/(3*b^3*d) - (2*e*(e*sin(c + d*x))^(5/2))/(5*b*d)]
    @test_int [(e*sin(c + d*x))^(5/2)/(a + b*cos(c + d*x)), x, 13, -(((-a^2 + b^2)^(3/4)*e^(5/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(b^(5/2)*d)) + ((-a^2 + b^2)^(3/4)*e^(5/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(b^(5/2)*d) - (a*(a^2 - b^2)*e^3*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(b^3*(b - sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) - (a*(a^2 - b^2)*e^3*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(b^3*(b + sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) + (2*a*e^2*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(b^2*d*sqrt(sin(c + d*x))) - (2*e*(e*sin(c + d*x))^(3/2))/(3*b*d)]
    @test_int [(e*sin(c + d*x))^(3/2)/(a + b*cos(c + d*x)), x, 13, ((-a^2 + b^2)^(1/4)*e^(3/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(b^(3/2)*d) + ((-a^2 + b^2)^(1/4)*e^(3/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(b^(3/2)*d) + (2*a*e^2*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(b^2*d*sqrt(e*sin(c + d*x))) - (a*(a^2 - b^2)*e^2*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(b^2*(a^2 - b*(b - sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) - (a*(a^2 - b^2)*e^2*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(b^2*(a^2 - b*(b + sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) - (2*e*sqrt(e*sin(c + d*x)))/(b*d)]
    @test_int [(e*sin(c + d*x))^(1/2)/(a + b*cos(c + d*x)), x, 9, -((sqrt(e)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(sqrt(b)*(-a^2 + b^2)^(1/4)*d)) + (sqrt(e)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(sqrt(b)*(-a^2 + b^2)^(1/4)*d) + (a*e*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(b*(b - sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) + (a*e*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(b*(b + sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x)))]
    @test_int [1/((a + b*cos(c + d*x))*(e*sin(c + d*x))^(1/2)), x, 9, (sqrt(b)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/((-a^2 + b^2)^(3/4)*d*sqrt(e)) + (sqrt(b)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/((-a^2 + b^2)^(3/4)*d*sqrt(e)) + (a*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/((a^2 - b*(b - sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) + (a*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/((a^2 - b*(b + sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x)))]
    @test_int [1/((a + b*cos(c + d*x))*(e*sin(c + d*x))^(3/2)), x, 13, -((b^(3/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/((-a^2 + b^2)^(5/4)*d*e^(3/2))) + (b^(3/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/((-a^2 + b^2)^(5/4)*d*e^(3/2)) + (2*(b - a*cos(c + d*x)))/((a^2 - b^2)*d*e*sqrt(e*sin(c + d*x))) - (a*b*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/((a^2 - b^2)*(b - sqrt(-a^2 + b^2))*d*e*sqrt(e*sin(c + d*x))) - (a*b*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/((a^2 - b^2)*(b + sqrt(-a^2 + b^2))*d*e*sqrt(e*sin(c + d*x))) - (2*a*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/((a^2 - b^2)*d*e^2*sqrt(sin(c + d*x)))]
    @test_int [1/((a + b*cos(c + d*x))*(e*sin(c + d*x))^(5/2)), x, 13, (b^(5/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/((-a^2 + b^2)^(7/4)*d*e^(5/2)) + (b^(5/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/((-a^2 + b^2)^(7/4)*d*e^(5/2)) + (2*(b - a*cos(c + d*x)))/(3*(a^2 - b^2)*d*e*(e*sin(c + d*x))^(3/2)) + (2*a*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(3*(a^2 - b^2)*d*e^2*sqrt(e*sin(c + d*x))) - (a*b^2*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/((a^2 - b^2)*(a^2 - b*(b - sqrt(-a^2 + b^2)))*d*e^2*sqrt(e*sin(c + d*x))) - (a*b^2*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/((a^2 - b^2)*(a^2 - b*(b + sqrt(-a^2 + b^2)))*d*e^2*sqrt(e*sin(c + d*x)))]
    @test_int [1/((a + b*cos(c + d*x))*(e*sin(c + d*x))^(7/2)), x, 14, -((b^(7/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/((-a^2 + b^2)^(9/4)*d*e^(7/2))) + (b^(7/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/((-a^2 + b^2)^(9/4)*d*e^(7/2)) + (2*(b - a*cos(c + d*x)))/(5*(a^2 - b^2)*d*e*(e*sin(c + d*x))^(5/2)) - (2*(5*b^3 + a*(3*a^2 - 8*b^2)*cos(c + d*x)))/(5*(a^2 - b^2)^2*d*e^3*sqrt(e*sin(c + d*x))) + (a*b^3*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/((a^2 - b^2)^2*(b - sqrt(-a^2 + b^2))*d*e^3*sqrt(e*sin(c + d*x))) + (a*b^3*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/((a^2 - b^2)^2*(b + sqrt(-a^2 + b^2))*d*e^3*sqrt(e*sin(c + d*x))) - (2*a*(3*a^2 - 8*b^2)*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(5*(a^2 - b^2)^2*d*e^4*sqrt(sin(c + d*x)))]

    @test_int [(e*sin(c + d*x))^(11/2)/(a + b*cos(c + d*x))^2, x, 15, (9*a*(-a^2 + b^2)^(5/4)*e^(11/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*b^(11/2)*d) + (9*a*(-a^2 + b^2)^(5/4)*e^(11/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*b^(11/2)*d) - (3*(21*a^4 - 28*a^2*b^2 + 5*b^4)*e^6*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(7*b^6*d*sqrt(e*sin(c + d*x))) + (9*a^2*(a^2 - b^2)^2*e^6*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*b^6*(a^2 - b*(b - sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) + (9*a^2*(a^2 - b^2)^2*e^6*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*b^6*(a^2 - b*(b + sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) + (3*e^5*(21*a*(a^2 - b^2) - b*(7*a^2 - 5*b^2)*cos(c + d*x))*sqrt(e*sin(c + d*x)))/(7*b^5*d) - (9*e^3*(7*a - 5*b*cos(c + d*x))*(e*sin(c + d*x))^(5/2))/(35*b^3*d) + (e*(e*sin(c + d*x))^(9/2))/(b*d*(a + b*cos(c + d*x)))]
    @test_int [(e*sin(c + d*x))^(9/2)/(a + b*cos(c + d*x))^2, x, 14, (-7*a*(-a^2 + b^2)^(3/4)*e^(9/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*b^(9/2)*d) + (7*a*(-a^2 + b^2)^(3/4)*e^(9/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*b^(9/2)*d) - (7*a^2*(a^2 - b^2)*e^5*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*b^5*(b - sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) - (7*a^2*(a^2 - b^2)*e^5*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*b^5*(b + sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) + (7*(5*a^2 - 3*b^2)*e^4*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(5*b^4*d*sqrt(sin(c + d*x))) - (7*e^3*(5*a - 3*b*cos(c + d*x))*(e*sin(c + d*x))^(3/2))/(15*b^3*d) + (e*(e*sin(c + d*x))^(7/2))/(b*d*(a + b*cos(c + d*x)))]
    @test_int [(e*sin(c + d*x))^(7/2)/(a + b*cos(c + d*x))^2, x, 14, (5*a*(-a^2 + b^2)^(1/4)*e^(7/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*b^(7/2)*d) + (5*a*(-a^2 + b^2)^(1/4)*e^(7/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*b^(7/2)*d) + (5*(3*a^2 - b^2)*e^4*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(3*b^4*d*sqrt(e*sin(c + d*x))) - (5*a^2*(a^2 - b^2)*e^4*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*b^4*(a^2 - b*(b - sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) - (5*a^2*(a^2 - b^2)*e^4*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*b^4*(a^2 - b*(b + sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) - (5*e^3*(3*a - b*cos(c + d*x))*sqrt(e*sin(c + d*x)))/(3*b^3*d) + (e*(e*sin(c + d*x))^(5/2))/(b*d*(a + b*cos(c + d*x)))]
    @test_int [(e*sin(c + d*x))^(5/2)/(a + b*cos(c + d*x))^2, x, 13, (-3*a*e^(5/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*b^(5/2)*(-a^2 + b^2)^(1/4)*d) + (3*a*e^(5/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*b^(5/2)*(-a^2 + b^2)^(1/4)*d) + (3*a^2*e^3*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*b^3*(b - sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) + (3*a^2*e^3*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*b^3*(b + sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) - (3*e^2*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(b^2*d*sqrt(sin(c + d*x))) + (e*(e*sin(c + d*x))^(3/2))/(b*d*(a + b*cos(c + d*x)))]
    @test_int [(e*sin(c + d*x))^(3/2)/(a + b*cos(c + d*x))^2, x, 13, (a*e^(3/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*b^(3/2)*(-a^2 + b^2)^(3/4)*d) + (a*e^(3/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*b^(3/2)*(-a^2 + b^2)^(3/4)*d) - (e^2*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(b^2*d*sqrt(e*sin(c + d*x))) + (a^2*e^2*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*b^2*(a^2 - b*(b - sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) + (a^2*e^2*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*b^2*(a^2 - b*(b + sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) + (e*sqrt(e*sin(c + d*x)))/(b*d*(a + b*cos(c + d*x)))]
    @test_int [(e*sin(c + d*x))^(1/2)/(a + b*cos(c + d*x))^2, x, 13, (a*sqrt(e)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*sqrt(b)*(-a^2 + b^2)^(5/4)*d) - (a*sqrt(e)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*sqrt(b)*(-a^2 + b^2)^(5/4)*d) + (a^2*e*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*b*(a^2 - b^2)*(b - sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) + (a^2*e*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*b*(a^2 - b^2)*(b + sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) + (Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/((a^2 - b^2)*d*sqrt(sin(c + d*x))) - (b*(e*sin(c + d*x))^(3/2))/((a^2 - b^2)*d*e*(a + b*cos(c + d*x)))]
    @test_int [1/((a + b*cos(c + d*x))^2*(e*sin(c + d*x))^(1/2)), x, 13, (-3*a*sqrt(b)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*(-a^2 + b^2)^(7/4)*d*sqrt(e)) - (3*a*sqrt(b)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*(-a^2 + b^2)^(7/4)*d*sqrt(e)) - (Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/((a^2 - b^2)*d*sqrt(e*sin(c + d*x))) + (3*a^2*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*(a^2 - b^2)*(a^2 - b*(b - sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) + (3*a^2*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*(a^2 - b^2)*(a^2 - b*(b + sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) - (b*sqrt(e*sin(c + d*x)))/((a^2 - b^2)*d*e*(a + b*cos(c + d*x)))]
    @test_int [1/((a + b*cos(c + d*x))^2*(e*sin(c + d*x))^(3/2)), x, 14, (5*a*b^(3/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*(-a^2 + b^2)^(9/4)*d*e^(3/2)) - (5*a*b^(3/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*(-a^2 + b^2)^(9/4)*d*e^(3/2)) - b/((a^2 - b^2)*d*e*(a + b*cos(c + d*x))*sqrt(e*sin(c + d*x))) + (5*a*b - (2*a^2 + 3*b^2)*cos(c + d*x))/((a^2 - b^2)^2*d*e*sqrt(e*sin(c + d*x))) - (5*a^2*b*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*(a^2 - b^2)^2*(b - sqrt(-a^2 + b^2))*d*e*sqrt(e*sin(c + d*x))) - (5*a^2*b*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*(a^2 - b^2)^2*(b + sqrt(-a^2 + b^2))*d*e*sqrt(e*sin(c + d*x))) - ((2*a^2 + 3*b^2)*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/((a^2 - b^2)^2*d*e^2*sqrt(sin(c + d*x)))]
    @test_int [1/((a + b*cos(c + d*x))^2*(e*sin(c + d*x))^(5/2)), x, 14, (-7*a*b^(5/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*(-a^2 + b^2)^(11/4)*d*e^(5/2)) - (7*a*b^(5/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*(-a^2 + b^2)^(11/4)*d*e^(5/2)) - b/((a^2 - b^2)*d*e*(a + b*cos(c + d*x))*(e*sin(c + d*x))^(3/2)) + (7*a*b - (2*a^2 + 5*b^2)*cos(c + d*x))/(3*(a^2 - b^2)^2*d*e*(e*sin(c + d*x))^(3/2)) + ((2*a^2 + 5*b^2)*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(3*(a^2 - b^2)^2*d*e^2*sqrt(e*sin(c + d*x))) - (7*a^2*b^2*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*(a^2 - b^2)^2*(a^2 - b*(b - sqrt(-a^2 + b^2)))*d*e^2*sqrt(e*sin(c + d*x))) - (7*a^2*b^2*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*(a^2 - b^2)^2*(a^2 - b*(b + sqrt(-a^2 + b^2)))*d*e^2*sqrt(e*sin(c + d*x)))]
    @test_int [1/((a + b*cos(c + d*x))^2*(e*sin(c + d*x))^(7/2)), x, 15, (9*a*b^(7/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*(-a^2 + b^2)^(13/4)*d*e^(7/2)) - (9*a*b^(7/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(2*(-a^2 + b^2)^(13/4)*d*e^(7/2)) - b/((a^2 - b^2)*d*e*(a + b*cos(c + d*x))*(e*sin(c + d*x))^(5/2)) + (9*a*b - (2*a^2 + 7*b^2)*cos(c + d*x))/(5*(a^2 - b^2)^2*d*e*(e*sin(c + d*x))^(5/2)) - (3*(15*a*b^3 + (2*a^4 - 10*a^2*b^2 - 7*b^4)*cos(c + d*x)))/(5*(a^2 - b^2)^3*d*e^3*sqrt(e*sin(c + d*x))) + (9*a^2*b^3*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*(a^2 - b^2)^3*(b - sqrt(-a^2 + b^2))*d*e^3*sqrt(e*sin(c + d*x))) + (9*a^2*b^3*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(2*(a^2 - b^2)^3*(b + sqrt(-a^2 + b^2))*d*e^3*sqrt(e*sin(c + d*x))) - (3*(2*a^4 - 10*a^2*b^2 - 7*b^4)*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(5*(a^2 - b^2)^3*d*e^4*sqrt(sin(c + d*x)))]

    @test_int [(e*sin(c + d*x))^(13/2)/(a + b*cos(c + d*x))^3, x, 15, (11*(9*a^4 - 11*a^2*b^2 + 2*b^4)*e^(13/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*b^(13/2)*(-a^2 + b^2)^(1/4)*d) - (11*(9*a^4 - 11*a^2*b^2 + 2*b^4)*e^(13/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*b^(13/2)*(-a^2 + b^2)^(1/4)*d) - (11*a*(9*a^4 - 11*a^2*b^2 + 2*b^4)*e^7*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*b^7*(b - sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) - (11*a*(9*a^4 - 11*a^2*b^2 + 2*b^4)*e^7*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*b^7*(b + sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) + (11*a*(45*a^2 - 37*b^2)*e^6*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(20*b^6*d*sqrt(sin(c + d*x))) - (11*e^5*(5*(9*a^2 - 2*b^2) - 27*a*b*cos(c + d*x))*(e*sin(c + d*x))^(3/2))/(60*b^5*d) + (11*e^3*(9*a + 2*b*cos(c + d*x))*(e*sin(c + d*x))^(7/2))/(28*b^3*d*(a + b*cos(c + d*x))) + (e*(e*sin(c + d*x))^(11/2))/(2*b*d*(a + b*cos(c + d*x))^2)]
    @test_int [(e*sin(c + d*x))^(11/2)/(a + b*cos(c + d*x))^3, x, 15, (-9*(7*a^4 - 9*a^2*b^2 + 2*b^4)*e^(11/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*b^(11/2)*(-a^2 + b^2)^(3/4)*d) - (9*(7*a^4 - 9*a^2*b^2 + 2*b^4)*e^(11/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*b^(11/2)*(-a^2 + b^2)^(3/4)*d) + (3*a*(21*a^2 - 13*b^2)*e^6*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(4*b^6*d*sqrt(e*sin(c + d*x))) - (9*a*(7*a^4 - 9*a^2*b^2 + 2*b^4)*e^6*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*b^6*(a^2 - b*(b - sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) - (9*a*(7*a^4 - 9*a^2*b^2 + 2*b^4)*e^6*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*b^6*(a^2 - b*(b + sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) - (3*e^5*(3*(7*a^2 - 2*b^2) - 7*a*b*cos(c + d*x))*sqrt(e*sin(c + d*x)))/(4*b^5*d) + (9*e^3*(7*a + 2*b*cos(c + d*x))*(e*sin(c + d*x))^(5/2))/(20*b^3*d*(a + b*cos(c + d*x))) + (e*(e*sin(c + d*x))^(9/2))/(2*b*d*(a + b*cos(c + d*x))^2)]
    @test_int [(e*sin(c + d*x))^(9/2)/(a + b*cos(c + d*x))^3, x, 14, (-7*(5*a^2 - 2*b^2)*e^(9/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*b^(9/2)*(-a^2 + b^2)^(1/4)*d) + (7*(5*a^2 - 2*b^2)*e^(9/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*b^(9/2)*(-a^2 + b^2)^(1/4)*d) + (7*a*(5*a^2 - 2*b^2)*e^5*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*b^5*(b - sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) + (7*a*(5*a^2 - 2*b^2)*e^5*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*b^5*(b + sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) - (35*a*e^4*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(4*b^4*d*sqrt(sin(c + d*x))) + (7*e^3*(5*a + 2*b*cos(c + d*x))*(e*sin(c + d*x))^(3/2))/(12*b^3*d*(a + b*cos(c + d*x))) + (e*(e*sin(c + d*x))^(7/2))/(2*b*d*(a + b*cos(c + d*x))^2)]
    @test_int [(e*sin(c + d*x))^(7/2)/(a + b*cos(c + d*x))^3, x, 14, (5*(3*a^2 - 2*b^2)*e^(7/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*b^(7/2)*(-a^2 + b^2)^(3/4)*d) + (5*(3*a^2 - 2*b^2)*e^(7/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*b^(7/2)*(-a^2 + b^2)^(3/4)*d) - (15*a*e^4*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(4*b^4*d*sqrt(e*sin(c + d*x))) + (5*a*(3*a^2 - 2*b^2)*e^4*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*b^4*(a^2 - b*(b - sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) + (5*a*(3*a^2 - 2*b^2)*e^4*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*b^4*(a^2 - b*(b + sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) + (5*e^3*(3*a + 2*b*cos(c + d*x))*sqrt(e*sin(c + d*x)))/(4*b^3*d*(a + b*cos(c + d*x))) + (e*(e*sin(c + d*x))^(5/2))/(2*b*d*(a + b*cos(c + d*x))^2)]
    @test_int [(e*sin(c + d*x))^(5/2)/(a + b*cos(c + d*x))^3, x, 14, (-3*(a^2 - 2*b^2)*e^(5/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*b^(5/2)*(-a^2 + b^2)^(5/4)*d) + (3*(a^2 - 2*b^2)*e^(5/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*b^(5/2)*(-a^2 + b^2)^(5/4)*d) - (3*a*(a^2 - 2*b^2)*e^3*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*b^3*(a^2 - b^2)*(b - sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) - (3*a*(a^2 - 2*b^2)*e^3*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*b^3*(a^2 - b^2)*(b + sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) + (3*a*e^2*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(4*b^2*(a^2 - b^2)*d*sqrt(sin(c + d*x))) + (e*(e*sin(c + d*x))^(3/2))/(2*b*d*(a + b*cos(c + d*x))^2) - (3*a*e*(e*sin(c + d*x))^(3/2))/(4*b*(a^2 - b^2)*d*(a + b*cos(c + d*x)))]
    @test_int [(e*sin(c + d*x))^(3/2)/(a + b*cos(c + d*x))^3, x, 14, -((a^2 + 2*b^2)*e^(3/2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*b^(3/2)*(-a^2 + b^2)^(7/4)*d) - ((a^2 + 2*b^2)*e^(3/2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*b^(3/2)*(-a^2 + b^2)^(7/4)*d) - (a*e^2*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(4*b^2*(a^2 - b^2)*d*sqrt(e*sin(c + d*x))) + (a*(a^2 + 2*b^2)*e^2*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*b^2*(a^2 - b^2)*(a^2 - b*(b - sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) + (a*(a^2 + 2*b^2)*e^2*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*b^2*(a^2 - b^2)*(a^2 - b*(b + sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) + (e*sqrt(e*sin(c + d*x)))/(2*b*d*(a + b*cos(c + d*x))^2) - (a*e*sqrt(e*sin(c + d*x)))/(4*b*(a^2 - b^2)*d*(a + b*cos(c + d*x)))]
    @test_int [(e*sin(c + d*x))^(1/2)/(a + b*cos(c + d*x))^3, x, 14, -((3*a^2 + 2*b^2)*sqrt(e)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*sqrt(b)*(-a^2 + b^2)^(9/4)*d) + ((3*a^2 + 2*b^2)*sqrt(e)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*sqrt(b)*(-a^2 + b^2)^(9/4)*d) + (a*(3*a^2 + 2*b^2)*e*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*b*(a^2 - b^2)^2*(b - sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) + (a*(3*a^2 + 2*b^2)*e*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*b*(a^2 - b^2)^2*(b + sqrt(-a^2 + b^2))*d*sqrt(e*sin(c + d*x))) + (5*a*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(4*(a^2 - b^2)^2*d*sqrt(sin(c + d*x))) - (b*(e*sin(c + d*x))^(3/2))/(2*(a^2 - b^2)*d*e*(a + b*cos(c + d*x))^2) - (5*a*b*(e*sin(c + d*x))^(3/2))/(4*(a^2 - b^2)^2*d*e*(a + b*cos(c + d*x)))]
    @test_int [1/((a + b*cos(c + d*x))^3*(e*sin(c + d*x))^(1/2)), x, 14, (3*sqrt(b)*(5*a^2 + 2*b^2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*(-a^2 + b^2)^(11/4)*d*sqrt(e)) + (3*sqrt(b)*(5*a^2 + 2*b^2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*(-a^2 + b^2)^(11/4)*d*sqrt(e)) - (7*a*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(4*(a^2 - b^2)^2*d*sqrt(e*sin(c + d*x))) + (3*a*(5*a^2 + 2*b^2)*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*(a^2 - b^2)^2*(a^2 - b*(b - sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) + (3*a*(5*a^2 + 2*b^2)*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*(a^2 - b^2)^2*(a^2 - b*(b + sqrt(-a^2 + b^2)))*d*sqrt(e*sin(c + d*x))) - (b*sqrt(e*sin(c + d*x)))/(2*(a^2 - b^2)*d*e*(a + b*cos(c + d*x))^2) - (7*a*b*sqrt(e*sin(c + d*x)))/(4*(a^2 - b^2)^2*d*e*(a + b*cos(c + d*x)))]
    @test_int [1/((a + b*cos(c + d*x))^3*(e*sin(c + d*x))^(3/2)), x, 15, (-5*b^(3/2)*(7*a^2 + 2*b^2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*(-a^2 + b^2)^(13/4)*d*e^(3/2)) + (5*b^(3/2)*(7*a^2 + 2*b^2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*(-a^2 + b^2)^(13/4)*d*e^(3/2)) - b/(2*(a^2 - b^2)*d*e*(a + b*cos(c + d*x))^2*sqrt(e*sin(c + d*x))) - (9*a*b)/(4*(a^2 - b^2)^2*d*e*(a + b*cos(c + d*x))*sqrt(e*sin(c + d*x))) + (5*b*(7*a^2 + 2*b^2) - a*(8*a^2 + 37*b^2)*cos(c + d*x))/(4*(a^2 - b^2)^3*d*e*sqrt(e*sin(c + d*x))) - (5*a*b*(7*a^2 + 2*b^2)*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*(a^2 - b^2)^3*(b - sqrt(-a^2 + b^2))*d*e*sqrt(e*sin(c + d*x))) - (5*a*b*(7*a^2 + 2*b^2)*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*(a^2 - b^2)^3*(b + sqrt(-a^2 + b^2))*d*e*sqrt(e*sin(c + d*x))) - (a*(8*a^2 + 37*b^2)*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(4*(a^2 - b^2)^3*d*e^2*sqrt(sin(c + d*x)))]
    @test_int [1/((a + b*cos(c + d*x))^3*(e*sin(c + d*x))^(5/2)), x, 15, (7*b^(5/2)*(9*a^2 + 2*b^2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*(-a^2 + b^2)^(15/4)*d*e^(5/2)) + (7*b^(5/2)*(9*a^2 + 2*b^2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*(-a^2 + b^2)^(15/4)*d*e^(5/2)) - b/(2*(a^2 - b^2)*d*e*(a + b*cos(c + d*x))^2*(e*sin(c + d*x))^(3/2)) - (11*a*b)/(4*(a^2 - b^2)^2*d*e*(a + b*cos(c + d*x))*(e*sin(c + d*x))^(3/2)) + (7*b*(9*a^2 + 2*b^2) - a*(8*a^2 + 69*b^2)*cos(c + d*x))/(12*(a^2 - b^2)^3*d*e*(e*sin(c + d*x))^(3/2)) + (a*(8*a^2 + 69*b^2)*Elliptic.F((c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(12*(a^2 - b^2)^3*d*e^2*sqrt(e*sin(c + d*x))) - (7*a*b^2*(9*a^2 + 2*b^2)*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*(a^2 - b^2)^3*(a^2 - b*(b - sqrt(-a^2 + b^2)))*d*e^2*sqrt(e*sin(c + d*x))) - (7*a*b^2*(9*a^2 + 2*b^2)*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*(a^2 - b^2)^3*(a^2 - b*(b + sqrt(-a^2 + b^2)))*d*e^2*sqrt(e*sin(c + d*x)))]
    @test_int [1/((a + b*cos(c + d*x))^3*(e*sin(c + d*x))^(7/2)), x, 16, (-9*b^(7/2)*(11*a^2 + 2*b^2)*atan((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*(-a^2 + b^2)^(17/4)*d*e^(7/2)) + (9*b^(7/2)*(11*a^2 + 2*b^2)*atanh((sqrt(b)*sqrt(e*sin(c + d*x)))/((-a^2 + b^2)^(1/4)*sqrt(e))))/(8*(-a^2 + b^2)^(17/4)*d*e^(7/2)) - b/(2*(a^2 - b^2)*d*e*(a + b*cos(c + d*x))^2*(e*sin(c + d*x))^(5/2)) - (13*a*b)/(4*(a^2 - b^2)^2*d*e*(a + b*cos(c + d*x))*(e*sin(c + d*x))^(5/2)) + (9*b*(11*a^2 + 2*b^2) - a*(8*a^2 + 109*b^2)*cos(c + d*x))/(20*(a^2 - b^2)^3*d*e*(e*sin(c + d*x))^(5/2)) - (3*(15*b^3*(11*a^2 + 2*b^2) + a*(8*a^4 - 64*a^2*b^2 - 139*b^4)*cos(c + d*x)))/(20*(a^2 - b^2)^4*d*e^3*sqrt(e*sin(c + d*x))) + (9*a*b^3*(11*a^2 + 2*b^2)*Elliptic.Pi((2*b)/(b - sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*(a^2 - b^2)^4*(b - sqrt(-a^2 + b^2))*d*e^3*sqrt(e*sin(c + d*x))) + (9*a*b^3*(11*a^2 + 2*b^2)*Elliptic.Pi((2*b)/(b + sqrt(-a^2 + b^2)), (c - pi/2 + d*x)/2, 2)*sqrt(sin(c + d*x)))/(8*(a^2 - b^2)^4*(b + sqrt(-a^2 + b^2))*d*e^3*sqrt(e*sin(c + d*x))) - (3*a*(8*a^4 - 64*a^2*b^2 - 139*b^4)*Elliptic.E((c - pi/2 + d*x)/2, 2)*sqrt(e*sin(c + d*x)))/(20*(a^2 - b^2)^4*d*e^4*sqrt(sin(c + d*x)))]

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

    #= ::Title:: =#
    #=Integrands*of*the*form*(g*csc(e+f*x))^p*(a+b*cos(e+f*x))^m=#

    #= ::Section:: =#
    #=Integrands*of*the*form*(g*csc(e+f*x))^p*(a+b*cos(e+f*x))^m*when*a^2-b^2=0=#

    #= ::Section:: =#
    #=Integrands*of*the*form*(g*csc(e+f*x))^p*(a+b*cos(e+f*x))^m=#
    end

max_stars_repo_path: test/4 Trig functions/4.2 Cosine/4.2.1.2 (g sin)^p (a+b cos)^m.jl


text:
# Autogenerated wrapper script for qwtw_jll for x86_64-linux-gnu-cxx03
export qwproc, qwtw

using Libglvnd_jll
using boost_jll
using Qt_jll
using qwt_jll
using CompilerSupportLibraries_jll
using FreeType2_jll
using OpenSSL_jll
using marble_jll
using MathGL_jll
JLLWrappers.@generate_wrapper_header("qwtw")
JLLWrappers.@declare_executable_product(qwproc)
JLLWrappers.@declare_library_product(qwtw, "libqwtw.so")
function __init__()
    JLLWrappers.@generate_init_header(Libglvnd_jll, boost_jll, Qt_jll, qwt_jll, CompilerSupportLibraries_jll, FreeType2_jll, OpenSSL_jll, marble_jll, MathGL_jll)
    JLLWrappers.@init_executable_product(
        qwproc,
        "bin/qwproc",
    )

    JLLWrappers.@init_library_product(
        qwtw,
        "lib64/libqwtw.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ b8293916-42e5-11eb-0a82-11d14d1b4ba1
using Pkg; Pkg.status()

# ╔═╡ d5ea666e-42e5-11eb-2d9b-432c08b84a98
using Camoforge

# ╔═╡ 7009f8c4-42e9-11eb-2775-eb9f6cce20a2
using Colors

# ╔═╡ 4f305f60-42eb-11eb-33aa-1548a90f866a
using StatsBase

# ╔═╡ 3aa78d38-42f1-11eb-1005-2bf725559cb3
using Images

# ╔═╡ 34165c74-42fb-11eb-2d74-512a600c05ba
using PyCall

# ╔═╡ e0c451bc-42e5-11eb-1e3e-417e400304ad
Camoforge.get_image("/home/brad/Desktop/test.jpg")

# ╔═╡ fd62b854-42e5-11eb-0880-cf34d4d189ae
test = Camoforge.get_image("/home/brad/Desktop/test.jpg")

# ╔═╡ 28456c56-42e6-11eb-31f3-dfadfd7bc8a1
test[end-9]

# ╔═╡ 4fd8f040-42ea-11eb-2548-bdb8e7cd15ec
md"### Info on arrays"

# ╔═╡ 5a3a1c3a-42ea-11eb-134a-514790026194
md"[https://en.wikibooks.org/wiki/Introducing_Julia/Arrays_and_tuples#Finding_out_about_an_array](https://en.wikibooks.org/wiki/Introducing_Julia/Arrays_and_tuples#Finding_out_about_an_array)"

# ╔═╡ 6b156e22-42e9-11eb-0e22-9d27615d22cc
Colors.hex(test[end])

# ╔═╡ d963bbd6-42e9-11eb-3c7d-3b2722e77bc2
ndims(test)

# ╔═╡ 37d03b70-42ea-11eb-1a05-ff87bd42f290
size(test)

# ╔═╡ 6f44c31e-42ea-11eb-3507-cbff33f47bc0
to_remove = mode(test)

# ╔═╡ c0b38aaa-42ea-11eb-06dc-6f53490d1ad3
md"Soo, could we find mode, then subtract that mean from image, find mode again, ad nauseaum?"

# ╔═╡ d22ada9a-42ea-11eb-380f-677bbeb7cf41
typeof(to_remove)

# ╔═╡ f8a4e1f6-42eb-11eb-0588-cd1265323630
splice!(collect(test), to_remove)

# ╔═╡ 46bf1bae-42ec-11eb-3dd6-3db6770bcf2b
get_outt = findall(x -> x == to_remove, test)

# ╔═╡ 9cabda98-42ec-11eb-3b7c-eda2202fc177
deleteat!(test, findall(x -> x == to_remove, collect(test)))

# ╔═╡ 060bc478-42ed-11eb-348a-9151d0e922f7
splice!(to_remove, test)

# ╔═╡ 59125ee8-42ed-11eb-3542-97ef27972eeb
to_remove

# ╔═╡ f078aa16-42ee-11eb-3a5e-b3e15fd996be
length(to_remove)

# ╔═╡ 30c5609e-42ef-11eb-2192-b55fbd3d9117
to_remove.b

# ╔═╡ 4136c37a-42ef-11eb-1cf3-730605a35461
to_remove.g

# ╔═╡ 446fd160-42ef-11eb-1c4a-6b3ef08da790
to_remove.r

# ╔═╡ 4c9365da-42ef-11eb-1267-bba560e658da
color_type(to_remove)

# ╔═╡ 7c035cf0-42ef-11eb-0e03-3deebb70be20
Gray.(test)

# ╔═╡ 2c5cf682-42f1-11eb-2df6-fd0f011892f3
floaty = float64.(test)

# ╔═╡ 43b0273c-42f1-11eb-1913-8fb68197b140
size(floaty)

# ╔═╡ 6a728abc-42f2-11eb-0ba4-35851b18884f
summary(floaty)

# ╔═╡ aeca591a-42f2-11eb-151d-6f91e731109c
floaty[end-5000:end]

# ╔═╡ dbdc8da6-42f2-11eb-3be6-6f2534c46a1c
mode(floaty)

# ╔═╡ dfcb378c-42f2-11eb-35d7-ef2bf7d5fb5b
b = floaty-mode(floaty)

# ╔═╡ eb5b08fc-42f2-11eb-2eef-e537c759b187
c = b - mode(b)

# ╔═╡ 0557a7c4-42f3-11eb-07c5-8939b379a026
d = c - mode(c)

# ╔═╡ 0bd7bca6-42f3-11eb-0465-31501db99fb6
e = d - mode(d)

# ╔═╡ 128ea79e-42f3-11eb-1ff6-1b3cc98007a5
floaty

# ╔═╡ caaf8810-42f9-11eb-1fdf-39cebeb0171d
function dewitt(img,n)
	counter = 0
	modes = []
	push!(modes, mode(img))

	while counter < n
		this_mode = mode((img-modes[1:end]))
		push!(modes, this_mode)
		counter += 1
	end

end

# ╔═╡ ee0082e6-42fa-11eb-2a7b-a96ad878ce9c
summary(floaty)

# ╔═╡ e663a5da-42f4-11eb-0f0c-e933e4de323e
onedim_floaty = reshape(floaty,(2073600))

# ╔═╡ 23ebde74-42f5-11eb-37ca-c3cbc43b9e98
ndims(onedim_floaty)

# ╔═╡ 311c3904-42f5-11eb-1240-1996737cce47
length(onedim_floaty)

# ╔═╡ d6a715b6-42f4-11eb-324d-5539c2ae9160
md"### Using python to find most common clusters
[https://towardsdatascience.com/finding-most-common-colors-in-python-47ea0767a06a](https://towardsdatascience.com/finding-most-common-colors-in-python-47ea0767a06a)"

# ╔═╡ bdd2d418-42fd-11eb-2227-d79ddf290602
# colorthief = pyimport("colorthief")

# ╔═╡ d7c255ba-42fd-11eb-308b-b1e8a03b4a2a
# color_thief = colorthief.ColorThief("/home/brad/Desktop/plains.jpg")

# ╔═╡ bde8c072-42fe-11eb-2a86-733ed6debe38
cols = color_thief.get_palette(color_count=8)

# ╔═╡ 799c8628-4301-11eb-3f7a-73ac9c544609
function convert_colors_to_julia(x::Array{Tuple{Int64,Int64,Int64}})

	color1 = RGB(cols[1][1]/255, cols[1][2]/255,  cols[1][3]/255)
	color2 = RGB(cols[2][1]/255, cols[2][2]/255,  cols[2][3]/255)
	color3 = RGB(cols[3][1]/255, cols[3][2]/255,  cols[3][3]/255)

	return (color1, color2, color3)
end

# ╔═╡ 983f1744-4301-11eb-3a99-43b1a6e47f37
convert_colors_to_julia(cols)

# ╔═╡ Cell order:
# ╠═b8293916-42e5-11eb-0a82-11d14d1b4ba1
# ╠═d5ea666e-42e5-11eb-2d9b-432c08b84a98
# ╠═e0c451bc-42e5-11eb-1e3e-417e400304ad
# ╠═fd62b854-42e5-11eb-0880-cf34d4d189ae
# ╠═28456c56-42e6-11eb-31f3-dfadfd7bc8a1
# ╟─4fd8f040-42ea-11eb-2548-bdb8e7cd15ec
# ╟─5a3a1c3a-42ea-11eb-134a-514790026194
# ╠═7009f8c4-42e9-11eb-2775-eb9f6cce20a2
# ╠═6b156e22-42e9-11eb-0e22-9d27615d22cc
# ╠═d963bbd6-42e9-11eb-3c7d-3b2722e77bc2
# ╠═37d03b70-42ea-11eb-1a05-ff87bd42f290
# ╠═4f305f60-42eb-11eb-33aa-1548a90f866a
# ╠═6f44c31e-42ea-11eb-3507-cbff33f47bc0
# ╟─c0b38aaa-42ea-11eb-06dc-6f53490d1ad3
# ╠═d22ada9a-42ea-11eb-380f-677bbeb7cf41
# ╠═f8a4e1f6-42eb-11eb-0588-cd1265323630
# ╠═46bf1bae-42ec-11eb-3dd6-3db6770bcf2b
# ╠═9cabda98-42ec-11eb-3b7c-eda2202fc177
# ╠═060bc478-42ed-11eb-348a-9151d0e922f7
# ╠═59125ee8-42ed-11eb-3542-97ef27972eeb
# ╠═f078aa16-42ee-11eb-3a5e-b3e15fd996be
# ╠═30c5609e-42ef-11eb-2192-b55fbd3d9117
# ╠═4136c37a-42ef-11eb-1cf3-730605a35461
# ╠═446fd160-42ef-11eb-1c4a-6b3ef08da790
# ╠═4c9365da-42ef-11eb-1267-bba560e658da
# ╠═7c035cf0-42ef-11eb-0e03-3deebb70be20
# ╠═3aa78d38-42f1-11eb-1005-2bf725559cb3
# ╠═2c5cf682-42f1-11eb-2df6-fd0f011892f3
# ╠═43b0273c-42f1-11eb-1913-8fb68197b140
# ╠═6a728abc-42f2-11eb-0ba4-35851b18884f
# ╠═aeca591a-42f2-11eb-151d-6f91e731109c
# ╠═dbdc8da6-42f2-11eb-3be6-6f2534c46a1c
# ╠═dfcb378c-42f2-11eb-35d7-ef2bf7d5fb5b
# ╠═eb5b08fc-42f2-11eb-2eef-e537c759b187
# ╠═0557a7c4-42f3-11eb-07c5-8939b379a026
# ╠═0bd7bca6-42f3-11eb-0465-31501db99fb6
# ╠═128ea79e-42f3-11eb-1ff6-1b3cc98007a5
# ╠═caaf8810-42f9-11eb-1fdf-39cebeb0171d
# ╠═ee0082e6-42fa-11eb-2a7b-a96ad878ce9c
# ╠═e663a5da-42f4-11eb-0f0c-e933e4de323e
# ╠═23ebde74-42f5-11eb-37ca-c3cbc43b9e98
# ╠═311c3904-42f5-11eb-1240-1996737cce47
# ╟─d6a715b6-42f4-11eb-324d-5539c2ae9160
# ╠═34165c74-42fb-11eb-2d74-512a600c05ba
# ╠═bdd2d418-42fd-11eb-2227-d79ddf290602
# ╠═d7c255ba-42fd-11eb-308b-b1e8a03b4a2a
# ╠═bde8c072-42fe-11eb-2a86-733ed6debe38
# ╠═799c8628-4301-11eb-3f7a-73ac9c544609
# ╠═983f1744-4301-11eb-3a99-43b1a6e47f37


text:
using Resample
using Test
using Dates
using DataFrames
using Statistics

times_long = DateTime(2019):Minute(5):DateTime(2019, 1, 31, 23, 55)
times_short = DateTime(2019):Hour(1):DateTime(2019, 1, 31, 23)

inds_long = eachindex(times_long)
inds_short = 1:12:length(times_long)

power_long = rand(length(times_long))
energy_long = power_long / 12

df = DataFrame(
    Index = inds_long,
    Time = times_long,
    Power = power_long,
    Energy = energy_long,
)

power_short_times = resample(power_long, times_long, times_short)
power_short_inds = resample(power_long, inds_long, inds_short)
power_short_step = resample(power_long, times_long, Hour(1))

energy_short_times = resample(energy_long, times_long, times_short, Sum())
energy_short_inds = resample(energy_long, inds_long, inds_short, Sum())
energy_short_step = resample(energy_long, times_long, Hour(1), Sum())

df_short_times = resample(df, :Time, times_short, [First(), Mean(), Sum()])
df_short_inds = resample(df, :Index, inds_short, [First(), Mean(), Sum()])
df_short_step = resample(df, :Time, Hour(1), [First(), Mean(), Sum()])

@testset "Resample" begin
    @test power_short_times == power_short_inds == power_short_step
    @test energy_short_times == energy_short_inds == energy_short_step
    @test df_short_times == df_short_inds == df_short_step

    @test mean(power_long) ≈ mean(power_short_times)
    @test sum(energy_long) ≈ sum(energy_short_times)
end


text:
const SUITS = ['♣', '♢', '♠', '♡']
const RANKS = ['A', Char.(50:57)..., 'T', 'J', 'Q', 'K']

const RANK_MASK = Dict(c => i for (i, c) in enumerate(RANKS))
const RANK_VALUE = Dict(c => i > 10 ? 10 : i for (i, c) in enumerate(RANKS))

struct Rank
    rank::Char
    value::Int64
    mask::Int64

    function Rank(rank::Char; ace_high::Bool=false)
        rank in RANKS || throw(ArgumentError("Incorrect rank: $rank"))
        value = RANK_VALUE[rank]
        mask = RANK_MASK[rank]
        if ace_high
            rank == 'A' || throw(ArgumentError("Incorrect rank: $(rank(x))"))
            mask = 14
        end
        new(rank, value, mask)
    end
end

rank(x::Rank) = x.rank
value(x::Rank) = x.value
mask(x::Rank) = x.mask

Base.string(x::Rank) = string(rank(x))
Base.show(io::IO, x::Rank) = print(io, string(x))

ace_high(x::Rank) = Rank(rank(x), ace_high=true)

struct Suit
    suit::Char

    function Suit(suit::Char)
        suit in SUITS || throw(ArgumentError("Incorrect suit: $suit"))
        new(suit)
    end
end

suit(x::Suit) = x.suit
Base.string(x::Suit) = string(suit(x))
Base.show(io::IO, x::Suit) = print(io, string(x))

struct Card
    rank::Rank
    suit::Suit
end

function Card(rank::Char, suit::Char; ace_high::Bool=false)
    Card(Rank(rank; ace_high=ace_high), Suit(suit))
end

function Card(str::AbstractString; ace_high::Bool=false)
    length(str) == 2 || throw(ArgumentError("`Card` string must be of length 2"))
    Card(str[1], str[2]; ace_high=ace_high)
end

rank(x::Card) = x.rank
value(x::Card) = value(rank(x))
mask(x::Card) = mask(rank(x))
suit(x::Card) = x.suit

ace_high(x::Card) = Card(ace_high(rank(x)), suit(x))

Base.isless(x::Card, y::Card) = mask(x) < mask(y)

Base.string(x::Card) = string(rank(x))*string(suit(x))
Base.show(io::IO, x::Card) = print(io, string(x))

# Methods to create `Card` instances by juxtaposing `Rank` and `Suit`, or `Int` and `Suit`

function Base.:*(n::Int, s::Suit)
    2 ≤ n ≤ 9 || error("")
    return Card(Rank(Char(48+n)), s)
end

# Exports ♡, ♠, ♢, ♣,
for s in SUITS
    name = Symbol("$s")
    @eval const $name = Suit($s)
end

# Exports A, T, J, Q, K,
for r in ['A', 'T', 'J', 'Q', 'K']
    name = Symbol("$r")
    @eval const $name = Rank($r)
end

# Exports
# A♡, A♠, A♢, A♣,
# T♡, T♠, T♢, T♣,
# J♡, J♠, J♢, J♣,
# Q♡, Q♠, Q♢, Q♣,
# K♡, K♠, K♢, K♣,
for s in SUITS, r in ['A', 'T', 'J', 'Q', 'K']
    name = Symbol("$(r)$(s)")
    @eval const $name = Card(Rank($r), Suit($s))
end

function deck()
    xs = sizehint!(Card[], 52)
    for s in SUITS, r in RANKS
        push!(xs, Card(r, s))
    end
    return xs
end

const DECK = deck()


text:
using Photometry.Aperture: bounds,
                           oblique_coefficients

@testset "Apertures" begin
    ap_ellipse = EllipticalAperture(0, 0, 20, 10, 0)
    @test bounds(ap_ellipse) == (-20, 20, -10, 10)
    @test size(ap_ellipse) == (21, 41)
    @test size(ap_ellipse, 1) == 21
    @test all(axes(ap_ellipse) .== (-10:10, -20:20))
    @test eachindex(ap_ellipse) == CartesianIndex(-10, -20):CartesianIndex(10, 20)

    @test EllipticalAperture([0, 0], 20, 10, 0) == ap_ellipse

    ap_ellipse = EllipticalAperture(0, 0, 2, 1, 45)
    @test bounds(ap_ellipse) == (-1, 1, -1, 1)
end

@testset "Elliptical Aperture" begin
    e = EllipticalAperture(0, 0, 20, 10, 0)
    @test sprint(show, e) == "EllipticalAperture(0, 0, a=20, b=10, θ=0°)"
end

@testset "oblique_coefficients" begin
    @test all(oblique_coefficients(2, 2, 0) .≈ (0.25, 0.25, 0.0))
    @test all(oblique_coefficients(2, 2, 90) .≈ (0.25, 0.25, 0.0))
    @test all(oblique_coefficients(2, 1, 30) .≈ (7 / 16, 13 / 16, -6sqrt(3) / 16))
end

@testset "Elliptical Annulus" begin
    e0 = EllipticalAnnulus(0, 0, 8, 16, 4, 45)
    @test sprint(show, e0) == "EllipticalAnnulus(0.0, 0.0, a_in=8.0, a_out=16.0, b_in=2.0, b_out=4.0, θ=45.0°)"
    @test EllipticalAnnulus([0, 0], 8, 16, 4, 45) == EllipticalAnnulus(0, 0, 8, 16, 4, 45)
end

@testset "Elliptical Annulus bounding box" begin
    e = EllipticalAnnulus(0, 0, 8, 16, 4, 0)

    @test bounds(e) == (-16, 16, -4, 4)
    @test size(e) == (9, 33)
    @test all(axes(e) .== (-4:4, -16:16))
    @test eachindex(e) == CartesianIndex(-4, -16):CartesianIndex(4, 16)
end

@testset "regression circular ellipse" begin
    # some weird bug where centered on-grid past a certain size (and rotated) would fail
    e = EllipticalAperture(5.5, 5.5, 4, 4, 20)
    @test sum(e) ≈ photometry(e, ones(9, 9)).aperture_sum # just a test that no errors occur
end


text:
import LinearAlgebra
import QuantumPropagators
import Base: -

@doc raw"""Extended generator for the standard dynamic gradient.

```julia
G̃ = TimeDependentGradGenerator(G)
```

contains the original time-dependent generator `G` (a Hamiltonian or
Liouvillian) in `G̃.G`, a vector of control derivatives ``∂G/∂ϵₗ(t)`` in
`G̃.control_derivs`, and the controls in `G̃.controls`.

For a generator ``G = Ĥ(t) = Ĥ₀ + ϵ₁(t) Ĥ₁ + … +  ϵₙ(t) Ĥₙ``, this extended
generator encodes the block-matrix

```math
G̃ = \begin{pmatrix}
         Ĥ(t)  &  0    &  \dots   &  0     &  Ĥ₁     \\
         0     &  Ĥ(t) &  \dots   &  0     &  Ĥ₂     \\
    \vdots     &       &  \ddots  &        &  \vdots \\
         0     &  0    &  \dots   &  Ĥ(t)  &  Ĥₙ     \\
         0     &  0    &  \dots   &  0     &  Ĥ(t)
\end{pmatrix}
```

Note that the ``∂G/∂ϵₗ(t)`` (``Ĥₗ`` in the above example) are functions, to
account for the possibility of non-linear control terms, see
[`getcontrolderiv`](@ref).
"""
struct TimeDependentGradGenerator{GT, CT}
    G :: GT
    control_derivs :: Vector{Function}
    controls :: Vector{CT}

    function TimeDependentGradGenerator(G::GT) where GT
        controls = collect(getcontrols(G))
        control_derivs = getcontrolderivs(G, controls)
        new{GT, eltype(controls)}(G, control_derivs, controls)
    end

end

"""Static generator for the standard dynamic gradient.

```julia
G::GradGenerator = evalcontrols(G_of_t::TimeDependentGradGenerator, vals_dict)
```

is the result of plugging in specific values for all controls in a
[`TimeDependentGradGenerator`](@ref). See [`evalcontrols`](@ref) and
[`evalcontrols!`](@ref).

The resulting object can be multiplied directly with a [`GradVector`](@ref),
e.g., in the process of evaluating a piecewise-constant time propagation.
"""
struct GradGenerator{GT, CGT}
    G :: GT
    control_derivs :: Vector{CGT}

    function GradGenerator(G_of_t::TimeDependentGradGenerator)
        dummy_vals = IdDict(control => 1.0  for control in G_of_t.controls)
        G = evalcontrols(G_of_t.G, dummy_vals)
        control_derivs = [μ(1.0) for μ in G_of_t.control_derivs]
        new{typeof(G), eltype(control_derivs)}(G, control_derivs)
    end
end

function evalcontrols!(G::GradGenerator, G_of_t::TimeDependentGradGenerator, vals_dict::AbstractDict)
    evalcontrols!(G.G, G_of_t.G, vals_dict)
    for (i, control) in enumerate(G_of_t.controls)
        G.control_derivs[i] = G_of_t.control_derivs[i](vals_dict[control])
        # In most cases (for linear controls), the above line will be a no-op.
        # Hence, we're not using `copyto!`.
    end
    return G
end

function evalcontrols(G_of_t::TimeDependentGradGenerator, vals_dict::AbstractDict)
    G = GradGenerator(G_of_t)
    evalcontrols!(G, G_of_t, vals_dict)
end

@doc raw"""Extended state-vector for the dynamic gradient.

```julia
Ψ̃ = GradVector(Ψ, num_controls)
```

for an initial state `Ψ` and `num_controls` control fields.

The `GradVector` conceptually corresponds to a direct-sum (block) column-vector
``Ψ̃ = (|Ψ̃₁⟩, |Ψ̃₂⟩, … |Ψ̃ₙ⟩, |Ψ⟩)^T``, where ``n`` is `num_controls`. With a
matching ``G̃`` as in the documentation of [`TimeDependentGradGenerator`](@ref),
we have

```math
G̃ Ψ̃ = \begin{pmatrix}
Ĥ |Ψ̃₁⟩ + Ĥ₁|Ψ⟩ \\
\vdots \\
Ĥ |Ψ̃ₙ⟩ + Ĥₙ|Ψ⟩ \\
Ĥ |Ψ⟩
\end{pmatrix}
```

and

```math
e^{-i G̃ dt} \begin{pmatrix} 0 \\ \vdots \\ 0 \\ |Ψ⟩ \end{pmatrix}
= \begin{pmatrix}
\frac{∂}{∂ϵ₁} e^{-i Ĥ dt} |Ψ⟩ \\
\vdots \\
\frac{∂}{∂ϵₙ} e^{-i Ĥ dt} |Ψ⟩ \\
e^{-i Ĥ dt} |Ψ⟩
\end{pmatrix}.
```
"""
struct GradVector{T}
    state :: T
    grad_states :: Vector{T}
end

function GradVector(Ψ::T, num_controls::Int64) where T
    grad_states = [similar(Ψ) for _ in 1:num_controls]
    for i = 1 : num_controls
        fill!(grad_states[i], 0.0)
    end
    GradVector{T}(copy(Ψ), grad_states)
end

"""Reset the given gradient vector for a new gradient evaluation.

```julia
resetgradvec!(Ψ̃::GradVector)
```

zeroes out `Ψ̃.grad_states` but leaves `Ψ̃.state` unaffected.

```julia
resetgradvec!(Ψ̃::GradVector, Ψ)
```

additionally sets `Ψ̃.state` to `Ψ`.
"""
function resetgradvec!(Ψ̃::GradVector{T}) where T
    for i = 1 : length(Ψ̃.grad_states)
        fill!(Ψ̃.grad_states[i], 0.0)
    end
end

function resetgradvec!(Ψ̃::GradVector{T}, Ψ::T) where T
    copyto!(Ψ̃.state, Ψ)
    resetgradvec!(Ψ̃)
end

function LinearAlgebra.mul!(Φ::GradVector, G::GradGenerator, Ψ::GradVector)
    LinearAlgebra.mul!(Φ.state, G.G, Ψ.state)
    for i = 1:length(Ψ.grad_states)
        LinearAlgebra.mul!(Φ.grad_states[i], G.G, Ψ.grad_states[i])
        LinearAlgebra.mul!(Φ.grad_states[i], G.control_derivs[i], Ψ.state, 1, 1)
    end
end

function LinearAlgebra.lmul!(c, Ψ::GradVector)
    LinearAlgebra.lmul!(c, Ψ.state)
    for i = 1:length(Ψ.grad_states)
        LinearAlgebra.lmul!(c, Ψ.grad_states[i])
    end
end

function LinearAlgebra.axpy!(a, X::GradVector, Y::GradVector)
    LinearAlgebra.axpy!(a, X.state, Y.state)
    for i = 1:length(X.grad_states)
        LinearAlgebra.axpy!(a, X.grad_states[i], Y.grad_states[i])
    end
end

function LinearAlgebra.norm(Ψ::GradVector)
    nrm = LinearAlgebra.norm(Ψ.state)
    for i = 1:length(Ψ.grad_states)
        nrm += LinearAlgebra.norm(Ψ.grad_states[i])
    end
    return nrm
end

function LinearAlgebra.dot(Ψ::GradVector, Φ::GradVector)
    c::ComplexF64 = LinearAlgebra.dot(Ψ.state, Φ.state)
    for i = 1:length(Ψ.grad_states)
        c += LinearAlgebra.dot(Ψ.grad_states[i], Φ.grad_states[i])
    end
    return c
end

function Base.copyto!(dest::GradVector, src::GradVector)
    copyto!(dest.state, src.state)
    for i = 1:length(src.grad_states)
        copyto!(dest.grad_states[i], src.grad_states[i])
    end
    return dest
end

function Base.copy(Ψ::GradVector)
    Φ = GradVector(Ψ.state, length(Ψ.grad_states))
    for i = 1:length(Ψ.grad_states)
        copyto!(Φ.grad_states[i], Ψ.grad_states[i])
    end
    return Φ
end

function Base.length(Ψ::GradVector)
    return length(Ψ.state) * (1 + length(Ψ.grad_states))
end

function Base.similar(Ψ::GradVector)
    return GradVector(similar(Ψ.state), [similar(ϕ) for ϕ ∈ Ψ.grad_states])
end

function Base.fill!(Ψ::GradVector, v)
    Base.fill!(Ψ.state, v)
    for i = 1:length(Ψ.grad_states)
        Base.fill!(Ψ.grad_states[i], v)
    end
end

function -(Ψ::GradVector, Φ::GradVector)
    res = copy(Ψ)
    LinearAlgebra.axpy!(-1, Φ.state, res.state)
    for i = 1:length(Ψ.grad_states)
        LinearAlgebra.axpy!(-1, Φ.grad_states[i], res.grad_states[i])
    end
    return res
end

struct DenseGradExpPropWrk{T}
    Ψ_full :: Vector{ComplexF64}
    Φ_full :: Vector{ComplexF64}
    G_full :: Matrix{T}
    function DenseGradExpPropWrk(
            Ψ̃::GradVector{Vector{ComplexF64}},
            G̃::GradGenerator{Matrix{T}, Matrix{T}}) where T
        Ψ_full = convert_gradvec_to_dense(Ψ̃)
        Φ_full = similar(Ψ_full)
        G_full = convert_gradgen_to_dense(G̃)
        new{T}(Ψ_full, Φ_full, G_full)
    end
end

function QuantumPropagators.initpropwrk(
        state::GradVector{Vector{ComplexF64}},
        tlist,
        method::Val{:expprop},
        generator::Vararg{GradGenerator{Matrix{T}, Matrix{T}}};
        kwargs...) where T
    return DenseGradExpPropWrk(state, generator[1])
end

@inline function convert_gradgen_to_dense(G)
    N = size(G.G)[1]
    L = length(G.control_derivs)
    G_full = zeros(eltype(G.G), N*(L+1), N*(L+1))
    convert_gradgen_to_dense!(G_full, G)
end

@inline function convert_gradgen_to_dense!(G_full, G)
    N = size(G.G)[1]
    L = length(G.control_derivs)
    @inbounds for i = 1:L+1
        G_full[(i-1)*N + 1 : i*N, (i-1)*N + 1 : i*N] .= G.G
    end
    # Set the control-derivatives in the last (block-)column
    @inbounds for i = 1:L
        G_full[(i-1)*N + 1 : i*N, L*N + 1 : (L+1)*N]  .= G.control_derivs[i]
    end
    return G_full
end

@inline function convert_gradvec_to_dense(Ψ)
    N = length(Ψ.state)
    L = length(Ψ.grad_states)
    Ψ_full = zeros(ComplexF64, N*(L+1))
    convert_gradvec_to_dense!(Ψ_full, Ψ)
end

@inline function convert_gradvec_to_dense!(Ψ_full, Ψ)
    N = length(Ψ.state)
    L = length(Ψ.grad_states)
    @inbounds for i = 1:L
        Ψ_full[(i-1)*N + 1 : i*N] .= Ψ.grad_states[i]
    end
    @inbounds Ψ_full[L*N + 1: (L+1)*N] .= Ψ.state
    return Ψ_full
end

@inline function convert_dense_to_gradvec!(Ψ, Ψ_full)
    N = length(Ψ.state)
    L = length(Ψ.grad_states)
    @inbounds for i = 1:L
        Ψ.grad_states[i] .= Ψ_full[(i-1)*N + 1 : i*N]
    end
    @inbounds Ψ.state .= Ψ_full[L*N + 1: (L+1)*N]
    return Ψ
end

function QuantumPropagators.propstep!(
        Ψ̃::GradVector{Vector{ComplexF64}},
        G̃::GradGenerator{Matrix{T}, Matrix{T}},
        dt::Float64,
        wrk::DenseGradExpPropWrk;
        kwargs...) where T
    func = get(kwargs, :func, H_dt -> exp(-1im * H_dt))
    convert_gradgen_to_dense!(wrk.G_full, G̃)
    convert_gradvec_to_dense!(wrk.Ψ_full, Ψ̃)
    U = func(wrk.G_full*dt)
    LinearAlgebra.mul!(wrk.Φ_full, U, wrk.Ψ_full)
    convert_dense_to_gradvec!(Ψ̃, wrk.Φ_full)
    return Ψ̃
end


text:
DEFAULT_HEADERS = ["ALLOW" => "POST, OPTIONS",
                   "accept" => "application/json",
                   "Content-Type" => "application/json",
                   "Content-Encoding" => "gzip",
                   "Access-Control-Allow-Origin" => "*",
                   "Access-Control-Allow-Methods" => "POST, OPTIONS",
                   "Access-Control-Allow-Headers" => "*",
                   "Accept-Encoding" => "text/plain", "gzip"]

function handle_http_request_method(request::HTTP.Request)
    request_method = getfield(request, :method)
    if request_method == "OPTIONS"
        return HTTP.Response(200, DEFAULT_HEADERS)
    elseif request_method != "POST"
        return HTTP.Response(400, DEFAULT_HEADERS)
    end
    return nothing
end

function check_content_type(request::HTTP.Request)
    content_type = HTTP.header(request, "Content-Type")
    if content_type != "application/json"
        throw(ArgumentError("Invalid Content-Type '$(content_type)'"))
    end
    return nothing
end

function parse_content_body(request::HTTP.Request)
    content_encoding = HTTP.header(request, "Content-Encoding")
    request_body = nothing
    if content_encoding == "gzip"
        try
            request_body = String(transcode(GzipDecompressor, request.body))
        catch e
            throw(ArgumentError("Cannot decompress with gzip."))
        end
    elseif content_encoding == ""
        request_body = String(request.body)
    else
        throw(ArgumentError("Invalid Content-Encoding '$(content_encoding)'"))
    end
    try
        return JSON.parse(request_body)
    catch e
        throw(ArgumentError("Cannot parse JSON."))
    end
    return nothing
end

function build_response(status, params)
    return HTTP.Response(status, DEFAULT_HEADERS;
                         body=transcode(GzipCompressor, JSON.json(Dict(:status => status, params...))))
end


text:
"""
A `WeightedDict` represents some function f(x) times an existing dictionary.
"""
struct WeightedDict{S,T} <: DerivedDict{S,T}
    superdict   ::  Dictionary{S}
    weightfun
end

WeightedDict(superdict::Dictionary{S,T},weightfun) where {S,T} = WeightedDict{S,T}(superdict,weightfun)

WeightedDict(superdict::Dictionary{S},weightfun,T) where S = WeightedDict{S,T}(superdict,weightfun)

const WeightedDict1d{S <: Number,T} = WeightedDict{S,T}
const WeightedDict2d{S <: Number,T} = WeightedDict{SVector{2,S},T}
const WeightedDict3d{S <: Number,T} = WeightedDict{SVector{3,S},T}
const WeightedDict4d{S <: Number,T} = WeightedDict{SVector{4,S},T}

==(dict1::WeightedDict,dict2::WeightedDict) =
    superdict(dict1)==superdict(dict2) && weightfunction(dict1)==weightfunction(dict2)

weightfunction(dict::WeightedDict) = dict.weightfun

similardictionary(dict1::WeightedDict, dict2::Dictionary) = WeightedDict(dict2, weightfunction(dict1))

name(dict::WeightedDict) = "Weighted " * name(superdict(dict))

isreal(dict::WeightedDict) = _isreal(dict, superdict(dict), weightfunction(dict))
_isreal(dict::WeightedDict, superdict, fun::Function) = isreal(superdict)

hasderivative(dict::WeightedDict) = hasderivative(superdict(dict))
isorthonormal(dict::WeightedDict) = false
isorthogonal(dict::WeightedDict) = false
# We can not compute antiderivatives in general.
hasantiderivative(dict::WeightedDict) = false

hasmeasure(dict::WeightedDict) = false

mapped_dict(d::WeightedDict, map) = MappedDict(d, map)

# We have to distinguish between 1d and higher-dimensional grids, since we
# have to splat the arguments to the weightfunction
eval_weight_on_grid(w, grid::AbstractGrid1d) = [w(x) for x in grid]

function eval_weight_on_grid(w, grid::AbstractGrid)
    # Perhaps the implementation here could be simpler, but [w(x...) for x in grid]
    # does not seem to respect the size of the grid, only its length
    a = zeros(prectype(grid), size(grid))
    for i in eachindex(grid)
        a[i] = w(grid[i]...)
    end
    a
end

# Evaluating basis functions: we multiply by the function of the dict
unsafe_eval_element(dict::WeightedDict, idx, x) = _unsafe_eval_element(dict, weightfunction(dict), idx, x)
_unsafe_eval_element(dict::WeightedDict1d, w, idx, x) = w(x) * unsafe_eval_element(superdict(dict), idx, x)
_unsafe_eval_element(dict::WeightedDict, w, idx, x) = w(x...) * unsafe_eval_element(superdict(dict), idx, x)

# Evaluate the derivative of 1d weighted sets
function unsafe_eval_element_derivative(dict::WeightedDict1d, idx, x, order)
	@assert order == 1
    diff(weightfunction(dict), x) * unsafe_eval_element(superdict(dict), idx, x) +
		weightfunction(dict)(x) * unsafe_eval_element_derivative(superdict(dict), idx, x, order)
end

# Special case for log, since it is not precise at 0.
diff(::typeof(log)) = x->1/x
diff(::typeof(log), x::T) where T<:Number= convert(T,1)/x

# Evaluate an expansion: same story
eval_expansion(dict::WeightedDict, coefficients, x) = _eval_expansion(dict, weightfunction(dict), coefficients, x)
# temporary, to remove an ambiguity
eval_expansion(dict::WeightedDict, coefficients, x::AbstractGrid) = _eval_expansion(dict, weightfunction(dict), coefficients, x)

_eval_expansion(dict::WeightedDict1d, w, coefficients, x::Number) = w(x) * eval_expansion(superdict(dict), coefficients, x)
_eval_expansion(dict::WeightedDict, w, coefficients, x) = w(x...) * eval_expansion(superdict(dict), coefficients, x)

_eval_expansion(dict::WeightedDict, w, coefficients, grid::AbstractGrid) =
    eval_weight_on_grid(w, grid) .* eval_expansion(superdict(dict), coefficients, grid)

# You can create a WeightedDict by multiplying a function with a dict, using
# left multiplication:
(*)(f::Function, dict::Dictionary) = WeightedDict(dict, f)

weightfun_scaling_operator(::Type{T}, gb::GridBasis, weightfunction) where {T} =
	weightfun_scaling_operator(T, gb, weightfunction, grid(gb))

weightfun_scaling_operator(::Type{T}, gb, weightfunction, grid::AbstractGrid1d) where {T} =
    DiagonalOperator(gb, gb, map(T,map(weightfunction, grid)))

function weightfun_scaling_operator(::Type{T}, gb, weightfunction, grid::AbstractGrid) where {T}
	A = map(T,map(x->weightfunction(x...), grid))
    DiagonalOperator(gb, gb, A[:])
end

function weightfun_scaling_operator(::Type{T}, gb, weightfunction, grid::ProductGrid) where {T}
	A = map(T, map(x->weightfunction(x...), grid))
    DiagonalOperator(gb, gb, A[:])
end

transform_to_grid(T, src::WeightedDict, dest::GridBasis, grid; options...) =
    weightfun_scaling_operator(T, dest, weightfunction(src)) * wrap_operator(src, dest, transform_to_grid(T, superdict(src), dest, grid; options...))

transform_from_grid(T, src::GridBasis, dest::WeightedDict, grid; options...) =
	wrap_operator(src, dest, transform_from_grid(T, src, superdict(dest), grid; options...)) * inv(weightfun_scaling_operator(T, src, weightfunction(dest)))

function derivative_dict(src::WeightedDict, order; options...)
	if order == 0
		src
	elseif order == 1
	    s = superdict(src)
	    f = weightfunction(src)
	    f_prime = diff(f)
	    s_prime = derivative_dict(s, order)
	    (f_prime * s) ⊕ (f * s_prime)
    else
        @assert order>1
        s = superdict(src)
        f = weightfunction(src)
        dfs = Any[]
        push!(dfs,f)
        for k in 1:order
            push!(dfs,diff(dfs[k]))
        end
        ⊕([dfs[k+1]*(derivative_dict(s, order-k))  for k in order:-1:0]...)
    end
end

# Assume order = 1...
function differentiation(::Type{T}, src::WeightedDict, dest::MultiDict, order; options...) where {T}
    @assert size(dest) == size(derivative_dict(src, order))
    if order==0
        return IdentityOperator{T}(src, component(dest, 1))
    end
    if order == 1
        I = IdentityOperator{T}(src, component(dest, 1))
        D = differentiation(T, superdict(src), order)
        DW = wrap_operator(src, component(dest, 2), D)
        block_column_operator([I,DW])
    elseif order > 1
        Ds = [wrap_operator(src,e,binomial(order,k)*differentiation(T, superdict(src), k))
            for (e,k) in zip(components(dest),0:order)]
        block_column_operator(Ds)
    else
        error("differentiation of order $order not implemented")
    end
end

function evaluation(::Type{T}, dict::WeightedDict, gb::GridBasis, grid::AbstractGrid; options...) where {T}
    A = evaluation(T, superdict(dict), gb, grid; options...)
    D = weightfun_scaling_operator(T, gb, weightfunction(dict))
    D * wrap_operator(dict, gb, A)
end

## Printing

Display.object_parentheses(d::WeightedDict) = true
Display.stencil_parentheses(d::WeightedDict) = true
Display.displaystencil(d::WeightedDict) = _stencil(d, superdict(d), weightfunction(d))
_stencil(d::WeightedDict, dict, weight) = [weight, " * ", dict]


text:
include("/Users/nirmitdeshpande/.julia/dev/impedance/src/plotting.jl")

plotNyquist()

max_stars_repo_name: nirmitdeshpande/EISModeling

max_issues_repo_name: nirmitdeshpande/EISModeling

max_forks_repo_name: nirmitdeshpande/EISModeling


text:
import DistributedBoosting

function main()
	if(ARGS[1] == "provide")
		if(length(ARGS) == 3)
			DistributedBoosting.provide(join(ARGS[2]),join(ARGS[3]))
		else
			DistributedBoosting.provide(join(ARGS[2]),join(ARGS[3]),join(ARGS[4]))
		end
	end
	if(ARGS[1] == "init")
		DistributedBoosting.boosting(parse(Int,ARGS[2]), parse(Int,ARGS[4]), parse(Int,ARGS[5]), join(ARGS[6]))
	end
	if(ARGS[1] == "boost")
		inputarray = Array{ASCIIString, 1}(length(ARGS)-3)
		for i = 1 : length(ARGS)-3
			inputarray[i] = join(ARGS[i+3])
		end
		DistributedBoosting.boosting(join(ARGS[2]), join(ARGS[3]), inputarray)
	end
end

max_stars_repo_name: danielazoeller/DistributedBoosting

max_issues_repo_name: danielazoeller/DistributedBoosting

max_forks_repo_name: danielazoeller/DistributedBoosting


text:
mutable struct PipsSolver <: AbstractPlasmoSolver
    options::Dict{Symbol,Any}
    partition_options::Dict{Symbol,Any}
    manager::MPI.MPIManager
    status
end

function PipsSolver(;n_workers = 1,master = 0,children = nothing,partitions = Vector{Vector{Int64}}(),master_partition = Vector{Int64}(),lift_link_constraints = false)
    solver = PipsSolver(
    Dict(:n_workers => n_workers,:master => master,:children => children,:lift_link_constraints => lift_link_constraints),
    Dict(:master_partition => master_partition,:sub_partitions => partitions),
    MPI.MPIManager(np = n_workers),
    nothing)
end

function Plasmo.solve(graph::ModelGraph,solver::PipsSolver)
    #TODO Might make more sense to directly look up the pips library environment variable
    !isempty(Libdl.find_library("libparpipsnlp")) || error("Could not find a PIPS-NLP installation")

    #If we partition
    #NOTE Might make more sense to move the creation of the PIPS structure to pipsnlpsolve
    if !isempty(solver.partition_options[:sub_partitions])
        println("Using partitions for PIPS-NLP")
        pips_graph = create_pips_tree(graph,solver.partition_options[:sub_partitions];master_partition = solver.partition_options[:master_partition], lift_link_constraints =
        solver.options[:lift_link_constraints])
        master = pips_graph.master_node_index
        children = pips_graph.sub_node_indices
    else #just use master and child indices from original graph
        println("Using graph structure for PIPS-NLP")
        master = solver.options[:master]
        children = solver.options[:children]
        pips_graph = graph
    end

    manager = solver.manager
    if length(manager.mpi2j) == 0
        addprocs(manager)
    end

    #TODO better structure checks
    #@assert length(children) + 1 == length(getnodes(graph))

    println("Preparing PIPS MPI environment")
    eval(quote @everywhere using Plasmo end)
    eval(quote @everywhere using PlasmoSolverInterface end)
    #eval(quote @everywhere using PlasmoSolverInterface.PlasmoPipsNlpInterface end)

    println("Finished loading Plasmo on workers")

    send_pips_data(manager,pips_graph,master,children)

    println("Solving with PIPS-NLP")
    #NOTE Might need to do status = pipsnlp_solve
    MPI.@mpi_do manager pipsnlp_solve(pips_graph,master,children)

    #Get solution
    rank_zero = manager.mpi2j[0]
    sol = fetch(@spawnat(rank_zero, getfield(Main, :pips_graph)))

    #Update the graph on the julia process if we used a PipsTree
    setsolution(sol,pips_graph)

    #Now move the pips_graph solution to the original model graph
    setsolution(pips_graph,graph)

    return nothing  #TODO retrieve solve status
end

function send_pips_data(manager::MPI.MPIManager,graph::AbstractModelGraph,master::Int,children::Vector{Int})
    julia_workers = collect(values(manager.mpi2j))
    r = RemoteChannel(1)
    @spawnat(1, put!(r, [graph,master,children]))
    @sync for to in julia_workers
        @spawnat(to, Core.eval(Main, Expr(:(=), :pips_graph, fetch(r)[1])))
        @spawnat(to, Core.eval(Main, Expr(:(=), :master, fetch(r)[2])))
        @spawnat(to, Core.eval(Main, Expr(:(=), :children, fetch(r)[3])))
    end
end

# #load PIPS-NLP if the library can be found
# function load_pips()
#     if  !isempty(Libdl.find_library("libparpipsnlp"))
#         #include("solver_interfaces/plasmoPipsNlpInterface3.jl")
#         #eval(quote using .PlasmoPipsNlpInterface3 end)
#         eval(macroexpand(quote @everywhere using .PlasmoPipsNlpInterface3 end))
#     else
#         pipsnlp_solve(Any...) = throw(error("Could not find a PIPS-NLP installation"))
#     end
# end

#load DSP if the library can be found
# function load_dsp()
#     if !isempty(Libdl.find_library("libDsp"))
#         include("solver_interfaces/plasmoDspInterface.jl")
#         using .PlasmoDspInterface
#     else
#         dsp_solve(Any...) = throw(error("Could not find a DSP installation"))
#     end
# end


text:
load("prettyshow/prettyshow.jl")

ex = :( function pprint(io::PrettyIO, ex::Expr)
    if contains([:(=), :(.), doublecolon], head) && nargs==2
        pprint(subblock(io), args[1], string(head), args[2])
    else
        pprint(io, ex.head, "(")
        let io=subblock(io)
            for (arg, k) in enumerate(ex.args)
                pprint(io, arg)
                if k < length(ex.args)
                    pprint(io, ", ")
                end
            end
        end
        return pprint(io, ")")
    end

    try
        a
    end
    try
        a
    catch x
        b
    end

    let
        x+=1
    end
    let x=1
        x+=1
    end
    let x=1,y=2
        x+=1
    end

    begin
        a=1
        b=2
    end

    abstract R<:S
    type T
        x::Int
    end
    typealias P T

    type Q{T}
        q::Q
    end
    while 1 < x < 5
        x += 1
    end
    local x::Int
    local x,y,z
    global y=4.5
    const z=3
    f=x->x^2
end )

pprint(ex)

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
export ifigure, addcanvas!, addaxes!, addcursor!
export HeatmapCanvas, ScatterCanvas, LineCanvas

const curviz = Ref{Union{Viz,Nothing}}(nothing)

const keymap = Dict(
  :pan_left => Keyboard.left,
  :pan_right => Keyboard.right,
  :pan_up => Keyboard.up,
  :pan_down => Keyboard.down,
  :zoom_x_out => Keyboard.left_bracket,
  :zoom_x_in => Keyboard.right_bracket,
  :zoom_y_out => Keyboard.minus,
  :zoom_y_in => Keyboard.equal,
  :brightness_down => Keyboard.comma,
  :brightness_up => Keyboard.period,
  :contrast_up => Keyboard.semicolon,
  :contrast_down => Keyboard.apostrophe,
  :reset => Keyboard._0
)

function ifigure(; width=1024, height=768, hold=false, show=true)
  hold && curviz[] !== nothing && return curviz[]
  curviz[] = Viz(
    scene = Scene(resolution=(width,height), camera=campixel!),
    ijrect = Node(ℛ(0, 0, width, height)),
    selrect = Node(ℛ(0, 0, width, height)),
    children = Canvas[]
  )
  show && display(curviz[].scene)
  on(curviz[].scene.events.window_area) do win
    curviz[].ijrect[] = ℛ(win.origin[1], win.origin[2], win.widths[1], win.widths[2])
  end
  on(curviz[].scene.events.window_open) do b
    b || (curviz[] = nothing)
  end
  on(select_rectangle(curviz[].scene)) do r
    curviz[].selrect[] = ℛ(
      round(Int, r.origin[1]),
      round(Int, r.origin[2]),
      round(Int, r.widths[1]),
      round(Int, r.widths[2]))
  end
  curviz[]
end

function ij2xy(i, j, c::Canvas)
  ijrect = c.ijhome[]
  xyrect = c.xyrect[]
  x = left(xyrect) + (i - left(ijrect)) * width(xyrect) / width(ijrect)
  y = bottom(xyrect) + (j - bottom(ijrect)) * height(xyrect) / height(ijrect)
  (x, y)
end

function ij2xy(p::Point2f0, c::Canvas)
  x, y = ij2xy(p[1], p[2], c)
  Point2f0(x, y)
end

function ij2xy(p::ℛ, c::Canvas)
  x1, y1 = ij2xy(left(p), bottom(p), c)
  x2, y2 = ij2xy(right(p), top(p), c)
  ℛ(x1, y1, x2-x1, y2-y1)
end

function xy2ij(x, y, c::Canvas)
  ijrect = c.ijhome[]
  xyrect = c.xyrect[]
  i = left(ijrect) + (x - left(xyrect)) * width(ijrect) / width(xyrect)
  j = bottom(ijrect) + (y - bottom(xyrect)) * height(ijrect) / height(xyrect)
  (i, j)
end

function xy2ij(p::Point2f0, c::Canvas)
  i, j = xy2ij(p[1], p[2], c)
  Point2f0(i, j)
end

function xy2ij(dtype::DataType, x, y, c::Canvas)
  i, j = xy2ij(x, y, c)
  (round(dtype, i), round(dtype, j))
end

function applyconstraints!(c::Canvas)
  # TODO: apply constraints xmin, xmax, ymin, ymax, cmin, cmax
  # TODO: apply constraints minwidth, maxwidth, minheight, maxheight
end

function realigncanvas!(c::Canvas)
  if c.ijhome[] != c.ijrect[]
    xyrect = c.xyrect[]
    xscale = width(xyrect) / width(c.ijrect)
    yscale = height(xyrect) / height(c.ijrect)
    c.xyrect[] = ℛ(
      left(xyrect) + xscale*(left(c.ijhome) - left(c.ijrect)),
      bottom(xyrect) + yscale*(bottom(c.ijhome) - bottom(c.ijrect)),
      xscale * width(c.ijhome),
      yscale * height(c.ijhome)
    )
    c.ijrect[] = c.ijhome[]
  end
end

function updatecanvas!(c::Canvas; invalidate=true, delay=0.1, first=false)
  invalidate && (c.dirty[] = true)
  isopen(c.parent.scene) || return
  c.dirty[] || return
  first || c.task[].state === :done || return
  c.task[] = @async begin
    sleep(delay)
    realigncanvas!(c)
    applyconstraints!(c)
    c.datasrc.generate!(c.buf[], c)
    c.buf[] = c.buf[]  # force observable triggers
    c.dirty[] = false
  end
end

function pancanvas!(c::Canvas, dx; rubberband=true)
  c.datasrc.xpan || c.datasrc.ypan || return
  c.datasrc.xpan || dx[1] == 0 || (dx = [0, dx[2]])
  c.datasrc.ypan || dx[2] == 0 || (dx = [dx[1], 0])
  ijrect = c.ijrect[]
  c.ijrect[] = ℛ(
    round(Int, left(ijrect) - dx[1]),
    round(Int, bottom(ijrect) + dx[2]),
    width(ijrect),
    height(ijrect)
  )
  (c.datasrc.rubberband && rubberband) || realigncanvas!(c)
  updatecanvas!(c)
end

function zoomcanvas!(c::Canvas, sx; rubberband=true)
  c.datasrc.xzoom || c.datasrc.yzoom || return
  c.datasrc.xzoom || sx[1] == 1 || (sx = [1, sx[2]])
  c.datasrc.yzoom || sx[2] == 1 || (sx = [sx[1], 1])
  ijrect = c.ijrect[]
  x0 = left(ijrect) + width(ijrect)/2
  y0 = bottom(ijrect) + height(ijrect)/2
  w = sx[1] * width(ijrect)
  h = sx[2] * height(ijrect)
  c.ijrect[] = ℛ(
    round(Int, x0 - w/2),
    round(Int, y0 - h/2),
    round(Int, w),
    round(Int, h)
  )
  (c.datasrc.rubberband && rubberband) || realigncanvas!(c)
  updatecanvas!(c)
end

function zoomcanvas!(c::Canvas, xyrect::ℛ)
  c.datasrc.xzoom || c.datasrc.yzoom || return
  c.datasrc.xpan || (xyrect = ℛ(xyrect; left=left(c.xyrect)))
  c.datasrc.ypan || (xyrect = ℛ(xyrect; bottom=bottom(c.xyrect)))
  c.datasrc.xzoom || (xyrect = ℛ(xyrect; width=width(c.xyrect)))
  c.datasrc.yzoom || (xyrect = ℛ(xyrect; height=height(c.xyrect)))
  if c.datasrc.xylock
    ratio = aspectratio(c.xyrect)
    if aspectratio(xyrect) > ratio
      xyrect = ℛ(xyrect; height=width(xyrect)/ratio)
    else
      xyrect = ℛ(xyrect; width=height(xyrect)*ratio)
    end
  end
  c.xyrect[] = xyrect
  updatecanvas!(c)
end

function brighten!(c::Canvas, dx)
  isdefined(c, :clim) || return
  clim = c.clim[]
  r = clim[2] - clim[1]
  c.clim[] = clim .- r*dx
end

function contrast!(c::Canvas, dx)
  isdefined(c, :clim) || return
  clim = c.clim[]
  r = clim[2] - clim[1]
  clim = clim .- [-r*dx, r*dx]
  clim[2] > clim[1] && (c.clim[] = (clim[1], clim[2]))
end

function resetcanvas!(c::Canvas)
  c.xyrect[] = c.datasrc.xyrect
  isdefined(c, :clim) && (c.clim[] = c.datasrc.clim)
  updatecanvas!(c; delay=0)
end

function bindevents!(c::Canvas)
  on(c.parent.scene.events.scroll) do dx
    pancanvas!(c, dx)
  end
  on(c.parent.scene.events.mousebuttons) do but
    ispressed(but, Mouse.right) && resetcanvas!(c)
  end
  on(c.parent.scene.events.keyboardbuttons) do but
    ispressed(but, keymap[:pan_left]) && pancanvas!(c, [-width(c.ijrect)/8, 0.0])
    ispressed(but, keymap[:pan_right]) && pancanvas!(c, [width(c.ijrect)/8, 0.0])
    ispressed(but, keymap[:pan_up]) && pancanvas!(c, [0.0, -height(c.ijrect)/8])
    ispressed(but, keymap[:pan_down]) && pancanvas!(c, [0.0, height(c.ijrect)/8])
    if c.datasrc.xylock
      ispressed(but, keymap[:zoom_x_out]) && zoomcanvas!(c, [1/1.2, 1/1.2])
      ispressed(but, keymap[:zoom_x_in]) && zoomcanvas!(c, [1.2, 1.2])
      ispressed(but, keymap[:zoom_y_out]) && zoomcanvas!(c,[1/1.2, 1/1.2])
      ispressed(but, keymap[:zoom_y_in]) && zoomcanvas!(c, [1.2, 1.2])
    else
      ispressed(but, keymap[:zoom_x_out]) && zoomcanvas!(c, [1/1.2, 1.0])
      ispressed(but, keymap[:zoom_x_in]) && zoomcanvas!(c, [1.2, 1.0])
      ispressed(but, keymap[:zoom_y_out]) && zoomcanvas!(c,[1.0, 1/1.2])
      ispressed(but, keymap[:zoom_y_in]) && zoomcanvas!(c, [1.0, 1.2])
    end
    ispressed(but, keymap[:brightness_down]) && brighten!(c, -0.1)
    ispressed(but, keymap[:brightness_up]) && brighten!(c, 0.1)
    ispressed(but, keymap[:contrast_down]) && contrast!(c, -0.1)
    ispressed(but, keymap[:contrast_up]) && contrast!(c, 0.1)
    ispressed(but, keymap[:reset]) && resetcanvas!(c)
  end
  on(c.parent.selrect) do r
    margin = 100
    ijrect = c.ijrect[]
    left(ijrect) - margin <= left(r) <= right(ijrect) + margin || return
    bottom(ijrect) - margin <= bottom(r) <= top(ijrect) + margin || return
    left(ijrect) - margin <= right(r) <= right(ijrect) + margin || return
    bottom(ijrect) - margin <= top(r) <= top(ijrect) + margin || return
    zoomcanvas!(c, ij2xy(r, c))
  end
end

function mkcanvas(::Type{HeatmapCanvas}, viz, ijhome, datasrc, kwargs)
  canvas = HeatmapCanvas(
    parent = viz,
    ijhome = ijhome,
    ijrect = Node(ijhome[]),
    xyrect = Node(datasrc.xyrect),
    clim = Node(datasrc.clim),
    buf = Node(zeros(Float32, width(ijhome), height(ijhome))),
    datasrc = datasrc,
    dirty = Node{Bool}(true),
    task = Ref{Task}()
  )
  x = lift(r -> left(r) : right(r), canvas.ijrect)
  y = lift(r -> bottom(r) : top(r), canvas.ijrect)
  heatmap!(viz.scene, x, y, canvas.buf; show_axis=false, colorrange=canvas.clim, kwargs...)
  on(ijhome) do r
    xscale = width(canvas.xyrect) / width(canvas.ijrect)
    yscale = height(canvas.xyrect) / height(canvas.ijrect)
    canvas.ijrect[] = r
    canvas.xyrect[] = ℛ(canvas.xyrect; width=xscale*width(canvas.ijrect), height=yscale*height(canvas.ijrect))
    canvas.buf[] = zeros(Float32, width(r), height(r))
    updatecanvas!(canvas)
  end
  canvas
end

function mkcanvas(::Type{ScatterCanvas}, viz, ijhome, datasrc, kwargs)
  canvas = ScatterCanvas(
    parent = viz,
    ijhome = ijhome,
    ijrect = Node(ijhome[]),
    xyrect = Node(datasrc.xyrect),
    buf = Node(Point2f0[]),
    datasrc = datasrc,
    dirty = Node{Bool}(true),
    task = Ref{Task}()
  )
  if :markersize ∈ keys(kwargs)
    scatter!(viz.scene, canvas.buf; show_axis=false, kwargs...)
  else
    scatter!(viz.scene, canvas.buf; show_axis=false, markersize=1, kwargs...)
  end
  on(ijhome) do r
    canvas.ijrect[] = r
    updatecanvas!(canvas)
  end
  canvas
end

function mkcanvas(::Type{LineCanvas}, viz, ijhome, datasrc, kwargs)
  canvas = LineCanvas(
    parent = viz,
    ijhome = ijhome,
    ijrect = Node(ijhome[]),
    xyrect = Node(datasrc.xyrect),
    buf = Node(Point2f0[(Inf32, Inf32)]),
    datasrc = datasrc,
    dirty = Node{Bool}(true),
    task = Ref{Task}()
  )
  lines!(viz.scene, canvas.buf; show_axis=false, kwargs...)
  on(ijhome) do r
    canvas.ijrect[] = r
    updatecanvas!(canvas)
  end
  canvas
end

function addcanvas!(ctype, viz::Viz, datasrc::DataSource; rect=nothing, kwargs...)
  if rect === nothing
    ijhome = lift(r -> r, viz.ijrect)
  elseif rect isa ℛ
    ijhome = Node(rect)
  elseif rect isa Node
    ijhome = rect
  else
    ijhome = Node(ℛ(rect[1], rect[2], rect[3], rect[4]))
  end
  canvas = mkcanvas(ctype, viz, ijhome, datasrc, kwargs)
  push!(viz.children, canvas)
  bindevents!(canvas)
  isopen(viz.scene) && updatecanvas!(canvas; first=true, delay=0)
  canvas
end

# TODO: add support for legend
# TODO: add support for colorbar

function addaxes!(c::Canvas; inset=0, color=:black, frame=false, grid=false, border=0, bordercolor=:white, xlabel=missing, ylabel=missing, xticks=5, yticks=5, ticksize=10, textsize=15.0)
  scene = c.parent.scene
  if border > 0
    poly!(scene, lift(r -> Point2f0[
      (left(r) - border, bottom(r)),
      (right(r) + border, bottom(r)),
      (right(r) + border, bottom(r) - border),
      (left(r) - border, bottom(r) - border)
    ], c.ijhome); show_axis=false, color=bordercolor)
    poly!(scene, lift(r -> Point2f0[
      (left(r) - border, top(r)),
      (right(r) + border, top(r)),
      (right(r) + border, top(r) + border),
      (left(r) - border, top(r) + border)
    ], c.ijhome); show_axis=false, color=bordercolor)
    poly!(scene, lift(r -> Point2f0[
      (left(r), bottom(r)),
      (left(r) - border, bottom(r)),
      (left(r) - border, top(r)),
      (left(r), top(r))
    ], c.ijhome); show_axis=false, color=bordercolor)
    poly!(scene, lift(r -> Point2f0[
      (right(r), bottom(r)),
      (right(r) + border, bottom(r)),
      (right(r) + border, top(r)),
      (right(r), top(r))
    ], c.ijhome); show_axis=false, color=bordercolor)
  end
  r = c.ijhome
  inset != 0 && (r = lift(r -> ℛ(
    left(r) + inset, bottom(r) + inset,
    width(r) - 2*inset - 1, height(r) - 2*inset - 1), c.ijhome))
  if frame
    lines!(scene, lift(r -> Point2f0[
      (right(r), bottom(r)),
      (left(r), bottom(r)),
      (left(r), top(r)),
      (right(r), top(r)),
      (right(r), bottom(r))
    ], r); show_axis=false, color=color)
  else
    lines!(scene, lift(r -> Point2f0[
      (right(r), bottom(r)),
      (left(r), bottom(r)),
      (left(r), top(r))
    ], r); show_axis=false, color=color)
  end
  xticktext(i, c) = formatnums(ij2xy(i, 0, c)[1])
  if xticks > 0
    linesegments!(scene, lift(r -> [
      Point2f0(left(r) + k * width(r) / xticks, bottom(r)) =>
      Point2f0(left(r) + k * width(r) / xticks, bottom(r) - ticksize)
      for k ∈ 0:xticks
    ], r); show_axis=false, color=color)
    grid && linesegments!(scene, lift(r -> [
      Point2f0(left(r) + k * width(r) / xticks, bottom(r)) =>
      Point2f0(left(r) + k * width(r) / xticks, top(r))
      for k ∈ 0:xticks
    ], r); show_axis=false, color=color, linestyle=:dot)
    for k ∈ 0:xticks
      text!(scene, lift((r, xy) -> xticktext(left(r) + k * width(r) / xticks, c), r, c.xyrect);
        position=lift(r -> (left(r) + k * width(r) / xticks, bottom(r) - ticksize - textsize/2), r),
        textsize=textsize, color=color,
        align=(:center, :center))
    end
  end
  yticktext(j, c) = formatnums(ij2xy(0, j, c)[2])
  if yticks > 0
    linesegments!(scene, lift(r -> [
      Point2f0(left(r), bottom(r) + k * height(r) / yticks) =>
      Point2f0(left(r) - ticksize, bottom(r) + k * height(r) / yticks)
      for k ∈ 0:yticks
    ], r); show_axis=false, color=color)
    grid && linesegments!(scene, lift(r -> [
      Point2f0(left(r), bottom(r) + k * height(r) / yticks) =>
      Point2f0(right(r), bottom(r) + k * height(r) / yticks)
      for k ∈ 0:yticks
    ], r); show_axis=false, color=color, linestyle=:dot)
    for k ∈ 0:yticks
      text!(scene, lift((r, xy) -> yticktext(bottom(r) + k * height(r) / yticks, c), r, c.xyrect);
        position=lift(r -> (float(left(r) - ticksize), bottom(r) + k * height(r) / yticks), r),
        textsize=textsize, color=color,
        align=(:right, :center))
    end
  end
  # TODO: remove hardcoded 50 and 75
  xlabel === missing || text!(scene, xlabel; position=lift(r ->
    (left(r) + width(r)/2 - 50, bottom(r) - 4*ticksize - textsize), r),
    textsize=textsize, color=color, align=(:center, :bottom))
  ylabel === missing || text!(scene, ylabel; position=lift(r ->
      (left(r) - 75, bottom(r) + height(r)÷2), r),
      textsize=textsize, color=color, align=(:center, :top), rotation=π/2)
  nothing
end

function addcursor!(c::Canvas; position=nothing, color=:black, textsize=15.0, align=(:center, :center))
  if position === nothing
    position = lift(r -> (right(r) - 20, top(r) - 20), c.ijhome)
    align=(:right, :top)
  end
  s = Node(" ")
  text!(c.parent.scene, s; position=position, align=align, textsize=textsize, color=color)
  on(c.parent.scene.events.mouseposition) do ij
    if left(c.ijhome) <= ij[1] <= right(c.ijhome) && bottom(c.ijhome) <= ij[2] <= top(c.ijhome)
      x, y = ij2xy(ij[1], ij[2], c)
      if c isa HeatmapCanvas
        try
          s[] = formatnums(x, y, c.buf[][round(Int, ij[1] - left(c.ijhome)) + 1, round(Int, ij[2] - bottom(c.ijhome)) + 1])
        catch
          s[] = " "
        end
      else
        s[] = formatnums(x, y)
      end
    else
      s[] = " "
    end
  end
  on(c.parent.scene.events.entered_window) do b
    b || (s[] = " ")
  end
  nothing
end

function formatnums(x...)
  s = ""
  for x1 ∈ x
    length(s) > 0 && (s *= ", ")
    s *= formatnum(x1)
  end
  " " * s * " "
end

function formatnum(x)
  x == 0 && return "0"
  abs(x) > 10000 && return @sprintf("%.0e", x)
  abs(x) >= 100 && return @sprintf("%.0f", x)
  abs(x) >= 10 && return @sprintf("%.1f", x)
  abs(x) >= 1 && return @sprintf("%.2f", x)
  abs(x) < 0.001 && return @sprintf("%.0e", x)
  @sprintf("%.3f", x)
end


text:
using AttractorNetworksBase ; const A=AttractorNetworksBase
using Test
using Calculus,LinearAlgebra,Statistics
using Random ; Random.seed!(1)
const _rtol = 1E-4

function test_Jgradient(utest::Vector{R},ntw::A.BaseNetwork) where R
    Jalloc = A.JAlloc(ntw)
    gradW=similar(ntw.weights)
    gradu=similar(ntw.weights)
    J=similar(gradW)
    # analytic gradient
    A.jacobian!(J,gradW,gradu,utest,ntw,Jalloc)
    # numerical gradient weights
    grad_numW = similar(gradW)
    ntw2=copy(ntw)
    for ij in eachindex(ntw2.weights)
        function gradfun(w)
            ntw2.weights[ij] = w
            A.jacobian!(J,nothing,nothing,utest,ntw2,Jalloc)
            return J[ij]
        end
        grad_numW[ij] = Calculus.gradient(gradfun,ntw2.weights[ij])
    end
    # numerical gradient u
    grad_numu = similar(gradu)
    for lj in CartesianIndices(gradu)
        (j,l)=Tuple(lj)
        function gradfun(u)
            u2 = copy(utest)
            u2[l]=u
            A.jacobian!(J,nothing,nothing,u2,ntw,Jalloc)
            return J[lj]
        end
        grad_numu[lj] = Calculus.gradient(gradfun,utest[l])
    end
    return gradW,grad_numW,gradu,grad_numu
end

##

@testset "Gain functions" begin
    # identity
    g = A.IOId{Float64}()
    x = randn(1_000)
    @test all(isapprox.(g.(x),x))
    @test all(isapprox.(A.ioinv.(x,g),x))
    xdest = A.ioprime.(x,g)
    @test all(isapprox.(xdest,1.0))
    xdest = A.ioprimeprime.(x,g)
    @test all(isapprox.(xdest,0.0))
    # quadratic
    x = 3 .* randn(1_000)
    g = A.IOQuad(1+rand())
    y = g.(x)
    @test all( y .> 0.0)
    xdest = g.(x)
    @test all( isapprox.(xdest,y))
    @test all(isapprox.(A.ioinv.(y,g),x ; rtol=_rtol))
    dgnum = Calculus.gradient.(_x->g(_x),x)
    xdest=A.ioprime.(x,g)
    @test all(isapprox.(xdest,dgnum; rtol=_rtol))
    ddgnum = Calculus.gradient.(_x->A.ioprime(_x,g),x)
    xdest=A.ioprimeprime.(x,g)
    @test all(isapprox.(xdest,ddgnum; rtol=_rtol))
end

@testset "Weight matrix" begin
    ne = 82
    ne_aux = 22
    ni = 64
    ni_aux = 17
    wmeanf = [ 0.333  -1.123
                12.122 -0.06432]
    W = A.make_wmat(ne,ni,wmeanf;noautapses=true)
    @test tr(W) == 0.0
    @test size(W,1) == size(W,2)
    @test begin
        n_tot = size(W,1)
        n_tot == ne+ni
    end
    # test if Dalian
    pe,pi = A.populations(ne,ni)
    @test all(W[:,pe] .>= 0.0)
    @test all(W[:,pi] .<= 0.0)
    @test  all(isapprox.(sum(W[pe,pe];dims=2), 0.333 ; rtol=1E-3))
    @test  all(isapprox.(sum(W[pi,pe];dims=2), 12.122 ; rtol=1E-3))
end

@testset "Network constructor" begin
    ne = 13
    ni = 10
    ntw = A.BaseNetwork(ne,ni ; gfun=A.IOQuad(0.123) )
    @test all( size(ntw.weights) .== (ne+ni) )
    @test A.n_neurons(ntw) == (ne+ni)
end

# Jacobian !
@testset "Jacobian matrix" begin
    ne,ni = 13,10
    ntot = ne+ni
    ntw = A.BaseNetwork(ne,ni ; gfun=A.IOQuad(0.123) )
    v_alloc = zeros(ntot)
    veli(u,i) = let v = A.velocity!(v_alloc,u,ntw.iofunction.(u),ntw) ; v[i]; end
    @info "building the Jacobian numerically"
    utest = randn(ntot)
    Jnum = Matrix{Float64}(undef,ntot,ntot)
    for i in 1:ntot
      Jnum[i,:] =  Calculus.gradient( u -> veli(u,i), utest )
    end
    Jan=A.jacobian(utest,ntw)
    @test all(isapprox.(Jnum,Jan;rtol=1E-4))
end

@testset "Jacobian gradients" begin
    ne,ni = 13,10
    ntot = ne+ni
    ntw = A.BaseNetwork(ne,ni ; gfun=A.IOQuad(0.123) )
    utest = randn(ntot)
    anaW,numW,anau,numu=test_Jgradient(utest,ntw)
    @test all(isapprox.(anaW,numW;rtol=1E-4))
    @test all(isapprox.(anau,numu;rtol=1E-4))
    ne,ni = 23,12
    ntot = ne+ni
    ntw = A.BaseNetwork(ne,ni ; gfun=A.IOId{Float64}() )
    utest = randn(ntot)
    # weights
    anaW,numW,anau,numu=test_Jgradient(utest,ntw)
    @test all(isapprox.(anaW,numW;rtol=1E-4))
    @test all(isapprox.(anau,numu;rtol=1E-4))
end

@testset "Dynamics" begin
    ne,ni = 20,23
    ntot = ne+ni
    ntw = A.BaseNetwork(ne,ni ; gfun=A.IOQuad(0.02))
    fill!(ntw.weights,0.0)
    u_start = randn(A.ndims(ntw))
    r_start = ntw.iofunction.(u_start)
    u_end,_=A.run_network_to_convergence(ntw,r_start)
    @test all(isapprox.(u_end,ntw.external_input;atol=0.05))
    ntw = A.BaseNetwork(ne,ni ; gfun=A.IOQuad(0.02))
    r_start = ntw.iofunction.(u_start)
    u_end,r_end=A.run_network_to_convergence(ntw,r_start)
    t,ut,rt=A.run_network(ntw,r_end,1.0)
    @test all(isapprox.(u_end,ut[:,end];atol=0.05))
end

@testset "Making attractors" begin
    ne,ni = 200,100
    natt=431
    ntw = A.AttractorNetwork(natt,ne,ni ;
        gfun=A.IOQuad(3.45),mu_attr=3.23,std_attr=2.2)
    attr_u = ntw.attractors_u
    attr_r = ntw.iofunction.(attr_u)
    @test isapprox(mean(attr_r),3.23;atol=0.05)
    @test isapprox(std(attr_r),2.2;atol=0.1)
    @test A.n_attractors(ntw) == natt
end