text:
@testset "4.2.1.1 (a+b cos)^n" begin
    (a, b, c, d, n, x, ) = @variables a b c d n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(a+b*cos(c+d*x))^n=#

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

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

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

    @test_int [(a + a*cos(c + d*x))^(7/2), x, 4, (256*a^4*sin(c + d*x))/(35*d*sqrt(a + a*cos(c + d*x))) + (64*a^3*sqrt(a + a*cos(c + d*x))*sin(c + d*x))/(35*d) + (24*a^2*(a + a*cos(c + d*x))^(3/2)*sin(c + d*x))/(35*d) + (2*a*(a + a*cos(c + d*x))^(5/2)*sin(c + d*x))/(7*d)]
    @test_int [(a + a*cos(c + d*x))^(5/2), x, 3, (64*a^3*sin(c + d*x))/(15*d*sqrt(a + a*cos(c + d*x))) + (16*a^2*sqrt(a + a*cos(c + d*x))*sin(c + d*x))/(15*d) + (2*a*(a + a*cos(c + d*x))^(3/2)*sin(c + d*x))/(5*d)]
    @test_int [(a + a*cos(c + d*x))^(3/2), x, 2, (8*a^2*sin(c + d*x))/(3*d*sqrt(a + a*cos(c + d*x))) + (2*a*sqrt(a + a*cos(c + d*x))*sin(c + d*x))/(3*d)]
    @test_int [(a + a*cos(c + d*x))^(1/2), x, 1, (2*a*sin(c + d*x))/(d*sqrt(a + a*cos(c + d*x)))]
    @test_int [1/(a + a*cos(c + d*x))^(1/2), x, 2, (sqrt(2)*atanh((sqrt(a)*sin(c + d*x))/(sqrt(2)*sqrt(a + a*cos(c + d*x)))))/(sqrt(a)*d)]
    @test_int [1/(a + a*cos(c + d*x))^(3/2), x, 3, atanh((sqrt(a)*sin(c + d*x))/(sqrt(2)*sqrt(a + a*cos(c + d*x))))/(2*sqrt(2)*a^(3/2)*d) + sin(c + d*x)/(2*d*(a + a*cos(c + d*x))^(3/2))]
    @test_int [1/(a + a*cos(c + d*x))^(5/2), x, 4, (3*atanh((sqrt(a)*sin(c + d*x))/(sqrt(2)*sqrt(a + a*cos(c + d*x)))))/(16*sqrt(2)*a^(5/2)*d) + sin(c + d*x)/(4*d*(a + a*cos(c + d*x))^(5/2)) + (3*sin(c + d*x))/(16*a*d*(a + a*cos(c + d*x))^(3/2))]

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

    @test_int [(a + a*cos(c + d*x))^(4/3), x, 2, (2*2^(5/6)*a*(a + a*cos(c + d*x))^(1/3)*HypergeometricFunctions._₂F₁(-(5/6), 1/2, 3/2, (1/2)*(1 - cos(c + d*x)))*sin(c + d*x))/(d*(1 + cos(c + d*x))^(5/6))]
    @test_int [(a + a*cos(c + d*x))^(2/3), x, 2, (2*2^(1/6)*(a + a*cos(c + d*x))^(2/3)*HypergeometricFunctions._₂F₁(-(1/6), 1/2, 3/2, (1/2)*(1 - cos(c + d*x)))*sin(c + d*x))/(d*(1 + cos(c + d*x))^(7/6))]
    @test_int [(a + a*cos(c + d*x))^(1/3), x, 2, (2^(5/6)*(a + a*cos(c + d*x))^(1/3)*HypergeometricFunctions._₂F₁(1/6, 1/2, 3/2, (1/2)*(1 - cos(c + d*x)))*sin(c + d*x))/(d*(1 + cos(c + d*x))^(5/6))]
    @test_int [1/(a + a*cos(c + d*x))^(1/3), x, 2, (2^(1/6)*HypergeometricFunctions._₂F₁(1/2, 5/6, 3/2, (1/2)*(1 - cos(c + d*x)))*sin(c + d*x))/(d*(1 + cos(c + d*x))^(1/6)*(a + a*cos(c + d*x))^(1/3))]
    @test_int [1/(a + a*cos(c + d*x))^(2/3), x, 2, ((1 + cos(c + d*x))^(1/6)*HypergeometricFunctions._₂F₁(1/2, 7/6, 3/2, (1/2)*(1 - cos(c + d*x)))*sin(c + d*x))/(2^(1/6)*d*(a + a*cos(c + d*x))^(2/3))]
    @test_int [1/(a + a*cos(c + d*x))^(4/3), x, 2, (HypergeometricFunctions._₂F₁(1/2, 11/6, 3/2, (1/2)*(1 - cos(c + d*x)))*sin(c + d*x))/(2^(5/6)*a*d*(1 + cos(c + d*x))^(1/6)*(a + a*cos(c + d*x))^(1/3))]

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

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

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

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

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

    @test_int [1/(5 + 3*cos(c + d*x)), x, 1, x/4 - atan(sin(c + d*x)/(3 + cos(c + d*x)))/(2*d)]
    @test_int [1/(5 + 3*cos(c + d*x))^2, x, 3, (5*x)/64 - (5*atan(sin(c + d*x)/(3 + cos(c + d*x))))/(32*d) - (3*sin(c + d*x))/(16*d*(5 + 3*cos(c + d*x)))]
    @test_int [1/(5 + 3*cos(c + d*x))^3, x, 4, (59*x)/2048 - (59*atan(sin(c + d*x)/(3 + cos(c + d*x))))/(1024*d) - (3*sin(c + d*x))/(32*d*(5 + 3*cos(c + d*x))^2) - (45*sin(c + d*x))/(512*d*(5 + 3*cos(c + d*x)))]
    @test_int [1/(5 + 3*cos(c + d*x))^4, x, 5, (385*x)/32768 - (385*atan(sin(c + d*x)/(3 + cos(c + d*x))))/(16384*d) - sin(c + d*x)/(16*d*(5 + 3*cos(c + d*x))^3) - (25*sin(c + d*x))/(512*d*(5 + 3*cos(c + d*x))^2) - (311*sin(c + d*x))/(8192*d*(5 + 3*cos(c + d*x)))]

    @test_int [1/(5 - 3*cos(c + d*x)), x, 1, x/4 + atan(sin(c + d*x)/(3 - cos(c + d*x)))/(2*d)]
    @test_int [1/(5 - 3*cos(c + d*x))^2, x, 3, (5*x)/64 + (5*atan(sin(c + d*x)/(3 - cos(c + d*x))))/(32*d) + (3*sin(c + d*x))/(16*d*(5 - 3*cos(c + d*x)))]
    @test_int [1/(5 - 3*cos(c + d*x))^3, x, 4, (59*x)/2048 + (59*atan(sin(c + d*x)/(3 - cos(c + d*x))))/(1024*d) + (3*sin(c + d*x))/(32*d*(5 - 3*cos(c + d*x))^2) + (45*sin(c + d*x))/(512*d*(5 - 3*cos(c + d*x)))]
    @test_int [1/(5 - 3*cos(c + d*x))^4, x, 5, (385*x)/32768 + (385*atan(sin(c + d*x)/(3 - cos(c + d*x))))/(16384*d) + sin(c + d*x)/(16*d*(5 - 3*cos(c + d*x))^3) + (25*sin(c + d*x))/(512*d*(5 - 3*cos(c + d*x))^2) + (311*sin(c + d*x))/(8192*d*(5 - 3*cos(c + d*x)))]

    @test_int [1/(-5 + 3*cos(c + d*x)), x, 1, -(x/4) - atan(sin(c + d*x)/(3 - cos(c + d*x)))/(2*d)]
    @test_int [1/(-5 + 3*cos(c + d*x))^2, x, 3, (5*x)/64 + (5*atan(sin(c + d*x)/(3 - cos(c + d*x))))/(32*d) + (3*sin(c + d*x))/(16*d*(5 - 3*cos(c + d*x)))]
    @test_int [1/(-5 + 3*cos(c + d*x))^3, x, 4, -((59*x)/2048) - (59*atan(sin(c + d*x)/(3 - cos(c + d*x))))/(1024*d) - (3*sin(c + d*x))/(32*d*(5 - 3*cos(c + d*x))^2) - (45*sin(c + d*x))/(512*d*(5 - 3*cos(c + d*x)))]
    @test_int [1/(-5 + 3*cos(c + d*x))^4, x, 5, (385*x)/32768 + (385*atan(sin(c + d*x)/(3 - cos(c + d*x))))/(16384*d) + sin(c + d*x)/(16*d*(5 - 3*cos(c + d*x))^3) + (25*sin(c + d*x))/(512*d*(5 - 3*cos(c + d*x))^2) + (311*sin(c + d*x))/(8192*d*(5 - 3*cos(c + d*x)))]

    @test_int [1/(-5 - 3*cos(c + d*x)), x, 1, -(x/4) + atan(sin(c + d*x)/(3 + cos(c + d*x)))/(2*d)]
    @test_int [1/(-5 - 3*cos(c + d*x))^2, x, 3, (5*x)/64 - (5*atan(sin(c + d*x)/(3 + cos(c + d*x))))/(32*d) - (3*sin(c + d*x))/(16*d*(5 + 3*cos(c + d*x)))]
    @test_int [1/(-5 - 3*cos(c + d*x))^3, x, 4, -((59*x)/2048) + (59*atan(sin(c + d*x)/(3 + cos(c + d*x))))/(1024*d) + (3*sin(c + d*x))/(32*d*(5 + 3*cos(c + d*x))^2) + (45*sin(c + d*x))/(512*d*(5 + 3*cos(c + d*x)))]
    @test_int [1/(-5 - 3*cos(c + d*x))^4, x, 5, (385*x)/32768 - (385*atan(sin(c + d*x)/(3 + cos(c + d*x))))/(16384*d) - sin(c + d*x)/(16*d*(5 + 3*cos(c + d*x))^3) - (25*sin(c + d*x))/(512*d*(5 + 3*cos(c + d*x))^2) - (311*sin(c + d*x))/(8192*d*(5 + 3*cos(c + d*x)))]

    @test_int [1/(3 + 5*cos(c + d*x)), x, 2, -(log(2*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x)))/(4*d)) + log(2*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x)))/(4*d)]
    @test_int [1/(3 + 5*cos(c + d*x))^2, x, 4, (3*log(2*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x))))/(64*d) - (3*log(2*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x))))/(64*d) + (5*sin(c + d*x))/(16*d*(3 + 5*cos(c + d*x)))]
    @test_int [1/(3 + 5*cos(c + d*x))^3, x, 5, -((43*log(2*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x))))/(2048*d)) + (43*log(2*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x))))/(2048*d) + (5*sin(c + d*x))/(32*d*(3 + 5*cos(c + d*x))^2) - (45*sin(c + d*x))/(512*d*(3 + 5*cos(c + d*x)))]
    @test_int [1/(3 + 5*cos(c + d*x))^4, x, 6, (279*log(2*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x))))/(32768*d) - (279*log(2*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x))))/(32768*d) + (5*sin(c + d*x))/(48*d*(3 + 5*cos(c + d*x))^3) - (25*sin(c + d*x))/(512*d*(3 + 5*cos(c + d*x))^2) + (995*sin(c + d*x))/(24576*d*(3 + 5*cos(c + d*x)))]

    @test_int [1/(3 - 5*cos(c + d*x)), x, 2, log(cos((1/2)*(c + d*x)) - 2*sin((1/2)*(c + d*x)))/(4*d) - log(cos((1/2)*(c + d*x)) + 2*sin((1/2)*(c + d*x)))/(4*d)]
    @test_int [1/(3 - 5*cos(c + d*x))^2, x, 4, -((3*log(cos((1/2)*(c + d*x)) - 2*sin((1/2)*(c + d*x))))/(64*d)) + (3*log(cos((1/2)*(c + d*x)) + 2*sin((1/2)*(c + d*x))))/(64*d) - (5*sin(c + d*x))/(16*d*(3 - 5*cos(c + d*x)))]
    @test_int [1/(3 - 5*cos(c + d*x))^3, x, 5, (43*log(cos((1/2)*(c + d*x)) - 2*sin((1/2)*(c + d*x))))/(2048*d) - (43*log(cos((1/2)*(c + d*x)) + 2*sin((1/2)*(c + d*x))))/(2048*d) - (5*sin(c + d*x))/(32*d*(3 - 5*cos(c + d*x))^2) + (45*sin(c + d*x))/(512*d*(3 - 5*cos(c + d*x)))]
    @test_int [1/(3 - 5*cos(c + d*x))^4, x, 6, -((279*log(cos((1/2)*(c + d*x)) - 2*sin((1/2)*(c + d*x))))/(32768*d)) + (279*log(cos((1/2)*(c + d*x)) + 2*sin((1/2)*(c + d*x))))/(32768*d) - (5*sin(c + d*x))/(48*d*(3 - 5*cos(c + d*x))^3) + (25*sin(c + d*x))/(512*d*(3 - 5*cos(c + d*x))^2) - (995*sin(c + d*x))/(24576*d*(3 - 5*cos(c + d*x)))]

    @test_int [1/(-3 + 5*cos(c + d*x)), x, 2, -(log(cos((1/2)*(c + d*x)) - 2*sin((1/2)*(c + d*x)))/(4*d)) + log(cos((1/2)*(c + d*x)) + 2*sin((1/2)*(c + d*x)))/(4*d)]
    @test_int [1/(-3 + 5*cos(c + d*x))^2, x, 4, -((3*log(cos((1/2)*(c + d*x)) - 2*sin((1/2)*(c + d*x))))/(64*d)) + (3*log(cos((1/2)*(c + d*x)) + 2*sin((1/2)*(c + d*x))))/(64*d) - (5*sin(c + d*x))/(16*d*(3 - 5*cos(c + d*x)))]
    @test_int [1/(-3 + 5*cos(c + d*x))^3, x, 5, -((43*log(cos((1/2)*(c + d*x)) - 2*sin((1/2)*(c + d*x))))/(2048*d)) + (43*log(cos((1/2)*(c + d*x)) + 2*sin((1/2)*(c + d*x))))/(2048*d) + (5*sin(c + d*x))/(32*d*(3 - 5*cos(c + d*x))^2) - (45*sin(c + d*x))/(512*d*(3 - 5*cos(c + d*x)))]
    @test_int [1/(-3 + 5*cos(c + d*x))^4, x, 6, -((279*log(cos((1/2)*(c + d*x)) - 2*sin((1/2)*(c + d*x))))/(32768*d)) + (279*log(cos((1/2)*(c + d*x)) + 2*sin((1/2)*(c + d*x))))/(32768*d) - (5*sin(c + d*x))/(48*d*(3 - 5*cos(c + d*x))^3) + (25*sin(c + d*x))/(512*d*(3 - 5*cos(c + d*x))^2) - (995*sin(c + d*x))/(24576*d*(3 - 5*cos(c + d*x)))]

    @test_int [1/(-3 - 5*cos(c + d*x)), x, 2, log(2*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x)))/(4*d) - log(2*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x)))/(4*d)]
    @test_int [1/(-3 - 5*cos(c + d*x))^2, x, 4, (3*log(2*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x))))/(64*d) - (3*log(2*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x))))/(64*d) + (5*sin(c + d*x))/(16*d*(3 + 5*cos(c + d*x)))]
    @test_int [1/(-3 - 5*cos(c + d*x))^3, x, 5, (43*log(2*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x))))/(2048*d) - (43*log(2*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x))))/(2048*d) - (5*sin(c + d*x))/(32*d*(3 + 5*cos(c + d*x))^2) + (45*sin(c + d*x))/(512*d*(3 + 5*cos(c + d*x)))]
    @test_int [1/(-3 - 5*cos(c + d*x))^4, x, 6, (279*log(2*cos((1/2)*(c + d*x)) - sin((1/2)*(c + d*x))))/(32768*d) - (279*log(2*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x))))/(32768*d) + (5*sin(c + d*x))/(48*d*(3 + 5*cos(c + d*x))^3) - (25*sin(c + d*x))/(512*d*(3 + 5*cos(c + d*x))^2) + (995*sin(c + d*x))/(24576*d*(3 + 5*cos(c + d*x)))]

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

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

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

    @test_int [(a + b*cos(c + d*x))^(4/3), x, 3, (sqrt(2)*(a + b)*AppellF1(1/2, 1/2, -(4/3), 3/2, (1/2)*(1 - cos(c + d*x)), (b*(1 - cos(c + d*x)))/(a + b))*(a + b*cos(c + d*x))^(1/3)*sin(c + d*x))/(d*sqrt(1 + cos(c + d*x))*((a + b*cos(c + d*x))/(a + b))^(1/3))]
    @test_int [(a + b*cos(c + d*x))^(2/3), x, 3, (sqrt(2)*AppellF1(1/2, 1/2, -(2/3), 3/2, (1/2)*(1 - cos(c + d*x)), (b*(1 - cos(c + d*x)))/(a + b))*(a + b*cos(c + d*x))^(2/3)*sin(c + d*x))/(d*sqrt(1 + cos(c + d*x))*((a + b*cos(c + d*x))/(a + b))^(2/3))]
    @test_int [(a + b*cos(c + d*x))^(1/3), x, 3, (sqrt(2)*AppellF1(1/2, 1/2, -(1/3), 3/2, (1/2)*(1 - cos(c + d*x)), (b*(1 - cos(c + d*x)))/(a + b))*(a + b*cos(c + d*x))^(1/3)*sin(c + d*x))/(d*sqrt(1 + cos(c + d*x))*((a + b*cos(c + d*x))/(a + b))^(1/3))]
    @test_int [1/(a + b*cos(c + d*x))^(1/3), x, 3, (sqrt(2)*AppellF1(1/2, 1/2, 1/3, 3/2, (1/2)*(1 - cos(c + d*x)), (b*(1 - cos(c + d*x)))/(a + b))*((a + b*cos(c + d*x))/(a + b))^(1/3)*sin(c + d*x))/(d*sqrt(1 + cos(c + d*x))*(a + b*cos(c + d*x))^(1/3))]
    @test_int [1/(a + b*cos(c + d*x))^(2/3), x, 3, (sqrt(2)*AppellF1(1/2, 1/2, 2/3, 3/2, (1/2)*(1 - cos(c + d*x)), (b*(1 - cos(c + d*x)))/(a + b))*((a + b*cos(c + d*x))/(a + b))^(2/3)*sin(c + d*x))/(d*sqrt(1 + cos(c + d*x))*(a + b*cos(c + d*x))^(2/3))]
    @test_int [1/(a + b*cos(c + d*x))^(4/3), x, 3, (sqrt(2)*AppellF1(1/2, 1/2, 4/3, 3/2, (1/2)*(1 - cos(c + d*x)), (b*(1 - cos(c + d*x)))/(a + b))*((a + b*cos(c + d*x))/(a + b))^(1/3)*sin(c + d*x))/((a + b)*d*sqrt(1 + cos(c + d*x))*(a + b*cos(c + d*x))^(1/3))]

    #= [(a + b*cos(c + d*x))^(4/3) - (4*a^2 + b^2 + 5*a*b*cos(c + d*x))/(4*(a + b*cos(c + d*x))^(2/3)), x, -11, (3*b*(a + b*cos(c + d*x))^(1/3)*sin(c + d*x))/(4*d)] =#

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

    @test_int [(a + b*cos(c + d*x))^n, x, 3, (sqrt(2)*AppellF1(1/2, 1/2, -n, 3/2, (1/2)*(1 - cos(c + d*x)), (b*(1 - cos(c + d*x)))/(a + b))*(a + b*cos(c + d*x))^n*sin(c + d*x))/(((a + b*cos(c + d*x))/(a + b))^n*(d*sqrt(1 + cos(c + d*x))))]
    end

max_stars_repo_path: test/4 Trig functions/4.2 Cosine/4.2.1.1 (a+b cos)^n.jl


text:
using MatrixPolynomials
using Test

include("testPolynomials.jl")
include("testMatrices.jl")
include("testEigenvalues.jl")

ext: jl


text:
# Repeated content
# ----------------

"""
    plan_repetition(cb::ChromosomeBlueprint, from::UnitRange{Int}, to::UnitRange{Int})

Plan a repetition in a chromosome, where the bases in the `from` region of the
chromosome, are guaranteed to occur again in the `to` region of the chromosome.

Every copy of the chromosome will have this same repetition.

Creates a new chromosome blueprint, based on the input blueprint `cb`.

!!! note
    In the new blueprint, the `to` region of the planned chromosome
    will be consumed, and cannot be used to plan any other subsequently added
    features.
"""
function plan_repetition(cb::ChromosomeBlueprint, from::UnitRange{Int}, to::UnitRange{Int})
    if length(from) != length(to)
        throw(ArgumentError("from and to must be the same length"))
    end
    checkavailability(cb, to)
    newspec = consume_regions(cb, to)
    newrep = RepetitionBlueprint(first(from), first(to), length(from))
    return ChromosomeBlueprint(newspec, push!(repetitions(newspec), newrep))
end

"""
    plan_repetition(cb::ChromosomeBlueprint, from::Int, to::Int, size::Int)

Plan a repetition in a chromosome, where the bases in the
`from`:(`from` + `size` - 1) region of the chromosome, are guaranteed to occur
again in the `to`:(`to` + `size` - 1) region of the chromosome.

Every copy of the chromosome will have this same repetition.

Creates a new chromosome blueprint, based on the input blueprint `cb`.
"""
function plan_repetition(cb::ChromosomeBlueprint, from::Int, to::Int, size::Int)
    rfrom = from:(from + size - 1)
    rto = to:(to + size - 1)
    return plan_repetition(cb, rfrom, rto)
end

"""
    plan_repetition(cb::ChromosomeBlueprint, intervals::Vector{UnitRange{Int}})

A conveinience method of `plan_repetition`.
Designed to ease the process of planing a series of repetitions in a chromosome.

Every copy of the chromosome will have these same repetitions.

Creates a new chromosome blueprint, based on the input blueprint `cb`.

!!! tip
    Use the [`suggest_regions`](@ref) function to help decide on a set of
    sites to make heterozygous.

!!! note
    The number of intervals provided must be an even number. This is because
    intervals 1 & 2  define the first repeat, intervals 3 & 4 define the second, and
    so on.

!!! note
    In the new blueprint, the regions the repetitions occupy, have been consumed,
    and cannot be used to plan any other subsequently added features.
"""
function plan_repetition(cb::ChromosomeBlueprint, intervals::Vector{UnitRange{Int}})
    @assert iseven(length(intervals))
    lastcb = cb
    for i in firstindex(intervals):2:lastindex(intervals)
        lastcb = plan_repetition(lastcb, intervals[i], intervals[i + 1])
    end
    return lastcb
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# n-queens (nqueens) solver, for nsquaresx-by-nsquaresy board

immutable Queen
    x::Int
    y::Int
end
hitshorz(queena, queenb) = queena.x == queenb.x
hitsvert(queena, queenb) = queena.y == queenb.y
hitsdiag(queena, queenb) = abs(queena.x - queenb.x) == abs(queena.y - queenb.y)
hitshvd(qa, qb) = hitshorz(qa, qb) || hitsvert(qa, qb) || hitsdiag(qa, qb)
hitsany(testqueen, queens) = any(q -> hitshvd(testqueen, q), queens)

function trysolve(nsquaresx, nsquaresy, nqueens, presqueens = ())
    nqueens == 0 && return presqueens
    for xsquare in 1:nsquaresx
        for ysquare in 1:nsquaresy
            testqueen = Queen(xsquare, ysquare)
            if !hitsany(testqueen, presqueens)
                tryqueens = (presqueens..., testqueen)
                maybesol = trysolve(nsquaresx, nsquaresy, nqueens - 1, tryqueens)
                maybesol !== nothing && return maybesol
            end
        end
    end
    return nothing
end


text:
const DEFAULT_DELIM = ','
const DEFAULT_QUOTE = '"'
const DEFAULT_NUM_ROWS = SizeType(2048)
const DEFAULT_LINE_LEN = DEFAULT_NUM_ROWS * SizeType(1024)
const DEFAULT_NUM_FIELDS = SizeType(16)
const DEFAULT_EAGER_PARSE_FIELDS = true

"""
A mutable struct for keeping a counter

This struct can be used as a field of another immutable struct
"""
mutable struct Counter
    v::Int
end

# Mutable File Position
const MutFilePos = Counter

# Immutable File Position
struct FilePos
    v::Int
end

Base.hash(x::FilePos, h::UInt) = hash(x.v)
Base.isless(x::FilePos, y::FilePos) = isless(x.v, y.v)
Base.isequal(x::FilePos, y::FilePos) = x.v == y.v

const DUMMY_STRING = ""
const DUMMY_WEAK_STRING = WeakRefString{UInt8}(pointer(DUMMY_STRING),0)

"""
This struct holds all the information in a record inside a CSV file

Note: it is mutable, as it'll be re-used multiple times (once per block of records)
"""
mutable struct Record
    is_valid::Bool
    pos::FilePos
    fields_buff::BufferedVector{WeakRefString{UInt8}}
    line::WeakRefString{UInt8}

    function Record(fields_buff_len::SizeType) where {T}
        fields = BufferedVector{WeakRefString{UInt8}}()
        resize!(fields, fields_buff_len)
        new(true, FilePos(0), fields, DUMMY_WEAK_STRING)
    end
end

struct File{IO_TYPE}
    input::IO_TYPE    # the underlying input IO to be parsed
    delim::UInt8      # the delimiter should be an ASCII character to fit in a single byte
    quotechar::UInt8  # the quote character should be an ASCII character to fit in a single byte
    escapechar::UInt8 # the escape character should be an ASCII character to fit in a single byte
    header_exists::Bool      # does this file have a header?
    eager_parse_fields::Bool # should we parse the fields eagerly, or just parsing a line
                             # (i.e., a record) without digging into its fields is enough?
                             # This is useful for skipping records.
    lines_buff::Vector{UInt8} # A buffer that contains the line for the current record.
                             # All the `WeakRefString` instances point to a location in this
                             # `Vector`.
    recs::BufferedVector{Record} # a buffer for the fields of the current record
    current_line::Counter        # The counter that keeps the index of current line
                                 # (i.e., record) and it's used for making sure that all
                                 # iterations over this file are done sequentially.
                                 # Note that the counter value is a local value and it's not
                                 # shared between several workers in concurrent data loading.
    current_byte_pos::MutFilePos # The byte position of the current row (i.e., record)
    function File(input::IO_TYPE, delim::UInt8, quotechar::UInt8, escapechar::UInt8,
                  header_exists::Bool, eager_parse_fields::Bool, lines_buff::Vector{UInt8},
                  recs::BufferedVector{Record},
                  start_file_pos::Int64) where {IO_TYPE}
        new{IO_TYPE}(input, delim, quotechar, escapechar, header_exists,
                     eager_parse_fields, lines_buff, recs, Counter(0),
                     MutFilePos(start_file_pos))
    end
end

function File(input::IO; delim::Char=DEFAULT_DELIM,
              eager_parse_fields::Bool=DEFAULT_EAGER_PARSE_FIELDS,
              lines_buff_len::SizeType=DEFAULT_LINE_LEN, fields_buff_len::SizeType=DEFAULT_NUM_FIELDS,
              quotechar::Char=DEFAULT_QUOTE, escapechar::Char=quotechar,
              header_exists::Bool=DEFAULT_HEADER_EXISTS, start_file_pos::Int64=0,
              num_rows::SizeType=DEFAULT_NUM_ROWS)
    buff = Vector{UInt8}()
    resize!(buff, lines_buff_len)

    recs = BufferedVector{Record}()
    resize!(recs, num_rows)
    for i in 1:num_rows
        recs[i] = Record(fields_buff_len)
    end
    empty!(recs)

    File(input, UInt8(delim), UInt8(quotechar), UInt8(escapechar), header_exists,
        eager_parse_fields, buff, recs, start_file_pos)
end

function Base.iterate(f::File, state::Int = 1)
    if f.current_line.v + 1 == state
        next_pos_line = read_csv_line!(f)
        if next_pos_line === nothing
            nothing
        else
            f.current_line.v += 1
            (next_pos_line, state+1)
        end
    else
        #TODO: shall we throw an error here?
        # we are expecting the user to always scan sequentially
        nothing
    end
end

result_line_collection(f::File) = Vector{String}()
result_collection(f::File) = Vector{Vector{String}}()
first_rec(f::File) = f.recs[1]
first_rec_fields(f::File) = first_rec(f).fields_buff

function materialize_line(f::File, pos::FilePos, line)
    line_vec = result_line_collection(f)
    for field in first_rec_fields(f)
        push!(line_vec, string(field))
    end
    line_vec
end

function materialize(f::File, vec)
    counter = 0
    while (pos_line = read_csv_line!(f)) != nothing
        (pos, line) = pos_line
        push!(vec, materialize_line(f, pos, line))
        counter += 1
    end
    f.current_line.v = counter
    vec
end

"""
A helper function for parsing and materializing the while `File` into a `Vector`
"""
function materialize(f::File)
    materialize(f, result_collection(f))
end

function num_fields_for_current_line(f::File)
    length(first_rec_fields(f))
end

"""
A helper function for counting the number of lines in a `File`
"""
function count_lines(csv_file::File)
    counter = UInt(0)
    for pos_line in csv_file
        counter += 1
    end
    counter
end

"""
A helper function for counting the number of fields in a `File`
"""
function count_fields(csv_file::File)
    counter = UInt(0)
    for pos_line in csv_file
        counter += num_fields_for_current_line(csv_file)
    end
    counter
end

"""
A helper function for converting a CSV field value to string and writing it into an output buffer
"""
function csv_field_string(buff::IO, csv_file::File, field, i, force_quote::Bool=false)
    quotechar = Char(csv_file.quotechar)
    escapechar = Char(csv_file.escapechar)
    delim = Char(csv_file.delim)
    has_delim = occursin(delim, field)
    has_quote = occursin(quotechar, field)
    if has_delim || has_quote || force_quote
        if has_quote
            field_escape_quotes = replace(field, quotechar => "$escapechar$quotechar")
            write(buff, "$quotechar$field_escape_quotes$quotechar")
        else
            write(buff, "$quotechar$field$quotechar")
        end
    else
        write(buff, field)
    end
    if i < length(first_rec_fields(csv_file))
        write(buff, csv_file.delim)
    end
end

export File, FilePos, num_fields_for_current_line, count_lines, count_fields


text:
const Prime = UInt64
const BasePrime = UInt32
const BasePrimesArray = Array{BasePrime,1}
const SieveBuffer = Array{UInt8,1}

# contains a lazy list of a secondary base primes arrays feed
# NOT thread safe; needs a Mutex gate to make it so...
abstract type BPAS end # stands in for BasePrimesArrays, not defined yet
mutable struct BasePrimesArrays <: BPAS
    thunk :: Union{Nothing,Function} # problem with efficiency - untyped function!!!!!!!!!
    value :: Union{Nothing,Tuple{BasePrimesArray, BPAS}}
    BasePrimesArrays(thunk::Function) = new(thunk)
end
Base.eltype(::Type{BasePrimesArrays}) = BasePrime
Base.IteratorSize(::Type{BasePrimesArrays}) = Base.SizeUnknown() # "infinite"...
function Base.iterate(BPAs::BasePrimesArrays, state::BasePrimesArrays = BPAs)
    if state.thunk !== nothing
        newvalue :: Union{Nothing,Tuple{BasePrimesArray, BasePrimesArrays}} =
            state.thunk() :: Union{Nothing,Tuple{BasePrimesArray
                                                 , BasePrimesArrays}}
        state.value = newvalue
        state.thunk = nothing
        return newvalue
    end
    state.value
end

# count the number of zero bits (primes) in a byte array,
# also works for part arrays/slices, best used as an `@view`...
function countComposites(cmpsts::AbstractArray{UInt8,1})
    foldl((a, b) -> a + count_zeros(b), cmpsts; init = 0)
end

# converts an entire sieved array of bytes into an array of UInt32 primes,
# to be used as a source of base primes...
function composites2BasePrimesArray(low::Prime, cmpsts::SieveBuffer)
    limiti = length(cmpsts) * 8
    len :: Int = countComposites(cmpsts)
    rslt :: BasePrimesArray = BasePrimesArray(undef, len)
    i :: Int = 0
    j :: Int = 1
    @inbounds(
    while i < limiti
        if cmpsts[i >>> 3 + 1] & (1 << (i & 7)) == 0
            rslt[j] = low + i + i
            j += 1
        end
        i += 1
    end)
    rslt
end

# sieving work done, based on low starting value for the given buffer and
# the given lazy list of base prime arrays...
function sieveComposites(low::Prime, buffer::Array{UInt8,1},
                                     bpas::BasePrimesArrays)
    lowi :: Int = (low - 3) ÷ 2
    len :: Int = length(buffer)
    limiti :: Int = len * 8 - 1
    nexti :: Int = lowi + limiti
    for bpa::BasePrimesArray in bpas
        for bp::BasePrime in bpa
            bpint :: Int = bp
            bpi :: Int = (bpint - 3) >>> 1
            starti :: Int = 2 * bpi * (bpi + 3) + 3
            starti >= nexti && return
            if starti >= lowi starti -= lowi
            else
                r :: Int = (lowi - starti) % bpint
                starti = r == 0 ? 0 : bpint - r
            end
            lmti :: Int = limiti - 40 * bpint
            @inbounds(
            if bpint <= (len >>> 2) starti <= lmti
                for i in 1:8
                    if starti > limiti break end
                    mask = convert(UInt8,1) << (starti & 7)
                    c = starti >>> 3 + 1
                    while c <= len
                        buffer[c] |= mask
                        c += bpint
                    end
                    starti += bpint
                end
            else
                c = starti
                while c <= limiti
                    buffer[c >>> 3 + 1] |= convert(UInt8,1) << (c & 7)
                    c += bpint
                end
            end)
        end
    end
    return
end

# starts the secondary base primes feed with minimum size in bits set to 4K...
# thus, for the first buffer primes up to 8293,
# the seeded primes easily cover it as 97 squared is 9409.
function makeBasePrimesArrays() :: BasePrimesArrays
    cmpsts :: SieveBuffer = Array{UInt8,1}(undef, 512)
    function nextelem(low::Prime, bpas::BasePrimesArrays) ::
                                    Tuple{BasePrimesArray, BasePrimesArrays}
        # calculate size so that the bit span is at least as big as the
        # maximum culling prime required, rounded up to minsizebits blocks...
        reqdsize :: Int = 2 + isqrt(1 + low)
        size :: Int = (reqdsize ÷ 4096 + 1) * 4096 ÷ 8 # size in bytes
        if size > length(cmpsts) cmpsts = Array{UInt8,1}(undef, size) end
        fill!(cmpsts, 0)
        sieveComposites(low, cmpsts, bpas)
        arr :: BasePrimesArray = composites2BasePrimesArray(low, cmpsts)
        next :: Prime = low + length(cmpsts) * 8 * 2
        arr, BasePrimesArrays(() -> nextelem(next, bpas))
    end
    # pre-seeding breaks recursive race,
    # as only known base primes used for first page...
    preseedarr :: BasePrimesArray = # pre-seed to 100, can sieve to 10,000...
        [ 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41
        , 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97
        ]
    nextfunc :: Function = () ->
        (nextelem(convert(Prime,101), makeBasePrimesArrays()))
    firstfunc :: Function = () -> (preseedarr, BasePrimesArrays(nextfunc))
    BasePrimesArrays(firstfunc)
end

# an iterator over successive sieved buffer composite arrays,
# returning a tuple of the value represented by the lowest possible prime
# in the sieved composites array and the array itself;
# the array has a 16 Kilobytes minimum size (CPU L1 cache), but
# will grow so that the bit span is larger than the
# maximum culling base prime required, possibly making it larger than
# the L1 cache for large ranges, but still reasonably efficient using
# the L2 cache: very efficient up to about 16e9 range;
# reasonably efficient to about 2.56e14 for two Megabyte L2 cache = > 1 week...
struct PrimesPages
    baseprimes :: BasePrimesArrays
    PrimesPages() = new(makeBasePrimesArrays())
end
Base.eltype(::Type{PrimesPages}) = SieveBuffer
Base.IteratorSize(::Type{PrimesPages}) = Base.SizeUnknown() # "infinite"...
function Base.iterate(PP::PrimesPages,
                      state :: Tuple{Prime,SieveBuffer} =
                            ( convert(Prime,3), Array{UInt8,1}(undef,16384) ))
    (low, cmpsts) = state
    # calculate size so that the bit span is at least as big as the
    # maximum culling prime required, rounded up to minsizebits blocks...
    reqdsize :: Int = 2 + isqrt(1 + low)
    size :: Int = (reqdsize ÷ 131072 + 1) * 131072 ÷ 8 # size in bytes
    if size > length(cmpsts) cmpsts = Array{UInt8,1}(undef, size) end
    fill!(cmpsts, 0)
    sieveComposites(low, cmpsts, PP.baseprimes)
    newlow :: Prime = low + length(cmpsts) * 8 * 2
    ( low, cmpsts ), ( newlow, cmpsts )
end

function countPrimesTo(range::Prime) :: Int64
    range < 3 && ((range < 2 && return 0) || return 1)
    count :: Int64 = 1
    for ( low, cmpsts ) in PrimesPages() # almost never exits!!!
        if low + length(cmpsts) * 8 * 2 > range
            lasti :: Int = (range - low) ÷ 2
            count += countComposites(@view cmpsts[1:lasti >>> 3])
            count += count_zeros(cmpsts[lasti >>> 3 + 1] |
                                 (0xFE << (lasti & 7)))
            return count
        end
        count += countComposites(cmpsts)
    end
    count
end

# iterator over primes from above page iterator;
# unless doing something special with individual primes, usually unnecessary;
# better to do manipulations based on the composites bit arrays...
# takes at least as long to enumerate the primes as sieve them...
mutable struct PrimesPaged
    primespages :: PrimesPages
    primespageiter :: Tuple{Tuple{Prime,SieveBuffer},Tuple{Prime,SieveBuffer}}
    PrimesPaged() = let PP = PrimesPages(); new(PP, Base.iterate(PP)) end
end
Base.eltype(::Type{PrimesPaged}) = Prime
Base.IteratorSize(::Type{PrimesPaged}) = Base.SizeUnknown() # "infinite"...
function Base.iterate(PP::PrimesPaged, state::Int = -1 )
    state < 0 && return Prime(2), 0
    (low, cmpsts) = PP.primespageiter[1]
    len = length(cmpsts) * 8
    @inbounds(
    while state < len && cmpsts[state >>> 3 + 1] &
                         (UInt8(1) << (state & 7)) != 0
        state += 1
    end)
    if state >= len
        PP.primespageiter = Base.iterate(PP.primespages, PP.primespageiter[2])
        return Base.iterate(PP, 0)
    end
    low + state + state, state + 1
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "Toroidal tests" begin

  R0 = rand()
  Z0 = rand() - 0.5
  majorradius(x) = sqrt(sum(x[1:2].^2))
  toroidalangle(x) = atan(x[2], x[1])
  lateral(x) = majorradius(x) - R0
  axial(x) = x[3] - Z0
  minorradius(x) = sqrt(lateral(x)^2 + axial(x)^2)
  poloidalangle(x) = atan(axial(x), lateral(x))
  f⁻¹ = FromCartesian(xyz -> [minorradius(xyz),
                              toroidalangle(xyz),
                              poloidalangle(xyz)], 3)

  function r2x(rϕθ)
    r, ϕ, θ = rϕθ
    R = R0 + r * cos(θ)
    x, y, z = R * cos(ϕ), R * sin(ϕ), Z0 + r * sin(θ)
    return [x, y, z]
  end
  f = ToCartesian(r2x, 3)

  @testset "Length of CoordinateTransform" begin
    @test length(f⁻¹) == 3
  end

  @testset "Gradient of CoordinateTransform" begin
    for i ∈ 1:10
      xyz = rand(3) .* 2 .- 1
      x, y, z = xyz
      r, ϕ, θ = f⁻¹(xyz)
      @assert r ≈ sqrt((sqrt(x^2 + y^2) - R0)^2 + (z - Z0)^2)
      drdx = x * (sqrt(x^2 + y^2) - R0) / sqrt(x^2 + y^2) / r
      drdy = y * (sqrt(x^2 + y^2) - R0) / sqrt(x^2 + y^2) / r
      drdz = (z - Z0) / r
      dtdx = - y / (x^2 + y^2)
      dtdy = x / (x^2 + y^2)
      dtdz = 0.0
      dpdx = - x * (z - Z0) / sqrt(x^2 + y^2) / r^2
      dpdy = - y * (z - Z0) / sqrt(x^2 + y^2) / r^2
      dpdz = (sqrt(x^2 + y^2) - R0) / r^2
      @test J⃗(f⁻¹, [r, ϕ, θ], Covariant) ≈ [drdx drdy drdz;
                                      dtdx dtdy dtdz;
                                      dpdx dpdy dpdz]
    end
  end

  @testset "Metric Tensors" begin
    @testset "Contravariant Metric Tensors" begin
      for i ∈ 1:10
        xyz = rand(3) .* 2 .- 1
        x, y, z = xyz
        rϕθ = f⁻¹(xyz)
        r, ϕ, θ = rϕθ
        R = sqrt(sum(xyz[1:2].^2))
        gij = Array{Float64}(I, 3, 3)
        gij[2, 2] = 1/R^2
        gij[3, 3] = 1/r^2
        @test gⁱʲ(f, rϕθ) ≈ gij
      end
    end
    @testset "Contravariant Metric Tensors" begin
      for i ∈ 1:10
        xyz = rand(3) .* 2 .- 1
        x, y, z = xyz
        r, ϕ, θ = f⁻¹(xyz)
        R = sqrt(sum(xyz[1:2].^2))
        g_ij = Array{Float64}(I, 3, 3)
        g_ij[2, 2] = R^2
        g_ij[3, 3] = r^2
        @test gᵢⱼ(f⁻¹, [r, ϕ, θ]) ≈ g_ij
      end
    end
  end

  @testset "Jacobian" begin
    for i ∈ 1:10
      xyz = rand(3) .* 2 .- 1
      x, y, z = xyz
      rϕθ = f⁻¹(xyz)
      r, ϕ, θ = rϕθ
      R = R0 + r * cos(θ)
      @test J(f, rϕθ) ≈ r * R
      @test J(f⁻¹, rϕθ) ≈ r * R
    end
  end

  @testset "invert" begin
    for i ∈ 1:10
      xyz = rand(3) .* 2 .- 1
      rϕθ = f⁻¹(xyz)
      @test f⁻¹(f(rϕθ)) ≈ rϕθ
      @test f(f⁻¹(xyz)) ≈ xyz
      @test invert(f⁻¹, rϕθ) ≈ xyz
      @test f⁻¹ \ rϕθ ≈ xyz
    end
  end

  @testset "divergence" begin
      a, b = rand(3), rand(2:5, 3, 3)
      fv = z -> [a .* prod(z .^ b[:, i]) for i in 1:3]
      r = [2, 2pi, 2] .* rand(3) .- [0, pi, 1]
      x = f(r)
      j = ForwardDiff.jacobian(fv, r)
      jr = ForwardDiff.jacobian(x->fv(x) * x[1], r)
      R = sqrt(x[1]^2 + x[2]^2)
      divv = jr[1, 1] / r[1] + j[2,2] / R + j[3,3] / r[1]

      v = CovariantVector(fv, f, true)
      @test div(v)(r) ≈ divv
      v = ContravariantVector(fv, f, true)
      @test div(v)(r) ≈ divv
  end

  @testset "curl" begin
      a, b = rand(3), rand(2:5, 3, 3)
      fv = z -> [a .* prod(z .^ b[:, i]) for i in 1:3]
      r = [2, 2pi, 2] .* rand(3) .- [0, pi, 1]
      x = f(r)
      j = ForwardDiff.jacobian(fv, r)
      jr = ForwardDiff.jacobian(x->fv(x) * x[1], r)
      R = sqrt(x[1]^2 + x[2]^2)
      curl1 = (j[3, 2] - j[2, 3]) / r[1]
      curl2 = (j[1, 3] - jr[3, 1]) / r[1]
      curl3 = j[2, 1] - j[1, 2] / R
      curlv = [curl1, curl2, curl3]

      v = CovariantVector(fv, f, true)
      @test curl(v)(r, true) ≈ curlv
      v = ContravariantVector(fv, f, true)
      @test curl(v)(r, true) ≈ curlv
  end
end

max_stars_repo_name: jwscook/DifferentialGeometries

max_issues_repo_name: jwscook/DifferentialGeometries

max_forks_repo_name: jwscook/DifferentialGeometries


text:
"""
	sphtriangulate(cmd0::String="", arg1=[], kwargs...)

Delaunay or Voronoi construction of spherical lon,lat data

Full option list at [`sphtriangulate`](http://gmt.soest.hawaii.edu/doc/latest/sphtriangulate .html)

Parameters
----------

- **A** : **area** : -- Bool or [] --

    Compute the area of the spherical triangles (Qd) or polygons (Qv) and write the areas in the output segment headers
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/sphtriangulate.html#a)
- **C** : **save_mem** : -- Bool or [] --

    For large data sets you can save some memory (at the expense of more processing).
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/sphdistance.html#c)
- **D** : **skip** : -- Bool or [] --

    Skip the last (repeated) input vertex at the end of a closed segment if it equals the first point in the segment.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/sphtriangulate.html#d)
- **L** : **unit** : -- Str --          Flags = e|f|k|m|n|u|d

    Specify the unit used for distance and area calculations.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/sphtriangulate.html#l)
- **N** : -- Str --         Flags = ``nfile``

    Write the information pertaining to each polygon to a separate file.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/sphtriangulate.html#n)
- **Q** : **voronoi** : -- Str --     Flags = d|v

    Append d for Delaunay triangles or v for Voronoi polygons [Delaunay].
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/sphtriangulate .html#q)
- **T** : -- Bool or Str --

    Write the unique arcs of the construction [Default writes fillable triangles or polygons].
    When used with -A we store arc length in the segment header in chosen unit.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/sphtriangulate.html#t)
- $(GMT.opt_V)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_swap_xy)
"""
function sphtriangulate(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("sphtriangulate ", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd = parse_V_params("", d)
	cmd, = parse_b(cmd, d)
	cmd, = parse_d(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

    cmd = add_opt(cmd, 'A', d, [:A :area])
    cmd = add_opt(cmd, 'C', d, [:C :save_mem])
	cmd = add_opt(cmd, 'D', d, [:D :skip])
	cmd = add_opt(cmd, 'L', d, [:L :unit])
	cmd = add_opt(cmd, 'N', d, [:N])
	cmd = add_opt(cmd, 'Q', d, [:Q :voronoi])
	cmd = add_opt(cmd, 'T', d, [:T])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "sphtriangulate", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
sphtriangulate(arg1=[], cmd0::String=""; kw...) = sphtriangulate(cmd0, arg1; kw...)


text:
module WhenCases
using MLStyle
using MLStyle.Sugars: Q
using MLStyle.AbstractPatterns: init_cfg

export @when, @otherwise, gen_when
@nospecialize
function split_case_and_block(stmts, first_bindings, first_source)
    blocks::Vector{Any} = []
    binding_seqs::Vector{Any} = [first_bindings]
    sources::Vector{LineNumberNode} = [first_source]
    current_block = []

    function make_block!()
        # avoid setting LineNumberNode in the end of block.
        block_size = length(current_block)
        take_size = block_size
        for i = block_size:-1:1
            take_size = i
            current_block[i] isa LineNumberNode || break
        end
        # push current_block(cache) to blocks, then clear cache
        push!(blocks, Expr(:block, view(current_block, 1:take_size)...))
        empty!(current_block)
        nothing
    end

    function append_block!(a)
        push!(current_block, a)
        nothing
    end

    for stmt in stmts
        @switch stmt begin
            @case :(@when $source begin
                $(bindings...)
            end) || Q[@when $source $elt] && let bindings = [elt]
            end

            push!(binding_seqs, bindings)
            push!(sources, source)
            make_block!()
            continue
            @case Q[@otherwise $source]
            push!(binding_seqs, [])
            push!(sources, source)
            make_block!()
            continue
            @case a
            append_block!(a)
            continue
        end
    end

    make_block!()
    # thus, the length of `sources`, `binding_seqs` and `blocks`
    # are guaranteed to be the same.
    collect(zip(sources, binding_seqs, blocks))
end

"""
Used in the `@when` block:

```julia
@when (a, b) = x begin
    a + b
@otherwise
    0
end
```

See also: [`@when`](@ref)
"""
macro otherwise()
    throw(
        SyntaxError(
            "@otherwise is only used inside @when block, as a token to indicate default case.",
        ),
    )
end

"""
Code generation for `@when`.
You should pass an `Expr(:let, ...)` as the first argument.
"""
function gen_when(let_expr, source::LineNumberNode, mod::Module)
    @switch let_expr begin
        @case Expr(
            :let,
            Expr(:block, bindings...) || a && let bindings = [a]
            end,
            Expr(:block, stmts...) || b && let stmts = [b]
            end,
        )
        sources_cases_blocks = split_case_and_block(stmts, bindings, source)
        return foldr(
            sources_cases_blocks,
            init = :nothing,
        ) do (source, bindings, block), last_block
            foldr(bindings, init = block) do each, last_ret
                @switch each begin
                    @case :($a = $b)
                    cbl = Expr(:block, source, :($a => $last_ret), :(_ => $last_block))
                    return gen_match(b, cbl, source, mod)

                    @case new_source::LineNumberNode
                    source = new_source
                    return last_ret

                    # match `cond.?` or `if cond end`
                    @case :(
                        if $a
                            $(_...)
                        end
                    ) || :($a.?)

                    return Expr(:block, source, :($a ? $last_ret : $last_block))

                    # match something like: `let a; a = 1; a end`
                    @case a
                    return Expr(:block, source, Expr(:let, a, last_ret))
                end
            end
        end

        @case a
        s = string(a)
        short_msg = SubString(s, 1, min(length(s), 20))
        throw(
            SyntaxError(
                "Expected a let expression, got a `$short_msg` at $(string(source)).",
            ),
        )
    end
end

"""
1. Allow destructuring in binding sequences of `let` syntax.

In binding sequences,
- For bindings in form of `a = b` and `f(x) = y`, it's regarded as pattern matching here.
- For others like `@inline f(x) = 1`, it's the same as the original let binding(not pattern matching).

```julia
    @when let (a, 1) = x,
          [b, c, 5] = y
        (a, b, c)
    end
```

2. For a regular assignment, like

```julia
    @when (a, 2) = x begin
        # dosomething
    end
```

It's nothing different from

```julia
    @match x begin
        (a, 2) => # dosomething
        _ => nothing
    end
```

3. For multiple branches, we can have

```julia
x = 1
@when let (_, _) = x
    :tuple
@when begin ::Float64 = x end
    :float
@when ::Int = x
    :int
@otherwise
    :unknown
end
```

4. Also, you can use predicates when matching, in the form
of `cond.?` or `if cond end`:

```julia
x = 1
y = (1, 2)
cond1 = true
cond2 = true
@when let cond1.?,
          (a, b) = x
    a + b
@when begin if cond2 end
            (a, b) = y
      end
    a + b
end
```
"""
macro when(let_expr)
    res = gen_when(let_expr, __source__, __module__)
    res = init_cfg(res)
    esc(res)
end

macro when(assignment, ret)
    @match assignment begin
        :($_ = $_) => let let_expr = Expr(:let, Expr(:block, assignment), ret)
            res = gen_when(let_expr, __source__, __module__)
            res = init_cfg(res)
            esc(res)
        end
        _ => throw(SyntaxError("Not match the form of `@when a = b expr`"))
    end
end
@specialize

end


text:
"""
Abstract distance type for placement.

API
---
* [`getdistance`](@ref)
* [`maxdistance`](@ref)

Basic Implementation
--------------------
* [`BasicDistance`](@ref)
"""
abstract type SADistance end

##################
# SADistance API #
##################

"""
    getdistance(A::SADistance, a, b)

Return the distance from `a` to `b` with respect to `A`. Both `a` and `b` will
implement [`getaddress`](@ref).

Method List
-----------
$(METHODLIST)
"""
function getdistance end

"""
    maxdistance(sa_struct, A::SADistance)

Return the maximum distance value of that occurs in `sa_struct` using the
distance metric imposed by `A`.

Method List
-----------
$(METHODLIST)
"""
function maxdistance end

################################################################################
# Basic Distance
################################################################################

"""
    $(SIGNATURES)

Basic implementation of [`SADistance`](@ref). Constructs a look up table of
distances between all address pairs. Addresses have a distance of 1 if there
exists even one link between components at those addresses.
"""
struct BasicDistance{D} <: SADistance
    """
    Simple look up table indexed by pairs of addresses. Returned value is the
    distance between the two addresses.

    If an architecture has dimension "N", then the dimension of `table` is `2N`.
    """
    table::Array{UInt8,D}
end

# Unwrap and
@propagate_inbounds @inline function getdistance(A::BasicDistance, a, b)
    # Constant propgation doesn't take this all the way for higher dimensions
    return A.table[Tuple(location(a))..., Tuple(location(b))...]
end
maxdistance(sa_struct, A::BasicDistance) = maximum(A.table)

################################################################################
# BFS Routines for building the distance look up table
################################################################################
function BasicDistance(toplevel::TopLevel, pathtable::PathTable)
    # The data type for the LUT
    element_type = UInt8
    # Pre-allocate the table. Get the size of the `toplevel` to figure out
    # how big to make the table.
    #
    # Since we want to get the distance between pairs of addresses, we need
    # to double size of `toplevel`.
    dims = size(toplevel)
    pretable = fill(typemin(element_type), dims..., dims...)

    # Get the adjacency information from toplevel.
    neighbors = neighbor_dict(toplevel)

    @debug "Building Distance Table"

    # Get all-pairs shorest path from the adjacency lists. For each pair
    # of addresses "a", "b", "table[a,b]" will be the distance from "a" to "b"
    for address in keys(neighbors)
        bfs!(pretable, address, neighbors)
    end

    table = expand(pretable, toplevel, pathtable)
    return BasicDistance(table)
end

expand(pretable::Array{T}, ::TopLevel{D}, ::PathTable{D}) where {T,D} = pretable

function expand(
    pretable::Array{T}, toplevel::TopLevel{D}, pathtable::PathTable{N}
) where {T,D,N}

    # Just Replicate the table.
    fullsize = size(pathtable)
    table = fill(typemin(T), fullsize..., fullsize...)
    table = map(CartesianIndices((fullsize..., fullsize...))) do idx
        a, b = half(Tuple(idx))
        return pretable[Base.tail(a)..., Base.tail(b)...]
    end
    return table
end
half(x::NTuple{N}) where {N} = x[1:(N >> 1)], x[((N >> 1) + 1):N]

#=
Simple data structure for keeping track of costs associated with addresses
Gets put the the queue for the BFS.
=#
struct CostAddress{U,D}
    cost::U
    address::Address{D}
end

# Implementation note:
# This function is only correct if the cost of each link is 1. If cost can vary,
# will have to code this using some kind of shortest path formulation.
function bfs!(distance::Array{U,N}, source::Address{D}, neighbors) where {U,N,D}
    # Create a queue for visiting addresses. Add source to get into the loop.
    q = Queue{CostAddress{U,D}}()
    enqueue!(q, CostAddress(zero(U), source))

    # Create a set of visited items to avoid visiting the same address twice
    # and mark the initial address as "seen".
    seen = Set{Address{D}}()
    push!(seen, source)

    # Basic BFS.
    while !isempty(q)
        u = dequeue!(q)
        distance[source, u.address] = u.cost
        for v in neighbors[u.address]
            in(v, seen) && continue

            enqueue!(q, CostAddress(u.cost + one(U), v))
            push!(seen, v)
        end
    end

    return nothing
end

function neighbor_dict(toplevel::TopLevel)
    @debug "Building Neighbor Table"
    # Get the connected component dictionary
    cc = MapperCore.connected_components(toplevel)

    # Addresses may be zero or negative, but since we're using an array as our
    # underlying data structure, we have to shift the whole array to addresses
    # start at 1.
    #
    # Get this offset amount and apply it to the array.
    #
    # Wrap the offset in a Ref to ensure the broadcast treats it as a scalar.
    offset = getoffset(toplevel)
    return Dict(a + offset => collect(s) .+ Ref(offset) for (a, s) in cc)
end


text:
@testset "problems" begin

include("1.two-sum.jl")
include("240.search-a-2d-matrix-ii.jl")

end


text:
module ode_dopri

export ode, ode_event

using Roots

const order = 5

# Runge-Kutta table
const a2 = Float64[1/5]
const a3 = Float64[3/40, 9/40]
const a4 = Float64[44/45, -56/15, 32/9]
const a5 = Float64[19372/6561, -25360/2187, 64448/6561, -212/729]
const a6 = Float64[9017/3168, -355/33, 46732/5247, 49/176, -5103/18656]
const a7 = Float64[35/384, 0., 500/1113, 125/192, -2187/6784, 11/84]
const a = Vector{Float64}[Float64[], a2, a3, a4, a5, a6, a7]

# Weights
const bs = Float64[5179/57600, 0., 7571/16695, 393/640, -92097/339200, 187/2100, 1/40]

# Nodes
const c = Float64[0., 1/5, 3/10, 4/5, 8/9, 1., 1.]

b7(θ) = (θ^2*(θ-1) + θ^2*(θ-1)^2*10*(7414447 - 829305*θ)/29380423)
b1(θ) = θ^2*(3-2*θ)*a7[1] + θ*(θ-1)^2 - θ^2*(θ-1)^2*5*(2558722523-31403016*θ)/11282082432
b3(θ) = θ^2*(3-2*θ)*a7[3] + θ^2*(θ-1)^2 * 100*(882725551 - 15701508*θ)/32700410799
b4(θ) = θ^2*(3-2*θ)*a7[4] - θ^2*(θ-1)^2 * 25*(443332067 - 31403016*θ)/1880347072
b5(θ) = θ^2*(3-2*θ)*a7[5] + θ^2*(θ-1)^2 * 32805*(23143187 - 3489224*θ)/199316789632
b6(θ) = θ^2*(3-2*θ)*a7[6] - θ^2*(θ-1)^2 * 55*(29972135 - 7076736*θ)/822651844

"""
Intermediate function value for the nth step.

Arguments:
----------
    * x: End value.
    * x0: Starting value.
    * h: Step size.
    * coeffs: Runge-kutta coefficents for the nth substep.
    * k: Function derivatives at the previous substep points.
"""
function substep{T}(x::Vector{T}, x0::Vector{T}, h::Float64, coeffs::Vector{Float64}, k::Vector{Vector{T}})
    @inbounds for m=1:length(x0)
        dx::T = 0.
        @inbounds for i=1:length(coeffs)
            dx += coeffs[i]::Float64 * k[i][m]
        end
        x[m] = x0[m] + h*dx
    end
    return nothing
end

"""
Perform one single Runge-Kutta step.

Arguments:
----------
    * F: Derivative function with signature F(t, y, dy) which writes the derivative into dy.
    * t0: Start time of the step.
    * h: Length of the step.
    * x0: Initial state.
    * xp: Final state after the Runge-Kutta step. (5th order)
    * xs: Final state after the Runge-Kutta step. (4th order)
    * k: States at the substeps.
"""
function step{T}(F::Function, t0::Float64, h::Float64,
                x0::Vector{T}, xp::Vector{T}, xs::Vector{T}, k::Vector{Vector{T}})
    for i=2:length(c)
        substep(xp, x0, h, a[i], k)
        F(t0 + h*c[i], xp, k[i])
    end
    substep(xs, x0, h, bs, k)
    return nothing
end

"""
Allocate necessary memory.
"""
function allocate_memory{T}(x::Vector{T})
    xp = zeros(T, length(x))
    xs = zeros(T, length(x))
    k = Vector{T}[]
    for i=1:7
        push!(k, zeros(T, length(x)))
    end
    return (xp, xs, k)
end

"""
Calculate state for intermediate steps (Dense output).

Arguments:
----------
    * t0: Start time of the step.
    * x0: Initial state.
    * h: Length of the step.
    * k: States at the substeps.
    * t: Time at which the state should be interpolated.
    * x: Resulting state at the given point of time.
"""
function interpolate{T}(t0::Float64, x0::Vector{T}, h::Float64, k::Vector{Vector{T}}, t::Float64, x::Vector{T})
    θ = (t-t0)/h
    b1_ = b1(θ); b3_ = b3(θ); b4_ = b4(θ); b5_ = b5(θ); b6_ = b6(θ); b7_ = b7(θ)
    for i=1:length(x0)
        dx = b1_*k[1][i] + b3_*k[3][i] + b4_*k[4][i]
        dx += b5_*k[5][i] + b6_*k[6][i] + b7_*k[7][i]
        x[i] = x0[i] + h*dx
    end
    nothing
end

"""
Estimate the error by comparing results of different order Runge-Kutta methods.
"""
function error_estimate(xp, xs, abstol, reltol)
    err::Float64 = 0.
    for i=1:length(xp)
        sc_i = abstol + reltol*max(abs(xp[i]), abs(xs[i]))
        err += abs2(xp[i]-xs[i])/sc_i^2
    end
    return sqrt(err/length(xp))
end

"""
Find a good guess for an initial time step.
"""
function initial_stepsize(F, t, x, k, abstol, reltol, tmp1, tmp2)
    d0 = 0.
    d1 = 0.
    for i=1:length(x)
        sc_i2 = (abstol + abs(x[i])*reltol)^2
        d0 += abs2(x[i])/sc_i2
        d1 += abs2(k[1][i])/sc_i2
    end
    d0 = sqrt(d0/length(x))
    d1 = sqrt(d1/length(x))
    h0 = ((d0<1e-5 || d1<1e-5) ? 1e-6 : 0.01*d0/d1)
    substep(tmp1, x, h0, [1.], k)
    F(t+h0, tmp1, tmp2)
    for i=1:length(tmp2)
        tmp2[i] -= k[1][i]
    end
    d2 = norm(tmp2, 2)
    if max(d1,d2)<1e-15
        h1 = max(1e-6, h0*1e-3)
    else
        h1 = (0.01/max(d1,d2))^(1./order)
    end
    return min(100*h0, h1)
end

"""
Determine the size of the next step.
"""
function stepsize_strategy(err, laststepaccepted, h, hmin, hmax)
    accept_step = err<1
    facmin = (laststepaccepted ? 5. : 1.)
    hnew = h*min(facmin, max(0.2, 0.9*(1./err)^(1./order)))
    hnew = min(hmax, hnew)
    if hnew<hmin
        error("Stepsize below hmin.")
    end
    return hnew, accept_step
end

"""
Display states at the given times which are in the current Runge-Kutta step.
"""
function display_steps{T}(fout::Function, tspan::Vector{Float64}, t::Float64, x::Vector{T}, h::Float64, k::Vector{Vector{T}}, xs::Vector{T})
    for tout=tspan
        if t<tout<=t+h
            interpolate(t, x, h, k, tout, xs)
            fout(tout, xs)
        end
    end
end

@enum CallbackCommand nojump jump stop

"""
Adaptive Runge-Kutta Dormand-Prince 4(5) solver with event handling and dense output.

Arguments:
----------
    * F: Derivative function with signature F(t, y, dy) which writes the derivative into dy.
    * tspan: Vector of times at which output should be displayd.
    * x0: Initial state.
    * event_locator: Function used to find events with signature
            event_locator(t, x) returning a real value. If the sign of the returned value
            changes the event_callback function is called.
    * event_callback: Function that is called when an event happens. Its
            signature is event_callback(t, x) and it should return a CallbackCommand.
            The possible CallBack commands are:
                - nojump: No changes in the dynamics. In this case x should not be changed
                          inside the callback function.
                - jump: the x vector has changed and time evolution continues from t_event.
                - stop: The ode solver stops ia the event time.

Optional Arguments:
-------------------
    * fout: Function called to display the state at the given points of time
            with signature fout(t, x). If no function is given, ode_event returns
            a vector with the states at all points in time given in tspan.
    * reltol: Relative error tolerance.
    * abstol: Absolute error tolerance.
    * h0: Initial guess for the size of the time step. If no number is given an
            initial timestep is chosen automatically.
    * hmin: If the automatic stepsize goes below this limit the ode solver stops
            with an error.
    * hmax: Stepsize is never increased above this limit.
    * display_initialvalue: Call fout function at tspan[1].
    * display_finalvalue: Call fout function at tspan[end].
    * display_intermediatesteps: Call fout function after every Runge-Kutta step.
    * display_beforeevent: Call fout function immediately before an event.
    * display_afterevent: Call fout function immediately after an event.
"""
function ode_event{T}(F, tspan::Vector{Float64}, x0::Vector{T}, fout::Function,
                    event_locator::Function, event_callback::Function;
                    reltol::Float64 = 1.0e-6,
                    abstol::Float64 = 1.0e-8,
                    h0::Float64 = NaN,
                    hmin::Float64 = (tspan[end]-tspan[1])/1e9,
                    hmax::Float64 = (tspan[end]-tspan[1]),
                    display_initialvalue::Bool = true,
                    display_finalvalue::Bool = true,
                    display_intermediatesteps::Bool = false,
                    display_beforeevent::Bool = false,
                    display_afterevent::Bool = false
                    )
    t, tfinal = tspan[1], tspan[end]
    display_initialvalue && fout(t, x0)

    # Allocate memory
    x = copy(x0) # Dont change initial state.
    xp, xs, k = allocate_memory(x0)

    # Initial derivative.
    F(t,x,k[1])

    # Find initial step size.
    h = (h0===NaN ? initial_stepsize(F, t, x, k, abstol, reltol, k[2], k[3]) : h0)
    h = max(hmin, h)
    h = min(hmax, h)

    accept_step = true
    while t < tfinal
        # Try Runge-Kutta steps until the error is below the tolerances.
        step(F, t, h, x, xp, xs, k)
        err = error_estimate(xp, xs, abstol, reltol)
        hnew, accept_step = stepsize_strategy(err, accept_step, h, hmin, hmax)
        if accept_step
            # Check if there is an event in the last time step.
            e1 = event_locator(t+hmin,x)
            e2 = event_locator(t+h,xp)
            if e2==0. || e1*e2 < 0. # Handle event case.
                t_event = fzero(t_->(interpolate(t, x, h, k, t_, xs); event_locator(t_, xs)), t, t+h)
                display_steps(fout, tspan[tspan.<t_event], t, x, h, k, xs)
                interpolate(t, x, h, k, t_event, xs)
                display_beforeevent && fout(t_event, xs)
                cmd = event_callback(t_event, xs)
                if typeof(cmd)!=CallbackCommand
                    error("Event callback function has to return a CallbackCommand.")
                end
                display_afterevent && fout(t_event, xs)
                if cmd == stop
                    return nothing
                elseif cmd == jump
                    F(t,x,k[1])
                    h = 0.9*hnew
                    xs, x = x, xs
                    t = t_event
                    continue
                elseif cmd == nojump
                    display_steps(fout, tspan[tspan.>t_event], t, x, h, k, xs)
                    display_intermediatesteps && fout(t, xp)
                    xp, x = x, xp
                    k[1], k[end] = k[end], k[1]
                    t = t + h
                else
                    error("Unrecognized event command.")
                end
            else # Handle no-event case.
                display_steps(fout, tspan, t, x, h, k, xs)
                if display_intermediatesteps && t+h<tfinal
                    fout(t+h, xp)
                end
                xp, x = x, xp
                k[1], k[end] = k[end], k[1]
                t = t + h
            end
        end
        h = hnew
    end
end

function ode_event{T}(F, tspan::Vector{Float64}, x0::Vector{T},
                    event_locator::Function, event_callback::Function;
                    args...)
    tout = Float64[]
    xout = Vector{T}[]
    fout = (t, x) -> (push!(tout, t); push!(xout, deepcopy(x)); nothing)
    ode_event(F, tspan, x0, fout, event_locator, event_callback; args...)
    return tout, xout
end

"""
Adaptive Runge-Kutta Dormand-Prince 4(5) solver with dense output.

Arguments:
----------
    * F: Derivative function with signature F(t, y, dy) which writes the derivative into dy.
    * tspan: Vector of times at which output should be displayd.
    * x0: Initial state.

Optional Arguments:
-------------------
    * fout: Function called to display the state at the given points of time
            with signature fout(t, x). If no function is given, ode_event returns
            a vector with the states at all points in time given in tspan.
    * reltol: Relative error tolerance.
    * abstol: Absolute error tolerance.
    * h0: Initial guess for the size of the time step. If no number is given an
            initial timestep is chosen automatically.
    * hmin: If the automatic stepsize goes below this limit the ode solver stops
            with an error.
    * hmax: Stepsize is never increased above this limit.
    * display_initialvalue: Call fout function at tspan[1].
    * display_finalvalue: Call fout function at tspan[end].
    * display_intermediatesteps: Call fout function after every Runge-Kutta step.
"""
function ode{T}(F, tspan::Vector{Float64}, x0::Vector{T}, fout::Function;
                    reltol::Float64 = 1.0e-6,
                    abstol::Float64 = 1.0e-8,
                    h0::Float64 = NaN,
                    hmin::Float64 = (tspan[end]-tspan[1])/1e9,
                    hmax::Float64 = (tspan[end]-tspan[1]),
                    display_initialvalue::Bool = true,
                    display_finalvalue::Bool = true,
                    display_intermediatesteps::Bool = false,
                    )
    t, tfinal = tspan[1], tspan[end]
    display_initialvalue && fout(t, x0)

    # Allocate memory
    x = copy(x0) # Dont change initial state.
    xp, xs, k = allocate_memory(x0)

    # Initial derivative.
    F(t,x,k[1])

    # Find initial step size.
    h = (h0===NaN ? initial_stepsize(F, t, x, k, abstol, reltol, k[2], k[3]) : h0)
    h = max(hmin, h)
    h = min(hmax, h)

    accept_step = true
    while t < tfinal
        # Try Runge-Kutta steps until the error is below the tolerances.
        step(F, t, h, x, xp, xs, k)
        err = error_estimate(xp, xs, abstol, reltol)
        hnew, accept_step = stepsize_strategy(err, accept_step, h, hmin, hmax)
        if accept_step
            display_steps(fout, tspan, t, x, h, k, xs)
            if display_intermediatesteps && t+h<tfinal
                fout(t+h, xp)
            end
            xp, x = x, xp
            k[1], k[end] = k[end], k[1]
            t = t + h
        end
        h = hnew
    end
end

function ode{T}(F, tspan::Vector{Float64}, x0::Vector{T}; args...)
    tout = Float64[]
    xout = Vector{T}[]
    fout = (t, x) -> (push!(tout, t); push!(xout, deepcopy(x)); nothing)
    ode(F, tspan, x0, fout; args...)
    return tout, xout
end

end # module


text:
include("hilbertcurve/hilbert_curve.jl")

function resample(s, w)
  # Samples new particles based on their weights. If you find algorithmic optimizations to this rou$
  N = length(w)
  bins = [0.; cumsum(exp(w))]
  j = zeros(Int64,N)
  bo = 1
  for i = 1:N
    for b = bo:N
      #WHY CAN I NOT SKIP THE RIGHT PART!?
      if bins[b] <= s[i] < bins[b+1]
        j[i] = b
        bo = b
        break
      end
    end
  end
  return j
end

function resample(w)
  N = length(w)
  bins = [0.; cumsum(exp(w))]
  s = collect((rand()/N+0):1/N:bins[end])
  resample(s,w)
end

""" `xhat = SQMC(f, g, ginv, x0, y, est, N=200, Nu=lenth(x0))`

Gives estimation of the states `x(:,1),...x(:,T)` given outputs `y(:,1)...y(:,T)` from system:

x(t) = f(x(t-1),t,u) \n
y = gn(x,v),\n
where `g(x,t)=E[gn(x,t,u)]`, `ginv(yhat, y, t) = p(y=gn(x,u)|yhat=g(x))`,
u is uniform noise on `[0,1)`, with `E[x(0)] = xhat0`,
and some estimator `est(xₚ,w,t)` that ouputs the estimate `xhat(:,t)`
given the particles `xₚ` and weights `w`
"""
function SQMC(f!, g!, ginv, randsf, xhat0, y, est; N = 256, Nu = length(xhat0), debug = false, xreal=xreal)
  aux= false
  normal = false
  if isa(randsf, Tuple)
    rands! = randsf[1]
    if randsf[2] == :aux
      aux = true
    end
    if length(randsf)>2
      if randsf[3] != :rand && randsf[3] != :randn
        error("Only allowing :rand or :randn as rands argument")
      elseif randsf[3] == :randn
        normal = true
      end
    end
  else
    rands! = randsf
  end
  if !isinteger(log(2,N))
    error("Number of particles must be a power of 2")
  end
  if debug
    pdata = Void
  end
  T = size(y,2)
  #Estimates at time 1:T
  xhat = Array{Float64,2}(length(xhat0), T)
  #Current predicted output for each particle
  yhat = Array{Float64,2}(size(y,1), N)
  x = Array{Float64,2}(length(xhat0), N)
  hilbertIndex = Array{Int64}(size(x,2))
  xtemp = similar(x)
  w = Array{Float64,1}(N)
  #x[:,1] = xhat0; #Hack to use the pfStep fuction. Works by letting `a` be ones
  x = repmat(xhat0,1,N);

  #Step (a) (Draw u and x), OBS u will be of size NxNu
  u, = sobol(Nu,N)
  if normal
    unit2Normal!(u)
  end
  if !aux
    rands!(u)
  else
    rands!(u,y,1)
  end

  if debug
    p = Void
    xOld = similar(x)
    xOld[:] = x[:]
  end
  pfStep!(f!,g!,ginv,x,xtemp,u,ones(Int,N),1:N,1,w,yhat,y,N, 1:Nu)
  xhat[:,1] = est(x,w,1)
  if debug
    #p = plotPoints(x, w, y, yhat, N, ones(Int,N), 1, xreal, xhat, xOld, p)
    p = pploti(x, exp(w)./sum(exp(w)), y, yhat, N, ones(Int,N), 1, xreal, xhat, xOld, p, density = true,leftOnly=true, xIndices = [1], yIndices = [],  slidef=1.0)
    xOld[:] = x[:]
  end
  #u = Array{Float64,2}(N,Nu+1)
  #nextseed, MeM = sobol!(u,Nu+1,N)
  u, nextseed, MeM = sobol(Nu+1,N)
  #The first column in u is for resample purpose
  Nus = 2:(Nu+1)
  for t = 2:T
    #nextseed = sobol!(u, nextseed, MeM)
    nextseed = sobol!(u, nextseed, MeM)
    if normal
      unit2Normal!(slice(u,:,Nus))
    end
    if !aux
      rands!(slice(u,:,Nus))
    else
      rands!(slice(u,:,Nus),y,t)
    end
    τ = sortperm(slice(u,:,1))
    σ = sortperm(slice(ψ!(x,hilbertIndex),:))
    a = resample(slice(u,τ,1), slice(w,σ))
    # Time update, no sigma here?
    #if randsf[4]
    #Make this faster
      x = x[:,σ]
      w = w[σ]
    #end
    pfStep!(f!,g!,ginv,x,xtemp,u,a,τ,t,w,yhat,y,N, Nus)
    xhat[:,t] = est(x,w,1)
    if debug
      #p = plotPoints(x, w, y, yhat, N, a, t, xreal, xhat, xOld, p)
      p = pploti(x, exp(w)./sum(exp(w)), y, yhat, N, a, t, xreal, xhat, xOld, p, density = true,leftOnly=true, xIndices = [1], yIndices = [],  slidef=1.0)
      xOld[:] = x[:]
    end
  end
  xhat
end

function unit2Normal!(u)
  for j = 1:size(u,2)
    for i = 1:size(u,1)
      u[i,j] = sqrt(2)*erfinv(2*u[i,j]-1);
    end
  end
end

#Should only be used for when 0<x<1
function ψUnit!(x,hbtIdx)
  #Figure out a resonable resolution (Number of indices has to be lower than max int)
  ncoords = size(x,1)
  N=floor(Int64,2^(60/ncoords))-1 #The resolution in each direction
  #Vector to keep the integer (scaled) value for each coordinate
  xInts = Array{Int64}(size(x,1))
  for j = 1:size(x,2)
    for i = 1:size(x,1)
      xInts[i] = fld(x[i,j]*N,1/N)
    end
    hbtIdx[j] = hilbert2int(xInts)
  end
  hbtIdx
end

function ψ!(x,hbtIdx)
  #Figure out a resonable resolution (Number of indices has to be lower than max int)
  ncoords = size(x,1)
  N=floor(Int64,2^(60/ncoords))-1 #The resolution in each direction
  #Find a reasonable scale
  maxVal = maximum(x,2)
  minVal = minimum(x,2)
  #Only use incices for which max!=min in sorting
  idx = find((v)->v!=0,maxVal-minVal)
  length(idx) == 0 && error("Error in sortnig, all values equal")
  #Vector to keep the integer (scaled) value for each coordinate
  xInts = Array{Int64}(length(idx))
  for j = 1:size(x,2)
    for (i,val) in enumerate(idx)
      xInts[i] = fld((x[val,j]-minVal[val])*N,maxVal[val]-minVal[val])
    end
    hbtIdx[j] = hilbert2int(xInts)
  end
  hbtIdx
end

#Randsf is either a function that fills the u matrix (only argument) with random numebrs or a tuple with a function and a symbol.
#If symbol is :aux then the function should accept the matrix as well as the measurement and time as arguments
function pf(f!, g!, ginv, randsf, xhat0, y, est; N = 200, Nu = length(xhat0), debug=false, xreal=xreal)
  aux= false
  if isa(randsf, Tuple)
    rands! = randsf[1]
    if randsf[2] == :aux
      aux = true
    end
  else
    rands! = randsf
  end
  T = size(y,2)

  #Estimates at time 1:T
  xhat = Array{Float64,2}(length(xhat0), T)
  #Current predicted output for each particle
  yhat = Array{Float64,2}(size(y,1), N)
  x = Array{Float64,2}(length(xhat0), N)
  xtemp = similar(x)
  w = Array{Float64,1}(N)
  τ = 1:N
  x = repmat(xhat0,1,N); #Hack to use the pfStep fuction. Works by letting `a` be ones

  u = Array{Float64,2}(N,Nu)
  if !aux
    rands!(u)
  else
    rands!(u,y,1)
  end

  if debug
    p = Void
    xOld = similar(x)
    xOld[:] = x[:]
  end
  pfStep!(f!,g!,ginv,x,xtemp,u,ones(Int,N),τ,1,w,yhat,y,N, 1:Nu)
  xhat[:,1] = est(x,w,1)
  if debug
    #p = plotPoints(x, w, y, yhat, N, ones(Int,N), 1, xreal, xhat, xOld, p)
    p = pploti(x, exp(w)./sum(exp(w)), y, yhat, N, ones(Int,N), 1, xreal, xhat, xOld, p, density = true, leftOnly=true, xIndices = [1], yIndices = [], slidef=1.0)
    xOld[:] = x[:]
  end
  for t = 2:T
    #Avoiding SQMC resample because of gaussian noise (change back?)
    a = resample(w)
    #a = resample(slice(u,τ,1), slice(w,:))
    #Rand
    if !aux
      rands!(u)
    else
      rands!(u,y,t)
    end
    # Time update, no sigma here?
    pfStep!(f!,g!,ginv,x,xtemp,u,a,τ,t,w,yhat,y,N,1:Nu)
    xhat[:,t] = est(x,w,1)
    if debug
      #p = plotPoints(x, w, y, yhat, N, a, t, xreal, xhat, xOld, p)
      p = pploti(x, exp(w)./sum(exp(w)), y, yhat, N, a, t, xreal, xhat, xOld, p, density = true,leftOnly=true, xIndices = [1], yIndices = [],  slidef=1.0)
      xOld[:] = x[:]
    end
  end
  xhat
end

#Time update and weighting. Chances `x`, `yhat` and `w`
function pfStep!(f!,g!,ginv,x,xtemp,u,a,τ,t,w,yhat,y,N,uRange)
  for i = 1:N
    xS = slice(x, :, a[i])
    uS = slice(u, τ[i], uRange)
    xtS = slice(xtemp,:,i)
    f!(xS, t, uS, xtS, i)
    g!(xtS, t, slice(yhat,:,i))
    #println(ginv(slice(yhat,:,i), slice(y,:,t), t))
    w[i] = ginv(slice(yhat,:,i), slice(y,:,t), t)
  end
  x[:] = xtemp[:]
  offset = maximum(w)
  mySum = sum(exp(w-offset))
  normConstant = log(mySum)+offset
  w[:] = w[:] - normConstant
end

function generateRealSequence(f!, gn!, x0, T, rands!, Ny, Nu)
  x = Array{Float64,2}(length(x0), T)
  xt = Array{Float64,2}(length(x0), 1)
  u = Array{Float64,1}(Nu)
  y = Array{Float64,2}(Ny, T)
  rands!(u, 1)
  f!(x0, 1, u,xt,1)
  x[:,1] = xt[:]
  gn!(slice(x,:,1), 1, u, slice(y,:,1))
  for t = 2:T
    rands!(u, t)
    f!(x[:,t-1], t, u,xt,t)
    x[:,t] = xt[:]
    gn!(slice(x,:,t), t, u, slice(y,:,t))
  end
  x, y
end

type FilterSettings
  f::Function
  g::Function
  gn::Function
  ginv::Function
  Ny::Int
  randsfGen
  randsfSim
  x0::Vector
  NuGen::Int
  NuSim::Int
  est
end

function runTest(N, method, s, T, debug = false)
  x, y = generateRealSequence(s.f, s.gn, s.x0, T, s.randsfGen, s.Ny, s.NuGen)
  xhat = method(s.f, s.g, s.ginv, s.randsfSim, s.x0, y, s.est,
  N = N, Nu = s.NuSim, debug = debug, xreal = x)
  x, xhat
end


text:
@elapsed a = rand(10 ^ 6)

@macroexpand @elapsed a = rand(10 ^ 6)

@allocated a = rand(10 ^ 6)

@macroexpand @allocated a = rand(10 ^ 6)

@time a = rand(10 ^ 6)

@macroexpand @time a = rand(10 ^ 6)

@timev a = rand(10 ^ 6)

@macroexpand @timev a = rand(10 ^ 6)

max_stars_repo_name: phrb/PDP

max_issues_repo_name: phrb/PDP

max_forks_repo_name: phrb/PDP


text:
function get_all_files(cur_item; package_name=nothing, is_testing=false, is_sorted=false)
  nested_files = get_nested_files(cur_item)
  clean_files_list(nested_files, package_name)

  sort!(nested_files)
  if is_testing && !is_sorted
    shuffle!(nested_files)
  end

  nested_files
end

function get_nested_files(cur_item)
  nested_files = []

  if isfile(cur_item)
    if endswith(cur_item, ".jl")
      push!(nested_files, cur_item)
    end

    return nested_files
  end

  for sub_item in readdir(cur_item)
    if startswith(sub_item, ".") ; continue ; end

    sub_nested_files = get_nested_files("$cur_item/$sub_item")
    append!(nested_files, sub_nested_files)
  end

  return nested_files
end

function clean_files_list(nested_files, package_name=nothing)
  delete_list = []

  if package_name == nothing
    split_limit = endswith(pwd(), "/test") ? 3 : 2
    package_name = replace(rsplit(pwd(), "/"; limit=split_limit)[2], ".jl", "")
  end

  push!(delete_list, "src/$package_name.jl")
  push!(delete_list, "src/Julz.jl")

  push!(delete_list, "test/runtests.jl")

  push!(delete_list, "/defaults.jl")

  push!(delete_list, "/input.jl")

  push!(delete_list, "src/functions/utils/flatten.jl")
  push!(delete_list, "src/functions/utils/delete_method.jl")
  push!(delete_list, "src/functions/utils/get_all_files.jl")
  push!(delete_list, "src/functions/utils/get_all_symbols.jl")
  push!(delete_list, "src/functions/utils/include_all_files.jl")
  push!(delete_list, "src/macros/export_all_files.jl")

  for deleted_file in delete_list
    for (index, nested_file) in enumerate(nested_files)
      if endswith(nested_file, deleted_file)
        deleteat!(nested_files, index)
        break
      end
    end
  end
end


text:
export iau2000a, nu2000k, fund_args

"""
    iau2000a(jd)

Compute the forced nutation of the non-rigid earth based on the IAU 2000A nutation model.

# Arguments
- `jd_high::Real`: High order part of the TT Julian date
- `jd_low::Real=0.0`: Low order part of the TT Julian date

# Return
`(dpsi,deps)` where
- `dpsi`: Nutation (luni-solar + planetary) in longitude, in radians
- `deps`: Nutation (luni-solar + planetary) in obliquity, in radians
"""
function iau2000a(jd_high::Real, jd_low::Real=0.0)
    # Interval between fundamental epoch J2000.0 and given date
    t = ((jd_high - T0) + jd_low) / 36525.0

    # Compute fundamental arguments in radians
    a = fund_args(t)

    Δϕ_ls = 0.0
    Δε_ls = 0.0

    @inbounds @simd for i in 1:size(nals_a)[1]
        arg = nals_a[i, 1] * a[1] +
              nals_a[i, 2] * a[2] +
              nals_a[i, 3] * a[3] +
              nals_a[i, 4] * a[4] +
              nals_a[i, 5] * a[5]
        sarg, carg = sincos(arg)
        Δϕ_ls += (cls_a[i, 1] + cls_a[i, 2] * t) * sarg + cls_a[i, 3] * carg
        Δε_ls += (cls_a[i, 4] + cls_a[i, 5] * t) * carg + cls_a[i, 6] * sarg
    end

    # Mean anomaly of the Moon.
    al = mod2pi(2.35555598 + 8328.6914269554 * t)
    # Mean anomaly of the Sun.
    alsu = mod2pi(6.24006013 + 628.301955 * t)
    #  Mean argument of the latitude of the Moon.
    af = mod2pi(1.627905234 + 8433.466158131 * t)
    #  Mean elongation of the Moon from the Sun.
    ad = mod2pi(5.198466741 + 7771.3771468121 * t)
    # Mean longitude of the ascending node of the Moon.
    aom = mod2pi(2.18243920 - 33.757045 * t)
    # General accumulated precession in longitude.
    apa = (0.02438175 + 0.00000538691 * t) * t
    # Planetary longitudes, Mercury through Neptune (Souchay et al. 1999).
    alme = mod2pi(4.402608842 + 2608.7903141574 * t)
    alve = mod2pi(3.176146697 + 1021.3285546211 * t)
    alea = mod2pi(1.753470314 + 628.3075849991 * t)
    alma = mod2pi(6.203480913 + 334.0612426700 * t)
    alju = mod2pi(0.599546497 + 52.9690962641 * t)
    alsa = mod2pi(0.874016757 + 21.3299104960 * t)
    alur = mod2pi(5.481293871 + 7.4781598567 * t)
    alne = mod2pi(5.321159000 + 3.8127774000 * t)

    Δϕ_pl = 0.0
    Δε_pl = 0.0

    @inbounds @simd for i in 1:size(napl_a)[1]
        arg = napl_a[i, 1] * al +
              napl_a[i, 2] * alsu +
              napl_a[i, 3] * af +
              napl_a[i, 4] * ad +
              napl_a[i, 5] * aom +
              napl_a[i, 6] * alme +
              napl_a[i, 7] * alve +
              napl_a[i, 8] * alea +
              napl_a[i, 9] * alma +
              napl_a[i, 10] * alju +
              napl_a[i, 11] * alsa +
              napl_a[i, 12] * alur +
              napl_a[i, 13] * alne +
              napl_a[i, 14] * apa
        sarg, carg = sincos(arg)
        Δϕ_pl += cpl_a[i, 1] * sarg + cpl_a[i, 2] * carg
        Δε_pl += cpl_a[i, 3] * sarg + cpl_a[i, 4] * carg
    end

    # Scale Results
    scaling_factor = 1.0e-7 * ASEC2RAD

    # Return result!
    return (Δϕ_ls + Δϕ_pl, Δε_ls + Δε_pl) .* scaling_factor
end

"""
    nu2000k(jd)

Compute the forced nutation of the non-rigid earth based on the NU2000K nutation model.

# Arguments
- `jd_high::Real`: High order part of the TT Julian date
- `jd_low::Real=0.0`: Low order part of the TT Julian date

# Return
`(dpsi,deps)` where
- `dpsi`: Nutation (luni-solar + planetary) in longitude, in radians
- `deps`: Nutation (luni-solar + planetary) in obliquity, in radians
"""
function nu2000k(jd_high::Real, jd_low::Real=0.0)
    # Interval between fundamental epoch J2000.0 and given date
    t = ((jd_high - T0) + jd_low) / 36525.0

    # Compute fundamental arguments in radians
    a = fund_args(t)

    Δϕ_ls = 0.0
    Δε_ls = 0.0

    @inbounds @simd for i in 1:size(nals_k)[1]
        arg = nals_k[i, 1] * a[1] +
              nals_k[i, 2] * a[2] +
              nals_k[i, 3] * a[3] +
              nals_k[i, 4] * a[4] +
              nals_k[i, 5] * a[5]
        sarg, carg = sincos(arg)
        Δϕ_ls += (cls_k[i, 1] + cls_k[i, 2] * t) * sarg + cls_k[i, 3] * carg
        Δε_ls += (cls_k[i, 4] + cls_k[i, 5] * t) * carg + cls_k[i, 6] * sarg
    end

    # Planetary longitudes, Mercury through Neptune, wrt mean dynamical
    # ecliptic and equinox of J2000, with high order terms omitted
    # (Simon et al. 1994, 5.8.1-5.8.8).
    alme = mod2pi(4.402608842461 + 2608.790314157421 * t)
    alve = mod2pi(3.176146696956 + 1021.328554621099 * t)
    alea = mod2pi(1.753470459496 + 628.307584999142 * t)
    alma = mod2pi(6.203476112911 + 334.061242669982 * t)
    alju = mod2pi(0.599547105074 + 52.969096264064 * t)
    alsa = mod2pi(0.874016284019 + 21.329910496032 * t)
    alur = mod2pi(5.481293871537 + 7.478159856729 * t)
    alne = mod2pi(5.311886286677 + 3.813303563778 * t)
    # General precession in longitude (Simon et al. 1994), equivalent to 5028.8200 arcsec/cy at J2000.
    apa = (0.024380407358 + 0.000005391235 * t) * t

    Δϕ_pl = 0.0
    Δε_pl = 0.0

    @inbounds @simd for i in 1:size(napl_k)[1]
        arg = napl_k[i, 1] * a[1] +
              napl_k[i, 2] * a[2] +
              napl_k[i, 3] * a[3] +
              napl_k[i, 4] * a[4] +
              napl_k[i, 5] * a[5] +
              napl_k[i, 6] * alme +
              napl_k[i, 7] * alve +
              napl_k[i, 8] * alea +
              napl_k[i, 9] * alma +
              napl_k[i, 10] * alju +
              napl_k[i, 11] * alsa +
              napl_k[i, 12] * alur +
              napl_k[i, 13] * alne +
              napl_k[i, 14] * apa
        sarg, carg = sincos(arg)
        Δϕ_pl += cpl_k[i, 1] * sarg + cpl_k[i, 2] * carg
        Δε_pl += cpl_k[i, 3] * sarg + cpl_k[i, 4] * carg
    end

    # Scale Results
    scaling_factor = 1.0e-7 * ASEC2RAD

    # Return result!
    return (Δϕ_ls + Δϕ_pl, Δε_ls + Δε_pl) .* scaling_factor
end


text:
using Chess, Chess.PGN

filename = "newspeed.pgn"
infile = open(filename)
@time gs = gamesfromstream(infile)
outf = open("newspeed.fen", "w")
for g in gs
        s = headervalue(g, "Result")
        r = 0
        try
                r = parse(Int, s[1])
        catch e
                println(position(infile))
        end
        if Bool(r) && s[2] == '/'
                r = 0.5
        end
        for b in boards(g)
                println(outf, fen(b), ',', r)
        end
end

close(infile)
close(outf)


text:
"""
Compute the inverse of the lattice. Takes special care of 1D or 2D cases.
"""
function compute_inverse_lattice(lattice::AbstractMatrix{T}) where {T}
    # Note: pinv pretty much does the same, but the implied SVD causes trouble
    #       with interval arithmetic and dual numbers, so we go for this version.
    n_dim = count(!iszero, eachcol(lattice))
    @assert 1 ≤ n_dim ≤ 3
    if n_dim == 3
        inv(lattice)
    else
        Mat3{T}([
            inv(lattice[1:n_dim, 1:n_dim])   zeros(T, n_dim, 3 - n_dim);
            zeros(T, 3 - n_dim, 3)
        ])
    end
end

"""
Compute the reciprocal lattice.
We use the convention that the reciprocal lattice is the set of G vectors such
that G ⋅ R ∈ 2π ℤ for all R in the lattice.
"""
function compute_recip_lattice(lattice::AbstractMatrix{T}) where {T}
    2T(π) * compute_inverse_lattice(lattice')
end

"""
Compute unit cell volume volume. In case of 1D or 2D case, the volume is the length/surface.
"""
function compute_unit_cell_volume(lattice)
    n_dim = count(!iszero, eachcol(lattice))
    abs(det(lattice[1:n_dim, 1:n_dim]))
end

"""Compute the diameter of the unit cell"""
function diameter(lattice::AbstractMatrix)
    # brute force search
    diam = zero(eltype(lattice))
    for vec in Vec3.(Iterators.product(-1:1, -1:1, -1:1))
        diam = max(diam, norm(lattice * vec))
    end
    diam
end


text:
module MatrixLM

using Distributed
using SharedArrays
using Statistics
using LinearAlgebra
using LinearAlgebra.BLAS
import LinearAlgebra.I, LinearAlgebra.mul!,
    LinearAlgebra.diag, LinearAlgebra.diagm
using DataFrames
using Random

export Response, Predictors, RawData, get_X, get_Z, get_Y, contr,
    kron_diag, add_intercept, remove_intercept, shuffle_rows, shuffle_cols,
    Mlm, mlm, t_stat, coef, predict, fitted, resid,
    perm_pvals, mlm_perms

# Data object types
include("data_types.jl")
# Contrasts
include("contr.jl")

# Predicted values
include("calc_preds.jl")
# Residuals
include("calc_resid.jl")

# Diagonal of kronecker product
include("kron_diag.jl")
# Variance shrinkage
include("shrink_sigma.jl")

# Miscellenous helper functions
include("misc_helpers.jl")
# Matrix linear model helper functions
include("mlm_helpers.jl")

# Matrix linear models
include("mlm.jl")
# Predictions and residuals
include("predict.jl")

# Permutations
include("perm_pvals.jl")
include("mlm_perms.jl")

end


text:
# ---
# title: 1125. Smallest Sufficient Team
# id: problem1125
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Bit Manipulation
# link: <https://leetcode.com/problems/smallest-sufficient-team/description/>
# hidden: true
# ---
#
# In a project, you have a list of required skills `req_skills`, and a list of
# `people`.  The i-th person `people[i]` contains a list of skills that person
# has.
#
# Consider a _sufficient team_ : a set of people such that for every required
# skill in `req_skills`, there is at least one person in the team who has that
# skill.  We can represent these teams by the index of each person: for example,
# `team = [0, 1, 3]` represents the people with skills `people[0]`, `people[1]`,
# and `people[3]`.
#
# Return **any**  sufficient team of the smallest possible size, represented by
# the index of each person.
#
# You may return the answer in any order.  It is guaranteed an answer exists.
#
#
#
# **Example 1:**
#
#
#
#     Input: req_skills = ["java","nodejs","reactjs"], people = [["java"],["nodejs"],["nodejs","reactjs"]]
#     Output: [0,2]
#
#
# **Example 2:**
#
#
#
#     Input: req_skills = ["algorithms","math","java","reactjs","csharp","aws"], people = [["algorithms","math","java"],["algorithms","math","reactjs"],["java","csharp","aws"],["reactjs","csharp"],["csharp","math"],["aws","java"]]
#     Output: [1,2]
#
#
#
#
# **Constraints:**
#
#   * `1 <= req_skills.length <= 16`
#   * `1 <= people.length <= 60`
#   * `1 <= people[i].length, req_skills[i].length, people[i][j].length <= 16`
#   * Elements of `req_skills` and `people[i]` are (respectively) distinct.
#   * `req_skills[i][j], people[i][j][k]` are lowercase English letters.
#   * Every skill in `people[i]` is a skill in `req_skills`.
#   * It is guaranteed a sufficient team exists.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# requires 11.6_produce.jl to have run first
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using PyPlot
using Statistics

function compute_phases(uu,vv,xx,yy) # also computes phase integral
    umean = mean(uu)
    vmean = mean(vv)
    ustd = std(uu)
    vstd = std(vv)
    N = size(uu, 1)

    phimat = atan.( (vv .- vmean) ./ vstd , (uu .- umean) ./ ustd)

    # compute path integrals
    kk = [0 1 1 0 0]
    ll = [0 0 1 1 0]

    xps = Float64[]    # coordinates of phase singularities (to be detected)
    yps = Float64[]
    for i=1:N-1
        for j = 1:N-1
            sum = 0
            for m = 1:4
                diffalpha = phimat[i+kk[m+1],j+ll[m+1]] - phimat[i+kk[m],j+ll[m]]
                diffalpha = mod( diffalpha + pi, 2 .* pi) - pi
                sum = sum + diffalpha
            end

            topcharge = sum / (2. * pi)

            if abs(topcharge) > 0.95   # phase singularity detected
                xpos = 0.5*(xx[i] + xx[i+1])
                ypos = 0.5*(yy[j] + yy[j+1])
                push!(xps, xpos)
                push!(yps, ypos)
                # xps = cat(xps, xpos, dims=1)
                # yps = cat(yps, ypos, dims=1)
                # println("xps = ",xps)
            end
        end
    end
    return phimat, xps, yps
end

data = wload(datadir("FitzHugh", "spiralwave.jld2"))
@unpack uout, tvec, params = data
@unpack a, b, d, ε, hsq6, N, L = params
h = L/N # spatial steps size
xx = h * (0:N-1)   # x coordinates of grid
yy = h * (0:N-1)   # y coordinates of grid

nt = length(tvec)
ures = uout[nt]

uu = ures[1,:,:]
vv = ures[2,:,:]

# compute udot and vdot (for contour lines).
# Requires FitzHugh_Nagumo_ODE to be in scope, from 11.6_produce.jl
du = copy(ures)
FitzHugh_Nagumo_ODE(du, ures, params, 0)
udot = du[1,:,:]
vdot = du[2,:,:]

phimat, xps, yps = compute_phases(uu,vv,xx,yy)

# %% plot
fig = plt.figure(figsize=(figx,figy))

ucmap = "gist_earth"
wcmap = "inferno"
pcmap = "hsv"
color_udot = COLORS[1]  # "w"   # 4 o.k.
color_vdot = COLORS[2]   #  6 o.k.

# plot u
axu = subplot(1,4,1)
axu.clear()

im1 = pcolormesh(xx, yy, uu; cmap = ucmap) #"inferno")
axu.contour(xx,yy,udot,levels=[0.],colors="0.5",linewidths=[3.], linestyles = "-")   # with cyan
axu.contour(xx,yy,udot,levels=[0.],colors=color_udot,linewidths=[2.], linestyles = "-")   # with cyan
im1.set_clim(-0.2, 1.)

axu.set_xticks([])
axu.set_yticks([])
axu.set_aspect("equal")

cbar = colorbar(im1; ax = axu, ticks=[-0.2, 1.], orientation="horizontal", pad=0.02)
cbar.set_label(L"u", labelpad = -30)
cbar.ax.set_xticklabels(["-0.2", "1"])  # horizontal colorbar

# plot v
axv = subplot(1,4,2)
im2 = axv.pcolormesh(xx,yy,vv,cmap = "inferno")
axv.contour(xx,yy,vdot,levels=[0.],colors=color_vdot,linewidths=[1.5])
im2.set_clim(0., 0.3)

axv.set_xticks([])
axv.set_yticks([])
axv.set_aspect("equal")

cbar = colorbar(im2; ax = axv, ticks=[0, 0.3], orientation="horizontal", pad=0.02)
cbar.set_label(L"w", labelpad = -30)
cbar.ax.set_xticklabels(["0", "0.3"])  # horizontal colorbar

# phases
axp = subplot(1,4,3)
axz = subplot(1,4,4)

im3 = nothing

for (j, ax) in enumerate((axp, axz))
    ax.clear()
    global im3 = ax.pcolormesh(xx,yy,phimat,cmap = pcmap) #  cyclic color map
    ax.contour(xx,yy,udot,levels=[0.],colors="0.5",linewidths=[3.], linestyles = "-")   # with cyan
    ax.contour(xx,yy,udot,levels=[0.],colors=color_udot,linewidths=[2.], linestyles = "-")   # with cyan
    ax.contour(xx,yy,vdot,levels=[0.],colors=color_vdot,linewidths=[3])
    im3.set_clim(-pi, pi)

    # add markers at location(s) of PS
    # --------------------------------
    for ips = 1:length(yps)
        ax.scatter(yps[ips],xps[ips]; marker = "o" , s=100,
        c="0.5", linewidths=2,  zorder = 2, edgecolors = "w")
    end

    ax.set_xticks([])
    ax.set_yticks([])
    ax.set_aspect("equal")

    cbar = colorbar(im3; ax=ax, ticks=[-pi, pi], orientation="horizontal", pad=0.02)
    cbar.ax.set_xticklabels([L"-\pi", L"\pi"])  # horizontal colorbar
    cbar.set_label(L"\alpha", labelpad=-30)

    if j == 1
        # add colorbar
    end
end

dxi = 40
xmin = xx[115] # 220 #
xmax = xx[115+dxi] # 300 # N = 600
ymin = yy[157] # 320 # N = 600
ymax = yy[157+dxi] # 400 # N = 600

zbox = ((xmin, xmax), (ymin, ymax))

axis_zoomin!(axz, axp, zbox, zbox, "k"; lw = 3)

add_identifiers!(fig, (axu, axv, axp, axz))

fig.subplots_adjust(wspace = 0.15, hspace = 0.51, left = 0.025,
right = 0.975, bottom = 0, top = 1.5)

wsave(plotsdir("11", "FitzHugh_Nagumo_PS"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
## 3D Panels
#==========================================================================================#

abstract type AbstractPanel3D <: AbstractPanel end

"""
    Panel3D(p1, p2, p3, p4)

Four Cartesian coordinates `p1, p2, p3, p4` representing corners of a panel in 3 dimensions. The following commutative diagram (math joke) depicts the order:

```
z → y
↓
x
        p1 —→— p4
        |       |
        ↓       ↓
        |       |
        p2 —→— p3
```
"""
struct Panel3D{T <: Real} <: AbstractPanel3D
    p1 :: SVector{3,T}
    p2 :: SVector{3,T}
    p3 :: SVector{3,T}
    p4 :: SVector{3,T}
end

Panel3D(p1, p2, p3, p4) = let T = promote_type(eltype(p1), eltype(p2), eltype(p3), eltype(p4)); Panel3D{T}(p1, p2, p3, p4) end

Panel3D((p1, p2, p3, p4)) = Panel3D(p1, p2, p3, p4)

Base.length(:: Panel3D) = 1

average_chord(panel :: Panel3D) = (p2(panel) - p1(panel) + p3(panel) - p4(panel)) / 2
average_width(panel :: Panel3D) = (p4(panel) - p1(panel) + p3(panel) - p2(panel)) / 2

"""
    panel_coordinates(panel :: Panel3D)

Compute the coordinates of a `Panel3D`.
"""
panel_coordinates(panel :: Panel3D) = [ p1(panel), p2(panel), p3(panel), p4(panel) ]

"""
    make_panels(xyzs)

Convert an array of coordinates corresponding to a wing, ordered from root to tip and leading-edge to trailing-edge, into panels.
"""
make_panels(xyzs :: AbstractArray{<: SVector{3,<: Real}}) = @views Panel3D.(xyzs[1:end-1,1:end-1], xyzs[2:end,1:end-1], xyzs[2:end,2:end], xyzs[1:end-1,2:end])

"""
    transform(panel :: Panel3D, rotation, translation)

Perform an affine transformation on the coordinates of a `Panel3D` given a rotation matrix and translation vector.
"""
function transform(panel :: Panel3D, rotation, translation)
    T = Translation(translation) ∘ LinearMap(rotation)
    Panel3D(T(panel.p1), T(panel.p2), T(panel.p3), T(panel.p4))
end

"""
    midpoint(panel :: Panel3D)

Compute the midpoint of a `Panel3D`.
"""
midpoint(panel :: Panel3D) = (p1(panel) + p2(panel) + p3(panel) + p4(panel)) / 4

"""
    panel_normal(panel :: Panel3D)

Compute the normal vector of a `Panel3D`.
"""
panel_normal(panel :: Panel3D) = let p31 = p3(panel) - p1(panel), p42 = p4(panel) - p2(panel); p31 × p42 end

"""
    transform_normal(panel :: Panel3D, h_l, g_l)

Transform the normal vector ``n̂₀`` of a `Panel3D` about the hinge axis ``ĥₗ`` by the control gain ``gₗ``.

The transformation is the following: ``n̂ₗ = gₗ ĥₗ × n̂₀`` (FVA, Drela, 6.36).
"""
transform_normal(panel :: Panel3D, h_l, g_l) = g_l * cross(h_l, panel_normal(panel))

"""
    panel_area(panel :: Panel3D)

Compute the (possibly non-planar, hence nonsensical) area of a `Panel3D`.
"""
panel_area(panel :: Panel3D) = 1/2 * norm(panel_normal(panel)) # (norm ∘ cross)(average_chord(panel), average_width(panel))

"""
    wetted_area(panels :: Array{Panel3D})

Compute the total wetted area by summing the areas of an array of `Panel3D`.
"""
wetted_area(panels) = sum(panel -> panel_area(panel), panels)

"""
    reflect_xz(panel :: Panel3D)

Reflect a Panel3D with respect to the ``x``-``z`` plane of its reference coordinate system.
"""
reflect_xz(panel :: Panel3D) = Panel3D((reflect_xz ∘ p1)(panel), (reflect_xz ∘ p2)(panel), (reflect_xz ∘ p3)(panel), (reflect_xz ∘ p4)(panel))

# Compute local axis coordinates
function local_coordinate_system(stream, normie)
    s_hat = normalize(stream)
    n_hat = normalize(normie)
    c_hat = s_hat × n_hat

    reduce(hcat, @SMatrix [ c_hat; s_hat; n_hat ])
end

# Compute local axis coordinates
local_coordinate_system(panel :: Panel3D) = local_coordinate_system((panel.p4 - panel.p1 + panel.p3 - panel.p2) / 2, panel_normal(panel))