text:
@testset "1.1.1.5 P(x) (a+b x)^m (c+d x)^n" begin
    (A, B, C, D, a, b, c, d, m, n, x, ) = @variables A B C D a b c d m n x

    #= ::Package:: =#

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

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

    #= ::Subsection:: =#
    #=n>0=#

    #= ::Subsection::Closed:: =#
    #=n<0=#

    @test_int [(a + b*x)^3*(A + B*x + C*x^2 + D*x^3)/sqrt(c + d*x), x, 2, -((2*(b*c - a*d)^3*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D)*sqrt(c + d*x))/d^7) - (2*(b*c - a*d)^2*(a*d*(2*c*C*d - B*d^2 - 3*c^2*D) - b*(5*c^2*C*d - 4*B*c*d^2 + 3*A*d^3 - 6*c^3*D))*(c + d*x)^(3/2))/(3*d^7) - (2*(b*c - a*d)*(a^2*d^2*(C*d - 3*c*D) - a*b*d*(8*c*C*d - 3*B*d^2 - 15*c^2*D) + b^2*(10*c^2*C*d - 6*B*c*d^2 + 3*A*d^3 - 15*c^3*D))*(c + d*x)^(5/2))/(5*d^7) + (2*(a^3*d^3*D + 3*a^2*b*d^2*(C*d - 4*c*D) - 3*a*b^2*d*(4*c*C*d - B*d^2 - 10*c^2*D) + b^3*(10*c^2*C*d - 4*B*c*d^2 + A*d^3 - 20*c^3*D))*(c + d*x)^(7/2))/(7*d^7) + (2*b*(3*a^2*d^2*D + 3*a*b*d*(C*d - 5*c*D) - b^2*(5*c*C*d - B*d^2 - 15*c^2*D))*(c + d*x)^(9/2))/(9*d^7) + (2*b^2*(b*C*d - 6*b*c*D + 3*a*d*D)*(c + d*x)^(11/2))/(11*d^7) + (2*b^3*D*(c + d*x)^(13/2))/(13*d^7)]
    @test_int [(a + b*x)^2*(A + B*x + C*x^2 + D*x^3)/sqrt(c + d*x), x, 2, (2*(b*c - a*d)^2*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D)*sqrt(c + d*x))/d^6 + (2*(b*c - a*d)*(a*d*(2*c*C*d - B*d^2 - 3*c^2*D) - b*(4*c^2*C*d - 3*B*c*d^2 + 2*A*d^3 - 5*c^3*D))*(c + d*x)^(3/2))/(3*d^6) + (2*(a^2*d^2*(C*d - 3*c*D) - 2*a*b*d*(3*c*C*d - B*d^2 - 6*c^2*D) + b^2*(6*c^2*C*d - 3*B*c*d^2 + A*d^3 - 10*c^3*D))*(c + d*x)^(5/2))/(5*d^6) + (2*(a^2*d^2*D + 2*a*b*d*(C*d - 4*c*D) - b^2*(4*c*C*d - B*d^2 - 10*c^2*D))*(c + d*x)^(7/2))/(7*d^6) + (2*b*(b*C*d - 5*b*c*D + 2*a*d*D)*(c + d*x)^(9/2))/(9*d^6) + (2*b^2*D*(c + d*x)^(11/2))/(11*d^6)]
    @test_int [(a + b*x)^1*(A + B*x + C*x^2 + D*x^3)/sqrt(c + d*x), x, 2, -((2*(b*c - a*d)*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D)*sqrt(c + d*x))/d^5) - (2*(a*d*(2*c*C*d - B*d^2 - 3*c^2*D) - b*(3*c^2*C*d - 2*B*c*d^2 + A*d^3 - 4*c^3*D))*(c + d*x)^(3/2))/(3*d^5) + (2*(a*d*(C*d - 3*c*D) - b*(3*c*C*d - B*d^2 - 6*c^2*D))*(c + d*x)^(5/2))/(5*d^5) + (2*(b*C*d - 4*b*c*D + a*d*D)*(c + d*x)^(7/2))/(7*d^5) + (2*b*D*(c + d*x)^(9/2))/(9*d^5)]
    @test_int [(a + b*x)^0*(A + B*x + C*x^2 + D*x^3)/sqrt(c + d*x), x, 2, (2*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D)*sqrt(c + d*x))/d^4 - (2*(2*c*C*d - B*d^2 - 3*c^2*D)*(c + d*x)^(3/2))/(3*d^4) + (2*(C*d - 3*c*D)*(c + d*x)^(5/2))/(5*d^4) + (2*D*(c + d*x)^(7/2))/(7*d^4)]
    @test_int [(A + B*x + C*x^2 + D*x^3)/((a + b*x)^1*sqrt(c + d*x)), x, 4, (2*(a^2*d^2*D - a*b*d*(C*d - c*D) - b^2*(c*C*d - B*d^2 - c^2*D))*sqrt(c + d*x))/(b^3*d^3) + (2*(b*C*d - 2*b*c*D - a*d*D)*(c + d*x)^(3/2))/(3*b^2*d^3) + (2*D*(c + d*x)^(5/2))/(5*b*d^3) - (2*(A*b^3 - a*(b^2*B - a*b*C + a^2*D))*atanh((sqrt(b)*sqrt(c + d*x))/sqrt(b*c - a*d)))/(b^(7/2)*sqrt(b*c - a*d))]
    @test_int [(A + B*x + C*x^2 + D*x^3)/((a + b*x)^2*sqrt(c + d*x)), x, 5, (2*(b*C*d - b*c*D - 2*a*d*D)*sqrt(c + d*x))/(b^3*d^2) - ((A - (a*(b^2*B - a*b*C + a^2*D))/b^3)*sqrt(c + d*x))/((b*c - a*d)*(a + b*x)) + (2*D*(c + d*x)^(3/2))/(3*b^2*d^2) - ((b^3*(2*B*c - A*d) - a*b^2*(4*c*C + B*d) - 5*a^3*d*D + 3*a^2*b*(C*d + 2*c*D))*atanh((sqrt(b)*sqrt(c + d*x))/sqrt(b*c - a*d)))/(b^(7/2)*(b*c - a*d)^(3/2))]
    @test_int [(A + B*x + C*x^2 + D*x^3)/((a + b*x)^3*sqrt(c + d*x)), x, 5, (2*D*sqrt(c + d*x))/(b^3*d) - ((A*b^3 - a*(b^2*B - a*b*C + a^2*D))*sqrt(c + d*x))/(2*b^3*(b*c - a*d)*(a + b*x)^2) - ((b^3*(4*B*c - 3*A*d) - a*b^2*(8*c*C + B*d) - 9*a^3*d*D + a^2*b*(5*C*d + 12*c*D))*sqrt(c + d*x))/(4*b^3*(b*c - a*d)^2*(a + b*x)) - ((b^3*(8*c^2*C - 4*B*c*d + 3*A*d^2) - 15*a^3*d^2*D + 3*a^2*b*d*(C*d + 12*c*D) - a*b^2*(8*c*C*d - B*d^2 + 24*c^2*D))*atanh((sqrt(b)*sqrt(c + d*x))/sqrt(b*c - a*d)))/(4*b^(7/2)*(b*c - a*d)^(5/2))]
    @test_int [(A + B*x + C*x^2 + D*x^3)/((a + b*x)^4*sqrt(c + d*x)), x, 5, -(((A*b^3 - a*(b^2*B - a*b*C + a^2*D))*sqrt(c + d*x))/(3*b^3*(b*c - a*d)*(a + b*x)^3)) - ((b^3*(6*B*c - 5*A*d) - a*b^2*(12*c*C + B*d) - 13*a^3*d*D + a^2*b*(7*C*d + 18*c*D))*sqrt(c + d*x))/(12*b^3*(b*c - a*d)^2*(a + b*x)^2) - ((b^3*(8*c^2*C - 6*B*c*d + 5*A*d^2) - 11*a^3*d^2*D + a^2*b*d*(C*d + 30*c*D) - a*b^2*(4*c*C*d - B*d^2 + 24*c^2*D))*sqrt(c + d*x))/(8*b^3*(b*c - a*d)^3*(a + b*x)) + ((5*a^3*d^3*D + a^2*b*d^2*(C*d - 18*c*D) - a*b^2*d*(4*c*C*d - B*d^2 - 24*c^2*D) + b^3*(8*c^2*C*d - 6*B*c*d^2 + 5*A*d^3 - 16*c^3*D))*atanh((sqrt(b)*sqrt(c + d*x))/sqrt(b*c - a*d)))/(8*b^(7/2)*(b*c - a*d)^(7/2))]
    @test_int [(A + B*x + C*x^2 + D*x^3)/((a + b*x)^5*sqrt(c + d*x)), x, 6, -(((A*b^3 - a*(b^2*B - a*b*C + a^2*D))*sqrt(c + d*x))/(4*b^3*(b*c - a*d)*(a + b*x)^4)) - ((b^3*(8*B*c - 7*A*d) - a*b^2*(16*c*C + B*d) - 17*a^3*d*D + 3*a^2*b*(3*C*d + 8*c*D))*sqrt(c + d*x))/(24*b^3*(b*c - a*d)^2*(a + b*x)^3) - ((b^3*(48*c^2*C - 40*B*c*d + 35*A*d^2) - 59*a^3*d^2*D + 3*a^2*b*d*(C*d + 56*c*D) - a*b^2*(16*c*C*d - 5*B*d^2 + 144*c^2*D))*sqrt(c + d*x))/(96*b^3*(b*c - a*d)^3*(a + b*x)^2) + ((5*a^3*d^3*D + 3*a^2*b*d^2*(C*d - 8*c*D) - a*b^2*d*(16*c*C*d - 5*B*d^2 - 48*c^2*D) + b^3*(48*c^2*C*d - 40*B*c*d^2 + 35*A*d^3 - 64*c^3*D))*sqrt(c + d*x))/(64*b^3*(b*c - a*d)^4*(a + b*x)) - (d*(5*a^3*d^3*D + 3*a^2*b*d^2*(C*d - 8*c*D) - a*b^2*d*(16*c*C*d - 5*B*d^2 - 48*c^2*D) + b^3*(48*c^2*C*d - 40*B*c*d^2 + 35*A*d^3 - 64*c^3*D))*atanh((sqrt(b)*sqrt(c + d*x))/sqrt(b*c - a*d)))/(64*b^(7/2)*(b*c - a*d)^(9/2))]

    @test_int [(a + b*x)^3*(A + B*x + C*x^2 + D*x^3)/(c + d*x)^(3/2), x, 2, (2*(b*c - a*d)^3*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D))/(d^7*sqrt(c + d*x)) - (2*(b*c - a*d)^2*(a*d*(2*c*C*d - B*d^2 - 3*c^2*D) - b*(5*c^2*C*d - 4*B*c*d^2 + 3*A*d^3 - 6*c^3*D))*sqrt(c + d*x))/d^7 - (2*(b*c - a*d)*(a^2*d^2*(C*d - 3*c*D) - a*b*d*(8*c*C*d - 3*B*d^2 - 15*c^2*D) + b^2*(10*c^2*C*d - 6*B*c*d^2 + 3*A*d^3 - 15*c^3*D))*(c + d*x)^(3/2))/(3*d^7) + (2*(a^3*d^3*D + 3*a^2*b*d^2*(C*d - 4*c*D) - 3*a*b^2*d*(4*c*C*d - B*d^2 - 10*c^2*D) + b^3*(10*c^2*C*d - 4*B*c*d^2 + A*d^3 - 20*c^3*D))*(c + d*x)^(5/2))/(5*d^7) + (2*b*(3*a^2*d^2*D + 3*a*b*d*(C*d - 5*c*D) - b^2*(5*c*C*d - B*d^2 - 15*c^2*D))*(c + d*x)^(7/2))/(7*d^7) + (2*b^2*(b*C*d - 6*b*c*D + 3*a*d*D)*(c + d*x)^(9/2))/(9*d^7) + (2*b^3*D*(c + d*x)^(11/2))/(11*d^7)]
    @test_int [(a + b*x)^2*(A + B*x + C*x^2 + D*x^3)/(c + d*x)^(3/2), x, 2, -((2*(b*c - a*d)^2*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D))/(d^6*sqrt(c + d*x))) + (2*(b*c - a*d)*(a*d*(2*c*C*d - B*d^2 - 3*c^2*D) - b*(4*c^2*C*d - 3*B*c*d^2 + 2*A*d^3 - 5*c^3*D))*sqrt(c + d*x))/d^6 + (2*(a^2*d^2*(C*d - 3*c*D) - 2*a*b*d*(3*c*C*d - B*d^2 - 6*c^2*D) + b^2*(6*c^2*C*d - 3*B*c*d^2 + A*d^3 - 10*c^3*D))*(c + d*x)^(3/2))/(3*d^6) + (2*(a^2*d^2*D + 2*a*b*d*(C*d - 4*c*D) - b^2*(4*c*C*d - B*d^2 - 10*c^2*D))*(c + d*x)^(5/2))/(5*d^6) + (2*b*(b*C*d - 5*b*c*D + 2*a*d*D)*(c + d*x)^(7/2))/(7*d^6) + (2*b^2*D*(c + d*x)^(9/2))/(9*d^6)]
    @test_int [(a + b*x)^1*(A + B*x + C*x^2 + D*x^3)/(c + d*x)^(3/2), x, 2, (2*(b*c - a*d)*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D))/(d^5*sqrt(c + d*x)) - (2*(a*d*(2*c*C*d - B*d^2 - 3*c^2*D) - b*(3*c^2*C*d - 2*B*c*d^2 + A*d^3 - 4*c^3*D))*sqrt(c + d*x))/d^5 + (2*(a*d*(C*d - 3*c*D) - b*(3*c*C*d - B*d^2 - 6*c^2*D))*(c + d*x)^(3/2))/(3*d^5) + (2*(b*C*d - 4*b*c*D + a*d*D)*(c + d*x)^(5/2))/(5*d^5) + (2*b*D*(c + d*x)^(7/2))/(7*d^5)]
    @test_int [(a + b*x)^0*(A + B*x + C*x^2 + D*x^3)/(c + d*x)^(3/2), x, 2, -((2*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D))/(d^4*sqrt(c + d*x))) - (2*(2*c*C*d - B*d^2 - 3*c^2*D)*sqrt(c + d*x))/d^4 + (2*(C*d - 3*c*D)*(c + d*x)^(3/2))/(3*d^4) + (2*D*(c + d*x)^(5/2))/(5*d^4)]
    @test_int [(A + B*x + C*x^2 + D*x^3)/((a + b*x)^1*(c + d*x)^(3/2)), x, 6, (2*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D))/(d^3*(b*c - a*d)*sqrt(c + d*x)) - (2*c*D*sqrt(c + d*x))/(b*d^3) + (2*(b*C*d - b*c*D - a*d*D)*sqrt(c + d*x))/(b^2*d^3) + (2*D*(c + d*x)^(3/2))/(3*b*d^3) - (2*(A*b^3 - a*(b^2*B - a*b*C + a^2*D))*atanh((sqrt(b)*sqrt(c + d*x))/sqrt(b*c - a*d)))/(b^(5/2)*(b*c - a*d)^(3/2))]
    @test_int [(A + B*x + C*x^2 + D*x^3)/((a + b*x)^2*(c + d*x)^(3/2)), x, 5, (a*b^2*B*d^3 - a^2*b*C*d^3 + a^3*d^3*D - b^3*(2*c^2*C*d - 2*B*c*d^2 + 3*A*d^3 - 2*c^3*D))/(b^3*d^2*(b*c - a*d)^2*sqrt(c + d*x)) - (A - (a*(b^2*B - a*b*C + a^2*D))/b^3)/((b*c - a*d)*(a + b*x)*sqrt(c + d*x)) + (2*D*sqrt(c + d*x))/(b^2*d^2) - ((b^3*(2*B*c - 3*A*d) - a*b^2*(4*c*C - B*d) - 3*a^3*d*D + a^2*b*(C*d + 6*c*D))*atanh((sqrt(b)*sqrt(c + d*x))/sqrt(b*c - a*d)))/(b^(5/2)*(b*c - a*d)^(5/2))]
    @test_int [(A + B*x + C*x^2 + D*x^3)/((a + b*x)^3*(c + d*x)^(3/2)), x, 5, -((a*b^2*B*d^3 - a^2*b*C*d^3 + a^3*d^3*D - b^3*(4*c^2*C*d - 4*B*c*d^2 + 5*A*d^3 - 4*c^3*D))/(2*b^3*d*(b*c - a*d)^3*sqrt(c + d*x))) - (A*b^3 - a*(b^2*B - a*b*C + a^2*D))/(2*b^3*(b*c - a*d)*(a + b*x)^2*sqrt(c + d*x)) - ((b^3*(4*B*c - 5*A*d) - a*b^2*(8*c*C - B*d) - 7*a^3*d*D + 3*a^2*b*(C*d + 4*c*D))*sqrt(c + d*x))/(4*b^2*(b*c - a*d)^3*(a + b*x)) - ((b^3*(8*c^2*C - 12*B*c*d + 15*A*d^2) - 3*a^3*d^2*D - a^2*b*d*(C*d - 12*c*D) + a*b^2*(8*c*C*d - 3*B*d^2 - 24*c^2*D))*atanh((sqrt(b)*sqrt(c + d*x))/sqrt(b*c - a*d)))/(4*b^(5/2)*(b*c - a*d)^(7/2))]
    @test_int [(A + B*x + C*x^2 + D*x^3)/((a + b*x)^4*(c + d*x)^(3/2)), x, 6, (a*b^2*B*d^3 - a^2*b*C*d^3 + a^3*d^3*D - b^3*(6*c^2*C*d - 6*B*c*d^2 + 7*A*d^3 - 6*c^3*D))/(3*b^3*(b*c - a*d)^4*sqrt(c + d*x)) - (A*b^3 - a*(b^2*B - a*b*C + a^2*D))/(3*b^3*(b*c - a*d)*(a + b*x)^3*sqrt(c + d*x)) - ((b^3*(6*B*c - 7*A*d) - a*b^2*(12*c*C - B*d) - 11*a^3*d*D + a^2*b*(5*C*d + 18*c*D))*sqrt(c + d*x))/(12*b^2*(b*c - a*d)^3*(a + b*x)^2) - ((b^3*(24*c^2*C - 42*B*c*d + 49*A*d^2) + 5*a^3*d^2*D - a^2*b*d*(11*C*d - 18*c*D) + a*b^2*(36*c*C*d - 7*B*d^2 - 72*c^2*D))*sqrt(c + d*x))/(24*b^2*(b*c - a*d)^4*(a + b*x)) - ((a^3*d^3*D + a^2*b*d^2*(C*d - 6*c*D) - a*b^2*d*(12*c*C*d - 5*B*d^2 - 24*c^2*D) - b^3*(24*c^2*C*d - 30*B*c*d^2 + 35*A*d^3 - 16*c^3*D))*atanh((sqrt(b)*sqrt(c + d*x))/sqrt(b*c - a*d)))/(8*b^(5/2)*(b*c - a*d)^(9/2))]

    @test_int [(a + b*x)^3*(A + B*x + C*x^2 + D*x^3)/(c + d*x)^(5/2), x, 2, (2*(b*c - a*d)^3*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D))/(3*d^7*(c + d*x)^(3/2)) + (2*(b*c - a*d)^2*(a*d*(2*c*C*d - B*d^2 - 3*c^2*D) - b*(5*c^2*C*d - 4*B*c*d^2 + 3*A*d^3 - 6*c^3*D)))/(d^7*sqrt(c + d*x)) - (2*(b*c - a*d)*(a^2*d^2*(C*d - 3*c*D) - a*b*d*(8*c*C*d - 3*B*d^2 - 15*c^2*D) + b^2*(10*c^2*C*d - 6*B*c*d^2 + 3*A*d^3 - 15*c^3*D))*sqrt(c + d*x))/d^7 + (2*(a^3*d^3*D + 3*a^2*b*d^2*(C*d - 4*c*D) - 3*a*b^2*d*(4*c*C*d - B*d^2 - 10*c^2*D) + b^3*(10*c^2*C*d - 4*B*c*d^2 + A*d^3 - 20*c^3*D))*(c + d*x)^(3/2))/(3*d^7) + (2*b*(3*a^2*d^2*D + 3*a*b*d*(C*d - 5*c*D) - b^2*(5*c*C*d - B*d^2 - 15*c^2*D))*(c + d*x)^(5/2))/(5*d^7) + (2*b^2*(b*C*d - 6*b*c*D + 3*a*d*D)*(c + d*x)^(7/2))/(7*d^7) + (2*b^3*D*(c + d*x)^(9/2))/(9*d^7)]
    @test_int [(a + b*x)^2*(A + B*x + C*x^2 + D*x^3)/(c + d*x)^(5/2), x, 2, -((2*(b*c - a*d)^2*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D))/(3*d^6*(c + d*x)^(3/2))) - (2*(b*c - a*d)*(a*d*(2*c*C*d - B*d^2 - 3*c^2*D) - b*(4*c^2*C*d - 3*B*c*d^2 + 2*A*d^3 - 5*c^3*D)))/(d^6*sqrt(c + d*x)) + (2*(a^2*d^2*(C*d - 3*c*D) - 2*a*b*d*(3*c*C*d - B*d^2 - 6*c^2*D) + b^2*(6*c^2*C*d - 3*B*c*d^2 + A*d^3 - 10*c^3*D))*sqrt(c + d*x))/d^6 + (2*(a^2*d^2*D + 2*a*b*d*(C*d - 4*c*D) - b^2*(4*c*C*d - B*d^2 - 10*c^2*D))*(c + d*x)^(3/2))/(3*d^6) + (2*b*(b*C*d - 5*b*c*D + 2*a*d*D)*(c + d*x)^(5/2))/(5*d^6) + (2*b^2*D*(c + d*x)^(7/2))/(7*d^6)]
    @test_int [(a + b*x)^1*(A + B*x + C*x^2 + D*x^3)/(c + d*x)^(5/2), x, 2, (2*(b*c - a*d)*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D))/(3*d^5*(c + d*x)^(3/2)) + (2*(a*d*(2*c*C*d - B*d^2 - 3*c^2*D) - b*(3*c^2*C*d - 2*B*c*d^2 + A*d^3 - 4*c^3*D)))/(d^5*sqrt(c + d*x)) + (2*(a*d*(C*d - 3*c*D) - b*(3*c*C*d - B*d^2 - 6*c^2*D))*sqrt(c + d*x))/d^5 + (2*(b*C*d - 4*b*c*D + a*d*D)*(c + d*x)^(3/2))/(3*d^5) + (2*b*D*(c + d*x)^(5/2))/(5*d^5)]
    @test_int [(a + b*x)^0*(A + B*x + C*x^2 + D*x^3)/(c + d*x)^(5/2), x, 2, -((2*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D))/(3*d^4*(c + d*x)^(3/2))) + (2*(2*c*C*d - B*d^2 - 3*c^2*D))/(d^4*sqrt(c + d*x)) + (2*(C*d - 3*c*D)*sqrt(c + d*x))/d^4 + (2*D*(c + d*x)^(3/2))/(3*d^4)]
    @test_int [(A + B*x + C*x^2 + D*x^3)/((a + b*x)^1*(c + d*x)^(5/2)), x, 4, (2*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D))/(3*d^3*(b*c - a*d)*(c + d*x)^(3/2)) + (2*(a*d*(2*c*C*d - B*d^2 - 3*c^2*D) - b*(c^2*C*d - A*d^3 - 2*c^3*D)))/(d^3*(b*c - a*d)^2*sqrt(c + d*x)) + (2*D*sqrt(c + d*x))/(b*d^3) - (2*(A*b^3 - a*(b^2*B - a*b*C + a^2*D))*atanh((sqrt(b)*sqrt(c + d*x))/sqrt(b*c - a*d)))/(b^(3/2)*(b*c - a*d)^(5/2))]
    @test_int [(A + B*x + C*x^2 + D*x^3)/((a + b*x)^2*(c + d*x)^(5/2)), x, 5, (3*a*b^2*B*d^3 - 3*a^2*b*C*d^3 + 3*a^3*d^3*D - b^3*(2*c^2*C*d - 2*B*c*d^2 + 5*A*d^3 - 2*c^3*D))/(3*b^3*d^2*(b*c - a*d)^2*(c + d*x)^(3/2)) - (A - (a*(b^2*B - a*b*C + a^2*D))/b^3)/((b*c - a*d)*(a + b*x)*(c + d*x)^(3/2)) - (a^2*b*C*d^3 - a^3*d^3*D + a*b^2*d*(4*c*C*d - 3*B*d^2 - 6*c^2*D) - b^3*(2*B*c*d^2 - 5*A*d^3 - 2*c^3*D))/(b^2*d^2*(b*c - a*d)^3*sqrt(c + d*x)) - ((b^3*(2*B*c - 5*A*d) - a*b^2*(4*c*C - 3*B*d) - a^3*d*D - a^2*b*(C*d - 6*c*D))*atanh((sqrt(b)*sqrt(c + d*x))/sqrt(b*c - a*d)))/(b^(3/2)*(b*c - a*d)^(7/2))]
    @test_int [(A + B*x + C*x^2 + D*x^3)/((a + b*x)^3*(c + d*x)^(5/2)), x, 6, -((3*a*b^2*B*d^3 - 3*a^2*b*C*d^3 + 3*a^3*d^3*D - b^3*(4*c^2*C*d - 4*B*c*d^2 + 7*A*d^3 - 4*c^3*D))/(6*b^3*d*(b*c - a*d)^3*(c + d*x)^(3/2))) - (A*b^3 - a*(b^2*B - a*b*C + a^2*D))/(2*b^3*(b*c - a*d)*(a + b*x)^2*(c + d*x)^(3/2)) + (a^2*b*C*d^2 + b^3*(2*c^2*C - 4*B*c*d + 7*A*d^2) - a^3*d^2*D + a*b^2*(4*c*C*d - 3*B*d^2 - 6*c^2*D))/(b^2*(b*c - a*d)^4*sqrt(c + d*x)) - ((b^3*(4*B*c - 7*A*d) - a*b^2*(8*c*C - 3*B*d) - 5*a^3*d*D + a^2*b*(C*d + 12*c*D))*sqrt(c + d*x))/(4*b*(b*c - a*d)^4*(a + b*x)) - ((b^3*(8*c^2*C - 20*B*c*d + 35*A*d^2) + a^3*d^2*D + 3*a^2*b*d*(C*d - 4*c*D) + 3*a*b^2*(8*c*C*d - 5*B*d^2 - 8*c^2*D))*atanh((sqrt(b)*sqrt(c + d*x))/sqrt(b*c - a*d)))/(4*b^(3/2)*(b*c - a*d)^(9/2))]

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

    @test_int [(a + b*x)^3*(c + d*x)^n*(A + B*x + C*x^2 + D*x^3), x, 2, -(((b*c - a*d)^3*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D)*(c + d*x)^(1 + n))/(d^7*(1 + n))) - ((b*c - a*d)^2*(a*d*(2*c*C*d - B*d^2 - 3*c^2*D) - b*(5*c^2*C*d - 4*B*c*d^2 + 3*A*d^3 - 6*c^3*D))*(c + d*x)^(2 + n))/(d^7*(2 + n)) - ((b*c - a*d)*(a^2*d^2*(C*d - 3*c*D) - a*b*d*(8*c*C*d - 3*B*d^2 - 15*c^2*D) + b^2*(10*c^2*C*d - 6*B*c*d^2 + 3*A*d^3 - 15*c^3*D))*(c + d*x)^(3 + n))/(d^7*(3 + n)) + ((a^3*d^3*D + 3*a^2*b*d^2*(C*d - 4*c*D) - 3*a*b^2*d*(4*c*C*d - B*d^2 - 10*c^2*D) + b^3*(10*c^2*C*d - 4*B*c*d^2 + A*d^3 - 20*c^3*D))*(c + d*x)^(4 + n))/(d^7*(4 + n)) + (b*(3*a^2*d^2*D + 3*a*b*d*(C*d - 5*c*D) - b^2*(5*c*C*d - B*d^2 - 15*c^2*D))*(c + d*x)^(5 + n))/(d^7*(5 + n)) + (b^2*(b*C*d - 6*b*c*D + 3*a*d*D)*(c + d*x)^(6 + n))/(d^7*(6 + n)) + (b^3*D*(c + d*x)^(7 + n))/(d^7*(7 + n))]
    @test_int [(a + b*x)^2*(c + d*x)^n*(A + B*x + C*x^2 + D*x^3), x, 2, ((b*c - a*d)^2*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D)*(c + d*x)^(1 + n))/(d^6*(1 + n)) + ((b*c - a*d)*(a*d*(2*c*C*d - B*d^2 - 3*c^2*D) - b*(4*c^2*C*d - 3*B*c*d^2 + 2*A*d^3 - 5*c^3*D))*(c + d*x)^(2 + n))/(d^6*(2 + n)) + ((a^2*d^2*(C*d - 3*c*D) - 2*a*b*d*(3*c*C*d - B*d^2 - 6*c^2*D) + b^2*(6*c^2*C*d - 3*B*c*d^2 + A*d^3 - 10*c^3*D))*(c + d*x)^(3 + n))/(d^6*(3 + n)) + ((a^2*d^2*D + 2*a*b*d*(C*d - 4*c*D) - b^2*(4*c*C*d - B*d^2 - 10*c^2*D))*(c + d*x)^(4 + n))/(d^6*(4 + n)) + (b*(b*C*d - 5*b*c*D + 2*a*d*D)*(c + d*x)^(5 + n))/(d^6*(5 + n)) + (b^2*D*(c + d*x)^(6 + n))/(d^6*(6 + n))]
    @test_int [(a + b*x)^1*(c + d*x)^n*(A + B*x + C*x^2 + D*x^3), x, 2, -(((b*c - a*d)*(c^2*C*d - B*c*d^2 + A*d^3 - c^3*D)*(c + d*x)^(1 + n))/(d^5*(1 + n))) - ((a*d*(2*c*C*d - B*d^2 - 3*c^2*D) - b*(3*c^2*C*d - 2*B*c*d^2 + A*d^3 - 4*c^3*D))*(c + d*x)^(2 + n))/(d^5*(2 + n)) + ((a*d*(C*d - 3*c*D) - b*(3*c*C*d - B*d^2 - 6*c^2*D))*(c + d*x)^(3 + n))/(d^5*(3 + n)) + ((b*C*d - 4*b*c*D + a*d*D)*(c + d*x)^(4 + n))/(d^5*(4 + n)) + (b*D*(c + d*x)^(5 + n))/(d^5*(5 + n))]
    @test_int [(a + b*x)^0*(c + d*x)^n*(A + B*x + C*x^2 + D*x^3), x, 2, ((c^2*C*d - B*c*d^2 + A*d^3 - c^3*D)*(c + d*x)^(1 + n))/(d^4*(1 + n)) - ((2*c*C*d - B*d^2 - 3*c^2*D)*(c + d*x)^(2 + n))/(d^4*(2 + n)) + ((C*d - 3*c*D)*(c + d*x)^(3 + n))/(d^4*(3 + n)) + (D*(c + d*x)^(4 + n))/(d^4*(4 + n))]
    @test_int [(c + d*x)^n*(A + B*x + C*x^2 + D*x^3)/(a + b*x)^1, x, 3, ((a^2*d^2*D - a*b*d*(C*d - c*D) - b^2*(c*C*d - B*d^2 - c^2*D))*(c + d*x)^(1 + n))/(b^3*d^3*(1 + n)) + ((b*C*d - 2*b*c*D - a*d*D)*(c + d*x)^(2 + n))/(b^2*d^3*(2 + n)) + (D*(c + d*x)^(3 + n))/(b*d^3*(3 + n)) - ((A*b^3 - a*(b^2*B - a*b*C + a^2*D))*(c + d*x)^(1 + n)*HypergeometricFunctions._₂F₁(1, 1 + n, 2 + n, (b*(c + d*x))/(b*c - a*d)))/(b^3*(b*c - a*d)*(1 + n))]
    @test_int [(c + d*x)^n*(A + B*x + C*x^2 + D*x^3)/(a + b*x)^2, x, 4, ((b*C*d - b*c*D - 2*a*d*D)*(c + d*x)^(1 + n))/(b^3*d^2*(1 + n)) - ((A - (a*(b^2*B - a*b*C + a^2*D))/b^3)*(c + d*x)^(1 + n))/((b*c - a*d)*(a + b*x)) + (D*(c + d*x)^(2 + n))/(b^2*d^2*(2 + n)) + ((a^3*d*D*(3 + n) - b^3*(B*c + A*d*n) + a*b^2*(2*c*C + B*d*(1 + n)) - a^2*b*(3*c*D + C*d*(2 + n)))*(c + d*x)^(1 + n)*HypergeometricFunctions._₂F₁(1, 1 + n, 2 + n, (b*(c + d*x))/(b*c - a*d)))/(b^3*(b*c - a*d)^2*(1 + n))]
    @test_int [(c + d*x)^n*(A + B*x + C*x^2 + D*x^3)/(a + b*x)^3, x, 4, (D*(c + d*x)^(1 + n))/(b^3*d*(1 + n)) - ((A*b^3 - a*(b^2*B - a*b*C + a^2*D))*(c + d*x)^(1 + n))/(2*b^3*(b*c - a*d)*(a + b*x)^2) - ((b^3*(2*B*c - A*d*(1 - n)) - a^3*d*D*(5 + n) - a*b^2*(4*c*C + B*d*(1 + n)) + a^2*b*(6*c*D + C*d*(3 + n)))*(c + d*x)^(1 + n))/(2*b^3*(b*c - a*d)^2*(a + b*x)) - ((b^3*(2*c^2*C + 2*B*c*d*n - A*d^2*(1 - n)*n) - a^3*d^2*D*(6 + 5*n + n^2) + a^2*b*d*(2 + n)*(6*c*D + C*d*(1 + n)) - a*b^2*(6*c^2*D + 4*c*C*d*(1 + n) + B*d^2*n*(1 + n)))*(c + d*x)^(1 + n)*HypergeometricFunctions._₂F₁(1, 1 + n, 2 + n, (b*(c + d*x))/(b*c - a*d)))/(2*b^3*(b*c - a*d)^3*(1 + n))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*P(x)*(a+b*x)^m*(c+d*x)^n*when*m*and*n*symbolic=#

    @test_int [(a + b*x)^m*(c + d*x)^n*(A + B*x), x, 3, (B*(a + b*x)^(1 + m)*(c + d*x)^(1 + n))/(b*d*(2 + m + n)) + ((A*b*d*(2 + m + n) - B*(b*c*(1 + m) + a*d*(1 + n)))*(a + b*x)^(1 + m)*(c + d*x)^n*HypergeometricFunctions._₂F₁(1 + m, -n, 2 + m, -((d*(a + b*x))/(b*c - a*d))))/(((b*(c + d*x))/(b*c - a*d))^n*(b^2*d*(1 + m)*(2 + m + n)))]
    @test_int [(a + b*x)^m*(c + d*x)^n*(A + B*x + C*x^2), x, 4, -(((a*C*d*(4 + m + 2*n) + b*(c*C*(2 + m) - B*d*(3 + m + n)))*(a + b*x)^(1 + m)*(c + d*x)^(1 + n))/(b^2*d^2*(2 + m + n)*(3 + m + n))) + (C*(a + b*x)^(2 + m)*(c + d*x)^(1 + n))/(b^2*d*(3 + m + n)) - (1/(b^3*d^2*(1 + m)*(2 + m + n)*(3 + m + n)))*(((d*(2 + m + n)*(a*b*c*C*(2 + m) + a^2*C*d*(1 + n) - A*b^2*d*(3 + m + n)) - (b*c*(1 + m) + a*d*(1 + n))*(a*C*d*(4 + m + 2*n) + b*(c*C*(2 + m) - B*d*(3 + m + n))))*(a + b*x)^(1 + m)*(c + d*x)^n*HypergeometricFunctions._₂F₁(1 + m, -n, 2 + m, -((d*(a + b*x))/(b*c - a*d))))/((b*(c + d*x))/(b*c - a*d))^n)]
    @test_int [(a + b*x)^m*(c + d*x)^n*(A + B*x + C*x^2 + D*x^3), x, 5, (1/(b^3*d^3*(2 + m + n)*(3 + m + n)*(4 + m + n)))*((a^2*d^2*D*(m^2 + m*(8 + 3*n) + 3*(6 + 5*n + n^2)) + b^2*(c^2*D*(6 + 5*m + m^2) - c*C*d*(2 + m)*(4 + m + n) + B*d^2*(12 + m^2 + 7*n + n^2 + m*(7 + 2*n))) + a*b*d*(c*D*(2 + m)*(6 + m + 3*n) - C*d*(m^2 + m*(8 + 3*n) + 2*(8 + 6*n + n^2))))*(a + b*x)^(1 + m)*(c + d*x)^(1 + n)) - ((a*d*D*(9 + 2*m + 3*n) + b*(c*D*(3 + m) - C*d*(4 + m + n)))*(a + b*x)^(2 + m)*(c + d*x)^(1 + n))/(b^3*d^2*(3 + m + n)*(4 + m + n)) + (D*(a + b*x)^(3 + m)*(c + d*x)^(1 + n))/(b^3*d*(4 + m + n)) + (1/(b^4*d^3*(1 + m)*(2 + m + n)*(3 + m + n)*(4 + m + n)))*(((d*(2 + m + n)*(a^3*d^2*D*(1 + n)*(6 + m + 2*n) + a*b^2*c*(2 + m)*(c*D*(3 + m) - C*d*(4 + m + n)) + A*b^3*d^2*(12 + m^2 + 7*n + n^2 + m*(7 + 2*n)) - a^2*b*d*(C*d*(1 + n)*(4 + m + n) - c*D*(2 + m)*(6 + m + 3*n))) - (b*c*(1 + m) + a*d*(1 + n))*(a^2*d^2*D*(m^2 + m*(8 + 3*n) + 3*(6 + 5*n + n^2)) + b^2*(c^2*D*(6 + 5*m + m^2) - c*C*d*(2 + m)*(4 + m + n) + B*d^2*(12 + m^2 + 7*n + n^2 + m*(7 + 2*n))) + a*b*d*(c*D*(2 + m)*(6 + m + 3*n) - C*d*(m^2 + m*(8 + 3*n) + 2*(8 + 6*n + n^2)))))*(a + b*x)^(1 + m)*(c + d*x)^n*HypergeometricFunctions._₂F₁(1 + m, -n, 2 + m, -((d*(a + b*x))/(b*c - a*d))))/((b*(c + d*x))/(b*c - a*d))^n)]
    #= [(a + b*x)^m*(c + d*x)^n*(A + B*x + C*x^2 + D*x^3 + F*x^4), x, 6, (1/(b^4*d^4*(2 + m + n)*(3 + m + n)*(4 + m + n)*(5 + m + n)))*((d*(3 + m + n)*(2*a*b^2*c*(3 + m)*(c*F*(4 + m) - d*D*(5 + m + n)) + b^3*B*d^2*(20 + m^2 + 9*n + n^2 + m*(9 + 2*n)) + 2*a^3*d^2*F*(m^2 + m*(11 + 4*n) + 4*(8 + 6*n + n^2)) + a^2*b*d*(4*c*F*(3 + m)*(6 + m + 2*n) - d*D*(m^2 + m*(11 + 4*n) + 3*(10 + 7*n + n^2)))) - b*c*(2 + m)*(a^2*d^2*F*(3*m^2 + m*(29 + 8*n) + 6*(12 + 7*n + n^2)) + b^2*(c*(3 + m)*(c*F*(4 + m) - d*D*(5 + m + n)) + C*d^2*(20 + m^2 + 9*n + n^2 + m*(9 + 2*n))) + a*b*d*(2*c*F*(3 + m)*(6 + m + 2*n) - d*D*(2*m^2 + m*(19 + 5*n) + 3*(15 + 8*n + n^2)))) - a*d*(2 + n)*(a^2*d^2*F*(3*m^2 + m*(29 + 8*n) + 6*(12 + 7*n + n^2)) + b^2*(c*(3 + m)*(c*F*(4 + m) - d*D*(5 + m + n)) + C*d^2*(20 + m^2 + 9*n + n^2 + m*(9 + 2*n))) + a*b*d*(2*c*F*(3 + m)*(6 + m + 2*n) - d*D*(2*m^2 + m*(19 + 5*n) + 3*(15 + 8*n + n^2)))))*(a + b*x)^(1 + m)*(c + d*x)^(1 + n)) + ((a^2*d^2*F*(3*m^2 + m*(29 + 8*n) + 6*(12 + 7*n + n^2)) + b^2*(c*(3 + m)*(c*F*(4 + m) - d*D*(5 + m + n)) + C*d^2*(20 + m^2 + 9*n + n^2 + m*(9 + 2*n))) + a*b*d*(2*c*F*(3 + m)*(6 + m + 2*n) - d*D*(2*m^2 + m*(19 + 5*n) + 3*(15 + 8*n + n^2))))*x*(a + b*x)^(1 + m)*(c + d*x)^(1 + n))/(b^3*d^3*(3 + m + n)*(4 + m + n)*(5 + m + n)) - ((a*d*F*(3*m + 4*(4 + n)) + b*(c*F*(4 + m) - d*D*(5 + m + n)))*(a + b*x)^(3 + m)*(c + d*x)^(1 + n))/(b^4*d^2*(4 + m + n)*(5 + m + n)) + (F*(a + b*x)^(4 + m)*(c + d*x)^(1 + n))/(b^4*d*(5 + m + n)) + (1/(b^5*d^4*(1 + m)*(2 + m + n)*(3 + m + n)*(4 + m + n)*(5 + m + n)))*(((d*(2 + m + n)*(d*(3 + m + n)*(A*b^4*d^2*(20 + m^2 + 9*n + n^2 + m*(9 + 2*n)) + a^2*(a^2*d^2*F*(1 + n)*(2*m + 3*(4 + n)) + b^2*c*(3 + m)*(c*F*(4 + m) - d*D*(5 + m + n)) - a*b*d*(d*D*(1 + n)*(5 + m + n) - 2*c*F*(3 + m)*(6 + m + 2*n)))) - a*b*c*(a^2*d^2*F*(3*m^2 + m*(29 + 8*n) + 6*(12 + 7*n + n^2)) + b^2*(c*(3 + m)*(c*F*(4 + m) - d*D*(5 + m + n)) + C*d^2*(20 + m^2 + 9*n + n^2 + m*(9 + 2*n))) + a*b*d*(2*c*F*(3 + m)*(6 + m + 2*n) - d*D*(2*m^2 + m*(19 + 5*n) + 3*(15 + 8*n + n^2))))) - (b*c*(1 + m) + a*d*(1 + n))*(d*(3 + m + n)*(2*a*b^2*c*(3 + m)*(c*F*(4 + m) - d*D*(5 + m + n)) + b^3*B*d^2*(20 + m^2 + 9*n + n^2 + m*(9 + 2*n)) + 2*a^3*d^2*F*(m^2 + m*(11 + 4*n) + 4*(8 + 6*n + n^2)) + a^2*b*d*(4*c*F*(3 + m)*(6 + m + 2*n) - d*D*(m^2 + m*(11 + 4*n) + 3*(10 + 7*n + n^2)))) - b*c*(2 + m)*(a^2*d^2*F*(3*m^2 + m*(29 + 8*n) + 6*(12 + 7*n + n^2)) + b^2*(c*(3 + m)*(c*F*(4 + m) - d*D*(5 + m + n)) + C*d^2*(20 + m^2 + 9*n + n^2 + m*(9 + 2*n))) + a*b*d*(2*c*F*(3 + m)*(6 + m + 2*n) - d*D*(2*m^2 + m*(19 + 5*n) + 3*(15 + 8*n + n^2)))) - a*d*(2 + n)*(a^2*d^2*F*(3*m^2 + m*(29 + 8*n) + 6*(12 + 7*n + n^2)) + b^2*(c*(3 + m)*(c*F*(4 + m) - d*D*(5 + m + n)) + C*d^2*(20 + m^2 + 9*n + n^2 + m*(9 + 2*n))) + a*b*d*(2*c*F*(3 + m)*(6 + m + 2*n) - d*D*(2*m^2 + m*(19 + 5*n) + 3*(15 + 8*n + n^2))))))*(a + b*x)^(1 + m)*(c + d*x)^n*HypergeometricFunctions._₂F₁(1 + m, -n, 2 + m, -((d*(a + b*x))/(b*c - a*d))))/((b*(c + d*x))/(b*c - a*d))^n)] =#
    end

max_stars_repo_path: test/1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.5 P(x) (a+b x)^m (c+d x)^n.jl


text:
"""
    FullFeaturedCPNN(;
        graphChain::Flux.Chain
        nodeChain::Flux.Chain
        globalChain::Flux.Chain
        outputChain::Flux.Dense
    )

This structure is here to provide a flexible way to create a nn model which respect this approach:
Making modification on the graph, then extract one node feature and modify it.
"""
Base.@kwdef struct FullFeaturedCPNN <: NNStructure
    graphChain::Flux.Chain = Flux.Chain()
    nodeChain::Flux.Chain = Flux.Chain()
    globalChain::Flux.Chain = Flux.Chain()
    outputChain::Union{Flux.Dense, Flux.Chain} = Flux.Chain()
end

# Enable the `|> gpu` syntax from Flux
Flux.@functor FullFeaturedCPNN

# TODO make it possible to use global features
function (nn::FullFeaturedCPNN)(states::BatchedDefaultTrajectoryState)
    @assert !isnothing(states.allValuesIdx)

    variableIdx = states.variableIdx
    batchSize = length(variableIdx)
    allValuesIdx = states.allValuesIdx
    actionSpaceSize = size(allValuesIdx, 1)

    # chain working on the graph(s) with the GNNs
    featuredGraph = nn.graphChain(states.fg)
    nodeFeatures = featuredGraph.nf
    globalFeatures = featuredGraph.gf

    # Extract the features corresponding to the varibales
    variableIndices = nothing
    Zygote.ignore() do
        variableIndices = Flux.unsqueeze(CartesianIndex.(variableIdx, 1:batchSize), 1)
    end
    variableFeatures = nodeFeatures[:, variableIndices]
    variableFeatures = reshape(nn.nodeChain(RL.flatten_batch(variableFeatures)), :, 1, batchSize)

    # Extract the features corresponding to the values
    valueIndices = nothing
    Zygote.ignore() do
        valueIndices = CartesianIndex.(allValuesIdx, repeat(transpose(1:batchSize); outer=(actionSpaceSize, 1)))
    end
    valueFeatures = nodeFeatures[:, valueIndices]
    valueFeatures = reshape(nn.nodeChain(RL.flatten_batch(valueFeatures)), :, actionSpaceSize, batchSize)

    finalFeatures = nothing
    if sizeof(globalFeatures) != 0

        # Extract the global features
        globalFeatures = reshape(nn.globalChain(globalFeatures), :, 1, batchSize)

        # Prepare the input of the outputChain
        finalFeatures = vcat(
            repeat(variableFeatures; outer=(1, actionSpaceSize, 1)),
            repeat(globalFeatures; outer=(1, actionSpaceSize, 1)),
            valueFeatures,
        )
        finalFeatures = RL.flatten_batch(finalFeatures)
    else
        # Prepare the input of the outputChain
        finalFeatures = vcat(
            repeat(variableFeatures; outer=(1, actionSpaceSize, 1)),
            valueFeatures,
        )
        finalFeatures = RL.flatten_batch(finalFeatures)
    end

    # output layer
    predictions = nn.outputChain(finalFeatures)
    output = reshape(predictions, actionSpaceSize, batchSize)

    return output
end

ext: jl


text:
module Elemental

using Compat
import Compat.String
using DistributedArrays

import Base: *, \, Ac_mul_B
import Base: convert, copy, copy!, countnz, dot, eltype, fill!, getindex, hcat, inv, length,
    logdet, pointer, print, setindex!, showarray, similar, size
import Base.LinAlg: A_mul_B!, Ac_mul_B!, axpy!, norm, scale!, svd, svdvals, svdvals!

const libEl = abspath(joinpath(dirname(@__FILE__), "..", "deps", "usr", "lib", "libEl"))

include("error.jl")

function Init()
    argc = Cint[0]
    argv = Vector{String}[String[""]]
    err = ccall((:ElInitialize, libEl), Cint,
                (Ptr{Cint}, Ptr{Ptr{Ptr{UInt8}}}), argc, pointer(argv))
    err == 0 || error("Error Initializing Elemental: $(ErrorString(err))")
    return nothing
end

function Initialized()
    active = Ref(zero(Cint))
    err = ccall((:ElInitialized, libEl), Cuint, (Ref{Cint},), active)
    err == 0 || throw(ElError(err))
    return active[] == 1
end

function Finalize()
    err = ccall((:ElFinalize, libEl), Cint, ())
    err == 0 || error("Error Finalizing Elemental: $(ErrorString(err))")
    return nothing
end

function __init__()
    Init()
    atexit() do
        Initialized() && Finalize()
    end
end

include("core/types.jl")
include("mpi.jl")
include("core/matrix.jl")
include("core/grid.jl")
include("core/sparsematrix.jl")
include("core/distmatrix.jl")
include("core/distsparsematrix.jl")
include("core/distmultivec.jl")
include("matrices.jl")
include("io.jl")
include("blas_like/level1.jl")
include("blas_like/level2.jl")
include("blas_like/level3.jl")
include("lapack_like/euclidean_min.jl")
include("lapack_like/factor.jl")
include("lapack_like/props.jl")
include("lapack_like/funcs.jl")
include("lapack_like/solve.jl")
include("lapack_like/spectral.jl")
include("optimization/solvers.jl")
include("optimization/models.jl")
include("julia/generic.jl")
include("julia/darray.jl")

end # module


text:
#===============================================================================

    amp_fit_model.jl

    function to fit amplification models
    this is a refactoring of the function
    previously called mod_bl_q() in amp.jl

    Author: Tom Price
    Date:   July 2019

===============================================================================#

## constant:
## a value that cannot be obtained by normal calculation of Ct
const AMP_CT_VAL_DOMAINERROR        = -99

#==============================================================================#

## amp_fit_model() definition for DFC models
"Model amplification data for a single well and channel."
function amp_fit_model(
    ::Type{Val{M}} where M <: DFCModel,
    ::Type{R},
    i                   ::AmpInput,
    fluos               ::Vector{Float_T},
    ## these parameters not used to fit DFC models
    baseline_cyc_bounds ::AbstractArray,
    cq_method           ::CqMethod,
    ct_fluo             ::Float_T,
) where R <: Union{AmpLongModelResults,AmpShortModelResults}
    bl() = coefs[1] .+
            i.amp_model in [MAK3, MAKERGAUL4] ?
                coefs[2] .* cycles :
                0.0
    d0() = coefs[findfirst(bl.bl_fit.coef_syms .== :d0)]

    debug(logger, "at amp_fit_model(::Type{Val{M}} where M <: DFCModel)")
    ## no fallback for baseline, because:
    ## (1) curve may fit well though :Error or :UserLimit
    ## (search step becomes very small but has not converge);
    ## (2) the guessed basedline (`start` of `fb`) is usually
    ## quite close to a sensible baseline.
    # const num_cycles = length(fluos)
    # const cycles = range(1.0, num_cycles)
    num_cycles = i.num_cycles
    cycles = Vector(i.cycles) ## allows non-contiguous sequences of cycles
    wts = ones(num_cycles)
    ## fit model
    bl_fit = fit(Val{M}, cycles, fluos, wts; solver = i.solver)
    coefs = bl_fit.coefs
    blsub_fluos = fluos .- bl()
    ## set output
    if isa(R, AmpShortModelResults)
        return AmpShortModelResults(
            fluos, ## rbbs_ary3,
            blsub_fluos,
            NaN_T, ## dr1_pred
            NaN_T, ## dr2_pred
            NaN_T, ## cq
            d0())
    else
        return AmpLongModelResults(
            fluos, ## rbbs_ary3,
            bl_fit,
            [string(amp_model)], ## bl_notes
            blsub_fluos,
            bl_fit, ## quant_fit
            bl_fit.status, ## quant_status
            coefs, ## quant_coefs
            d0(),
            pred_from_cycs(Val{M}, cycles, coefs...), ## quant_fluos
            NaN_T, ## dr1_pred,
            NaN_T, ## dr2_pred,
            Inf_T, ## max_dr1
            Inf_T, ## max_dr2
            OrderedDict(), ## cyc_vals_4cq
            OrderedDict(), ## eff_vals_4cq
            NaN_T, ## cq_raw
            NaN_T, ## cq
            NaN_T, ## eff
            NaN_T) ## cq_fluo
    end ## if Q
end ## fit_amplification_model() where DFCModel

#==============================================================================#

## fit_amplification_model() definition for SFC models
function amp_fit_model(
    ::Type{Val{SFCModel}},
    ::Type{R},
    i                   ::AmpInput,
    fluos               ::Vector{Float_T},
    baseline_cyc_bounds ::AbstractArray,
    cq_method           ::CqMethod,
    ct_fluo             ::Float_T,
) where R <: AmpModelResults
    "Calculates weights used to estimate the baseline by SFC model."
    @inline function SFC_wts()
        local wts
        debug(logger, "at SFC_wts")
        if i.bl_method == lin_1ft || i.bl_method == lin_2ft ## linear models
            wts = zeros(num_cycles)
            wts[bl_cycs()] .= 1
            return wts
        else
            ## some kind of sigmoid model is used to estimate amplification curve
            ## issue: why are `baseline_cyc_bounds` not baked into the weights as per above ???
            if num_cycles >= last_cyc_wt0
                return vcat(zeros(last_cyc_wt0), ones(num_cycles - last_cyc_wt0))
            else
                return zeros(num_cycles)
            end
        end
    end ## SFC_wts()

    "Check that the SFC model used to calculate the baseline terminated appropriately."
    @inline function good_status()
        if bl_status == :Optimal || bl_status == :UserLimit
            (min_bfd, max_bfd) = extrema(blsub_fluos) ## `bfd` - blsub_fluos_draft
            if !(max_bfd - min_bfd <= abs(min_bfd))
                push!(bl_notes, "model-derived baseline")
                return true
            else
                push!(bl_notes, "fallback")
                push!(bl_notes, "max_bfd ($max_bfd) - min_bfd ($min_bfd) " *
                    "== $(max_bfd - min_bfd) <= abs(min_bfd)")
                return false
            end
        end ## if bl_fit.status
        push!(bl_notes, "fallback")
        # (bl_status == :Error) && return false
        ## other status codes include
        ## :Infeasible, :Unbounded, :DualityFailure, and possibly others
        ## https://mathprogbasejl.readthedocs.io/en/latest/solverinterface.html
        ## My suggestion is to treat them all the same as :Error (TP Jan 2019):
        return false
        ## Alternatively, an error could be raised:
        # error(logger, ExceptionError("Baseline estimation returned " *
        #     "unrecognized termination status $(bl_status)"))
    end ## good_status()

    @inline function bl_cycs()
        len_bcb = length(baseline_cyc_bounds)
        debug(logger, "at bl_cycs")
        if !(len_bcb in [0, 2])
            throw(ArgumentError("length of `baseline_cyc_bounds` must be 0 or 2"))
        elseif len_bcb == 2
            push!(bl_notes, "user-defined")
            return colon(baseline_cyc_bounds...)
        elseif len_bcb == 0
            return auto_choose_bl_cycs()
        end
        ## fallthrough
        throw(ErrorException("too few cycles to estimate baseline"))
    end ## bl_cycs()

    ## uses `fluos`, `last_cyc_wt0`; updates `bl_notes` using push!()
    ## `last_cyc_wt0 == floor(i.min_reliable_cyc) - 1`
    "Automatically choose baseline cycles as the flat part of the curve."
    @inline function auto_choose_bl_cycs()
        debug(logger, "at auto_choose_bl_cycs")
        (min_fluo, min_fluo_cyc) = findmin(fluos)
        ## finite diff is used to estimate derivative because
        ## `Dierckx.Spline1D` resulted in all `NaN` in some cases
        dr2_cfd = finite_diff(cycles, fluos; nu = 2)

        dr2_cfd_left  = dr2_cfd[1:min_fluo_cyc]
        dr2_cfd_right = dr2_cfd[min_fluo_cyc:end]
        (max_dr2_left_cyc, max_dr2_right_cyc) =
            (dr2_cfd_left, dr2_cfd_right) |> mold(indmax)
        if max_dr2_right_cyc <= last_cyc_wt0
            ## fluo on fitted spline may not be close to raw fluo
            ## at `max_dr2_left_cyc` and `max_dr2_right_cyc`
            # push!(bl_notes, "max_dr2_right_cyc ($max_dr2_right_cyc) <= " *
            #     "last_cyc_wt0 ($last_cyc_wt0), bl_cycs = $(last_cyc_wt0+1):$num_cycles")
            return colon(last_cyc_wt0 + 1, num_cycles)
        end
        ## max_dr2_right_cyc > last_cyc_wt0
        bl_cyc_start = max(last_cyc_wt0 + 1, max_dr2_left_cyc)
        # push!(bl_notes, "max_dr2_right_cyc ($max_dr2_right_cyc) > " *
        #     "last_cyc_wt0 ($last_cyc_wt0), bl_cyc_start = $bl_cyc_start " *
        #     "(max(last_cyc_wt0+1, max_dr2_left_cyc), i.e. " *
        #     "max($(last_cyc_wt0+1), $max_dr2_left_cyc))")
        if max_dr2_right_cyc - bl_cyc_start <= 1
            # push!(bl_notes, "max_dr2_right_cyc ($max_dr2_right_cyc) - " *
            #     "bl_cyc_start ($bl_cyc_start) <= 1")
            if (max_dr2_right_cyc < num_cycles)
                (max_dr2_right_2, max_dr2_right_cyc_2_shifted) =
                    findmax(dr2_cfd[max_dr2_right_cyc + 1:end])
            else
                max_dr2_right_cyc_2_shifted = 0
            end
            max_dr2_right_cyc_2 = max_dr2_right_cyc_2_shifted + max_dr2_right_cyc
            if max_dr2_right_cyc_2 - max_dr2_right_cyc <= 1
                bl_cyc_end = num_cycles
                # push!(bl_notes, "max_dr2_right_cyc_2 ($max_dr2_right_cyc_2) - " *
                #     "max_dr2_right_cyc ($max_dr2_right_cyc) == 1")
            else
                ## max_dr2_right_cyc_2 - max_dr2_right_cyc > 1
                # push!(bl_notes, "max_dr2_right_cyc_2 ($max_dr2_right_cyc_2) - " *
                #     "max_dr2_right_cyc ($max_dr2_right_cyc) != 1")
                bl_cyc_end = max_dr2_right_cyc_2
            end ## if
        else
            ## max_dr2_right_cyc - bl_cyc_start > 1
            # push!(bl_notes, "max_dr2_right_cyc ($max_dr2_right_cyc) - " *
            #     "bl_cyc_start ($bl_cyc_start) > 1")
            bl_cyc_end = max_dr2_right_cyc
        end ## if
        # push!(bl_notes, "bl_cyc_end = $bl_cyc_end")
        # push!(bl_notes, "bl_cycles = $bl_cyc_start:$bl_cyc_end")
        return bl_cyc_start:bl_cyc_end
    end ## auto_choose_bl_cycs()

    denser_len(i , n :: Int_T) =
        i.denser_factor * (n - 1) + 1
    #
    interpolated_cycles() =
        range(1.0, 1.0/i.denser_factor, len_denser)

    ## functions used to calculate cq / cq_raw / cq_fluo >>

    dr1_pred_() =
        funcs_pred[:dr1](cycles_denser, quant_coefs...)
    #
    dr2_pred_() =
        funcs_pred[:dr2](cycles_denser, quant_coefs...)
    #
    cyc_max_dr1_(dr1_pred ::Vector{Float_T}) =
        cycles_denser[indmax(dr1_pred)]
    #
    cyc_max_dr2_(dr2_pred ::Vector{Float_T}) =
        cycles_denser[indmax(dr2_pred)]
    #
    cy0_(max_dr1 ::Float_T, idx_max_dr1 ::Int_T) =
        cy0_(max_dr1, idx_max_dr1, cycles_denser[idx_max_dr1])
    cy0_(max_dr1 ::Float_T, idx_max_dr1 ::Int_T, cyc_max_dr1 ::Float_T) =
        cyc_max_dr1 - funcs_pred[:f](cyc_max_dr1, quant_coefs...) / max_dr1
    #
    ct_() = try
        funcs_pred[:inv](ct_fluo, quant_coefs...)
    catch err
        isa(err, DomainError) ?
            AMP_CT_VAL_DOMAINERROR :
            rethrow()
    end ## ct_()
    #
    max_eff() = cycles_denser[indmax(map(func_pred_eff, cycles_denser))]
    #
    nonpos2NaN(x ::Float_T) =
        x <= zero(x) ? NaN_T : x
    #
    calc_cq_fluo(cq_raw ::Float_T) =
        funcs_pred[:f](nonpos2NaN(cq_raw), quant_coefs...)
    #
    @inline function calc_cq_raw()
        (cq_method == ct)      && return ct_()
        (cq_method == max_eff) && return max_eff_()
        calc_cq_raw(dr1_pred_())
    end
    @inline function calc_cq_raw(dr1_pred ::Vector{Float_T})
        (cq_method == cp_dr1)  && return cyc_max_dr1_(dr1_pred)
        (cq_method == Cy0)     && return cy0_(findmax(dr1_pred)...)
        calc_cq_raw(dr1_pred, dr2_pred_())
    end
    @inline function calc_cq_raw(dr1_pred ::Vector{Float_T}, dr2_pred ::Vector{Float_T})
        (cq_method == cp_dr2)  && return cyc_max_dr2_(dr2_pred)
        (cq_method == cp_dr1)  && return cyc_max_dr1_(dr1_pred)
        (cq_method == Cy0)     && return cy0_(findmax(dr1_pred)...)
        (cq_method == ct)      && return ct_()
        (cq_method == max_eff) && return max_eff_()
        ## fallthrough
        throw(ArgumentError("`cq_method` $cq_method not implemented"))
    end

    ## function needed because `Cy0` may not be in `cycles_denser`
    @inline function func_pred_eff(cyc)
        try
            -(map([0.5, -0.5]) do epsilon
                log2(funcs_pred[:f](cyc + epsilon, quant_coefs...))
            end...)
        catch err
            isa(err, DomainError) ?
                NaN_T :
                rethrow()
        end ## try
    end ## func_pred_eff()

    ## << end of function definitions nested within fit_baseline_model() ## SFC

    debug(logger, "at amp_fit_model(Val{SFCModel})")
    #
    ## fit baseline model
    # const num_cycles = length(fluos)
    # const cycles = range(1.0, num_cycles)
    ## the following allows non-contiguous sequences of cycles
    ## NB. Ipopt solver hangs unless cycles converted from SVector to Vector
    num_cycles = i.num_cycles
    cycles = Vector(i.cycles)
    ## to determine weights (`wts`) for sigmoid fitting per `i.min_reliable_cyc`
    last_cyc_wt0 = 0
    if i.bl_method in SFC_MODEL_NAMES
        ## fit model to find baseline
        wts = SFC_wts()
        len_denser = denser_len(i, num_cycles)
        cycles_denser = interpolated_cycles()
        bl_fit = i.SFC_model_def_func(i.bl_method).func_fit(
            cycles[i.min_reliable_cyc:end], fluos[i.min_reliable_cyc:end], wts; solver = i.solver)
        bl_status = bl_fit.status
        len_bcb = length(i.baseline_cyc_bounds)
        if len_bcb==0
            if num_cycles<=i.min_reliable_cyc
                auto_cycs=1:num_cycles
            else
                if (bl_status == :Optimal || bl_status == :UserLimit) #&& (abs(maximum(fluos)-minimum(fluos))>20000)
                    dr2_pred_tm=i.SFC_model_def_func(i.bl_method).funcs_pred[:dr2](
                        cycles_denser, bl_fit.coefs...)
                    min_fluo_cyc_=findmin(fluos)[2]
                    if (min_fluo_cyc_<i.min_reliable_cyc)
                        fluos_muted=vcat(fill(findmax(fluos)[1],i.min_reliable_cyc-1),fluos[i.min_reliable_cyc:end])
                        min_fluo_cyc_=findmin(fluos_muted)[2]
                    end
                    approx_min_dr2=findmin(abs.(cycles_denser-min_fluo_cyc_))[2]
                    basecyc_1st=i.min_reliable_cyc
                    if !(min_fluo_cyc_==i.min_reliable_cyc)
                        dr2_left2mf=dr2_pred_tm[1:approx_min_dr2-1]
                        if !(all(dr2_left2mf.<dr2_pred_tm[approx_min_dr2]))
                            basecyc_1st=Int(floor(cycles_denser[findmax(dr2_left2mf)[2]]))
                            basecyc_1st=i.min_reliable_cyc<=basecyc_1st?basecyc_1st:i.min_reliable_cyc
                        end
                    end

                    if (min_fluo_cyc_==i.num_cycles || length(cycles_denser)==approx_min_dr2)
                        basecyc_last=i.num_cycles
                    else
                        cycmax_dr2=findmax(dr2_pred_tm)[2]
                        if (cycmax_dr2 == approx_min_dr2)
                            basecyc_last=Int(floor(cycles_denser[cycmax_dr2]))
                        else
                            dr2_right2mf=dr2_pred_tm[approx_min_dr2+1:end]
                            basecyc_last=Int(floor(cycles_denser[findmax(dr2_right2mf)[2]+approx_min_dr2]))
                            basecyc_last=i.min_reliable_cyc==basecyc_last?i.num_cycles:basecyc_last
                        end
                    end
                    auto_cycs=basecyc_1st:basecyc_last

                    # my method, also worked:)
                    # cyc_idx_=Int(floor(cycles_denser[findmax(dr2_pred_tm)[2]]))-1
                    # println(cyc_idx_)
                    # if cyc_idx_ <= i.min_reliable_cyc
                    #     auto_cycs=cyc_idx_==0?bl_cycs():(1:cyc_idx_)
                    #     baseline = i.SFC_model_def_func(i.bl_method).funcs_pred[:bl](
                    #                          cycles_denser, bl_fit.coefs...)[1]
                    # else
                    #     auto_cycs=i.min_reliable_cyc:cyc_idx_
                    #     baseline = median(fluos[auto_cycs])
                    # end
                else
                    auto_cycs=i.min_reliable_cyc:i.num_cycles
                end
            end
            baseline = median(fluos[auto_cycs])
        else
            auto_cycs=i.baseline_cyc_bounds[1]:i.baseline_cyc_bounds[2]
            baseline = median(fluos[auto_cycs])
        end

        blsub_fluos= fluos .-baseline
        blsub_fluos_flb=blsub_fluos

        # debug(logger, "bl notes: $bl_notes, good_results: $have_good_results")
    else
        ## bl_method == take_the_median
        ## do not fit model to find baseline
        wts = ones(num_cycles)
        # bl_fit = FIT[amp_model]() ## empty model fit
        # have_good_results = false
        bl_notes = ["no bl_status", "no fallback"]
        bl_status = ""
        bl_func = median
        baseline = bl_func(fluos[bl_cycs()])
        blsub_fluos= fluos .- baseline
        blsub_fluos_flb=blsub_fluos
    end ## if bl_method
    # if !have_good_results
    #     blsub_fluos = blsub_fluos_flb
    # end
    #
    ## fit quantitation model
    qm = i.SFC_model_def_func(i.quant_method)
    quant_fit = qm.func_fit(cycles[auto_cycs[1]:end], blsub_fluos[auto_cycs[1]:end], wts; solver = i.solver)
    quant_coefs = quant_fit.coefs
    len_denser = denser_len(i, num_cycles)
    cycles_denser = interpolated_cycles()
    funcs_pred = qm.funcs_pred
	#if quant_fit.status==Symbol("Optimal")

    dr1_pred = dr1_pred_()
    dr2_pred = dr2_pred_()

	#end
    #
    ## results by R
    if R == AmpCqFluoModelResults
        ## the following method for calculating cq_raw
        ## is efficient for DEFAULT_AMP_CQ_FLUO_METHOD = cp_dr1
        cq_raw = calc_cq_raw(dr1_pred_())
        return AmpCqFluoModelResults(
            quant_fit.status, ## quant_status
            calc_cq_fluo(cq_raw)) ## cq_fluo
    end ## if
    #

    function discrete_dr1_dr2(fluos)
        len=length(fluos)
        dr1_val=zeros(len)
        dr2_val=zeros(len)
        dr2_val_v2=zeros(len)
        for el=1:len
            if el == 1
                dr1_val[el]=(fluos[el+1]-fluos[el])
            elseif el == len
                dr1_val[el]=(fluos[el]-fluos[el-1])
            else
                dr1_val[el]=(fluos[el+1]-fluos[el-1])/2
            end
        end
        for el=1:len
            if el == 1
                dr2_val[el]=(dr1_val[el+1]-dr1_val[el])
                dr2_val_v2[el]=dr2_val[el]
            elseif el == len
                dr2_val[el]=(dr1_val[el]-dr1_val[el-1])
                dr2_val_v2[el]=dr2_val[el]
            else
                dr2_val[el]=(fluos[el+1]+fluos[el-1]-2*fluos[el])
                dr2_val_v2[el]=(dr1_val[el+1]-dr1_val[el-1])/2
            end
        end
        return (dr1_val,dr2_val,dr2_val_v2)
    end
    #criteria_norm_err

    if R == AmpShortModelResults
        raw_cycs_index = colon(1, i.denser_factor, len_denser)
        dr1_values1,dr2_values1,dr2_values1_v2=discrete_dr1_dr2(fluos)
        return AmpShortModelResults(
            # fluos, ## rbbs_ary3,
            blsub_fluos,
            blsub_fluos_flb,
            dr1_pred[raw_cycs_index],
            dr2_pred[raw_cycs_index],
            dr1_values1,
            dr2_values1_v2,
            quant_fit.status,
            quant_coefs,
            nonpos2NaN(calc_cq_raw(dr1_pred, dr2_pred)), ## cq
            NaN_T) ## d0
    end ## R
    #
    ## R == AmpLongModelResults
    (max_dr1, idx_max_dr1) = findmax(dr1_pred)
    cyc_max_dr1 = cycles_denser[idx_max_dr1]
    (max_dr2, idx_max_dr2) = findmax(dr2_pred)
    cyc_max_dr2 = cycles_denser[idx_max_dr2]
    eff_pred = map(func_pred_eff, cycles_denser)
    (eff_max, idx_max_eff) = findmax(eff_pred)
    cyc_vals_4cq = OrderedDict(
        :cp_dr1  => cyc_max_dr1,
        :cp_dr2  => cyc_max_dr2,
        :Cy0     => cy0_(max_dr1, idx_max_dr1, cyc_max_dr1),
        :ct      => ct_(),
        :max_eff => cycles_denser[idx_max_eff])
    cq_raw = cyc_vals_4cq[Symbol(cq_method)]
    eff_vals_4cq =
        OrderedDict(
            map(keys(cyc_vals_4cq)) do key
                key => (key == :max_eff) ?
                    eff_max :
                    func_pred_eff(cyc_vals_4cq[key])
            end)
    return AmpLongModelResults(
        fluos, ## rbbs_ary3
        bl_fit,
        bl_notes,
        blsub_fluos,
        blsub_fluos_flb,
        quant_fit,
        quant_fit.status, ## quant_status
        quant_coefs,
        NaN_T, ## d0
        funcs_pred[:f](cycles, quant_coefs...), ## blsub_fitted
        dr1_pred[raw_cycs_index],
        dr2_pred[raw_cycs_index],
        dr1_values1,
        dr2_values1_v2,
        max_dr1,
        max_dr2,
        cyc_vals_4cq,
        eff_vals_4cq,
        cq_raw,
        nonpos2NaN(cq_raw), ## cq
        copy(eff_vals_4cq[Symbol(cq_method)]), ## eff
        calc_cq_fluo(cq_raw)) ## cq_fluo
end ## fit_amplification_model(Val{SFCModel})

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
struct BIO
    dict::Dict
    tags::Vector
end
BIO() = BIO(Dict("O"=>1), ["O"])

Base.length(tagset::BIO) = length(tagset.tags)

function encode(tagset::BIO, tags::Vector{String})
    basetag = ""
    map(tags) do tag
        tag == "O" && return 1
        if tag == "_"
            tag = "I-" * basetag
        else
            basetag = tag
            tag = "B-" * basetag
        end
        get!(tagset.dict, tag) do
            id = length(tagset.dict) + 1
            push!(tagset.tags, tag)
            id
        end
    end
end

function decode(tagset::BIO, ids::Vector{Int})
    spans = Tuple{Int,Int,String}[]
    bpos = 0
    for i = 1:length(ids)
        tag = tagset.tags[ids[i]]
        tag == "O" && continue
        startswith(tag, "B-") && (bpos = i)
        nexttag = i == length(ids) ? "O" : tagset.tags[ids[i+1]]
        if !startswith(nexttag,"I-") && bpos > 0
            basetag = tagset.tags[ids[bpos]][3:end]
            push!(spans, (bpos,i,basetag))
            bpos = 0
        end
    end
    spans
end

struct BIO2
    dict::Dict
    tags::Vector
end
BIO2() = BIO2(Dict(), String[])

Base.length(tagset::BIO2) = length(tagset.tags)

function encode(tagset::BIO2, tags::Vector{String})
    basetag = ""
    ids = Int[]
    for tag in tags
        if tag == "O"
        elseif tag == "_"
            tag = "I-" * basetag
        else
            basetag = tag
            tag = "B-" * basetag
        end
        for suffix in ("-1","-2")
            t = tag * suffix
            id = get!(tagset.dict, t) do
                id = length(tagset.dict) + 1
                push!(tagset.tags, t)
                id
            end
            suffix == "-1" && push!(ids,id)
        end
    end
    ids
end

function decode(tagset::BIO2, ids::Vector{Int})
    spans = Tuple{Int,Int,String}[]
    bpos = 0
    for i = 1:length(ids)
        tag = tagset.tags[ids[i]]
        startswith(tag, "O-") && continue
        startswith(tag, "B-") && (bpos = i)
        nexttag = i == length(ids) ? "O" : tagset.tags[ids[i+1]]
        if !startswith(nexttag,"I-") && bpos > 0
            strs = Vector{String}(split(tagset.tags[ids[bpos]], "-"))
            basetag = strs[2]
            push!(spans, (bpos,i,basetag))
            bpos = 0
        end
    end
    spans
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


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

using Markdown
using InteractiveUtils

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

# ╔═╡ c7383ff4-ba4c-11eb-1977-b31b330b20d0
begin
	import Pkg
	Pkg.activate(@__DIR__)
	Pkg.instantiate()
	using Plots
	using HDF5
	using Unitful
	using UnitfulRecipes
	using PlutoUI
end

# ╔═╡ 3e85f1c4-32ca-4e8e-ab9f-2fefaaabffb7
md"## Load HDF5 output file"

# ╔═╡ 72c703c9-3c32-45ab-b910-25a3c6652bdc
filename = "particle_da.h5"

# ╔═╡ 7cbbad5a-3767-422c-9d6b-f32524d4bf04
fh = h5open(filename, "r")

# ╔═╡ 7c316637-defe-4f60-a7e5-2e10511f7044
md"## Set these parameters to choose what to plot"

# ╔═╡ e1acfe11-0ead-4ada-946e-eff674e6d44e
begin
	timestamps = keys(fh["data_syn"])
	md"""
	Select the timestamp
	$(@bind timestamp_idx Slider(1:length(timestamps)))
	"""
end

# ╔═╡ 775fbb1e-a760-4862-829d-455051942255
timestamp = timestamps[timestamp_idx]

# ╔═╡ d442ef52-1566-484d-af31-ba3565307502
md"""
Select the field
$(@bind field Select([f => f for f in keys(fh["data_syn"]["t0000"])]))
"""

# ╔═╡ 384116c1-ed19-4a8c-85a7-6c0dfd9c164f
md"## Contour plots of surface height"

# ╔═╡ 8f67e2b3-9a01-42a3-a56a-11b84776a5e1
md"## Scatter plot of particle weights"

# ╔═╡ a35f5895-5066-4e1d-b252-6172888aa92d
begin
	weights = read(fh["weights"][timestamp])

	p1 = scatter(weights, marker=:star)
	p2 = scatter(weights, marker=:star, yscale=:log10)

	for plt in (p1, p2)
	    plot!(plt; xlabel="Particle ID", ylabel="Weight")
	end

	plot(p1, p2, label="")
end

# ╔═╡ e5335e14-bdb6-432c-94ab-c666c304efc6
md"## Time series of Estimated Sample Size"

# ╔═╡ 343a1d50-38f8-4457-81dc-5d962a2acb4a
plot([1 / sum(read(w) .^ 2) for w in fh["weights"]];
     label="", marker=:o, xlabel="Time step", ylabel="Estimated Sample Size (1 / sum(weight^2))")

# ╔═╡ a52a3f7e-1d8e-4153-b0c2-2cb47584c447
md"## Animation"

# ╔═╡ 8520dcbb-0bd8-4020-aea3-009e24df2099
md"## Collect data from the output file"

# ╔═╡ cff1a64f-03ba-4150-9501-fa4803901808
# All time-independent quantities
begin
	field_unit = read(fh["data_syn"]["t0000"][field]["Unit"])
	var_unit = read(fh["data_var"]["t0000"][field]["Unit"])

	field_desc = read(fh["data_syn"]["t0000"][field]["Description"])

	x_unit = read(fh["grid"]["x"]["Unit"])
	y_unit = read(fh["grid"]["y"]["Unit"])
	x_st_unit = read(fh["stations"]["x"]["Unit"])
	y_st_unit = read(fh["stations"]["y"]["Unit"])

	x = read(fh["grid"]["x"]) .* uparse(x_unit) .|> u"km"
	y = read(fh["grid"]["y"]) .* uparse(y_unit) .|> u"km"

	x_st = read(fh["stations"]["x"]) .* uparse(x_st_unit) .|> u"km"
	y_st = read(fh["stations"]["y"]) .* uparse(y_st_unit) .|> u"km"
end

# ╔═╡ e17da3a7-39e4-4326-aa2e-08f70b574878
function plot_data(x, y, z_t, z_avg, z_std, field_desc)
    n_contours = 100
    zmax = max(maximum(z_t), maximum(z_avg))
    zmin = min(minimum(z_t), minimum(z_avg))
    levels = range(zmin, zmax; length=n_contours)

    # Note that for heatmaps we need to permute the dimensions of the z matrix
    p1 = heatmap(x, y, z_t'; title="True $(lowercase(field_desc))")
    p2 = heatmap(x, y, z_avg'; title="Assimilated $(lowercase(field_desc))")
    p3 = heatmap(x, y, z_std'; title="Std of assimilated $(lowercase(field_desc))")

    for (i, plt) in enumerate((p1, p2, p3))
        # Set labels
        plot!(plt; xlabel="x", ylabel="y")
        # Set range of color bar for first two plots
        i ∈ (1, 2) && plot!(plt; clims=(ustrip(zmin), ustrip(zmax)))
        # Add the positions of the stations
        scatter!(plt, x_st, y_st, color=:red, marker=:star, label="")
    end

    plot(p1, p2, p3; titlefontsize=8, guidefontsize=8)
end

# ╔═╡ ba33b9a1-7d73-4247-b298-ccf30acc8859
function animate_data(fh, field, field_unit, var_unit, x, y)
	animation = @animate for timestamp ∈ keys(fh["data_syn"])
	    z_t = read(fh["data_syn"][timestamp][field]) .* uparse(field_unit)
	    z_avg = read(fh["data_avg"][timestamp][field]) .* uparse(field_unit)
	    z_std = sqrt.(read(fh["data_var"][timestamp][field]) .* uparse(var_unit))

	    plot_data(x, y, z_t, z_avg, z_std, field_desc)
	end

	return mp4(animation, "animation_jl.mp4"; fps=5)
end

# ╔═╡ a9343779-de40-4d33-8487-27d53ec095c0
animate_data(fh, field, field_unit, var_unit, x, y)

# ╔═╡ da1315e0-71de-4df6-9d74-259979571e1e
# Quantities specific to the current timestamp
begin
	z_t = read(fh["data_syn"][timestamp][field]) .* uparse(field_unit)
	z_avg = read(fh["data_avg"][timestamp][field]) .* uparse(field_unit)
	z_std = sqrt.(read(fh["data_var"][timestamp][field]) .* uparse(var_unit))
end

# ╔═╡ 1d230245-5f29-4895-b0cb-4e49f6c125ff
plot_data(x, y, z_t, z_avg, z_std, field_desc)

# ╔═╡ Cell order:
# ╠═c7383ff4-ba4c-11eb-1977-b31b330b20d0
# ╟─3e85f1c4-32ca-4e8e-ab9f-2fefaaabffb7
# ╠═72c703c9-3c32-45ab-b910-25a3c6652bdc
# ╟─7cbbad5a-3767-422c-9d6b-f32524d4bf04
# ╟─7c316637-defe-4f60-a7e5-2e10511f7044
# ╟─e1acfe11-0ead-4ada-946e-eff674e6d44e
# ╟─775fbb1e-a760-4862-829d-455051942255
# ╟─d442ef52-1566-484d-af31-ba3565307502
# ╟─384116c1-ed19-4a8c-85a7-6c0dfd9c164f
# ╟─e17da3a7-39e4-4326-aa2e-08f70b574878
# ╟─1d230245-5f29-4895-b0cb-4e49f6c125ff
# ╟─8f67e2b3-9a01-42a3-a56a-11b84776a5e1
# ╟─a35f5895-5066-4e1d-b252-6172888aa92d
# ╟─e5335e14-bdb6-432c-94ab-c666c304efc6
# ╟─343a1d50-38f8-4457-81dc-5d962a2acb4a
# ╟─a52a3f7e-1d8e-4153-b0c2-2cb47584c447
# ╟─ba33b9a1-7d73-4247-b298-ccf30acc8859
# ╟─a9343779-de40-4d33-8487-27d53ec095c0
# ╟─8520dcbb-0bd8-4020-aea3-009e24df2099
# ╟─cff1a64f-03ba-4150-9501-fa4803901808
# ╟─da1315e0-71de-4df6-9d74-259979571e1e


text:
@testset "custom nnz (number of non-zeros)" begin

    using Calculus

    L = 3
    N = 2*L

    X1 = [n*l for n=1:N, l=1:L]
    X2 = [n*l for n=1:N, l=1:L]
    β1 = range(-1., stop=1., length=L)
    β2 = ones(L)

    # equation depends on β2 only if β1 is positive, so jacobian will have varying number of elements depending on β1
    y = X1*β1 + X2 * ((β1.>0) .* β2)

    function resid(par)
        return y - (X1*par[1:L] + X2 * ((par[1:L].>0) .* par[L+1:end]))
    end
    lower = -500*ones(L*2)
    upper = 500*ones(L*2)
    guess = -300*ones(L*2)

    # number of non-zero elements in Jacobian for negative/positive guess
    nnz1 = PATHSolver.count_nonzeros(Calculus.jacobian(resid, -ones(L*2), :central))
    nnz2 = PATHSolver.count_nonzeros(Calculus.jacobian(resid, ones(L*2), :central))
    @test nnz2>nnz1

    jac(par) = Calculus.jacobian(resid, par, :central)
    # Calculus jacobian used
    status, zero, f_zero = solveMCP(resid, jac, lower, upper, guess; nnz=length(lower)^2)
    @test status == :Solved

    # ForwardDiff jacobian
    status, zero, f_zero = solveMCP(resid, lower, upper, guess; nnz=length(lower)^2)
    @test status == :Solved

end


text:
using HydrogenDimerMRCC

r0 = 0.776777
R  = 6.0
theta_values = range(0.0, stop = pi, length = 7)

dimer1 = PolarCoords(r0, 0.0, 0.0)
for (i, theta) in enumerate(theta_values)
	dimer2 = PolarCoords(r0, theta, 0.0)
	dparam = DimerParameters(R, dimer1, dimer2)
	make_mrcc_input_file(dparam, "MRCCinputfiles/example2/MINP$i")
end


text:
import JSON

function read_response(
    path = "./data/sample_response_direction.json"
    )
    response = String(read(path))
    response_json = JSON.parse(response)
    return response_json
end

function predict(clustor_lat, clustor_lng, lat, lng, time_frame)
    virus_factor = rand(1:100) / 100
    response = read_response()
    route = first(response["routes"])
    leg = first(route["legs"])
    eta = leg["duration"]["value"]
    return eta * virus_factor / time_frame
end

max_issues_repo_name: MLH-Fellowship/CodeVidLive


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 876e6928-abdc-11eb-0fea-d51de2d4b971
md""" ### Quick Julia tutorial...
"""

# ╔═╡ 9e95da34-fc21-4e7e-99a9-aea83edacd19
# global variable
example = 9.5

# ╔═╡ 34282b28-b73f-4aa7-96be-31dd6700084f
areaofcircle(r) = pi * r^2

# ╔═╡ 451499e9-76d9-4088-af2f-f35a9222afca
#= this updates
~automatically~
if you change the first cell
(and here's how you do multi-line commenting!) =#

bep = areaofcircle(example)

# ╔═╡ 1a3a323e-3f6b-4d65-bcbe-7c1c98aeb150
md"""
##### Some words on scope...
In Julia, like all other languages, there are globally and locally scoped variables. Locally scoped can be considered hard or soft; difference being that soft local variables unassigned locally and defined globally reassign the global variable, whereas hard scope remains local.

Julia uses lexical scoping, meaning that a function's scope refers to the one in which it was defined. For example, calling the module "Bleep" from a global scope will use the scope independent to Bleep, NOT the caller's global scope:
"""

# ╔═╡ 3f2ba373-77c6-4462-b30b-94fc8602298e
module Bleep
	example = 1
	foo() = example
end

# ╔═╡ 6315a551-f39b-40be-8c84-d3fbf34de4d0
Bleep.foo()

# ╔═╡ fa32bf73-e41a-4db5-ab2a-3e5a9e383aef
md"""
And this new module's independent scope is separate from the caller's global scope:
"""

# ╔═╡ 0cd019ae-3521-4323-9371-020190603204
module Deeee
	beh = example
end

# ╔═╡ 92cfd5a5-09ef-401d-a4a6-53d2576b314f
md"""Also, changing a global variable with a function will result in that variable remaining modified - called 'pass-by-sharing'. Conventionally, these functions will have '!' after their names.
"""

# ╔═╡ 92713261-8080-4797-8e48-0f0c2ce0e391
function changer!(bees)
    bees .+= 1
    return
end

# ╔═╡ 997f3fe0-d920-4055-b81e-b259fce128a3
ty = [1.0,2.0,3.0]

# ╔═╡ b6fa822c-7a55-46fe-b3a3-f191fe3ebd3d
begin
	changer!(ty)
	ty
end

# ╔═╡ 883334fb-f81b-4c24-9527-571a26ce7ce8
md"""
Additionally, if one variable is defined as another, any modification to the new variable will result in the same change in the original.
"""

# ╔═╡ dc423ed7-c9ee-416a-a0b2-a3c1408221cd
begin
	blip = [1,2,3]
	behhh = blip
end

# ╔═╡ b9aab8cc-c738-4f3c-a845-52f5078d1478
behhh[2] = 4

# ╔═╡ 0b04f67a-e413-4c9a-b298-610a294561b8
blip

# ╔═╡ 251313ac-b425-47c5-aae8-d78e5167b052
md"""
Other fun things include binding with a variety of inputs (copied from their sample notebook)... click the eyeball on the left to see the cell!
"""

# ╔═╡ 9f9fb7dd-b0f6-4bb8-b523-6053b5bc39dc
md"""
`zoot = ` $(@bind zoot html"<input type=range min=5 max=35>")

`beef = ` $(@bind beef html"<input type=text >")

`crr = ` $(@bind crr html"<input type=button value='Click'>")

`drrt = ` $(@bind drrt html"<input type=checkbox >")

`efff = ` $(@bind efff html"<select><option value='one'>First</option><option value='two'>Second</option></select>")

`fooop = ` $(@bind fooop html"<input type=color >")

"""

# ╔═╡ b03cc059-d281-4973-a22a-fabeb5de06e0
(zoot, beef, crr, drrt, efff, fooop)

# ╔═╡ c423ec0a-dbf7-4400-99f9-8357ac59edd1
md"""
#### Working with arrays, concatenation, and plotting...
"""

# ╔═╡ d5e187ee-c6b0-4490-80e0-12a47e809d27
begin
	years = collect(2001:2005)
	oranges = [3, 7, 12, 89, 20]
	pears = zeros(5,1)
	apples = rand(1:100, 5)
	bananas = rand(1:100, 5)
end

# ╔═╡ d50d1e34-9ef2-41de-8dac-bb763685fd5f
fruits = hcat(apples, oranges, pears, bananas)

# ╔═╡ 2f248615-aaca-44ce-a0d9-0ac31e12ba52
begin
	using Plots
	plot(years, fruits, layout=(4, 1), legend=false)
end

# ╔═╡ 6658d1dc-aefd-449b-8d81-f2b9f6a6fe3c
sort(fruits, dims=2)

# ╔═╡ 3bc93ea9-0718-47b9-ae4c-d1286641c9f2
all(i->(4<=i<=6), [4,5,6])

# ╔═╡ 28108af6-a838-417b-85a5-5e3fff0ae49b
md"""
#### Fun with structs!

Julia structs do not require inheritance like OOP languages do - of which Julia is not considered to be, due to the fact that you cannot attach a method (a behavior of a function) to a type. Dispatch, or the choice of method, depends on the first argument in OOP, while in Julia, it chooses which of a function's methods to call depending on the types of all the arguments. This is called multiple dispatch, and is the primary differentiator from other languages.

Types can be either system-defined or user-defined, and allow for structs to define a "supertype" which can be called by any number of functions or structs, as follows.

"""

# ╔═╡ cca3129b-8ccb-4747-b75d-cb34f3ff5b23
struct Limbs
	hands :: Float64
	leg :: Integer
end

# ╔═╡ b1ce5d0c-c3f6-411c-8f04-c44854fa89b1
struct boy
	body_parts :: Limbs #declaring body_parts as a Composite Type, defined by Limbs
	# only inner constructors have access to this "new" function
	boy(hands, leg) = new(Limbs(hands, leg))
end

# ╔═╡ 26730fa4-254c-4bef-81e5-39803ebcaa2e
begin
	joseph = boy(2.,1)
	joseph.body_parts.hands
end

# ╔═╡ 921051ca-54dc-4563-a116-0756cca43557
md"""
This will come into play with user-defined type "Cell" that can be acted on by any number of functions.
"""

# ╔═╡ 8a3ee669-bad2-4086-b90a-1744b85dcf76
md"""
##### Differences from other languages that may not be immediately obvious...

https://docs.julialang.org/en/v1/manual/noteworthy-differences/

https://erik-engheim.medium.com/defining-custom-units-in-julia-and-python-513c34a4c971
"""

# ╔═╡ 39d444ec-f84a-4d1b-9ba4-19af636bf3dd

# ╔═╡ Cell order:
# ╟─876e6928-abdc-11eb-0fea-d51de2d4b971
# ╠═9e95da34-fc21-4e7e-99a9-aea83edacd19
# ╠═34282b28-b73f-4aa7-96be-31dd6700084f
# ╠═451499e9-76d9-4088-af2f-f35a9222afca
# ╟─1a3a323e-3f6b-4d65-bcbe-7c1c98aeb150
# ╠═3f2ba373-77c6-4462-b30b-94fc8602298e
# ╠═6315a551-f39b-40be-8c84-d3fbf34de4d0
# ╟─fa32bf73-e41a-4db5-ab2a-3e5a9e383aef
# ╠═0cd019ae-3521-4323-9371-020190603204
# ╟─92cfd5a5-09ef-401d-a4a6-53d2576b314f
# ╠═92713261-8080-4797-8e48-0f0c2ce0e391
# ╠═997f3fe0-d920-4055-b81e-b259fce128a3
# ╠═b6fa822c-7a55-46fe-b3a3-f191fe3ebd3d
# ╟─883334fb-f81b-4c24-9527-571a26ce7ce8
# ╠═dc423ed7-c9ee-416a-a0b2-a3c1408221cd
# ╠═b9aab8cc-c738-4f3c-a845-52f5078d1478
# ╠═0b04f67a-e413-4c9a-b298-610a294561b8
# ╟─251313ac-b425-47c5-aae8-d78e5167b052
# ╟─9f9fb7dd-b0f6-4bb8-b523-6053b5bc39dc
# ╠═b03cc059-d281-4973-a22a-fabeb5de06e0
# ╟─c423ec0a-dbf7-4400-99f9-8357ac59edd1
# ╠═d5e187ee-c6b0-4490-80e0-12a47e809d27
# ╠═d50d1e34-9ef2-41de-8dac-bb763685fd5f
# ╠═6658d1dc-aefd-449b-8d81-f2b9f6a6fe3c
# ╠═2f248615-aaca-44ce-a0d9-0ac31e12ba52
# ╠═3bc93ea9-0718-47b9-ae4c-d1286641c9f2
# ╟─28108af6-a838-417b-85a5-5e3fff0ae49b
# ╠═cca3129b-8ccb-4747-b75d-cb34f3ff5b23
# ╠═b1ce5d0c-c3f6-411c-8f04-c44854fa89b1
# ╠═26730fa4-254c-4bef-81e5-39803ebcaa2e
# ╟─921051ca-54dc-4563-a116-0756cca43557
# ╟─8a3ee669-bad2-4086-b90a-1744b85dcf76
# ╠═39d444ec-f84a-4d1b-9ba4-19af636bf3dd

max_stars_repo_name: andi-bergeson/simgym

max_issues_repo_name: andi-bergeson/simgym

max_forks_repo_name: andi-bergeson/simgym