text:
@testset "8.5 Hyperbolic integral functions" begin
    (a, b, c, d, e, m, n, x, ) = @variables a b c d e m n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integration*Problems*Involving*Hyperbolic*Integral*Functions=#

    #= ::Section::Closed:: =#
    #=Hyperbolic*sine*integral*function=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*SinhIntegral(b*x)^n=#

    @test_int [x^m*SinhIntegral(b*x), x, 5, -((x^m*Gamma(1 + m, (-b)*x))/(((-b)*x)^m*(2*b*(1 + m)))) - (x^m*Gamma(1 + m, b*x))/((b*x)^m*(2*b*(1 + m))) + (x^(1 + m)*SinhIntegral(b*x))/(1 + m)]

    @test_int [x^3*SinhIntegral(b*x), x, 6, -((3*x*cosh(b*x))/(2*b^3)) - (x^3*cosh(b*x))/(4*b) + (3*sinh(b*x))/(2*b^4) + (3*x^2*sinh(b*x))/(4*b^2) + (1/4)*x^4*SinhIntegral(b*x)]
    @test_int [x^2*SinhIntegral(b*x), x, 5, -((2*cosh(b*x))/(3*b^3)) - (x^2*cosh(b*x))/(3*b) + (2*x*sinh(b*x))/(3*b^2) + (1/3)*x^3*SinhIntegral(b*x)]
    @test_int [x^1*SinhIntegral(b*x), x, 4, -((x*cosh(b*x))/(2*b)) + sinh(b*x)/(2*b^2) + (1/2)*x^2*SinhIntegral(b*x)]
    @test_int [x^0*SinhIntegral(b*x), x, 1, -(cosh(b*x)/b) + x*SinhIntegral(b*x)]
    @test_int [SinhIntegral(b*x)/x^1, x, 1, (1/2)*b*x*HypergeometricPFQ([1, 1, 1], [2, 2, 2], (-b)*x) + (1/2)*b*x*HypergeometricPFQ([1, 1, 1], [2, 2, 2], b*x)]
    @test_int [SinhIntegral(b*x)/x^2, x, 4, b*CoshIntegral(b*x) - sinh(b*x)/x - SinhIntegral(b*x)/x]
    @test_int [SinhIntegral(b*x)/x^3, x, 5, -((b*cosh(b*x))/(4*x)) - sinh(b*x)/(4*x^2) + (1/4)*b^2*SinhIntegral(b*x) - SinhIntegral(b*x)/(2*x^2)]

    @test_int [x^m*SinhIntegral(b*x)^2, x, 0, CannotIntegrate(x^m*SinhIntegral(b*x)^2, x)]

    @test_int [x^3*SinhIntegral(b*x)^2, x, 19, x^2/(2*b^2) - (3*CoshIntegral(2*b*x))/(2*b^4) + (3*log(x))/(2*b^4) - (x*cosh(b*x)*sinh(b*x))/b^3 + (2*sinh(b*x)^2)/b^4 + (x^2*sinh(b*x)^2)/(4*b^2) - (3*x*cosh(b*x)*SinhIntegral(b*x))/b^3 - (x^3*cosh(b*x)*SinhIntegral(b*x))/(2*b) + (3*sinh(b*x)*SinhIntegral(b*x))/b^4 + (3*x^2*sinh(b*x)*SinhIntegral(b*x))/(2*b^2) + (1/4)*x^4*SinhIntegral(b*x)^2]
    @test_int [x^2*SinhIntegral(b*x)^2, x, 15, (5*x)/(6*b^2) - (5*cosh(b*x)*sinh(b*x))/(6*b^3) + (x*sinh(b*x)^2)/(3*b^2) - (4*cosh(b*x)*SinhIntegral(b*x))/(3*b^3) - (2*x^2*cosh(b*x)*SinhIntegral(b*x))/(3*b) + (4*x*sinh(b*x)*SinhIntegral(b*x))/(3*b^2) + (1/3)*x^3*SinhIntegral(b*x)^2 + (2*SinhIntegral(2*b*x))/(3*b^3)]
    @test_int [x^1*SinhIntegral(b*x)^2, x, 10, -(CoshIntegral(2*b*x)/(2*b^2)) + log(x)/(2*b^2) + sinh(b*x)^2/(2*b^2) - (x*cosh(b*x)*SinhIntegral(b*x))/b + (sinh(b*x)*SinhIntegral(b*x))/b^2 + (1/2)*x^2*SinhIntegral(b*x)^2]
    @test_int [x^0*SinhIntegral(b*x)^2, x, 6, -((2*cosh(b*x)*SinhIntegral(b*x))/b) + x*SinhIntegral(b*x)^2 + SinhIntegral(2*b*x)/b]
    @test_int [SinhIntegral(b*x)^2/x^1, x, 0, CannotIntegrate(SinhIntegral(b*x)^2/x, x)]
    @test_int [SinhIntegral(b*x)^2/x^2, x, 0, CannotIntegrate(SinhIntegral(b*x)^2/x^2, x)]
    @test_int [SinhIntegral(b*x)^2/x^3, x, 0, CannotIntegrate(SinhIntegral(b*x)^2/x^3, x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*SinhIntegral(a+b*x)^n=#

    @test_int [x^m*SinhIntegral(a + b*x), x, 1, -((b*CannotIntegrate((x^(1 + m)*sinh(a + b*x))/(a + b*x), x))/(1 + m)) + (x^(1 + m)*SinhIntegral(a + b*x))/(1 + m)]

    @test_int [x^3*SinhIntegral(a + b*x), x, 14, (a*cosh(a + b*x))/(2*b^4) + (a^3*cosh(a + b*x))/(4*b^4) - (3*x*cosh(a + b*x))/(2*b^3) - (a^2*x*cosh(a + b*x))/(4*b^3) + (a*x^2*cosh(a + b*x))/(4*b^2) - (x^3*cosh(a + b*x))/(4*b) + (3*sinh(a + b*x))/(2*b^4) + (a^2*sinh(a + b*x))/(4*b^4) - (a*x*sinh(a + b*x))/(2*b^3) + (3*x^2*sinh(a + b*x))/(4*b^2) - (a^4*SinhIntegral(a + b*x))/(4*b^4) + (1/4)*x^4*SinhIntegral(a + b*x)]
    @test_int [x^2*SinhIntegral(a + b*x), x, 10, -((2*cosh(a + b*x))/(3*b^3)) - (a^2*cosh(a + b*x))/(3*b^3) + (a*x*cosh(a + b*x))/(3*b^2) - (x^2*cosh(a + b*x))/(3*b) - (a*sinh(a + b*x))/(3*b^3) + (2*x*sinh(a + b*x))/(3*b^2) + (a^3*SinhIntegral(a + b*x))/(3*b^3) + (1/3)*x^3*SinhIntegral(a + b*x)]
    @test_int [x^1*SinhIntegral(a + b*x), x, 7, (a*cosh(a + b*x))/(2*b^2) - (x*cosh(a + b*x))/(2*b) + sinh(a + b*x)/(2*b^2) - (a^2*SinhIntegral(a + b*x))/(2*b^2) + (1/2)*x^2*SinhIntegral(a + b*x)]
    @test_int [x^0*SinhIntegral(a + b*x), x, 1, -(cosh(a + b*x)/b) + ((a + b*x)*SinhIntegral(a + b*x))/b]
    @test_int [SinhIntegral(a + b*x)/x^1, x, 0, CannotIntegrate(SinhIntegral(a + b*x)/x, x)]
    @test_int [SinhIntegral(a + b*x)/x^2, x, 7, (b*CoshIntegral(b*x)*sinh(a))/a + (b*cosh(a)*SinhIntegral(b*x))/a - (b*SinhIntegral(a + b*x))/a - SinhIntegral(a + b*x)/x]
    @test_int [SinhIntegral(a + b*x)/x^3, x, 11, (b^2*cosh(a)*CoshIntegral(b*x))/(2*a) - (b^2*CoshIntegral(b*x)*sinh(a))/(2*a^2) - (b*sinh(a + b*x))/(2*a*x) - (b^2*cosh(a)*SinhIntegral(b*x))/(2*a^2) + (b^2*sinh(a)*SinhIntegral(b*x))/(2*a) + (b^2*SinhIntegral(a + b*x))/(2*a^2) - SinhIntegral(a + b*x)/(2*x^2)]

    @test_int [x^m*SinhIntegral(a + b*x)^2, x, 0, CannotIntegrate(x^m*SinhIntegral(a + b*x)^2, x)]

    @test_int [x^2*SinhIntegral(a + b*x)^2, x, 39, (2*x)/(3*b^2) - (a*cosh(2*a + 2*b*x))/(3*b^3) + (x*cosh(2*a + 2*b*x))/(6*b^2) + (a*CoshIntegral(2*a + 2*b*x))/b^3 - (a*log(a + b*x))/b^3 - (2*cosh(a + b*x)*sinh(a + b*x))/(3*b^3) - sinh(2*a + 2*b*x)/(12*b^3) - (4*cosh(a + b*x)*SinhIntegral(a + b*x))/(3*b^3) - (2*a^2*cosh(a + b*x)*SinhIntegral(a + b*x))/(3*b^3) + (2*a*x*cosh(a + b*x)*SinhIntegral(a + b*x))/(3*b^2) - (2*x^2*cosh(a + b*x)*SinhIntegral(a + b*x))/(3*b) - (2*a*sinh(a + b*x)*SinhIntegral(a + b*x))/(3*b^3) + (4*x*sinh(a + b*x)*SinhIntegral(a + b*x))/(3*b^2) + (a^2*(a + b*x)*SinhIntegral(a + b*x)^2)/(3*b^3) - (a*x*(a + b*x)*SinhIntegral(a + b*x)^2)/(3*b^2) + (x^2*(a + b*x)*SinhIntegral(a + b*x)^2)/(3*b) + (2*SinhIntegral(2*a + 2*b*x))/(3*b^3) + (a^2*SinhIntegral(2*a + 2*b*x))/b^3]
    @test_int [x^1*SinhIntegral(a + b*x)^2, x, 17, cosh(2*a + 2*b*x)/(4*b^2) - CoshIntegral(2*a + 2*b*x)/(2*b^2) + log(a + b*x)/(2*b^2) + (a*cosh(a + b*x)*SinhIntegral(a + b*x))/b^2 - (x*cosh(a + b*x)*SinhIntegral(a + b*x))/b + (sinh(a + b*x)*SinhIntegral(a + b*x))/b^2 - (a*(a + b*x)*SinhIntegral(a + b*x)^2)/(2*b^2) + (x*(a + b*x)*SinhIntegral(a + b*x)^2)/(2*b) - (a*SinhIntegral(2*a + 2*b*x))/b^2]
    @test_int [x^0*SinhIntegral(a + b*x)^2, x, 5, -((2*cosh(a + b*x)*SinhIntegral(a + b*x))/b) + ((a + b*x)*SinhIntegral(a + b*x)^2)/b + SinhIntegral(2*a + 2*b*x)/b]
    @test_int [SinhIntegral(a + b*x)^2/x^1, x, 0, CannotIntegrate(SinhIntegral(a + b*x)^2/x, x)]
    @test_int [SinhIntegral(a + b*x)^2/x^2, x, 0, CannotIntegrate(SinhIntegral(a + b*x)^2/x^2, x)]
    @test_int [SinhIntegral(a + b*x)^2/x^3, x, 0, CannotIntegrate(SinhIntegral(a + b*x)^2/x^3, x)]

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

    @test_int [x^2*SinhIntegral(d*(a + b*log(c*x^n))), x, 7, (x^3*ExpIntegralEi(((3 - b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(6*ℯ^((3*a)/(b*n))*(c*x^n)^(3/n)) - (x^3*ExpIntegralEi(((3 + b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(6*ℯ^((3*a)/(b*n))*(c*x^n)^(3/n)) + (x^3*SinhIntegral(d*(a + b*log(c*x^n))))/3]
    @test_int [x^1*SinhIntegral(d*(a + b*log(c*x^n))), x, 7, (x^2*ExpIntegralEi(((2 - b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(4*ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)) - (x^2*ExpIntegralEi(((2 + b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(4*ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)) + (x^2*SinhIntegral(d*(a + b*log(c*x^n))))/2]
    @test_int [x^0*SinhIntegral(d*(a + b*log(c*x^n))), x, 7, (x*ExpIntegralEi(((1 - b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(2*ℯ^(a/(b*n))*(c*x^n)^n^(-1)) - (x*ExpIntegralEi(((1 + b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(2*ℯ^(a/(b*n))*(c*x^n)^n^(-1)) + x*SinhIntegral(d*(a + b*log(c*x^n)))]
    @test_int [SinhIntegral(d*(a + b*log(c*x^n)))/x^1, x, 3, -(cosh(d*(a + b*log(c*x^n)))/(b*d*n)) + ((a + b*log(c*x^n))*SinhIntegral(d*(a + b*log(c*x^n))))/(b*n)]
    @test_int [SinhIntegral(d*(a + b*log(c*x^n)))/x^2, x, 7, (ℯ^(a/(b*n))*(c*x^n)^n^(-1)*ExpIntegralEi(-(((1 - b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(2*x) - (ℯ^(a/(b*n))*(c*x^n)^n^(-1)*ExpIntegralEi(-(((1 + b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(2*x) - SinhIntegral(d*(a + b*log(c*x^n)))/x]
    @test_int [SinhIntegral(d*(a + b*log(c*x^n)))/x^3, x, 7, (ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)*ExpIntegralEi(-(((2 - b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(4*x^2) - (ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)*ExpIntegralEi(-(((2 + b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(4*x^2) - SinhIntegral(d*(a + b*log(c*x^n)))/(2*x^2)]

    @test_int [(e*x)^m*SinhIntegral(d*(a + b*log(c*x^n))), x, 7, (x*(e*x)^m*ExpIntegralEi(((1 + m - b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(2*ℯ^((a*(1 + m))/(b*n))*(1 + m)*(c*x^n)^((1 + m)/n)) - (x*(e*x)^m*ExpIntegralEi(((1 + m + b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(2*ℯ^((a*(1 + m))/(b*n))*(1 + m)*(c*x^n)^((1 + m)/n)) + ((e*x)^(1 + m)*SinhIntegral(d*(a + b*log(c*x^n))))/(e*(1 + m))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Hyper(b*x)*SinhIntegral(b*x)^n=#

    @test_int [sinh(b*x)*SinhIntegral(b*x)/x^3, x, 14, b^2*CoshIntegral(2*b*x) - (b*cosh(b*x)*sinh(b*x))/(2*x) - sinh(b*x)^2/(4*x^2) - (b*sinh(2*b*x))/(4*x) - (b*cosh(b*x)*SinhIntegral(b*x))/(2*x) - (sinh(b*x)*SinhIntegral(b*x))/(2*x^2) + (1/4)*b^2*SinhIntegral(b*x)^2]
    @test_int [sinh(b*x)*SinhIntegral(b*x)/x^2, x, 5, b*CannotIntegrate((cosh(b*x)*SinhIntegral(b*x))/x, x) - sinh(b*x)^2/x - (sinh(b*x)*SinhIntegral(b*x))/x + b*SinhIntegral(2*b*x)]
    @test_int [sinh(b*x)*SinhIntegral(b*x)/x, x, 1, (1/2)*SinhIntegral(b*x)^2]
    @test_int [sinh(b*x)*SinhIntegral(b*x), x, 5, (cosh(b*x)*SinhIntegral(b*x))/b - SinhIntegral(2*b*x)/(2*b)]
    @test_int [x*sinh(b*x)*SinhIntegral(b*x), x, 9, CoshIntegral(2*b*x)/(2*b^2) - log(x)/(2*b^2) - sinh(b*x)^2/(2*b^2) + (x*cosh(b*x)*SinhIntegral(b*x))/b - (sinh(b*x)*SinhIntegral(b*x))/b^2]
    @test_int [x^2*sinh(b*x)*SinhIntegral(b*x), x, 14, -((5*x)/(4*b^2)) + (5*cosh(b*x)*sinh(b*x))/(4*b^3) - (x*sinh(b*x)^2)/(2*b^2) + (2*cosh(b*x)*SinhIntegral(b*x))/b^3 + (x^2*cosh(b*x)*SinhIntegral(b*x))/b - (2*x*sinh(b*x)*SinhIntegral(b*x))/b^2 - SinhIntegral(2*b*x)/b^3]
    @test_int [x^3*sinh(b*x)*SinhIntegral(b*x), x, 18, -(x^2/b^2) + (3*CoshIntegral(2*b*x))/b^4 - (3*log(x))/b^4 + (2*x*cosh(b*x)*sinh(b*x))/b^3 - (4*sinh(b*x)^2)/b^4 - (x^2*sinh(b*x)^2)/(2*b^2) + (6*x*cosh(b*x)*SinhIntegral(b*x))/b^3 + (x^3*cosh(b*x)*SinhIntegral(b*x))/b - (6*sinh(b*x)*SinhIntegral(b*x))/b^4 - (3*x^2*sinh(b*x)*SinhIntegral(b*x))/b^2]

    @test_int [cosh(b*x)*SinhIntegral(b*x)/x^3, x, 12, -((b*cosh(2*b*x))/(4*x)) + (1/2)*b^2*CannotIntegrate((cosh(b*x)*SinhIntegral(b*x))/x, x) - (b*sinh(b*x)^2)/(2*x) - sinh(2*b*x)/(8*x^2) - (cosh(b*x)*SinhIntegral(b*x))/(2*x^2) - (b*sinh(b*x)*SinhIntegral(b*x))/(2*x) + b^2*SinhIntegral(2*b*x)]
    @test_int [cosh(b*x)*SinhIntegral(b*x)/x^2, x, 7, b*CoshIntegral(2*b*x) - sinh(2*b*x)/(2*x) - (cosh(b*x)*SinhIntegral(b*x))/x + (1/2)*b*SinhIntegral(b*x)^2]
    @test_int [cosh(b*x)*SinhIntegral(b*x)/x, x, 0, CannotIntegrate((cosh(b*x)*SinhIntegral(b*x))/x, x)]
    @test_int [cosh(b*x)*SinhIntegral(b*x), x, 5, -(CoshIntegral(2*b*x)/(2*b)) + log(x)/(2*b) + (sinh(b*x)*SinhIntegral(b*x))/b]
    @test_int [x*cosh(b*x)*SinhIntegral(b*x), x, 9, x/(2*b) - (cosh(b*x)*sinh(b*x))/(2*b^2) - (cosh(b*x)*SinhIntegral(b*x))/b^2 + (x*sinh(b*x)*SinhIntegral(b*x))/b + SinhIntegral(2*b*x)/(2*b^2)]
    @test_int [x^2*cosh(b*x)*SinhIntegral(b*x), x, 13, x^2/(4*b) - CoshIntegral(2*b*x)/b^3 + log(x)/b^3 - (x*cosh(b*x)*sinh(b*x))/(2*b^2) + (5*sinh(b*x)^2)/(4*b^3) - (2*x*cosh(b*x)*SinhIntegral(b*x))/b^2 + (2*sinh(b*x)*SinhIntegral(b*x))/b^3 + (x^2*sinh(b*x)*SinhIntegral(b*x))/b]
    @test_int [x^3*cosh(b*x)*SinhIntegral(b*x), x, 20, (4*x)/b^3 + x^3/(6*b) - (4*cosh(b*x)*sinh(b*x))/b^4 - (x^2*cosh(b*x)*sinh(b*x))/(2*b^2) + (2*x*sinh(b*x)^2)/b^3 - (6*cosh(b*x)*SinhIntegral(b*x))/b^4 - (3*x^2*cosh(b*x)*SinhIntegral(b*x))/b^2 + (6*x*sinh(b*x)*SinhIntegral(b*x))/b^3 + (x^3*sinh(b*x)*SinhIntegral(b*x))/b + (3*SinhIntegral(2*b*x))/b^4]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Hyper(b*x)*SinhIntegral(d*x)^n=#

    @test_int [sinh(5*x)*SinhIntegral(2*x), x, 6, (1/5)*cosh(5*x)*SinhIntegral(2*x) + (1/10)*SinhIntegral(3*x) - (1/10)*SinhIntegral(7*x)]

    @test_int [cosh(5*x)*SinhIntegral(2*x), x, 6, (1/10)*CoshIntegral(3*x) - (1/10)*CoshIntegral(7*x) + (1/5)*sinh(5*x)*SinhIntegral(2*x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Hyper(a+b*x)*SinhIntegral(a+b*x)^n=#

    #= [x^3*sinh(a + b*x)*SinhIntegral(a + b*x), x, 32, (7*a*x)/(4*b^3) - x^2/b^2 + (3*cosh(a + b*x)^2)/b^4 - (a^2*cosh(a + b*x)^2)/(2*b^4) - (3*CosIntegral(2*a + 2*b*x))/b^4 + (3*a^2*CosIntegral(2*a + 2*b*x))/(2*b^4) + (3*log(a + b*x))/b^4 - (3*a^2*log(a + b*x))/(2*b^4) - (7*a*cosh(a + b*x)*sinh(a + b*x))/(4*b^4) + (2*x*cosh(a + b*x)*sinh(a + b*x))/b^3 - sinh(a + b*x)^2/b^4 - (a*x*sinh(a + b*x)^2)/(2*b^3) + (x^2*sinh(a + b*x)^2)/(2*b^2) + (6*x*cosh(a + b*x)*SinhIntegral(a + b*x))/b^3 - (x^3*cosh(a + b*x)*SinhIntegral(a + b*x))/b - (6*sinh(a + b*x)*SinhIntegral(a + b*x))/b^4 + (3*x^2*sinh(a + b*x)*SinhIntegral(a + b*x))/b^2 + (3*a*SinhIntegral(2*a + 2*b*x))/b^4 - (a^3*SinhIntegral(2*a + 2*b*x))/(2*b^4)] =#
    @test_int [x^2*sinh(a + b*x)*SinhIntegral(a + b*x), x, 21, -(x/b^2) + (a*cosh(2*a + 2*b*x))/(4*b^3) - (x*cosh(2*a + 2*b*x))/(4*b^2) - (a*CoshIntegral(2*a + 2*b*x))/b^3 + (a*log(a + b*x))/b^3 + (cosh(a + b*x)*sinh(a + b*x))/b^3 + sinh(2*a + 2*b*x)/(8*b^3) + (2*cosh(a + b*x)*SinhIntegral(a + b*x))/b^3 + (x^2*cosh(a + b*x)*SinhIntegral(a + b*x))/b - (2*x*sinh(a + b*x)*SinhIntegral(a + b*x))/b^2 - SinhIntegral(2*a + 2*b*x)/b^3 - (a^2*SinhIntegral(2*a + 2*b*x))/(2*b^3)]
    @test_int [x^1*sinh(a + b*x)*SinhIntegral(a + b*x), x, 11, -(cosh(2*a + 2*b*x)/(4*b^2)) + CoshIntegral(2*a + 2*b*x)/(2*b^2) - log(a + b*x)/(2*b^2) + (x*cosh(a + b*x)*SinhIntegral(a + b*x))/b - (sinh(a + b*x)*SinhIntegral(a + b*x))/b^2 + (a*SinhIntegral(2*a + 2*b*x))/(2*b^2)]
    @test_int [x^0*sinh(a + b*x)*SinhIntegral(a + b*x), x, 4, (cosh(a + b*x)*SinhIntegral(a + b*x))/b - SinhIntegral(2*a + 2*b*x)/(2*b)]
    @test_int [sinh(a + b*x)*SinhIntegral(a + b*x)/x^1, x, 0, CannotIntegrate((sinh(a + b*x)*SinhIntegral(a + b*x))/x, x)]

    #= [x^3*cosh(a + b*x)*SinhIntegral(a + b*x), x, 32, (4*x)/b^3 - (a^2*x)/(2*b^3) + (a*x^2)/(4*b^2) - x^3/(6*b) - (3*a*cosh(a + b*x)^2)/(2*b^4) + (3*a*CosIntegral(2*a + 2*b*x))/b^4 - (a^3*CosIntegral(2*a + 2*b*x))/(2*b^4) - (3*a*log(a + b*x))/b^4 + (a^3*log(a + b*x))/(2*b^4) - (4*cosh(a + b*x)*sinh(a + b*x))/b^4 + (a^2*cosh(a + b*x)*sinh(a + b*x))/(2*b^4) - (a*x*cosh(a + b*x)*sinh(a + b*x))/(2*b^3) + (x^2*cosh(a + b*x)*sinh(a + b*x))/(2*b^2) + (a*sinh(a + b*x)^2)/(4*b^4) - (2*x*sinh(a + b*x)^2)/b^3 - (6*cosh(a + b*x)*SinhIntegral(a + b*x))/b^4 + (3*x^2*cosh(a + b*x)*SinhIntegral(a + b*x))/b^2 - (6*x*sinh(a + b*x)*SinhIntegral(a + b*x))/b^3 + (x^3*sinh(a + b*x)*SinhIntegral(a + b*x))/b + (3*SinhIntegral(2*a + 2*b*x))/b^4 - (3*a^2*SinhIntegral(2*a + 2*b*x))/(2*b^4)] =#
    @test_int [x^2*cosh(a + b*x)*SinhIntegral(a + b*x), x, 21, -((a*x)/(2*b^2)) + x^2/(4*b) + cosh(2*a + 2*b*x)/(2*b^3) - CoshIntegral(2*a + 2*b*x)/b^3 - (a^2*CoshIntegral(2*a + 2*b*x))/(2*b^3) + log(a + b*x)/b^3 + (a^2*log(a + b*x))/(2*b^3) + (a*cosh(a + b*x)*sinh(a + b*x))/(2*b^3) - (x*cosh(a + b*x)*sinh(a + b*x))/(2*b^2) + sinh(a + b*x)^2/(4*b^3) - (2*x*cosh(a + b*x)*SinhIntegral(a + b*x))/b^2 + (2*sinh(a + b*x)*SinhIntegral(a + b*x))/b^3 + (x^2*sinh(a + b*x)*SinhIntegral(a + b*x))/b - (a*SinhIntegral(2*a + 2*b*x))/b^3]
    @test_int [x^1*cosh(a + b*x)*SinhIntegral(a + b*x), x, 12, x/(2*b) + (a*CoshIntegral(2*a + 2*b*x))/(2*b^2) - (a*log(a + b*x))/(2*b^2) - (cosh(a + b*x)*sinh(a + b*x))/(2*b^2) - (cosh(a + b*x)*SinhIntegral(a + b*x))/b^2 + (x*sinh(a + b*x)*SinhIntegral(a + b*x))/b + SinhIntegral(2*a + 2*b*x)/(2*b^2)]
    @test_int [x^0*cosh(a + b*x)*SinhIntegral(a + b*x), x, 4, -(CoshIntegral(2*a + 2*b*x)/(2*b)) + log(a + b*x)/(2*b) + (sinh(a + b*x)*SinhIntegral(a + b*x))/b]
    @test_int [cosh(a + b*x)*SinhIntegral(a + b*x)/x^1, x, 0, CannotIntegrate((cosh(a + b*x)*SinhIntegral(a + b*x))/x, x)]

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

    #= [x^2*sinh(a + b*x)*SinhIntegral(c + d*x), x, 46, -((c*cosh(a - c + (b - d)*x))/(2*b*(b - d)*d)) + (x*cosh(a - c + (b - d)*x))/(2*b*(b - d)) + (c*cosh(a + c + (b + d)*x))/(2*b*d*(b + d)) - (x*cosh(a + c + (b + d)*x))/(2*b*(b + d)) + (CoshIntegral(((b - d)*(c + d*x))/d)*(2*b*c*d*cosh(a - (b*c)/d) + (b^2*c^2 + 2*d^2)*sinh(a - (b*c)/d)))/(2*b^3*d^2) - (CoshIntegral(((b + d)*(c + d*x))/d)*(2*b*c*d*cosh(a - (b*c)/d) + (b^2*c^2 + 2*d^2)*sinh(a - (b*c)/d)))/(2*b^3*d^2) - sinh(a - c + (b - d)*x)/(2*b*(b - d)^2) - sinh(a - c + (b - d)*x)/(b^2*(b - d)) + sinh(a + c + (b + d)*x)/(2*b*(b + d)^2) + sinh(a + c + (b + d)*x)/(b^2*(b + d)) + (((2 + b^2*x^2)*cosh(a + b*x) - 2*b*x*sinh(a + b*x))*SinhIntegral(c + d*x))/b^3 + (((b^2*c^2 + 2*d^2)*cosh(a - (b*c)/d) + 2*b*c*d*sinh(a - (b*c)/d))*SinhIntegral(((b - d)*(c + d*x))/d))/(2*b^3*d^2) - (((b^2*c^2 + 2*d^2)*cosh(a - (b*c)/d) + 2*b*c*d*sinh(a - (b*c)/d))*SinhIntegral(((b + d)*(c + d*x))/d))/(2*b^3*d^2)] =#
    @test_int [x^1*sinh(a + b*x)*SinhIntegral(c + d*x), x, 24, cosh(a - c + (b - d)*x)/(2*b*(b - d)) - cosh(a + c + (b + d)*x)/(2*b*(b + d)) - (cosh(a - (b*c)/d)*CoshIntegral((c*(b - d))/d + (b - d)*x))/(2*b^2) + (cosh(a - (b*c)/d)*CoshIntegral((c*(b + d))/d + (b + d)*x))/(2*b^2) - (c*CoshIntegral((c*(b - d))/d + (b - d)*x)*sinh(a - (b*c)/d))/(2*b*d) + (c*CoshIntegral((c*(b + d))/d + (b + d)*x)*sinh(a - (b*c)/d))/(2*b*d) - (c*cosh(a - (b*c)/d)*SinhIntegral((c*(b - d))/d + (b - d)*x))/(2*b*d) - (sinh(a - (b*c)/d)*SinhIntegral((c*(b - d))/d + (b - d)*x))/(2*b^2) + (x*cosh(a + b*x)*SinhIntegral(c + d*x))/b - (sinh(a + b*x)*SinhIntegral(c + d*x))/b^2 + (c*cosh(a - (b*c)/d)*SinhIntegral((c*(b + d))/d + (b + d)*x))/(2*b*d) + (sinh(a - (b*c)/d)*SinhIntegral((c*(b + d))/d + (b + d)*x))/(2*b^2)]
    @test_int [x^0*sinh(a + b*x)*SinhIntegral(c + d*x), x, 9, (CoshIntegral((c*(b - d))/d + (b - d)*x)*sinh(a - (b*c)/d))/(2*b) - (CoshIntegral((c*(b + d))/d + (b + d)*x)*sinh(a - (b*c)/d))/(2*b) + (cosh(a - (b*c)/d)*SinhIntegral((c*(b - d))/d + (b - d)*x))/(2*b) + (cosh(a + b*x)*SinhIntegral(c + d*x))/b - (cosh(a - (b*c)/d)*SinhIntegral((c*(b + d))/d + (b + d)*x))/(2*b)]
    @test_int [sinh(a + b*x)*SinhIntegral(c + d*x)/x^1, x, 0, CannotIntegrate((sinh(a + b*x)*SinhIntegral(c + d*x))/x, x)]

    #= [x^2*cosh(a + b*x)*SinhIntegral(c + d*x), x, 46, -(cosh(a - c + (b - d)*x)/(2*b*(b - d)^2)) - cosh(a - c + (b - d)*x)/(b^2*(b - d)) + cosh(a + c + (b + d)*x)/(2*b*(b + d)^2) + cosh(a + c + (b + d)*x)/(b^2*(b + d)) + (CoshIntegral(((b - d)*(c + d*x))/d)*((b^2*c^2 + 2*d^2)*cosh(a - (b*c)/d) + 2*b*c*d*sinh(a - (b*c)/d)))/(2*b^3*d^2) - (CoshIntegral(((b + d)*(c + d*x))/d)*((b^2*c^2 + 2*d^2)*cosh(a - (b*c)/d) + 2*b*c*d*sinh(a - (b*c)/d)))/(2*b^3*d^2) - (c*sinh(a - c + (b - d)*x))/(2*b*(b - d)*d) + (x*sinh(a - c + (b - d)*x))/(2*b*(b - d)) + (c*sinh(a + c + (b + d)*x))/(2*b*d*(b + d)) - (x*sinh(a + c + (b + d)*x))/(2*b*(b + d)) - ((2*b*x*cosh(a + b*x) - (2 + b^2*x^2)*sinh(a + b*x))*SinhIntegral(c + d*x))/b^3 + ((2*b*c*d*cosh(a - (b*c)/d) + (b^2*c^2 + 2*d^2)*sinh(a - (b*c)/d))*SinhIntegral(((b - d)*(c + d*x))/d))/(2*b^3*d^2) - ((2*b*c*d*cosh(a - (b*c)/d) + (b^2*c^2 + 2*d^2)*sinh(a - (b*c)/d))*SinhIntegral(((b + d)*(c + d*x))/d))/(2*b^3*d^2)] =#
    @test_int [x^1*cosh(a + b*x)*SinhIntegral(c + d*x), x, 24, -((c*cosh(a - (b*c)/d)*CoshIntegral((c*(b - d))/d + (b - d)*x))/(2*b*d)) + (c*cosh(a - (b*c)/d)*CoshIntegral((c*(b + d))/d + (b + d)*x))/(2*b*d) - (CoshIntegral((c*(b - d))/d + (b - d)*x)*sinh(a - (b*c)/d))/(2*b^2) + (CoshIntegral((c*(b + d))/d + (b + d)*x)*sinh(a - (b*c)/d))/(2*b^2) + sinh(a - c + (b - d)*x)/(2*b*(b - d)) - sinh(a + c + (b + d)*x)/(2*b*(b + d)) - (cosh(a - (b*c)/d)*SinhIntegral((c*(b - d))/d + (b - d)*x))/(2*b^2) - (c*sinh(a - (b*c)/d)*SinhIntegral((c*(b - d))/d + (b - d)*x))/(2*b*d) - (cosh(a + b*x)*SinhIntegral(c + d*x))/b^2 + (x*sinh(a + b*x)*SinhIntegral(c + d*x))/b + (cosh(a - (b*c)/d)*SinhIntegral((c*(b + d))/d + (b + d)*x))/(2*b^2) + (c*sinh(a - (b*c)/d)*SinhIntegral((c*(b + d))/d + (b + d)*x))/(2*b*d)]
    @test_int [x^0*cosh(a + b*x)*SinhIntegral(c + d*x), x, 9, (cosh(a - (b*c)/d)*CoshIntegral((c*(b - d))/d + (b - d)*x))/(2*b) - (cosh(a - (b*c)/d)*CoshIntegral((c*(b + d))/d + (b + d)*x))/(2*b) + (sinh(a - (b*c)/d)*SinhIntegral((c*(b - d))/d + (b - d)*x))/(2*b) + (sinh(a + b*x)*SinhIntegral(c + d*x))/b - (sinh(a - (b*c)/d)*SinhIntegral((c*(b + d))/d + (b + d)*x))/(2*b)]
    @test_int [cosh(a + b*x)*SinhIntegral(c + d*x)/x^1, x, 0, CannotIntegrate((cosh(a + b*x)*SinhIntegral(c + d*x))/x, x)]

    #= ::Section::Closed:: =#
    #=Hyperbolic*cosine*integral*function=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*CoshIntegral(b*x)^n=#

    @test_int [x^m*CoshIntegral(b*x), x, 5, (x^(1 + m)*CoshIntegral(b*x))/(1 + m) - (x^m*Gamma(1 + m, (-b)*x))/(((-b)*x)^m*(2*b*(1 + m))) + (x^m*Gamma(1 + m, b*x))/((b*x)^m*(2*b*(1 + m)))]

    @test_int [x^3*CoshIntegral(b*x), x, 6, (3*cosh(b*x))/(2*b^4) + (3*x^2*cosh(b*x))/(4*b^2) + (1/4)*x^4*CoshIntegral(b*x) - (3*x*sinh(b*x))/(2*b^3) - (x^3*sinh(b*x))/(4*b)]
    @test_int [x^2*CoshIntegral(b*x), x, 5, (2*x*cosh(b*x))/(3*b^2) + (1/3)*x^3*CoshIntegral(b*x) - (2*sinh(b*x))/(3*b^3) - (x^2*sinh(b*x))/(3*b)]
    @test_int [x^1*CoshIntegral(b*x), x, 4, cosh(b*x)/(2*b^2) + (1/2)*x^2*CoshIntegral(b*x) - (x*sinh(b*x))/(2*b)]
    @test_int [x^0*CoshIntegral(b*x), x, 1, x*CoshIntegral(b*x) - sinh(b*x)/b]
    @test_int [CoshIntegral(b*x)/x^1, x, 1, (-(1/2))*b*x*HypergeometricPFQ([1, 1, 1], [2, 2, 2], (-b)*x) + (1/2)*b*x*HypergeometricPFQ([1, 1, 1], [2, 2, 2], b*x) + EulerGamma*log(x) + (1/2)*log(b*x)^2]
    @test_int [CoshIntegral(b*x)/x^2, x, 4, -(cosh(b*x)/x) - CoshIntegral(b*x)/x + b*SinhIntegral(b*x)]
    @test_int [CoshIntegral(b*x)/x^3, x, 5, -(cosh(b*x)/(4*x^2)) + (1/4)*b^2*CoshIntegral(b*x) - CoshIntegral(b*x)/(2*x^2) - (b*sinh(b*x))/(4*x)]

    @test_int [x^m*CoshIntegral(b*x)^2, x, 0, CannotIntegrate(x^m*CoshIntegral(b*x)^2, x)]

    @test_int [x^3*CoshIntegral(b*x)^2, x, 19, -(x^2/(4*b^2)) + (3*cosh(b*x)^2)/(8*b^4) + (3*cosh(b*x)*CoshIntegral(b*x))/b^4 + (3*x^2*cosh(b*x)*CoshIntegral(b*x))/(2*b^2) + (1/4)*x^4*CoshIntegral(b*x)^2 - (3*CoshIntegral(2*b*x))/(2*b^4) - (3*log(x))/(2*b^4) - (x*cosh(b*x)*sinh(b*x))/b^3 - (3*x*CoshIntegral(b*x)*sinh(b*x))/b^3 - (x^3*CoshIntegral(b*x)*sinh(b*x))/(2*b) + (13*sinh(b*x)^2)/(8*b^4) + (x^2*sinh(b*x)^2)/(4*b^2)]
    @test_int [x^2*CoshIntegral(b*x)^2, x, 15, -(x/(2*b^2)) + (4*x*cosh(b*x)*CoshIntegral(b*x))/(3*b^2) + (1/3)*x^3*CoshIntegral(b*x)^2 - (5*cosh(b*x)*sinh(b*x))/(6*b^3) - (4*CoshIntegral(b*x)*sinh(b*x))/(3*b^3) - (2*x^2*CoshIntegral(b*x)*sinh(b*x))/(3*b) + (x*sinh(b*x)^2)/(3*b^2) + (2*SinhIntegral(2*b*x))/(3*b^3)]
    @test_int [x^1*CoshIntegral(b*x)^2, x, 10, (cosh(b*x)*CoshIntegral(b*x))/b^2 + (1/2)*x^2*CoshIntegral(b*x)^2 - CoshIntegral(2*b*x)/(2*b^2) - log(x)/(2*b^2) - (x*CoshIntegral(b*x)*sinh(b*x))/b + sinh(b*x)^2/(2*b^2)]
    @test_int [x^0*CoshIntegral(b*x)^2, x, 6, x*CoshIntegral(b*x)^2 - (2*CoshIntegral(b*x)*sinh(b*x))/b + SinhIntegral(2*b*x)/b]
    @test_int [CoshIntegral(b*x)^2/x^1, x, 0, CannotIntegrate(CoshIntegral(b*x)^2/x, x)]
    @test_int [CoshIntegral(b*x)^2/x^2, x, 0, CannotIntegrate(CoshIntegral(b*x)^2/x^2, x)]
    @test_int [CoshIntegral(b*x)^2/x^3, x, 0, CannotIntegrate(CoshIntegral(b*x)^2/x^3, x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*CoshIntegral(a+b*x)^n=#

    @test_int [x^m*CoshIntegral(a + b*x), x, 1, (x^(1 + m)*CoshIntegral(a + b*x))/(1 + m) - (b*CannotIntegrate((x^(1 + m)*cosh(a + b*x))/(a + b*x), x))/(1 + m)]

    @test_int [x^3*CoshIntegral(a + b*x), x, 14, (3*cosh(a + b*x))/(2*b^4) + (a^2*cosh(a + b*x))/(4*b^4) - (a*x*cosh(a + b*x))/(2*b^3) + (3*x^2*cosh(a + b*x))/(4*b^2) - (a^4*CoshIntegral(a + b*x))/(4*b^4) + (1/4)*x^4*CoshIntegral(a + b*x) + (a*sinh(a + b*x))/(2*b^4) + (a^3*sinh(a + b*x))/(4*b^4) - (3*x*sinh(a + b*x))/(2*b^3) - (a^2*x*sinh(a + b*x))/(4*b^3) + (a*x^2*sinh(a + b*x))/(4*b^2) - (x^3*sinh(a + b*x))/(4*b)]
    @test_int [x^2*CoshIntegral(a + b*x), x, 10, -((a*cosh(a + b*x))/(3*b^3)) + (2*x*cosh(a + b*x))/(3*b^2) + (a^3*CoshIntegral(a + b*x))/(3*b^3) + (1/3)*x^3*CoshIntegral(a + b*x) - (2*sinh(a + b*x))/(3*b^3) - (a^2*sinh(a + b*x))/(3*b^3) + (a*x*sinh(a + b*x))/(3*b^2) - (x^2*sinh(a + b*x))/(3*b)]
    @test_int [x^1*CoshIntegral(a + b*x), x, 7, cosh(a + b*x)/(2*b^2) - (a^2*CoshIntegral(a + b*x))/(2*b^2) + (1/2)*x^2*CoshIntegral(a + b*x) + (a*sinh(a + b*x))/(2*b^2) - (x*sinh(a + b*x))/(2*b)]
    @test_int [x^0*CoshIntegral(a + b*x), x, 1, ((a + b*x)*CoshIntegral(a + b*x))/b - sinh(a + b*x)/b]
    @test_int [CoshIntegral(a + b*x)/x^1, x, 0, CannotIntegrate(CoshIntegral(a + b*x)/x, x)]
    @test_int [CoshIntegral(a + b*x)/x^2, x, 7, (b*cosh(a)*CoshIntegral(b*x))/a - (b*CoshIntegral(a + b*x))/a - CoshIntegral(a + b*x)/x + (b*sinh(a)*SinhIntegral(b*x))/a]
    @test_int [CoshIntegral(a + b*x)/x^3, x, 11, -((b*cosh(a + b*x))/(2*a*x)) - (b^2*cosh(a)*CoshIntegral(b*x))/(2*a^2) + (b^2*CoshIntegral(a + b*x))/(2*a^2) - CoshIntegral(a + b*x)/(2*x^2) + (b^2*CoshIntegral(b*x)*sinh(a))/(2*a) + (b^2*cosh(a)*SinhIntegral(b*x))/(2*a) - (b^2*sinh(a)*SinhIntegral(b*x))/(2*a^2)]

    @test_int [x^m*CoshIntegral(a + b*x)^2, x, 0, CannotIntegrate(x^m*CoshIntegral(a + b*x)^2, x)]

    @test_int [x^2*CoshIntegral(a + b*x)^2, x, 39, -((2*x)/(3*b^2)) - (a*cosh(2*a + 2*b*x))/(3*b^3) + (x*cosh(2*a + 2*b*x))/(6*b^2) - (2*a*cosh(a + b*x)*CoshIntegral(a + b*x))/(3*b^3) + (4*x*cosh(a + b*x)*CoshIntegral(a + b*x))/(3*b^2) + (a^2*(a + b*x)*CoshIntegral(a + b*x)^2)/(3*b^3) - (a*x*(a + b*x)*CoshIntegral(a + b*x)^2)/(3*b^2) + (x^2*(a + b*x)*CoshIntegral(a + b*x)^2)/(3*b) + (a*CoshIntegral(2*a + 2*b*x))/b^3 + (a*log(a + b*x))/b^3 - (2*cosh(a + b*x)*sinh(a + b*x))/(3*b^3) - (4*CoshIntegral(a + b*x)*sinh(a + b*x))/(3*b^3) - (2*a^2*CoshIntegral(a + b*x)*sinh(a + b*x))/(3*b^3) + (2*a*x*CoshIntegral(a + b*x)*sinh(a + b*x))/(3*b^2) - (2*x^2*CoshIntegral(a + b*x)*sinh(a + b*x))/(3*b) - sinh(2*a + 2*b*x)/(12*b^3) + (2*SinhIntegral(2*a + 2*b*x))/(3*b^3) + (a^2*SinhIntegral(2*a + 2*b*x))/b^3]
    @test_int [x^1*CoshIntegral(a + b*x)^2, x, 17, cosh(2*a + 2*b*x)/(4*b^2) + (cosh(a + b*x)*CoshIntegral(a + b*x))/b^2 - (a*(a + b*x)*CoshIntegral(a + b*x)^2)/(2*b^2) + (x*(a + b*x)*CoshIntegral(a + b*x)^2)/(2*b) - CoshIntegral(2*a + 2*b*x)/(2*b^2) - log(a + b*x)/(2*b^2) + (a*CoshIntegral(a + b*x)*sinh(a + b*x))/b^2 - (x*CoshIntegral(a + b*x)*sinh(a + b*x))/b - (a*SinhIntegral(2*a + 2*b*x))/b^2]
    @test_int [x^0*CoshIntegral(a + b*x)^2, x, 5, ((a + b*x)*CoshIntegral(a + b*x)^2)/b - (2*CoshIntegral(a + b*x)*sinh(a + b*x))/b + SinhIntegral(2*a + 2*b*x)/b]
    @test_int [CoshIntegral(a + b*x)^2/x^1, x, 0, CannotIntegrate(CoshIntegral(a + b*x)^2/x, x)]
    @test_int [CoshIntegral(a + b*x)^2/x^2, x, 0, CannotIntegrate(CoshIntegral(a + b*x)^2/x^2, x)]
    @test_int [CoshIntegral(a + b*x)^2/x^3, x, 0, CannotIntegrate(CoshIntegral(a + b*x)^2/x^3, x)]

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

    @test_int [x^2*CoshIntegral(d*(a + b*log(c*x^n))), x, 7, (x^3*CoshIntegral(d*(a + b*log(c*x^n))))/3 - (x^3*ExpIntegralEi(((3 - b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(6*ℯ^((3*a)/(b*n))*(c*x^n)^(3/n)) - (x^3*ExpIntegralEi(((3 + b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(6*ℯ^((3*a)/(b*n))*(c*x^n)^(3/n))]
    @test_int [x^1*CoshIntegral(d*(a + b*log(c*x^n))), x, 7, (x^2*CoshIntegral(d*(a + b*log(c*x^n))))/2 - (x^2*ExpIntegralEi(((2 - b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(4*ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)) - (x^2*ExpIntegralEi(((2 + b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(4*ℯ^((2*a)/(b*n))*(c*x^n)^(2/n))]
    @test_int [x^0*CoshIntegral(d*(a + b*log(c*x^n))), x, 7, x*CoshIntegral(d*(a + b*log(c*x^n))) - (x*ExpIntegralEi(((1 - b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(2*ℯ^(a/(b*n))*(c*x^n)^n^(-1)) - (x*ExpIntegralEi(((1 + b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(2*ℯ^(a/(b*n))*(c*x^n)^n^(-1))]
    @test_int [CoshIntegral(d*(a + b*log(c*x^n)))/x^1, x, 3, (CoshIntegral(d*(a + b*log(c*x^n)))*(a + b*log(c*x^n)))/(b*n) - sinh(d*(a + b*log(c*x^n)))/(b*d*n)]
    @test_int [CoshIntegral(d*(a + b*log(c*x^n)))/x^2, x, 7, -(CoshIntegral(d*(a + b*log(c*x^n)))/x) + (ℯ^(a/(b*n))*(c*x^n)^n^(-1)*ExpIntegralEi(-(((1 - b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(2*x) + (ℯ^(a/(b*n))*(c*x^n)^n^(-1)*ExpIntegralEi(-(((1 + b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(2*x)]
    @test_int [CoshIntegral(d*(a + b*log(c*x^n)))/x^3, x, 7, -CoshIntegral(d*(a + b*log(c*x^n)))/(2*x^2) + (ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)*ExpIntegralEi(-(((2 - b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(4*x^2) + (ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)*ExpIntegralEi(-(((2 + b*d*n)*(a + b*log(c*x^n)))/(b*n))))/(4*x^2)]

    @test_int [(e*x)^m*CoshIntegral(d*(a + b*log(c*x^n))), x, 7, ((e*x)^(1 + m)*CoshIntegral(d*(a + b*log(c*x^n))))/(e*(1 + m)) - (x*(e*x)^m*ExpIntegralEi(((1 + m - b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(2*ℯ^((a*(1 + m))/(b*n))*(1 + m)*(c*x^n)^((1 + m)/n)) - (x*(e*x)^m*ExpIntegralEi(((1 + m + b*d*n)*(a + b*log(c*x^n)))/(b*n)))/(2*ℯ^((a*(1 + m))/(b*n))*(1 + m)*(c*x^n)^((1 + m)/n))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Hyper(b*x)*CoshIntegral(b*x)^n=#

    @test_int [cosh(b*x)*CoshIntegral(b*x)/x^3, x, 14, -(cosh(b*x)^2/(4*x^2)) - (cosh(b*x)*CoshIntegral(b*x))/(2*x^2) + (1/4)*b^2*CoshIntegral(b*x)^2 + b^2*CoshIntegral(2*b*x) - (b*cosh(b*x)*sinh(b*x))/(2*x) - (b*CoshIntegral(b*x)*sinh(b*x))/(2*x) - (b*sinh(2*b*x))/(4*x)]
    @test_int [cosh(b*x)*CoshIntegral(b*x)/x^2, x, 5, -(cosh(b*x)^2/x) - (cosh(b*x)*CoshIntegral(b*x))/x + b*CannotIntegrate((CoshIntegral(b*x)*sinh(b*x))/x, x) + b*SinhIntegral(2*b*x)]
    @test_int [cosh(b*x)*CoshIntegral(b*x)/x, x, 1, (1/2)*CoshIntegral(b*x)^2]
    @test_int [cosh(b*x)*CoshIntegral(b*x), x, 5, (CoshIntegral(b*x)*sinh(b*x))/b - SinhIntegral(2*b*x)/(2*b)]
    @test_int [x*cosh(b*x)*CoshIntegral(b*x), x, 9, -((cosh(b*x)*CoshIntegral(b*x))/b^2) + CoshIntegral(2*b*x)/(2*b^2) + log(x)/(2*b^2) + (x*CoshIntegral(b*x)*sinh(b*x))/b - sinh(b*x)^2/(2*b^2)]
    @test_int [x^2*cosh(b*x)*CoshIntegral(b*x), x, 14, (3*x)/(4*b^2) - (2*x*cosh(b*x)*CoshIntegral(b*x))/b^2 + (5*cosh(b*x)*sinh(b*x))/(4*b^3) + (2*CoshIntegral(b*x)*sinh(b*x))/b^3 + (x^2*CoshIntegral(b*x)*sinh(b*x))/b - (x*sinh(b*x)^2)/(2*b^2) - SinhIntegral(2*b*x)/b^3]
    @test_int [x^3*cosh(b*x)*CoshIntegral(b*x), x, 18, x^2/(2*b^2) - (3*cosh(b*x)^2)/(4*b^4) - (6*cosh(b*x)*CoshIntegral(b*x))/b^4 - (3*x^2*cosh(b*x)*CoshIntegral(b*x))/b^2 + (3*CoshIntegral(2*b*x))/b^4 + (3*log(x))/b^4 + (2*x*cosh(b*x)*sinh(b*x))/b^3 + (6*x*CoshIntegral(b*x)*sinh(b*x))/b^3 + (x^3*CoshIntegral(b*x)*sinh(b*x))/b - (13*sinh(b*x)^2)/(4*b^4) - (x^2*sinh(b*x)^2)/(2*b^2)]

    @test_int [sinh(b*x)*CoshIntegral(b*x)/x^3, x, 12, -((b*cosh(b*x)^2)/(2*x)) - (b*cosh(2*b*x))/(4*x) - (b*cosh(b*x)*CoshIntegral(b*x))/(2*x) + (1/2)*b^2*CannotIntegrate((CoshIntegral(b*x)*sinh(b*x))/x, x) - (CoshIntegral(b*x)*sinh(b*x))/(2*x^2) - sinh(2*b*x)/(8*x^2) + b^2*SinhIntegral(2*b*x)]
    @test_int [sinh(b*x)*CoshIntegral(b*x)/x^2, x, 7, (1/2)*b*CoshIntegral(b*x)^2 + b*CoshIntegral(2*b*x) - (CoshIntegral(b*x)*sinh(b*x))/x - sinh(2*b*x)/(2*x)]
    @test_int [sinh(b*x)*CoshIntegral(b*x)/x, x, 0, CannotIntegrate((CoshIntegral(b*x)*sinh(b*x))/x, x)]
    @test_int [sinh(b*x)*CoshIntegral(b*x), x, 5, (cosh(b*x)*CoshIntegral(b*x))/b - CoshIntegral(2*b*x)/(2*b) - log(x)/(2*b)]
    @test_int [x*sinh(b*x)*CoshIntegral(b*x), x, 9, -(x/(2*b)) + (x*cosh(b*x)*CoshIntegral(b*x))/b - (cosh(b*x)*sinh(b*x))/(2*b^2) - (CoshIntegral(b*x)*sinh(b*x))/b^2 + SinhIntegral(2*b*x)/(2*b^2)]
    @test_int [x^2*sinh(b*x)*CoshIntegral(b*x), x, 13, -(x^2/(4*b)) + cosh(b*x)^2/(4*b^3) + (2*cosh(b*x)*CoshIntegral(b*x))/b^3 + (x^2*cosh(b*x)*CoshIntegral(b*x))/b - CoshIntegral(2*b*x)/b^3 - log(x)/b^3 - (x*cosh(b*x)*sinh(b*x))/(2*b^2) - (2*x*CoshIntegral(b*x)*sinh(b*x))/b^2 + sinh(b*x)^2/b^3]
    @test_int [x^3*sinh(b*x)*CoshIntegral(b*x), x, 20, -((5*x)/(2*b^3)) - x^3/(6*b) + (x*cosh(b*x)^2)/(2*b^3) + (6*x*cosh(b*x)*CoshIntegral(b*x))/b^3 + (x^3*cosh(b*x)*CoshIntegral(b*x))/b - (4*cosh(b*x)*sinh(b*x))/b^4 - (x^2*cosh(b*x)*sinh(b*x))/(2*b^2) - (6*CoshIntegral(b*x)*sinh(b*x))/b^4 - (3*x^2*CoshIntegral(b*x)*sinh(b*x))/b^2 + (3*x*sinh(b*x)^2)/(2*b^3) + (3*SinhIntegral(2*b*x))/b^4]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Hyper(b*x)*CoshIntegral(d*x)^n=#

    @test_int [sinh(5*x)*CoshIntegral(2*x), x, 6, (1/5)*cosh(5*x)*CoshIntegral(2*x) - (1/10)*CoshIntegral(3*x) - (1/10)*CoshIntegral(7*x)]

    @test_int [cosh(5*x)*CoshIntegral(2*x), x, 6, (1/5)*CoshIntegral(2*x)*sinh(5*x) - (1/10)*SinhIntegral(3*x) - (1/10)*SinhIntegral(7*x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Hyper(a+b*x)*CoshIntegral(a+b*x)^n=#

    #= [x^3*sinh(a + b*x)*CoshIntegral(a + b*x), x, 32, -((5*x)/(2*b^3)) - (a^2*x)/(2*b^3) + (a*x^2)/(4*b^2) - x^3/(6*b) - (7*a*cosh(a + b*x)^2)/(4*b^4) + (x*cosh(a + b*x)^2)/(2*b^3) + (6*x*cosh(a + b*x)*CoshIntegral(a + b*x))/b^3 + (x^3*cosh(a + b*x)*CoshIntegral(a + b*x))/b + (3*a*CoshIntegral(2*a + 2*b*x))/b^4 + (a^3*CoshIntegral(2*a + 2*b*x))/(2*b^4) + (3*a*log(a + b*x))/b^4 + (a^3*log(a + b*x))/(2*b^4) - (4*cosh(a + b*x)*sinh(a + b*x))/b^4 - (a^2*cosh(a + b*x)*sinh(a + b*x))/(2*b^4) + (a*x*cosh(a + b*x)*sinh(a + b*x))/(2*b^3) - (x^2*cosh(a + b*x)*sinh(a + b*x))/(2*b^2) - (6*CoshIntegral(a + b*x)*sinh(a + b*x))/b^4 - (3*x^2*CoshIntegral(a + b*x)*sinh(a + b*x))/b^2 + (3*x*sinh(a + b*x)^2)/(2*b^3) + (3*SinhIntegral(2*a + 2*b*x))/b^4 + (3*a^2*SinhIntegral(2*a + 2*b*x))/(2*b^4)] =#
    @test_int [x^2*sinh(a + b*x)*CoshIntegral(a + b*x), x, 21, (a*x)/(2*b^2) - x^2/(4*b) + cosh(a + b*x)^2/(4*b^3) + cosh(2*a + 2*b*x)/(2*b^3) + (2*cosh(a + b*x)*CoshIntegral(a + b*x))/b^3 + (x^2*cosh(a + b*x)*CoshIntegral(a + b*x))/b - CoshIntegral(2*a + 2*b*x)/b^3 - (a^2*CoshIntegral(2*a + 2*b*x))/(2*b^3) - log(a + b*x)/b^3 - (a^2*log(a + b*x))/(2*b^3) + (a*cosh(a + b*x)*sinh(a + b*x))/(2*b^3) - (x*cosh(a + b*x)*sinh(a + b*x))/(2*b^2) - (2*x*CoshIntegral(a + b*x)*sinh(a + b*x))/b^2 - (a*SinhIntegral(2*a + 2*b*x))/b^3]
    @test_int [x^1*sinh(a + b*x)*CoshIntegral(a + b*x), x, 12, -(x/(2*b)) + (x*cosh(a + b*x)*CoshIntegral(a + b*x))/b + (a*CoshIntegral(2*a + 2*b*x))/(2*b^2) + (a*log(a + b*x))/(2*b^2) - (cosh(a + b*x)*sinh(a + b*x))/(2*b^2) - (CoshIntegral(a + b*x)*sinh(a + b*x))/b^2 + SinhIntegral(2*a + 2*b*x)/(2*b^2)]
    @test_int [x^0*sinh(a + b*x)*CoshIntegral(a + b*x), x, 4, (cosh(a + b*x)*CoshIntegral(a + b*x))/b - CoshIntegral(2*a + 2*b*x)/(2*b) - log(a + b*x)/(2*b)]
    @test_int [sinh(a + b*x)*CoshIntegral(a + b*x)/x^1, x, 0, CannotIntegrate((sinh(a + b*x)*CoshIntegral(a + b*x))/x, x)]

    #= [x^3*cosh(a + b*x)*CoshIntegral(a + b*x), x, 32, -((5*a*x)/(4*b^3)) + x^2/(2*b^2) - (15*cosh(a + b*x)^2)/(4*b^4) - (a^2*cosh(a + b*x)^2)/(2*b^4) - (6*cosh(a + b*x)*CoshIntegral(a + b*x))/b^4 - (3*x^2*cosh(a + b*x)*CoshIntegral(a + b*x))/b^2 + (3*CoshIntegral(2*a + 2*b*x))/b^4 + (3*a^2*CoshIntegral(2*a + 2*b*x))/(2*b^4) + (3*log(a + b*x))/b^4 + (3*a^2*log(a + b*x))/(2*b^4) - (7*a*cosh(a + b*x)*sinh(a + b*x))/(4*b^4) + (2*x*cosh(a + b*x)*sinh(a + b*x))/b^3 + (6*x*CoshIntegral(a + b*x)*sinh(a + b*x))/b^3 + (x^3*CoshIntegral(a + b*x)*sinh(a + b*x))/b - sinh(a + b*x)^2/(4*b^4) + (a*x*sinh(a + b*x)^2)/(2*b^3) - (x^2*sinh(a + b*x)^2)/(2*b^2) + (3*a*SinhIntegral(2*a + 2*b*x))/b^4 + (a^3*SinhIntegral(2*a + 2*b*x))/(2*b^4)] =#
    @test_int [x^2*cosh(a + b*x)*CoshIntegral(a + b*x), x, 21, x/b^2 + (a*cosh(2*a + 2*b*x))/(4*b^3) - (x*cosh(2*a + 2*b*x))/(4*b^2) - (2*x*cosh(a + b*x)*CoshIntegral(a + b*x))/b^2 - (a*CoshIntegral(2*a + 2*b*x))/b^3 - (a*log(a + b*x))/b^3 + (cosh(a + b*x)*sinh(a + b*x))/b^3 + (2*CoshIntegral(a + b*x)*sinh(a + b*x))/b^3 + (x^2*CoshIntegral(a + b*x)*sinh(a + b*x))/b + sinh(2*a + 2*b*x)/(8*b^3) - SinhIntegral(2*a + 2*b*x)/b^3 - (a^2*SinhIntegral(2*a + 2*b*x))/(2*b^3)]
    @test_int [x^1*cosh(a + b*x)*CoshIntegral(a + b*x), x, 11, -(cosh(2*a + 2*b*x)/(4*b^2)) - (cosh(a + b*x)*CoshIntegral(a + b*x))/b^2 + CoshIntegral(2*a + 2*b*x)/(2*b^2) + log(a + b*x)/(2*b^2) + (x*CoshIntegral(a + b*x)*sinh(a + b*x))/b + (a*SinhIntegral(2*a + 2*b*x))/(2*b^2)]
    @test_int [x^0*cosh(a + b*x)*CoshIntegral(a + b*x), x, 4, (CoshIntegral(a + b*x)*sinh(a + b*x))/b - SinhIntegral(2*a + 2*b*x)/(2*b)]
    @test_int [cosh(a + b*x)*CoshIntegral(a + b*x)/x^1, x, 0, CannotIntegrate((cosh(a + b*x)*CoshIntegral(a + b*x))/x, x)]

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

    #= [x^2*sinh(a + b*x)*CoshIntegral(c + d*x), x, 46, cosh(a - c + (b - d)*x)/(2*b*(b - d)^2) + cosh(a - c + (b - d)*x)/(b^2*(b - d)) + cosh(a + c + (b + d)*x)/(2*b*(b + d)^2) + cosh(a + c + (b + d)*x)/(b^2*(b + d)) - (CoshIntegral(((b - d)*(c + d*x))/d)*((b^2*c^2 + 2*d^2)*cosh(a - (b*c)/d) + 2*b*c*d*sinh(a - (b*c)/d)))/(2*b^3*d^2) - (CoshIntegral(((b + d)*(c + d*x))/d)*((b^2*c^2 + 2*d^2)*cosh(a - (b*c)/d) + 2*b*c*d*sinh(a - (b*c)/d)))/(2*b^3*d^2) + (CoshIntegral(c + d*x)*((2 + b^2*x^2)*cosh(a + b*x) - 2*b*x*sinh(a + b*x)))/b^3 + (c*sinh(a - c + (b - d)*x))/(2*b*(b - d)*d) - (x*sinh(a - c + (b - d)*x))/(2*b*(b - d)) + (c*sinh(a + c + (b + d)*x))/(2*b*d*(b + d)) - (x*sinh(a + c + (b + d)*x))/(2*b*(b + d)) - ((2*b*c*d*cosh(a - (b*c)/d) + (b^2*c^2 + 2*d^2)*sinh(a - (b*c)/d))*SinhIntegral(((b - d)*(c + d*x))/d))/(2*b^3*d^2) - ((2*b*c*d*cosh(a - (b*c)/d) + (b^2*c^2 + 2*d^2)*sinh(a - (b*c)/d))*SinhIntegral(((b + d)*(c + d*x))/d))/(2*b^3*d^2)] =#
    @test_int [x^1*sinh(a + b*x)*CoshIntegral(c + d*x), x, 24, (c*cosh(a - (b*c)/d)*CoshIntegral((c*(b - d))/d + (b - d)*x))/(2*b*d) + (x*cosh(a + b*x)*CoshIntegral(c + d*x))/b + (c*cosh(a - (b*c)/d)*CoshIntegral((c*(b + d))/d + (b + d)*x))/(2*b*d) + (CoshIntegral((c*(b - d))/d + (b - d)*x)*sinh(a - (b*c)/d))/(2*b^2) + (CoshIntegral((c*(b + d))/d + (b + d)*x)*sinh(a - (b*c)/d))/(2*b^2) - (CoshIntegral(c + d*x)*sinh(a + b*x))/b^2 - sinh(a - c + (b - d)*x)/(2*b*(b - d)) - sinh(a + c + (b + d)*x)/(2*b*(b + d)) + (cosh(a - (b*c)/d)*SinhIntegral((c*(b - d))/d + (b - d)*x))/(2*b^2) + (c*sinh(a - (b*c)/d)*SinhIntegral((c*(b - d))/d + (b - d)*x))/(2*b*d) + (cosh(a - (b*c)/d)*SinhIntegral((c*(b + d))/d + (b + d)*x))/(2*b^2) + (c*sinh(a - (b*c)/d)*SinhIntegral((c*(b + d))/d + (b + d)*x))/(2*b*d)]
    @test_int [x^0*sinh(a + b*x)*CoshIntegral(c + d*x), x, 9, -((cosh(a - (b*c)/d)*CoshIntegral((c*(b - d))/d + (b - d)*x))/(2*b)) + (cosh(a + b*x)*CoshIntegral(c + d*x))/b - (cosh(a - (b*c)/d)*CoshIntegral((c*(b + d))/d + (b + d)*x))/(2*b) - (sinh(a - (b*c)/d)*SinhIntegral((c*(b - d))/d + (b - d)*x))/(2*b) - (sinh(a - (b*c)/d)*SinhIntegral((c*(b + d))/d + (b + d)*x))/(2*b)]
    @test_int [sinh(a + b*x)*CoshIntegral(c + d*x)/x^1, x, 0, CannotIntegrate((CoshIntegral(c + d*x)*sinh(a + b*x))/x, x)]

    #= [x^2*cosh(a + b*x)*CoshIntegral(c + d*x), x, 46, (c*cosh(a - c + (b - d)*x))/(2*b*(b - d)*d) - (x*cosh(a - c + (b - d)*x))/(2*b*(b - d)) + (c*cosh(a + c + (b + d)*x))/(2*b*d*(b + d)) - (x*cosh(a + c + (b + d)*x))/(2*b*(b + d)) - (CoshIntegral(((b - d)*(c + d*x))/d)*(2*b*c*d*cosh(a - (b*c)/d) + (b^2*c^2 + 2*d^2)*sinh(a - (b*c)/d)))/(2*b^3*d^2) - (CoshIntegral(((b + d)*(c + d*x))/d)*(2*b*c*d*cosh(a - (b*c)/d) + (b^2*c^2 + 2*d^2)*sinh(a - (b*c)/d)))/(2*b^3*d^2) - (CoshIntegral(c + d*x)*(2*b*x*cosh(a + b*x) - (2 + b^2*x^2)*sinh(a + b*x)))/b^3 + sinh(a - c + (b - d)*x)/(2*b*(b - d)^2) + sinh(a - c + (b - d)*x)/(b^2*(b - d)) + sinh(a + c + (b + d)*x)/(2*b*(b + d)^2) + sinh(a + c + (b + d)*x)/(b^2*(b + d)) - (((b^2*c^2 + 2*d^2)*cosh(a - (b*c)/d) + 2*b*c*d*sinh(a - (b*c)/d))*SinhIntegral(((b - d)*(c + d*x))/d))/(2*b^3*d^2) - (((b^2*c^2 + 2*d^2)*cosh(a - (b*c)/d) + 2*b*c*d*sinh(a - (b*c)/d))*SinhIntegral(((b + d)*(c + d*x))/d))/(2*b^3*d^2)] =#
    @test_int [x^1*cosh(a + b*x)*CoshIntegral(c + d*x), x, 24, -(cosh(a - c + (b - d)*x)/(2*b*(b - d))) - cosh(a + c + (b + d)*x)/(2*b*(b + d)) + (cosh(a - (b*c)/d)*CoshIntegral((c*(b - d))/d + (b - d)*x))/(2*b^2) - (cosh(a + b*x)*CoshIntegral(c + d*x))/b^2 + (cosh(a - (b*c)/d)*CoshIntegral((c*(b + d))/d + (b + d)*x))/(2*b^2) + (c*CoshIntegral((c*(b - d))/d + (b - d)*x)*sinh(a - (b*c)/d))/(2*b*d) + (c*CoshIntegral((c*(b + d))/d + (b + d)*x)*sinh(a - (b*c)/d))/(2*b*d) + (x*CoshIntegral(c + d*x)*sinh(a + b*x))/b + (c*cosh(a - (b*c)/d)*SinhIntegral((c*(b - d))/d + (b - d)*x))/(2*b*d) + (sinh(a - (b*c)/d)*SinhIntegral((c*(b - d))/d + (b - d)*x))/(2*b^2) + (c*cosh(a - (b*c)/d)*SinhIntegral((c*(b + d))/d + (b + d)*x))/(2*b*d) + (sinh(a - (b*c)/d)*SinhIntegral((c*(b + d))/d + (b + d)*x))/(2*b^2)]
    @test_int [x^0*cosh(a + b*x)*CoshIntegral(c + d*x), x, 9, -((CoshIntegral((c*(b - d))/d + (b - d)*x)*sinh(a - (b*c)/d))/(2*b)) - (CoshIntegral((c*(b + d))/d + (b + d)*x)*sinh(a - (b*c)/d))/(2*b) + (CoshIntegral(c + d*x)*sinh(a + b*x))/b - (cosh(a - (b*c)/d)*SinhIntegral((c*(b - d))/d + (b - d)*x))/(2*b) - (cosh(a - (b*c)/d)*SinhIntegral((c*(b + d))/d + (b + d)*x))/(2*b)]
    @test_int [cosh(a + b*x)*CoshIntegral(c + d*x)/x^1, x, 0, CannotIntegrate((cosh(a + b*x)*CoshIntegral(c + d*x))/x, x)]
    end

max_stars_repo_path: test/8 Special functions/8.5 Hyperbolic integral functions.jl


text:
#ubound-comparison.jl

#runs comparisons on the ubound type

import Base: ==, <, >

if options[:devmode]
  macro matchsign_devmode(a, b)
    esc(quote
      if options[:devmode]
        (@signof a) == (@signof b) || throw(ArgumentError("sign of arguments is mismatched"))
      end
    end)
  end
else
  macro matchsign_devmode(a, b); nothing; end
end

@universal function cmp_inner_bound(a::Unum, b::Unum)
  @matchsign_devmode a b
  resolve_degenerates!(a)
  resolve_degenerates!(b)
  decode_exp(a) == decode_exp(b) || return false
  return a.fraction == b.fraction
end

#TODO:  Make this comparison not depend on making new Unums.
@universal function cmp_outer_bound(a::Unum, b::Unum)
  @matchsign_devmode a b
  resolve_degenerates!(a)
  resolve_degenerates!(b)
  decode_exp(a) == decode_exp(b) || return false
  return outer_exact(a) == outer_exact(b)
end

@universal function cmp_lower_bound(a::Unum, b::Unum)
  @matchsign_devmode a b
  is_positive(a) ? cmp_inner_bound(a, b) : cmp_outer_bound(a, b)
end
@universal function cmp_upper_bound(a::Unum, b::Unum)
  @matchsign_devmode a b
  is_positive(a) ? cmp_outer_bound(a, b) : cmp_inner_bound(a, b)
end
#==============================================================================#
#equality comparison
@universal function ==(a::Ubound, b::Ubound)
  #first, check to make sure that the left and right sides of the unum have the
  #same exact vs. inexact character.  This is a quick first-pass check to make
  low_exact = is_exact(a.lower)
  high_exact = is_exact(a.upper)
  (low_exact != is_exact(b.lower)) && return false
  (high_exact != is_exact(b.upper)) && return false
  #also check signs.
  (@signof(a.lower) != @signof(b.lower)) && (!is_zero(a.lower)) && (!is_zero(b.lower)) && return false
  (@signof(a.upper) != @signof(b.upper)) && (!is_zero(a.upper)) && (!is_zero(b.upper)) && return false

  #in the case they're exact then checking the end bounds is straightforward equality.
  if low_exact
    (a.lower != b.lower) && return false
  else
    glb(a.lower) == glb(b.lower) || return false
  end

  if high_exact
    (a.upper != b.upper) && return false
  else
    lub(a.upper) == lub(b.upper) || return false
  end

  return true
end

function =={ESS,FSS}(a::Ubound{ESS,FSS}, b::Unum{ESS,FSS})
  #resolve the ubound then check against the unum value.  For now, this returns
  #false.  A correct implementation will do a more detailed check.
  is_exact(b) && return false
  is_exact(a.lower) && return false
  is_exact(a.upper) && return false

  @signof(a.lower) == @signof(b) || return false
  @signof(a.upper) == @signof(b) || return false

  glb(a.lower) == glb(b) || return false
  lub(a.upper) == lub(b) || return false
end

#just flip the previous function to make things easier.
=={ESS,FSS}(a::Unum{ESS,FSS}, b::Ubound{ESS,FSS}) = (b == a)

#repeat the process, except for isequal.
function Base.isequal{ESS,FSS}(a::Ubound{ESS,FSS}, b::Unum{ESS,FSS})
  return false
end

Base.isequal{ESS,FSS}(a::Unum{ESS,FSS}, b::Ubound{ESS,FSS}) = isequal(b, a)

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

<{ESS,FSS}(a::Ubound{ESS,FSS}, b::Ubound{ESS,FSS}) = a.upper < b.lower
<{ESS,FSS}(a::Unum{ESS,FSS}, b::Ubound{ESS,FSS}) = a < b.lower
<{ESS,FSS}(a::Ubound{ESS,FSS}, b::Unum{ESS,FSS}) = a.upper < b

>{ESS,FSS}(a::Ubound{ESS,FSS}, b::Ubound{ESS,FSS}) = a.lower > b.upper
>{ESS,FSS}(a::Unum{ESS,FSS}, b::Ubound{ESS,FSS}) = a > b.upper
>{ESS,FSS}(a::Ubound{ESS,FSS}, b::Unum{ESS,FSS}) = a.lower > b

@universal ≸(a::Ubound, b::Unum) = notgtr(a.lower, b) && notless(a.upper, b)
@universal ≸(a::Unum, b::Ubound) = b ≸ a
@universal function ≸(a::Ubound, b::Ubound)
  (notgtr(b.lower, a.lower) && notgtr(a.lower, b.upper)) && return true
  (notgtr(a.lower, b.lower) && notgtr(b.lower, a.upper)) && return true
  (notgtr(a.lower, b.lower) && notgtr(b.upper, a.upper)) && return true
  (notgtr(b.lower, a.lower) && notgtr(a.upper, b.upper)) && return true
  return false
end

≹(a::Ubound, b::Unum) = ≸(a, b)
≹(a::Unum, b::Ubound) = ≸(a, b)
≹(a::Ubound, b::Ubound) = ≸(a, b)

ext: jl

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
lookup = Dict(
    :Normal => (
        pmcall=:(pm.Normal()),
        pmargs=[:mu,:sd])
)

function pymc3(expr)
    result = []
    if @capture(expr, v_ ~ dist_)
        if @capture(dist, dname_(args__))
            lhs = v
            dinfo = lookup[dname]
            rhs = dinfo.pmcall
            push!(rhs.args, "$v")
            for (jlarg,pmarg) in zip(args,dinfo.pmargs)
                push!(rhs.args, :($pmarg = $jlarg))
            end
        :($lhs = $rhs)
        end
    else expr
    end
end

ex = :(x ~ Normal(μ, σ))
pymc3(ex)


text:
function calculate_error(predictions::T₁,
                         targets::T₂) where {T₁ <: AbstractVector,
                                             T₂ <: AbstractVector}
    y = predictions
    t = targets
    N = length(y)

    sum(sign.(y) .!= sign.(t)) / N
end


text:
include("./execute.jl")

mutable struct BinaryCFile
	filename::String
end

@doc raw"""
Check if the filename exists, if yes read the file and
execute the code, else throw an exception
"""
function is_valid_file(file::BinaryCFile):Nothing
	if isfile(file.filename)
		file  = open(file.filename, "r")

		file_read_data = read(file, String)
		compiled_code = execute_binary_c(file_read_data)

		write_to_file(compiled_code)

		close(file)
	else
		error("$(file.filename) does not exist")
	end
	return Nothing
end

max_stars_repo_name: Technetium-Inc/BinaryC

max_issues_repo_name: Technetium-Inc/BinaryC

max_forks_repo_name: Technetium-Inc/BinaryC


text:
module NodeJS

using Pkg.Artifacts
import Base: SHA1

export nodejs_cmd, npm_cmd

function conditional_nodejs_load()
    artifact_info = artifact_meta("nodejs_app", joinpath(@__DIR__, "..", "Artifacts.toml"))

    artifact_info === nothing && return nothing

    return artifact_path(SHA1(artifact_info["git-tree-sha1"]))
end

const nodejs_path = conditional_nodejs_load()

const node_exe_name = Sys.iswindows() ? "node.exe" : "node"
const npm_exe_name = Sys.iswindows() ? "npm.cmd" : "npm"

const node_executable_path = nodejs_path === nothing ? node_exe_name : Sys.iswindows() ? joinpath(nodejs_path, node_exe_name) : joinpath(nodejs_path, "bin", node_exe_name)
const npm_executable_path = nodejs_path === nothing ? npm_exe_name : Sys.iswindows() ? joinpath(nodejs_path, npm_exe_name) : joinpath(nodejs_path, "bin", npm_exe_name)

"""
Return the full path of the node command.
"""
function nodejs_cmd()
    return `$node_executable_path`
end

"""
Return the full path of the npm command.
"""
function npm_cmd()
    if nodejs_path === nothing
        return `npm`
    elseif Sys.iswindows()
        return `$npm_executable_path`
    else
        return `$node_executable_path $npm_executable_path`
    end
end

end # module


text:
struct SecantRelaxation <: AbstractRelaxation end

function reformulate_unary_quadratic_term!(
    model::MOI.ModelLike,
    reformulation::Reformulation{SecantRelaxation},
    _x::VI,
    _y::VI,
    l::Float64,
    u::Float64,
)
    x, y = _map_to_interval_if_needed(model, _x, _y, l, u)
    # Secant cut: y <= x
    MOI.add_constraint(model, 1.0SV(x) - SV(y), MOI.GreaterThan(0.0))
    return nothing
end


text:
using Distributions
using LazySets: sample, dim

# TODO generalize to affine dynamics

function sample_points_backward(As::AbstractVector, Bs::AbstractVector,
                                x_last::AbstractVector, P_last::LazySet,
                                times::AbstractVector, ε::Number,
                                nsamples::Integer;
                                log::Bool=false)
    samples_sat = Vector{Vector{Float64}}()
    samples_unsat = Vector{Vector{Float64}}()
    box = box_approximation(P_last)
    if isflat(box)
        samples = [center(box)]
    else
        samples = sample(P_last, nsamples)
    end
    l = length(As)
    if log
        println("sampling $(length(samples)) points")
    end
    for y0 in samples
        samples = samples_sat
        x = x_last
        y = y0
        for k in l:-1:1
            A = -As[k]
            B = -Bs[k]
            T = times[k]
            T_prev = (k == 1) ? 0 : times[k-1]
            time = (T_prev, T)
            max_diff = check_sufficient(A, B, x, y, time)
            if max_diff > ε
                samples = samples_unsat
                break
            end
            T = time[2] - time[1]
            x = exp(A * T) * x
            y = exp(B * T) * y
        end
        push!(samples, y0)
    end
    if log
        println("$(length(samples_sat)) sat samples and " *
                "$(length(samples_unsat)) unsat samples")
    end
    return samples_sat, samples_unsat
end

function sample_points(As::AbstractVector, Bs::AbstractVector,
                       xs::AbstractVector, P_last::LazySet,
                       times::AbstractVector, ε::Number,
                       nsamples::Integer;
                       log::Bool=false)
    return sample_points_backward(As, Bs, xs[end], P_last, times, ε, nsamples;
                                  log=log)
end

function sample_points_init(As::AbstractVector, Bs::AbstractVector,
                            x0::AbstractVector, P_last::LazySet,
                            times::AbstractVector, ε::Number,
                            nsamples::Integer;
                            log::Bool=false)
    x_last = x0
    T = 0.0
    for k in 1:length(As)
        A = As[k]
        T = times[k] - T
        x_last = exp(A * T) * x_last
    end
    return sample_points_backward(As, Bs, x_last, P_last, times, ε, nsamples;
                                  log=log)
end

max_stars_repo_name: HySynth/HySynth

max_issues_repo_name: HySynth/HySynth

max_forks_repo_name: HySynth/HySynth


text:
###
# Maps
###

"""
A subtype of `Map` is used as a one-to-one map between two domains
via `view`. The domain of the map `m` is `axes(m,1)` and the range
is `union(m)`.

Maps must also overload `invmap` to give the inverse of the map, which
is equivalent to `invmap(m)[x] == findfirst(isequal(x), m)`.
"""

abstract type Map{T} <: AbstractQuasiVector{T} end

invmap(M::Map) = error("Overload invmap(::$(typeof(M)))")

Base.in(x, m::Map) = x in union(m)
Base.issubset(d::Map, b) = union(d) ⊆ b
Base.union(d::Map) = axes(invmap(d),1)

for find in (:findfirst, :findlast)
    @eval function $find(f::Base.Fix2{typeof(isequal)}, d::Map)
        f.x in d || return nothing
        $find(isequal(invmap(d)[f.x]), union(d))
    end
end

@eval function findall(f::Base.Fix2{typeof(isequal)}, d::Map)
    f.x in d || return eltype(axes(d,1))[]
    findall(isequal(invmap(d)[f.x]), union(d))
end

function Base.getindex(d::Map, x::Inclusion)
    x == axes(d,1) || throw(BoundsError(d, x))
    d
end

# Affine map represents A*x .+ b
abstract type AbstractAffineQuasiVector{T,AA,X,B} <: Map{T} end

summary(io::IO, a::AbstractAffineQuasiVector) = print(io, "$(a.A) * $(a.x) .+ ($(a.b))")

MemoryLayout(::Type{<:AbstractAffineQuasiVector}) = PolynomialLayout()

struct AffineQuasiVector{T,AA,X,B} <: AbstractAffineQuasiVector{T,AA,X,B}
    A::AA
    x::X
    b::B
end

AffineQuasiVector(A::AA, x::X, b::B) where {AA,X,B} =
    AffineQuasiVector{promote_type(eltype(AA), eltype(X), eltype(B)),AA,X,B}(A,x,b)

AffineQuasiVector(A, x) = AffineQuasiVector(A, x, zero(promote_type(eltype(A),eltype(x))))
AffineQuasiVector(x) = AffineQuasiVector(one(eltype(x)), x)

AffineQuasiVector(A, x::AffineQuasiVector, b) = AffineQuasiVector(A*x.A, x.x, A*x.b .+ b)

axes(A::AbstractAffineQuasiVector) = axes(A.x)

affine_getindex(A, k) = A.A*A.x[k] .+ A.b
Base.unsafe_getindex(A::AbstractAffineQuasiVector, k) = A.A*Base.unsafe_getindex(A.x,k) .+ A.b
getindex(A::AbstractAffineQuasiVector, k::Number) = affine_getindex(A, k)
function getindex(A::AbstractAffineQuasiVector, k::Inclusion)
    @boundscheck A.x[k] # throws bounds error if k ≠ x
    A
end

getindex(A::AbstractAffineQuasiVector, ::Colon) = copy(A)

copy(A::AbstractAffineQuasiVector) = A

inbounds_getindex(A::AbstractAffineQuasiVector{<:Any,<:Any,<:Inclusion}, k::Number) = A.A*k .+ A.b
isempty(A::AbstractAffineQuasiVector) = isempty(A.x)
==(a::AbstractAffineQuasiVector, b::AbstractAffineQuasiVector) = a.A == b.A && a.x == b.x && a.b == b.b
==(a::AbstractAffineQuasiVector, b::Inclusion) =  a == affine(b,b)
==(a::Inclusion, b::AbstractAffineQuasiVector) =  affine(a,a) == b

BroadcastStyle(::Type{<:AbstractAffineQuasiVector}) = LazyQuasiArrayStyle{1}()

for op in(:*, :\, :+, :-)
    @eval broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), a::Number, x::Inclusion) = broadcast($op, a, AffineQuasiVector(x))
end
for op in(:/, :+, :-)
    @eval broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), x::Inclusion, a::Number) = broadcast($op, AffineQuasiVector(x), a)
end

broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(*), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(a, x)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(\), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(inv(a), x)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(/), x::AbstractAffineQuasiVector, a::Number) = AffineQuasiVector(inv(a), x)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(+), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(one(eltype(x)), x, a)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(+), x::AbstractAffineQuasiVector, b::Number) = AffineQuasiVector(one(eltype(x)), x, b)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), a::Number, x::AbstractAffineQuasiVector) = AffineQuasiVector(-one(eltype(x)), x, a)
broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(-), x::AbstractAffineQuasiVector, b::Number) = AffineQuasiVector(one(eltype(x)), x, -b)

Base.@propagate_inbounds _affine_checkindex(inds, r) = isempty(r) | (checkindex(Bool, inds, Base.to_indices(inds, (first(r),))...) & checkindex(Bool, inds,  Base.to_indices(inds, (last(r),))...))
Base.@propagate_inbounds checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::Inclusion{<:Any,<:AbstractInterval}) = _affine_checkindex(inds, r)
Base.@propagate_inbounds checkindex(::Type{Bool}, inds::Inclusion{<:Any,<:AbstractInterval}, r::AbstractAffineQuasiVector) = _affine_checkindex(inds, r)

minimum(d::AbstractAffineQuasiVector) = signbit(d.A) ? last(d) : first(d)
maximum(d::AbstractAffineQuasiVector) = signbit(d.A) ? first(d) : last(d)

union(d::AbstractAffineQuasiVector) = Inclusion(minimum(d)..maximum(d))
invmap(d::AbstractAffineQuasiVector) = affine(union(d), axes(d,1))

struct AffineMap{T,D,R} <: AbstractAffineQuasiVector{T,T,D,T}
    domain::D
    range::R
end

AffineMap(domain::AbstractQuasiVector{T}, range::AbstractQuasiVector{V}) where {T,V} =
    AffineMap{promote_type(T,V), typeof(domain),typeof(range)}(domain,range)

measure(x::Inclusion{<:Any,<:AbstractInterval}) = last(x)-first(x)

function getproperty(A::AffineMap, d::Symbol)
    domain, range = getfield(A, :domain), getfield(A, :range)
    d == :x && return domain
    d == :A && return measure(range)/measure(domain)
    d == :b && return (last(domain)*first(range) - first(domain)*last(range))/measure(domain)
    getfield(A, d)
end

function getindex(A::AffineMap{T}, k::Number)  where T
    # ensure we exactly hit range
    k == first(A.domain) && return convert(T, first(A.range))::T
    k == last(A.domain) && return convert(T, last(A.range))::T
    convert(T, affine_getindex(A, k))::T
end

first(A::AffineMap{T}) where T = convert(T, first(A.range))::T
last(A::AffineMap{T}) where T = convert(T, last(A.range))::T

affine(a::AbstractQuasiVector, b::AbstractQuasiVector) = AffineMap(a, b)
affine(a, b::AbstractQuasiVector) = affine(Inclusion(a), b)
affine(a::AbstractQuasiVector, b) = affine(a, Inclusion(b))
affine(a, b) = affine(Inclusion(a), Inclusion(b))

# mapped vectors
const AffineMappedQuasiVector = SubQuasiArray{<:Any, 1, <:Any, <:Tuple{AbstractAffineQuasiVector}}
const AffineMappedQuasiMatrix = SubQuasiArray{<:Any, 2, <:Any, <:Tuple{AbstractAffineQuasiVector,Slice}}

==(a::AffineMappedQuasiVector, b::AffineMappedQuasiVector) = parentindices(a) == parentindices(b) && parent(a) == parent(b)

_sum(V::AffineMappedQuasiVector, ::Colon) = parentindices(V)[1].A \ sum(parent(V))

# pretty print for bases
summary(io::IO, P::AffineMappedQuasiMatrix) = print(io, "$(parent(P)) affine mapped to $(parentindices(P)[1].x.domain)")
summary(io::IO, P::AffineMappedQuasiVector) = print(io, "$(parent(P)) affine mapped to $(parentindices(P)[1].x.domain)")


text:
cue() = CressieRead(1)
hd()  = CressieRead(-1/2)

################################################################################
## Cressie Read Divergence
################################################################################

#=---------------
Evaluate
---------------=#
function evaluate(div::CressieRead, a::T, b::T) where T <: AbstractFloat
    α = div.α
    u = a/b
    if u > 0
        u = ( (u^(1+α)-1)/(α*(α+1)) - u/α + 1/α )*b
    elseif u == 0
        u = b/(1+α)
    else
        u = oftype(a, Inf)
    end
    u
end

function evaluate(div::CressieRead, a::AbstractVector{T}) where T <: AbstractFloat
    α = div.α
    r = zero(T)
    n = length(a)::Int64
    @inbounds for i = 1:n
        u = a[i]
        r += evaluate(div, u, one(T))
    end
    return r
end

function evaluate(div::CressieRead, a::AbstractVector{T}, b::AbstractVector{T}) where T <: AbstractFloat
    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    r = zero(T)
    @inbounds for i = eachindex(a, b)
        ai = a[i]
        bi = b[i]
        r += evaluate(div, ai, bi)
    end
    return r
end

#=---------------
gradient
---------------=#
function gradient(div::CressieRead, a::T, b::T) where T <: AbstractFloat
    α = div.α
    if a >= 0 && b > 0
        u = ((a/b)^α-1.0)/α
    elseif a == 0 && b == 0
        u = zero(T)
    else
        u = oftype(a, Inf)
    end
    return u
end

function gradient(div::CressieRead, a::T) where T <: AbstractFloat
    return gradient(div, a, one(T))
end

#=---------------
hessian
---------------=#
function hessian(div::CressieRead, a::T, b::T) where T <: AbstractFloat
    α    = div.α
    if a > 0 && b > 0
        u = (a/b)^α/a
    elseif a == 0 && b > 0
        if α >= 0
            u = zero(T)
        else
            u = oftype(a, Inf)
        end
    elseif a == 0 && b == 0
        u = zero(T)
    elseif a > 0 && b == 0
        u = oftype(a, Inf)
    end
    return u
end

function hessian(div::CressieRead, a::T) where T <: AbstractFloat
    return hessian(div, a, one(T))
end


text:
mutable struct Tensor{T, N} <: AbstractArray{T, N}
    type::Type
    ndims::Int64
    pointer::Ptr
    data::Union{Nothing,Array{T, N}}

    function Tensor{T, N}(p::Ptr, data) where {T, N}
        if !haskey(TYPE_MAP, T)
            error("Type $T is not supported.")
        end
        ret = new(T, N, p, data)
        finalizer(ret) do t
            ccall((:tensor_destroy, libtorch_capi),
                  Cvoid, (Ptr{Cvoid},),
                  t.pointer)
        end
        return ret
    end

end

Base.setproperty!(::Tensor, ::Symbol, x) = error("Can't change field of Tensor.")

function Tensor{T}(array::Array{U, N};
                   detach=false, requires_grad=false) where {T, U, N}
    if !haskey(TYPE_MAP, T)
        error("Type $T is not supported.")
    end

    dims = collect(size(array))
    stri = collect(strides(array))
    if T != U
        array = convert.(T, array)
    end
    grad = requires_grad ? 1 : 0
    copy_data = detach ? 1 : 0

    ptr = ccall((:tensor_from_data, libtorch_capi),
                Ptr{Cvoid},
                (Ptr{Cvoid}, Csize_t, Cchar,
                 Ptr{Clonglong}, Ptr{Clonglong}, Csize_t, Cint, Cint),
                array, sizeof(array), TYPE_MAP[T], dims, stri, N, copy_data, grad)
    Tensor{T, N}(ptr, detach ? nothing : array)
end

function Tensor(array::Array{T, N}; detach=false, requires_grad=false) where {T, N}
    Tensor{T}(array, detach=detach, requires_grad=requires_grad)
end

# 0-dim Tensor
function Tensor(s::Int64; requires_grad=false)
    grad = requires_grad ? 1 : 0
    ptr = ccall((:tensor_int64_0dim, libtorch_capi),
                Ptr{Cvoid},
                (Clonglong, Cint), s, grad)
    Tensor{Int64, 0}(ptr, nothing)
end

function Tensor(s::T; requires_grad=false) where {T <: TorchNumber}
    data = T[s]
    grad = requires_grad ? 1 : 0
    ptr = ccall((:tensor_from_data, libtorch_capi),
                Ptr{Cvoid},
                (Ptr{Cvoid}, Csize_t, Cchar,
                 Ptr{Clonglong}, Ptr{Clonglong}, Csize_t, Cint, Cint),
                data, sizeof(T), TYPE_MAP[T], C_NULL, C_NULL, 0, 1, grad)
    Tensor{T, 0}(ptr, nothing)
end

function Tensor(a0::Array{T, 0}; requires_grad=false) where {T <: TorchNumber}
    Tensor(a0[], requires_grad=requires_grad)
end

function tensor_from_ptr(p::Ptr)
    n_dims = ccall((:tensor_method_ndimension, libtorch_capi),
                   Clonglong, (Ptr{Cvoid},),
                   p)
    # sizes = zeros(Int64, n_dims)
    # ccall((:tensor_method_sizes, libtorch_capi),
    #       Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}),
    #       p, sizes)
    dtype = ccall((:tensor_method_dtype, libtorch_capi),
                  Cchar, (Ptr{Cvoid},),
                  p)
    Tensor{REVERSE_TYPE_MAP[dtype], n_dims}(p, nothing)
end

function Base.convert(::Type{Array}, t::Tensor{T, N}) where {T, N}
    if t.data != nothing
        return t.data
    end
    dims = size(t)
    ret = Array{T, N}(undef, reverse(dims))
    ccall((:tensor_method_data_copy, libtorch_capi),
          Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
          t.pointer, ret, sizeof(T) * prod(dims))
    if strides(t)[1] != 1
        return permutedims(ret, collect(N:-1:1))
    end
    return reshape(ret, dims)
end
Base.convert(::Type{T}, x::Tensor{T, 0}) where T = x[]

function Base.string(t::Tensor)
    str = ccall((:tensor_to_string, libtorch_capi),
                Ptr{UInt8}, (Ptr{Cvoid},),
                t.pointer)
    ret = unsafe_string(str)
    ccall(:free, Cvoid, (Ptr{Cvoid},), str)
    return ret
end

function Base.show(io::IO, t::Tensor{T, N}) where {T, N}
    write(io, "PyTorch.Tensor{$T, $N}:\n")
    write(io, string(t))
    write(io, "\n")
end

function Base.display(t::Tensor)
    show(stdout, t)
end

# array interface
Base.eltype(::Type{Tensor{T}}) where {T} = Tensor{T, 0}
Base.ndims(t::Tensor{T, N}) where {T, N} = N

eltype_id(::Tensor{T}) where {T} = Int(TYPE_MAP[T])
eltype_id(::Type{T}) where {T <: TorchNumber} = Int(TYPE_MAP[T])

function Base.strides(t::Tensor)
    n_dims = ccall((:tensor_method_ndimension, libtorch_capi),
                   Clonglong, (Ptr{Cvoid},),
                   t.pointer)
    strides = zeros(Int64, n_dims)
    ccall((:tensor_method_strides, libtorch_capi),
          Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}),
          t.pointer, strides)
    strides
end

function Base.size(t::Tensor{T, N}) where {T, N}
    n_dims = ccall((:tensor_method_ndimension, libtorch_capi),
                   Clonglong, (Ptr{Cvoid},),
                   t.pointer)
    @assert N == n_dims "Dimension mismatch!"
    sizes = zeros(Int64, n_dims)
    ccall((:tensor_method_sizes, libtorch_capi),
          Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}),
          t.pointer, sizes)
    tuple(sizes...)
end

function _tensor_indices(t::Tensor, I)
    indices = collect.(to_indices(t, I))
    shape = collect(Iterators.flatten(size.(indices)))
    indices = map(x -> x .- 1, indices)
    collect(indices), shape
end

_to_dim_0(t::Tensor) = ThC.opt_reshape(t, Int64[])
_to_dim_1_1(t::Tensor) = ThC.opt_reshape(t, [1, 1])

function Base.getindex(t::Tensor, I...)
    ts, shape = _tensor_indices(t, I)
    ret = t
    for i in 1:length(ts)
        ret = ThC.opt_index_select(ret, i - 1, ts[i])
    end
    all(x -> x == 1, size(ret)) && shape == Union{}[] && return _to_dim_0(ret)
    ThC.opt_reshape(ret, shape)
end
Base.getindex(t::Tensor{T}) where T = item(t)
Base.getindex(t::Tensor, i::Int64) = t[eachindex(t)[i]]
Base.getindex(t::Tensor{T, 1}, i::Int64) where T =
    ThC.opt_index_select(t, 0, (i - 1)) |> _to_dim_0
function Base.getindex(t::Tensor, I::UnitRange{Int64})
    t = vcat(map(i->_to_dim_1_1(t[i]), eachindex(t)[I])...)
    ThC.opt_reshape(t, [length(t)])
end

function Base.setindex!(t::Tensor{T}, v::Tensor{T}, I...) where T
    @assert length(I) > 0  "no indices given"
    @assert(!any(i -> i isa StepRange, I),
            "StepRange indices are not supported in Tensor assignment")
    ts, _1 = _tensor_indices(t, I)
    ret = t
    for i in 1:(length(ts) - 1)
        ret = ThC.narrow(ret, i - 1, ts[i][1], length(ts[i]))
    end
    dshape = length.(ts)
    ThC.index_copy!(ret, length(ts) - 1, Tensor(ts[end]), reshape(v, dshape))
    v
end
Base.setindex!(t::Tensor{T}, v::Array, I...) where T =
    setindex!(t, Tensor{T}(v), I...)
Base.setindex!(t::Tensor{T}, v::TorchNumber, i::Int64) where T =
    setindex!(t, Tensor{T}([v]), (eachindex(t)[i].I)...)
Base.setindex!(t::Tensor{T, 1}, v::TorchNumber, i::Int64) where T =
    setindex!(t, Tensor{T}([v]), i)

function Base.setindex!(t::Tensor{T}, v::Array, I::UnitRange{Int64}) where T
    indices = eachindex(t)[I]
    @assert length(v) == length(indices)
    for idx in 1:length(v)
        setindex!(t, Tensor{T}(v[[idx]]), (indices[idx].I)...)
    end
end

function Base.iterate(t::Tensor, state=(eachindex(t),))
    y = iterate(state...)
    y === nothing && return nothing
    t[y[1]], (state[1], Base.tail(y)...)
end

Base.cat(I::Vararg{Tensor}; dims) = cat(collect(I), dims)
Base.vcat(I::Vararg{Tensor}) = cat(collect(I), 0)
Base.hcat(I::Vararg{Tensor}) = cat(collect(I), 1)
function Base.hvcat(rows::Tuple{Vararg{Int}}, I::Vararg{Tensor,N}) where N
    ts = Iterators.Stateful(I)
    hs = map(n -> collect(Iterators.take(ts, n)), rows)
    hs = [hcat(t...) for t in hs]
    vcat(hs...)
end

# methods
function item(t::Tensor{T,N}) where {T,N}
    @assert(N == 0 || prod(size(t)) == 1,
            "The Tensor must contain only one element.")
    data = T[zero(T)]
    ccall((:tensor_method_item, libtorch_capi),
          Cvoid, (Ptr{Cvoid}, Cchar, Ptr{Cvoid}),
          t.pointer, TYPE_MAP[T], data)
    return data[1]
end

# devices

abstract type Device end
struct CPU <: Device end
struct CUDA <: Device
    index::Int
end

Base.convert(::Type{Int}, ::CPU) = -1
Base.convert(::Type{Int}, d::CUDA) = d.index

to(t::Tensor, d::Device) = ThC.to(t, convert(Int, d))
to(t::Tensor, ::Type{T}) where T <: TorchNumber = ThC.to2(t, eltype_id(T), 0, 0)
to(t::Tensor, ::Type{T}, d::Device) where T <: TorchNumber = to(t, d, T)
to(t::Tensor, d::Device, ::Type{T}) where T <: TorchNumber =
    ThC.to4(t, convert(Int, d), eltype_id(T), 0, 0)

function on(t::Tensor)
    data = Int64[0, 0]
    ccall((:tensor_method_device, libtorch_capi),
          Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), t.pointer, data)

    data[1] == -1 && return CPU()
    return CUDA(data[2])
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libnode"
version = v"14.17.3"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://nodejs.org/dist/v14.17.3/node-v14.17.3.tar.gz", "dcbd156506ee79ee48439257626ca0a6db3d7eab8cb0208db6979125ae7d3a95"),
]

host_platform = HostPlatform()

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd node-*

export CC_host=$HOSTCC
export CXX_host=$HOSTCXX
# Build & install libnode
if [[ $target == $MACHTYPE ]]
then
    ./configure --prefix=${prefix} --shared --without-inspector --no-cross-compiling
else
    DEST_CPU=x86_64
    if [[ $target == *aarch64* ]]; then DEST_CPU=arm64; fi
    DEST_OS=linux
    ./configure --prefix=${prefix} --shared --without-inspector --without-intl --cross-compiling --dest-cpu=$DEST_CPU --dest-os=$DEST_OS
fi

export LDFLAGS=-lrt
export CPPFLAGS=-D__STDC_FORMAT_MACROS
make -j`nproc`
make install
cp ./out/Release/node ${bindir}
install_license ./LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# Only 64-bit Linux platforms.
# Requires XCode and MSVC to compile on Mac OS and Windows, respectively.
# TODO: support FreeBSD/Mac/Windows
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),

    Platform("x86_64", "linux"; libc="musl"),
    Platform("aarch64", "linux"; libc="musl"),

    # Platform("x86_64", "macos"),
    # Platform("aarch64", "macos"),

    # Platform("x86_64", "windows"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libnode", :libnode),
    ExecutableProduct("node", :node),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

init_block = raw"""
"""

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(
    ARGS, name, version, sources, script,
    platforms, products, dependencies;
    julia_compat="1.6",
    preferred_gcc_version = v"7",
    init_block = init_block
)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Test

function _timed_include(file)
    println("Testing ", file)
    start = time()
    include(file)
    run_time = round(time() - start, digits = 1)
    return println("    Took $(run_time) seconds")
end

@testset "$(dir)" for dir in [".", "Variable", "Constraint", "Objective"]
    @testset "$(file)" for file in readdir(joinpath(@__DIR__, dir))
        if !endswith(file, ".jl")
            continue
        elseif (dir == ".") && (file == "Bridges.jl" || file == "utilities.jl")
            continue
        end
        _timed_include(joinpath(@__DIR__, dir, file))
    end
end


text:
# Import libraries.
using Turing, Flux, Plots, Random, ReverseDiff

# Hide sampling progress.
Turing.setprogress!(false);

# Use reverse_diff due to the number of parameters in neural networks.
Turing.setadbackend(:reversediff)

# Number of points to generate.
N = 80
M = round(Int, N / 4)
Random.seed!(1234)

# Generate artificial data.
x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
xt1s = Array([[x1s[i] + 0.5; x2s[i] + 0.5] for i = 1:M])
x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
append!(xt1s, Array([[x1s[i] - 5; x2s[i] - 5] for i = 1:M]))

x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
xt0s = Array([[x1s[i] + 0.5; x2s[i] - 5] for i = 1:M])
x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
append!(xt0s, Array([[x1s[i] - 5; x2s[i] + 0.5] for i = 1:M]))

# Store all the data for later.
xs = [xt1s; xt0s]
ts = [ones(2*M); zeros(2*M)]

# Plot data points.
function plot_data()
    x1 = map(e -> e[1], xt1s)
    y1 = map(e -> e[2], xt1s)
    x2 = map(e -> e[1], xt0s)
    y2 = map(e -> e[2], xt0s)

    Plots.scatter(x1,y1, color="red", clim = (0,1))
    Plots.scatter!(x2, y2, color="blue", clim = (0,1))
end

plot_data()

# Turn a vector into a set of weights and biases.
function unpack(nn_params::AbstractVector)
    W₁ = reshape(nn_params[1:6], 3, 2);
    b₁ = reshape(nn_params[7:9], 3)

    W₂ = reshape(nn_params[10:15], 2, 3);
    b₂ = reshape(nn_params[16:17], 2)

    Wₒ = reshape(nn_params[18:19], 1, 2);
    bₒ = reshape(nn_params[20:20], 1)
    return W₁, b₁, W₂, b₂, Wₒ, bₒ
end

# Construct a neural network using Flux and return a predicted value.
function nn_forward(xs, nn_params::AbstractVector)
    W₁, b₁, W₂, b₂, Wₒ, bₒ = unpack(nn_params)
    nn = Chain(Dense(W₁, b₁, tanh),
               Dense(W₂, b₂, tanh),
               Dense(Wₒ, bₒ, σ))
    return nn(xs)
end;

# Create a regularization term and a Gaussain prior variance term.
alpha = 0.09
sig = sqrt(1.0 / alpha)

# Specify the probabalistic model.
@model function bayes_nn(xs, ts)
    # Create the weight and bias vector.
    nn_params ~ MvNormal(zeros(20), sig .* ones(20))

    # Calculate predictions for the inputs given the weights
    # and biases in theta.
    preds = nn_forward(xs, nn_params)

    # Observe each prediction.
    for i = 1:length(ts)
        ts[i] ~ Bernoulli(preds[i])
    end
end;

# Perform inference.
N = 5000
ch = sample(bayes_nn(hcat(xs...), ts), HMC(0.05, 4), N);

# Extract all weight and bias parameters.
theta = MCMCChains.group(ch, :nn_params).value;

# Plot the data we have.
plot_data()

# Find the index that provided the highest log posterior in the chain.
_, i = findmax(ch[:lp])

# Extract the max row value from i.
i = i.I[1]

# Plot the posterior distribution with a contour plot.
x_range = collect(range(-6,stop=6,length=25))
y_range = collect(range(-6,stop=6,length=25))
Z = [nn_forward([x, y], theta[i, :])[1] for x=x_range, y=y_range]
contour!(x_range, y_range, Z)

# Return the average predicted value across
# multiple weights.
function nn_predict(x, theta, num)
    mean([nn_forward(x, theta[i,:])[1] for i in 1:10:num])
end;

# Plot the average prediction.
plot_data()

n_end = 1500
x_range = collect(range(-6,stop=6,length=25))
y_range = collect(range(-6,stop=6,length=25))
Z = [nn_predict([x, y], theta, n_end)[1] for x=x_range, y=y_range]
contour!(x_range, y_range, Z)

# Number of iterations to plot.
n_end = 500

anim = @gif for i=1:n_end
    plot_data()
    Z = [nn_forward([x, y], theta[i,:])[1] for x=x_range, y=y_range]
    contour!(x_range, y_range, Z, title="Iteration $i", clim = (0,1))
end every 5

using Bijectors
using Turing: Variational
using AdvancedVI

m = bayes_nn(hcat(xs...), ts);

q = Variational.meanfield(m)

μ = randn(length(q))
ω = -1 .* ones(length(q))

q = AdvancedVI.update(q, μ, exp.(ω));

advi = ADVI(10, 5_000)
q_hat = vi(m, advi, q);

samples = transpose(rand(q_hat, 5000))
ch_vi = Chains(reshape(samples, size(samples)..., 1), string.(MCMCChains.namesingroup(ch, :nn_params)));

# Extract all weight and bias parameters.
theta = MCMCChains.group(ch_vi, :nn_params).value;

# Plot the average prediction.
plot_data()

n_end = 1500
x_range = collect(range(-6,stop=6,length=25))
y_range = collect(range(-6,stop=6,length=25))
Z = [nn_predict([x, y], theta, n_end)[1] for x=x_range, y=y_range]
contour!(x_range, y_range, Z)

# Specify the network architecture.
network_shape = [
    (3,2, :tanh),
    (2,3, :tanh),
    (1,2, :σ)]

# Regularization, parameter variance, and total number of
# parameters.
alpha = 0.09
sig = sqrt(1.0 / alpha)
num_params = sum([i * o + i for (i, o, _) in network_shape])

# This modification of the unpack function generates a series of vectors
# given a network shape.
function unpack(θ::AbstractVector, network_shape::AbstractVector)
    index = 1
    weights = []
    biases = []
    for layer in network_shape
        rows, cols, _ = layer
        size = rows * cols
        last_index_w = size + index - 1
        last_index_b = last_index_w + rows
        push!(weights, reshape(θ[index:last_index_w], rows, cols))
        push!(biases, reshape(θ[last_index_w+1:last_index_b], rows))
        index = last_index_b + 1
    end
    return weights, biases
end

# Generate an abstract neural network given a shape,
# and return a prediction.
function nn_forward(x, θ::AbstractVector, network_shape::AbstractVector)
    weights, biases = unpack(θ, network_shape)
    layers = []
    for i in eachindex(network_shape)
        push!(layers, Dense(weights[i],
            biases[i],
            eval(network_shape[i][3])))
    end
    nn = Chain(layers...)
    return nn(x)
end

# General Turing specification for a BNN model.
@model bayes_nn_general(xs, ts, network_shape, num_params) = begin
    θ ~ MvNormal(zeros(num_params), sig .* ones(num_params))
    preds = nn_forward(xs, θ, network_shape)
    for i = 1:length(ts)
        ts[i] ~ Bernoulli(preds[i])
    end
end

# Perform inference.
num_samples = 500
ch2 = sample(bayes_nn_general(hcat(xs...), ts, network_shape, num_params), NUTS(0.65), num_samples);

# This function makes predictions based on network shape.
function nn_predict(x, theta, num, network_shape)
    mean([nn_forward(x, theta[i,:], network_shape)[1] for i in 1:10:num])
end;

# Extract the θ parameters from the sampled chain.
params2 = MCMCChains.group(ch2, :θ).value

plot_data()

x_range = collect(range(-6,stop=6,length=25))
y_range = collect(range(-6,stop=6,length=25))
Z = [nn_predict([x, y], params2, length(ch2), network_shape)[1] for x=x_range, y=y_range]
contour!(x_range, y_range, Z)

if isdefined(Main, :TuringTutorials)
    Main.TuringTutorials.tutorial_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
end

max_stars_repo_name: leachim/TuringTutorials

max_issues_repo_name: leachim/TuringTutorials

max_forks_repo_name: leachim/TuringTutorials


text:
@testset "interface" begin
    @testset "single-valued" begin
        x_seed = [0.0, 0.0]
        g_seed = [0.0, 0.0]
        h_seed = [0.0 0.0; 0.0 0.0]
        f_x_seed = 8157.682077608529
        g_x_seed = [-218.39260013257694, -48618.50356545231]
        h_x_seed = [982.7667005965963 0.0; 0.0 307917.1892478646]

        x_alt = [1.0, 1.0]
        f_x_alt = 57.316431861603284
        g_x_alt = [-5.43656365691809, -218.39260013257694]
        h_x_alt = [16.30969097075427 0.; 0. 982.7667005965963]

        # Construct instances
        nd = NonDifferentiable(exponential, x_seed)
        od = OnceDifferentiable(exponential, exponential_gradient!, exponential_value_gradient!, x_seed, 0.0, g_seed)
        td = TwiceDifferentiable(exponential, exponential_gradient!, exponential_value_gradient!, exponential_hessian!, x_seed, 0.0, g_seed, h_seed)

        # Force evaluation
        value!!(nd, x_seed)
        value_gradient!!(od, x_seed)
        value_gradient!!(td, x_seed)
        hessian!!(td, x_seed)

        # Test that values are the same, and that values match those
        # calculated by the value(obj, x) methods
        @test value(nd) == value(od) == value(td) == f_x_seed
        @test value(nd, x_seed) == value(od, x_seed) == value(td, x_seed)

        # Test that the gradients match the intended values
        @test gradient(od) == gradient(td) == g_x_seed
        # Test that the Hessian matches the intended value
        @test hessian(td) == h_x_seed

        # Test that the call counters got incremented
        @test nd.f_calls == od.f_calls == td.f_calls == [2]
        @test od.df_calls == td.df_calls == [1]
        @test td.h_calls == [1]

        # Test that the call counters do not get incremented
        # with single-"bang" methods...
        value!(nd, x_seed)
        value_gradient!(od, x_seed)
        value_gradient!(td, x_seed)
        hessian!(td, x_seed)

        @test nd.f_calls == od.f_calls == td.f_calls == [2]
        @test od.df_calls == td.df_calls == [1]
        @test td.h_calls == [1]

        # ... and that they do with double-"bang" methods
        value!!(nd, x_seed)
        value_gradient!!(od, x_seed)
        value_gradient!!(td, x_seed)
        hessian!!(td, x_seed)

        @test nd.f_calls == od.f_calls == td.f_calls == [3]
        @test od.df_calls == td.df_calls == [2]
        @test td.h_calls == [2]

        # Test that gradient doesn't work for NonDifferentiable, but does otherwise
        @test_throws ErrorException gradient!(nd, x_alt)
        gradient!(od, x_alt)
        gradient!(td, x_alt)

        @test value(nd) == value(od) == value(td) == f_x_seed
        @test gradient(td) == g_x_alt
        @test gradient(td) == [gradient(td, i) for i = 1:length(x_seed)]
        @test hessian(td) == h_x_seed
        @test nd.f_calls == od.f_calls == td.f_calls == [3]
        @test od.df_calls == td.df_calls == [3]
        @test td.h_calls == [2]

        @test_throws ErrorException hessian!(nd, x_alt)
        @test_throws ErrorException hessian!(od, x_alt)
        hessian!(td, x_alt)

        @test value(nd) == value(od) == value(td) == f_x_seed
        @test gradient(td) == g_x_alt
        @test hessian(td) == h_x_alt
        @test nd.f_calls == od.f_calls == td.f_calls == [3]
        @test od.df_calls == td.df_calls == [3]
        @test td.h_calls == [3]

        value!(nd, x_alt)
        value!(od, x_alt)
        value!(td, x_alt)
        @test value(nd) == value(od) == value(td) == f_x_alt
        @test gradient(td) == g_x_alt
        @test hessian(td) == h_x_alt
        @test nd.f_calls == od.f_calls == td.f_calls == [4]
        @test od.df_calls == td.df_calls == [3]
        @test td.h_calls == [3]

        @test_throws ErrorException value_gradient!(nd, x_seed)
        value_gradient!(od, x_seed)
        value_gradient!(td, x_seed)
        @test value(od) == value(td) == f_x_seed
        # change x_f manually to test branch
        od.x_f .*= 0
        td.x_f .*= 0
        value_gradient!(od, x_seed)
        value_gradient!(td, x_seed)
        @test value(od) == value(td) == f_x_seed
        # change x_df manually to test branch
        # only df is meant to be re-calculated as
        # d.x_f == x_seed
        # and that the hessian counter doesn't increment
        od.x_df .*= 0
        td.x_df .*= 0
        value_gradient!(od, x_seed)
        value_gradient!(td, x_seed)
        @test value(od) == value(td) == f_x_seed
        @test gradient(td) == g_x_seed
        @test hessian(td) == h_x_alt
        @test od.f_calls == td.f_calls == [5]
        @test od.df_calls == td.df_calls == [4]
        @test td.h_calls == [3]

        # test the non-mutating gradient() function
        gradient!(od, fill(1.0, 2))
        gradient!(td, fill(1.0, 2))
        od_df_old = copy(gradient(od))
        td_df_old = copy(gradient(td))
        gradient!(od, fill(-1.0, 2))
        gradient!(td, fill(-1.0, 2))
        @test od_df_old == gradient(od, fill(1.0, 2))
        @test td_df_old == gradient(td, fill(1.0, 2))
        # Mutate obj.DF directly to check that the gradient
        # is not recalculated if the same x is reused.
        # We obviously cannot just rerun for the same x twice
        # as we would then not be able to tell if it was just
        # calculated again, or the cache was simply returned
        # as intended.
        od.DF .= 0.0
        td.DF .= 0.0
        @test gradient(od) == zeros(2)
        @test gradient(td) == zeros(2)

        clear!(nd)
        clear!(od)
        clear!(td)
        @test nd.f_calls == [0,]
        @test od.f_calls == [0,]
        @test td.f_calls == [0,]
        @test od.df_calls == [0,]
        @test td.df_calls == [0,]
        @test td.h_calls == [0,]

        @testset "TwiceDifferentiableHV" begin
            for (name, prob) in MultivariateProblems.UnconstrainedProblems.examples
                if prob.istwicedifferentiable
                    hv!(storage::Vector, x::Vector, v::Vector) = begin
                        n = length(x)
                        H = Matrix{Float64}(undef, n, n)
                        MVP.hessian(prob)(H, x)
                        storage .= H * v
                    end
                    fg!(G::Vector, x::Vector) = begin
                        MVP.gradient(prob)(G, x)
                        MVP.objective(prob)(x), G
                    end
                    ddf = TwiceDifferentiableHV(MVP.objective(prob), fg!, hv!, prob.initial_x)
                    x = rand(prob.initial_x, length(prob.initial_x))
                    v = rand(prob.initial_x, length(prob.initial_x))
                    G = NLSolversBase.alloc_DF(x, 0.0)
                    H = NLSolversBase.alloc_H(x)
                    MVP.hessian(prob)(H, x)
                    @test hv_product!(ddf, x, v) == H*v
                    @test hv_product(ddf) == H*v
                    @test hv_product(ddf) == ddf.Hv
                    F, G = fg!(G, x)
                    @test gradient!(ddf, x) == G
                    @test value!(ddf, x) == F
                    @test f_calls(ddf) == 1
                    @test g_calls(ddf) == 1
                    @test h_calls(ddf) == 0
                    @test hv_calls(ddf) == 1
                    clear!(ddf)
                    @test f_calls(ddf) == 0
                    @test g_calls(ddf) == 0
                    @test h_calls(ddf) == 0
                    @test hv_calls(ddf) == 0
                end
            end
        end
    end
    @testset "multivalued" begin
        # Test example: Rosenbrock MINPACK
        function f!(F::Vector, x::Vector)
            F[1] = 1 - x[1]
            F[2] = 10(x[2]-x[1]^2)
            F
        end
        function j!(J::Matrix, x::Vector)
            J[1,1] = -1
            J[1,2] = 0
            J[2,1] = -20x[1]
            J[2,2] = 10
            J
        end

        x_seed = [0.0, 0.0]
        F_seed = [0.0, 0.0]
        J_seed = [0.0 0.0; 0.0 0.0]
        F_x_seed = [1.0, 0.0]
        J_x_seed = [-1.0 0.0; -0.0 10.0]

        x_alt = [0.5, 0.5]
        F_x_alt = [0.5, 2.5]
        J_x_alt = [-1.0 0.0; -10.0 10.0]

        # Construct instances
        nd = NonDifferentiable(f!, x_seed, F_seed)
        od = OnceDifferentiable(f!, j!, x_seed, F_seed, J_seed)

        # Force evaluation
        value!!(nd, nd.F, x_seed)
        value_jacobian!!(od, od.F, od.DF, x_seed)

        # Test that values are the same, and that values match those
        # calculated by the value(obj, x) methods
        @test value(nd) == value(od) == F_x_seed
        @test value(nd, x_seed) == value(od, x_seed)
        @test value(nd, x_alt) == value(od, x_alt)

        # Test that the Jacobians match the intended values
        @test jacobian(od) == J_x_seed

        # Test that the call counters got incremented
        @test nd.f_calls == od.f_calls == [3]
        @test od.df_calls == [1]

        # Test that the call counters do not get incremented
        # with single-"bang" methods...
        value!(nd, x_seed)
        value_jacobian!(od, x_seed)

        @test nd.f_calls == od.f_calls == [3]
        @test od.df_calls == [1]

        # ... and that they do with double-"bang" methods
        value!!(nd, x_seed)
        value_jacobian!!(od, x_seed)

        @test nd.f_calls == od.f_calls == [4]
        @test od.df_calls == [2]

        # Test that jacobian doesn't work for NonDifferentiable, but does otherwise
        @test_throws ErrorException jacobian!(nd, x_alt)
        jacobian!(od, x_alt)

        @test value(nd) == value(od) == F_x_seed
        @test jacobian(od) == J_x_alt
        @test nd.f_calls == od.f_calls == [4]
        @test od.df_calls == [3]

        @test value(nd) == value(od) == F_x_seed
        @test jacobian(od) == J_x_alt
        @test nd.f_calls == od.f_calls == [4]
        @test od.df_calls == [3]

        value!(nd, x_alt)
        value!(od, x_alt)
        @test value(nd) == value(od) == F_x_alt
        @test jacobian(od) == J_x_alt
        @test nd.f_calls == od.f_calls == [5]
        @test od.df_calls == [3]

        @test_throws ErrorException value_jacobian!(nd, x_seed)
        value_jacobian!(od, x_seed)
        @test value(od) == F_x_seed
        # change x_f manually to test branch
        od.x_f .*= 0
        value_jacobian!(od, x_seed)
        @test value(od) == F_x_seed
        # change x_df manually to test branch
        # only df is meant to be re-calculated as
        # d.x_f == x_seed
        od.x_df .*= 0
        value_jacobian!(od, x_seed)
        @test value(od) == F_x_seed
        @test jacobian(od) == J_x_seed
        @test od.f_calls == [6]
        @test od.df_calls == [4]

        clear!(nd)
        clear!(od)
        @test nd.f_calls == [0,]
        @test od.f_calls == [0,]
        @test od.df_calls == [0,]

        # Test that the correct branch gets called if jacobian hasn't been
        # calculated yet
        xx = rand(2)
        @test jacobian!(od, xx) == j!(NLSolversBase.alloc_DF(xx, xx), xx)
        @test od.x_df == xx

        # Test the branching in value_jacobian! works as expected
        xxx = rand(2)
        value_jacobian!(od, xxx)
        @test xxx == od.x_f == od.x_df
        xxx2 = rand(2)
        value!(od, xxx2)
        @test xxx2 == od.x_f
        @test xxx == od.x_df
        xxx3 = rand(2)
        jacobian!(od, xxx3)
        @test xxx2 == od.x_f
        @test xxx3 == od.x_df

    end
end


text:
# Code used to generate smallest circle around multiple ellipsoids problems. Save in Testproblems/SmallestCircle with an accurante
# solution from MOSEK

# Original problem has the following format:
# min_{t,xc,γ,τ1,...,τk)   t
# s.t.    [I -xc; xc' t+γ] ⪳ τi *  [Ai bi; bi' ci], i+1,...,k
#         [I, xc; xc' t+γ] ⪴ 0
#         τi ≥ 0, i+1,...,k

workspace()
include("../../OSSDP/Code/src/Solver.jl")
using JuMP, Mosek, JLD

rng = MersenneTwister(12345)
dirPath = "../DataFiles/Julia/SmallestCircle/"
!ispath(dirPath) && mkdir(dirPath)
nn = 25

 for iii =1:1:nn
  # set dimension of ellipsoids
  n = 2
  # choose number of elipsoids
  ne = rand(rng,3:15)
  # create random ellipsoid data
  Ae = []
  be = []
  ce = []
  Pe = []
  xce = []
  for kkk=1:ne
    # create mapping matrix P that distorts unit circle by scaling principal axis and rotating
    # define center of ellisoid xc
    Peigs = rand(rng,2)*(8-0.25)+0.25
    theta = rand(rng)*2*pi-pi
    P = diagm(Peigs)
    R = [cos(theta) -sin(theta); sin(theta) cos(theta)]
    P = R*P*R'
    xc = rand(rng,2,1)*(10.-(-10))-10

    # determine corresponding matrices for sublevel representation f(x)=x'Ax+2x'b+c
    Ainv = P\eye(2)
    A = Ainv*Ainv
    b = -A*xc
    c = (xc'A*xc-1)[1]
    push!(Ae,A)
    push!(be,b)
    push!(ce,c)
    push!(Pe,P)
    push!(xce,xc)
  end

  # put problem into solver format -> Pa,qa,Aa,ba,K
  # Define order of decision variables x=[t;xc;γ;τ1;...;τk]
  dimLMI = n+1

  # 1) constraint 1: τi ≥ 0, i+1,...,k
  A1 = [spzeros(ne,2+n) -speye(ne)]
  Aa = A1
  ba = zeros(ne,1)

  # 2) constraint 2: [I, xc; xc' t+γ] ⪴ 0
  # NOTE: HERE I ASSUME DIM(xc) = 2
  Mt = zeros(dimLMI,dimLMI)
  Mt[dimLMI,dimLMI] = 1
  Mγ = Mt
  Mxc1 = zeros(dimLMI,dimLMI)
  Mxc1[1,3] = 1
  Mxc1[3,1] = 1
  Mxc2 = zeros(dimLMI,dimLMI)
  Mxc2[2,3] = 1
  Mxc2[3,2] = 1
  M0 = zeros(dimLMI,dimLMI)
  M0[1:n,1:n] = eye(n)
  A2 = [vec(Mt) vec(Mxc1) vec(Mxc2) vec(Mγ) zeros(dimLMI^2,ne)]

  # append LMI constraint to augmented Data Matrix Aa and ba
  Aa = [Aa;-A2]
  ba = [ba;vec(M0)]

  # 3) constraint 3: [I -xc; xc' t+γ]  ⪳ τi *  [Ai bi; bi' ci], i+1,...,k
  Mt = zeros(dimLMI,dimLMI)
  Mγ = zeros(dimLMI,dimLMI)
  Mγ[dimLMI,dimLMI] = 1
  Mxc1 = zeros(dimLMI,dimLMI)
  Mxc1[1,3] = -1
  Mxc1[3,1] = -1
  Mxc2 = zeros(dimLMI,dimLMI)
  Mxc2[2,3] = -1
  Mxc2[3,2] = -1
  M0 = zeros(dimLMI,dimLMI)
  M0[1:n,1:n] = eye(n)
  for kkk=1:ne
    Mτ = zeros(dimLMI,dimLMI)
    Mτ[1:n,1:n] = -Ae[kkk]
    Mτ[1:n,n+1] = -be[kkk]
    Mτ[n+1,1:n] = -be[kkk]'
    Mτ[dimLMI,dimLMI] = -ce[kkk]

    # append LMI constraint to Augmented Data Matrix Aa and ba
    A3 = [vec(Mt) vec(Mxc1) vec(Mxc2) vec(Mγ) zeros(dimLMI^2,kkk-1) vec(Mτ) zeros(dimLMI^2,ne-kkk)]
    Aa = [Aa;A3]
    ba = [ba;-vec(M0)]
  end
  Pa = spzeros(size(Aa,2),size(Aa,2))
  qa = [1;zeros(size(Aa,2)-1)]
  ra = 0.

  # define cone dimensions
  Kf = 0
  Kl = ne
  Kq = []
  Ks = dimLMI^2*ones(ne+1)

  # solve accurately once with mosek
  model = Model(solver=MosekSolver())
  @variable(model, xc[1:2])
  @variable(model, t)
  @variable(model, g)
  @variable(model, tau[1:ne] >= 0)
  @objective(model, Min, t)
  @SDconstraint(model,[eye(2) xc;xc' t+g] >= 0)
  for jjj=1:ne
    @SDconstraint(model,[eye(2) -xc;-xc' g] <= tau[jjj]*[Ae[jjj] be[jjj]; be[jjj]' ce[jjj]])
  end
  status = JuMP.solve(model)

  @test status == :Optimal
  objTrue = getobjectivevalue(model)
  xc = getvalue(xc)
  g = getvalue(g)
  radius = sqrt(xc'*xc-g)

  # # try mosek with qocs solver format
  # model2 = Model(solver=MosekSolver())
  # @variable(model2, S0[1:ne] >= 0)
  # @variable(model2, S1[1:dimLMI,1:dimLMI],SDP)
  # @variable(model2, S2[1:dimLMI,1:dimLMI,1:ne])
  # @variable(model2, xc2[1:2])
  # @variable(model2, t2)
  # @variable(model2, g2)
  # @variable(model2, tau2[1:ne] >= 0)

  # x = [t2;xc2;g2;tau2]ss
  # s = [S0;vec(S1)]
  # for jjj=1:ne
  #   s = [s;vec(S2[:,:,jjj])]
  #   @SDconstraint(model2,S2[:,:,jjj] >= 0)
  # end
  # @objective(model2, Min, qa'*x)
  # @constraint(model2, Aa*x+s.==ba)
  # status2 = JuMP.solve(model2)

  # objTrue2 = getobjectivevalue(model2)
  # xc2 = getvalue(xc2)
  # g2 = getvalue(g2)
  # radius2 = sqrt(xc2'*xc2-g2)
  nr = "$(iii)"
  if iii < 10
    nr = "0$(iii)"
  end
  problemType = "Smallest Circle around ellipsoids"
  problemName = "SmallestCircle"*nr
  fn = "SmallestCircle"*nr*".jld"
  JLD.save(dirPath*fn,"n",size(Aa,2),"m",size(Aa,1),"A",Aa,"b",ba,"P",Pa,"q",qa,"r",ra,"objTrue",objTrue,"radius",radius,"xc",xc,"problemType",problemType,"problemName",problemName,"Kf",Kf,"Kl",Kl,"Kq",Kq,"Ks",Ks)
  println("$(iii)/$(nn) completed!")
end

end


text:
# ## Data

using DataFrames
using Turing
using TuringModels

import CSV

data_path = joinpath(TuringModels.project_root, "data", "rugged.csv")
df = CSV.read(data_path, DataFrame)

df.log_gdp = log.(df.rgdppc_2000)
dropmissing!(df)

df = select(df, :log_gdp, :rugged, :cont_africa);

df.log_gdp_std = df.log_gdp ./ mean(df.log_gdp)
df.rugged_std = df.rugged ./ maximum(df.rugged)

first(df, 8)

# ## Model

@model function model_fn(log_gdp_std, rugged_std, mean_rugged)
    α ~ Normal(1, 0.1)
    β ~ Normal(0, 0.3)
    σ ~ Exponential(1)

    μ = α .+ β * (rugged_std .- mean_rugged)
    log_gdp_std ~ MvNormal(μ, σ)
end

model = model_fn(df.log_gdp, df.rugged_std, mean(df.rugged_std));

# ## Output

chns = sample(model, NUTS(), MCMCThreads(), 1000, 3)

# \defaultoutput{}

#


text:
julia> sum([1,2,3,4,5].^2)
55

julia> sum([x^2 for x in [1,2,3,4,5]])
55

julia> mapreduce(x->x^2,+,[1:5])
55

julia> sum([x^2 for x in []])
0

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode