text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((a+b*x)*(c+d*x)*(e+f*x)*(g+h*x), x) == :(x^2*((1/2)*a*(c*e*h+c*f*g+d*e*g)+(1/2)*b*c*e*g)+x^3*(1//3*a*(c*f*h+d*e*h+d*f*g)+1//3*b*(c*e*h+c*f*g+d*e*g))+x^4*(1//4*b*(c*f*h+d*e*h+d*f*g)+1//4*a*d*f*h)+a*c*e*g*x+1//5*b*d*f*h*x^5)
@test integrate((g+h*x)^-1*(a+b*x)*(c+d*x)*(e+f*x), x) == :(x*h^-3*(b*(d*g+-1*c*h)*(f*g+-1*e*h)+-1*a*h*(d*f*g+-1*c*f*h+-1*d*e*h))+(1/2)*h^-2*x^2*(-1*b*(d*f*g+-1*c*f*h+-1*d*e*h)+a*d*f*h)+-1*h^-4*(b*g+-1*a*h)*(d*g+-1*c*h)*(f*g+-1*e*h)*log(g+h*x)+1//3*b*d*f*h^-1*x^3)
@test integrate((e+f*x)^-1*(g+h*x)^-1*(a+b*x)*(c+d*x), x) == :(b*d*x*f^-1*h^-1+f^-2*(f*g+-1*e*h)^-1*(b*e+-1*a*f)*(d*e+-1*c*f)*log(e+f*x)+-1*h^-2*(f*g+-1*e*h)^-1*(b*g+-1*a*h)*(d*g+-1*c*h)*log(g+h*x))
@test integrate((c+d*x)^-1*(e+f*x)^-1*(g+h*x)^-1*(a+b*x), x) == :((d*e+-1*c*f)^-1*(f*g+-1*e*h)^-1*(b*e+-1*a*f)*log(e+f*x)+-1*(d*e+-1*c*f)^-1*(d*g+-1*c*h)^-1*(b*c+-1*a*d)*log(c+d*x)+-1*(d*g+-1*c*h)^-1*(f*g+-1*e*h)^-1*(b*g+-1*a*h)*log(g+h*x))
@test integrate((a+b*x)^-1*(c+d*x)^-1*(e+f*x)^-1*(g+h*x)^-1, x) == :(b^2*(b*c+-1*a*d)^-1*(b*e+-1*a*f)^-1*(b*g+-1*a*h)^-1*log(a+b*x)+f^2*(b*e+-1*a*f)^-1*(d*e+-1*c*f)^-1*(f*g+-1*e*h)^-1*log(e+f*x)+-1*d^2*(b*c+-1*a*d)^-1*(d*e+-1*c*f)^-1*(d*g+-1*c*h)^-1*log(c+d*x)+-1*h^2*(b*g+-1*a*h)^-1*(d*g+-1*c*h)^-1*(f*g+-1*e*h)^-1*log(g+h*x))
@test integrate(x*(1+x)^-1*(2+x)^-1*(3+x)^-1, x) == :(2*log(2+x)+-3//2*log(3+x)+-1//2*log(1+x))
@test integrate((-6+x)^-1*(3+5x)^-3*(x^3+-1*x^2), x) == :(201*(45375+75625x)^-1+-12//1375*(3+5x)^-2+20//3993*log(6+-1x)+1493//499125*log(3+5x))
@test integrate(x^-1*(a+b*x)^3*(c+d*x)^(1/2)*(e+f*x), x) == :(2*e*a^3*(c+d*x)^(1/2)+2//315*d^-4*(c+d*x)^3//2*(8*b^3*c^2*(-2*c*f+3*d*e)+40*f*a^3*d^3+3*b*d*x*(-1*(-4*a*d+4*b*c)*(-2*b*c*f+2*a*d*f+3*b*d*e)+21*a*b*e*d^2)+6*b*a^2*d^2*(-16*c*f+45*d*e)+-18*a*c*d*b^2*(-4*c*f+7*d*e))+-2*e*a^3*c^(1/2)*arctanh(c^-1//2*(c+d*x)^(1/2))+1//21*d^-2*(a+b*x)^2*(c+d*x)^3//2*(-4*b*c*f+4*a*d*f+6*b*d*e)+2//9*f*d^-1*(a+b*x)^3*(c+d*x)^3//2)
@test integrate(x^-1*(a+b*x)^2*(c+d*x)^(1/2)*(e+f*x), x) == :(2*e*a^2*(c+d*x)^(1/2)+2//105*d^-3*(c+d*x)^3//2*(-2*c*b^2*(-4*c*f+7*d*e)+20*f*a^2*d^2+3*b*d*x*(-4*b*c*f+4*a*d*f+7*b*d*e)+14*a*b*d*(-2*c*f+5*d*e))+-2*e*a^2*c^(1/2)*arctanh(c^-1//2*(c+d*x)^(1/2))+2//7*f*d^-1*(a+b*x)^2*(c+d*x)^3//2)
@test integrate(x^-1*(c+d*x)^(1/2)*(a+b*x)*(e+f*x), x) == :(2*a*e*(c+d*x)^(1/2)+-2//15*d^-2*(c+d*x)^3//2*(-5*d*(a*f+b*e)+2*b*c*f+-3*b*d*f*x)+-2*a*e*c^(1/2)*arctanh(c^-1//2*(c+d*x)^(1/2)))
@test integrate(x^-1*(c+d*x)^(1/2)*(e+f*x), x) == :(2*e*(c+d*x)^(1/2)+-2*e*c^(1/2)*arctanh(c^-1//2*(c+d*x)^(1/2))+2//3*f*d^-1*(c+d*x)^3//2)
@test integrate(x^-1*(a+b*x)^-1*(c+d*x)^(1/2)*(e+f*x), x) == :(2*f*b^-1*(c+d*x)^(1/2)+-2*e*a^-1*c^(1/2)*arctanh(c^-1//2*(c+d*x)^(1/2))+2*a^-1*b^-3//2*(b*c+-1*a*d)^(1/2)*(b*e+-1*a*f)*arctanh(b^(1/2)*(c+d*x)^(1/2)*(b*c+-1*a*d)^-1//2))
@test integrate(x^-1*(a+b*x)^-2*(c+d*x)^(1/2)*(e+f*x), x) == :(-2*e*a^-2*c^(1/2)*arctanh(c^-1//2*(c+d*x)^(1/2))+a^-1*b^-1*(a+b*x)^-1*(c+d*x)^(1/2)*(b*e+-1*a*f)+a^-2*b^-3//2*(b*c+-1*a*d)^-1//2*(-1*a*d*(a*f+b*e)+2*c*e*b^2)*arctanh(b^(1/2)*(c+d*x)^(1/2)*(b*c+-1*a*d)^-1//2))
@test integrate(x^-1*(a+b*x)^-3*(c+d*x)^(1/2)*(e+f*x), x) == :(-2*e*a^-3*c^(1/2)*arctanh(c^-1//2*(c+d*x)^(1/2))+(1/2)*a^-1*b^-1*(a+b*x)^-2*(c+d*x)^(1/2)*(b*e+-1*a*f)+1//4*a^-3*b^-3//2*(b*c+-1*a*d)^-3//2*(f*a^3*d^2+8*e*b^3*c^2+3*b*e*a^2*d^2+-12*a*c*d*e*b^2)*arctanh(b^(1/2)*(c+d*x)^(1/2)*(b*c+-1*a*d)^-1//2)+1//4*a^-2*b^-1*(a+b*x)^-1*(c+d*x)^(1/2)*(b*c+-1*a*d)^-1*(-1*d*f*a^2+4*c*e*b^2+-3*a*b*d*e))
@test integrate(x^-1*(a+b*x)^(1/2)*(c+d*x)^3*(e+f*x), x) == :(2*e*c^3*(a+b*x)^(1/2)+-2//315*b^-4*(a+b*x)^3//2*(-10*b^3*c^2*(4*c*f+27*d*e)+16*f*a^3*d^3+-24*b*a^2*d^2*(d*e+3*c*f)+-3*b*d*x*((-4*a*d+4*b*c)*(-2*a*d*f+2*b*c*f+3*b*d*e)+21*c*d*e*b^2)+6*a*c*d*b^2*(16*c*f+21*d*e))+-2*e*a^(1/2)*c^3*arctanh(a^-1//2*(a+b*x)^(1/2))+1//21*b^-2*(a+b*x)^3//2*(c+d*x)^2*(-4*a*d*f+4*b*c*f+6*b*d*e)+2//9*f*b^-1*(a+b*x)^3//2*(c+d*x)^3)
@test integrate(x^-1*(a+b*x)^(1/2)*(c+d*x)^2*(e+f*x), x) == :(2*e*c^2*(a+b*x)^(1/2)+2//105*b^-3*(a+b*x)^3//2*(8*f*a^2*d^2+10*c*b^2*(2*c*f+7*d*e)+-14*a*b*d*(d*e+2*c*f)+3*b*d*x*(-4*a*d*f+4*b*c*f+7*b*d*e))+-2*e*a^(1/2)*c^2*arctanh(a^-1//2*(a+b*x)^(1/2))+2//7*f*b^-1*(a+b*x)^3//2*(c+d*x)^2)
@test integrate(x^-1*(a+b*x)^(1/2)*(c+d*x)*(e+f*x), x) == :(2*c*e*(a+b*x)^(1/2)+-2//15*b^-2*(a+b*x)^3//2*(-5*b*(c*f+d*e)+2*a*d*f+-3*b*d*f*x)+-2*c*e*a^(1/2)*arctanh(a^-1//2*(a+b*x)^(1/2)))
@test integrate(x^-1*(a+b*x)^(1/2)*(e+f*x), x) == :(2*e*(a+b*x)^(1/2)+-2*e*a^(1/2)*arctanh(a^-1//2*(a+b*x)^(1/2))+2//3*f*b^-1*(a+b*x)^3//2)
@test integrate(x^-1*(a+b*x)^(1/2)*(c+d*x)^-1*(e+f*x), x) == :(2*f*d^-1*(a+b*x)^(1/2)+-2*e*a^(1/2)*c^-1*arctanh(a^-1//2*(a+b*x)^(1/2))+2*c^-1*d^-3//2*(b*c+-1*a*d)^(1/2)*(d*e+-1*c*f)*arctan(d^(1/2)*(a+b*x)^(1/2)*(b*c+-1*a*d)^-1//2))
@test integrate(x^-1*(a+b*x)^(1/2)*(c+d*x)^-2*(e+f*x), x) == :(-2*e*a^(1/2)*c^-2*arctanh(a^-1//2*(a+b*x)^(1/2))+c^-1*d^-1*(a+b*x)^(1/2)*(c+d*x)^-1*(d*e+-1*c*f)+-1*c^-2*d^-3//2*(b*c+-1*a*d)^-1//2*(-1*b*c*(c*f+d*e)+2*a*e*d^2)*arctan(d^(1/2)*(a+b*x)^(1/2)*(b*c+-1*a*d)^-1//2))
@test integrate(x^-1*(a+b*x)^(1/2)*(c+d*x)^-3*(e+f*x), x) == :(-2*e*a^(1/2)*c^-3*arctanh(a^-1//2*(a+b*x)^(1/2))+(1/2)*c^-1*d^-1*(a+b*x)^(1/2)*(c+d*x)^-2*(d*e+-1*c*f)+-1//4*c^-3*d^-3//2*(b*c+-1*a*d)^-3//2*(-1*b^2*c^2*(c*f+3*d*e)+-8*e*a^2*d^3+12*a*b*c*e*d^2)*arctan(d^(1/2)*(a+b*x)^(1/2)*(b*c+-1*a*d)^-1//2)+-1//4*c^-2*d^-1*(a+b*x)^(1/2)*(c+d*x)^-1*(b*c+-1*a*d)^-1*(-1*b*c*(c*f+3*d*e)+4*a*e*d^2))
@test integrate(x^3*(a*x)^-1//2*(1+-1*a*x)^-1//2*(1+a*x), x) == :(-75//128*a^-4*arcsin(1+-2*a*x)+-75//64*a^-4*(a*x)^(1/2)*(1+-1*a*x)^(1/2)+-25//32*a^-4*(a*x)^3//2*(1+-1*a*x)^(1/2)+-5//8*a^-4*(a*x)^5//2*(1+-1*a*x)^(1/2)+-1//4*a^-4*(a*x)^7//2*(1+-1*a*x)^(1/2))
@test integrate(x^2*(a*x)^-1//2*(1+-1*a*x)^-1//2*(1+a*x), x) == :(-11//16*a^-3*arcsin(1+-2*a*x)+-11//8*a^-3*(a*x)^(1/2)*(1+-1*a*x)^(1/2)+-11//12*a^-3*(a*x)^3//2*(1+-1*a*x)^(1/2)+-1//3*a^-3*(a*x)^5//2*(1+-1*a*x)^(1/2))
@test integrate(x*(a*x)^-1//2*(1+-1*a*x)^-1//2*(1+a*x), x) == :(-7//8*a^-2*arcsin(1+-2*a*x)+-7//4*a^-2*(a*x)^(1/2)*(1+-1*a*x)^(1/2)+-1//2*a^-2*(a*x)^3//2*(1+-1*a*x)^(1/2))
@test integrate((a*x)^-1//2*(1+-1*a*x)^-1//2*(1+a*x), x) == :(-3//2*a^-1*arcsin(1+-2*a*x)+-1*a^-1*(a*x)^(1/2)*(1+-1*a*x)^(1/2))
@test integrate(x^-1*(a*x)^-1//2*(1+-1*a*x)^-1//2*(1+a*x), x) == :(-1*arcsin(1+-2*a*x)+-2*(a*x)^-1//2*(1+-1*a*x)^(1/2))
@test integrate(x^-2*(a*x)^-1//2*(1+-1*a*x)^-1//2*(1+a*x), x) == :(-10//3*a*(a*x)^-1//2*(1+-1*a*x)^(1/2)+-2//3*a*(a*x)^-3//2*(1+-1*a*x)^(1/2))
@test integrate(x^-3*(a*x)^-1//2*(1+-1*a*x)^-1//2*(1+a*x), x) == :(-12//5*a^2*(a*x)^-1//2*(1+-1*a*x)^(1/2)+-6//5*a^2*(a*x)^-3//2*(1+-1*a*x)^(1/2)+-2//5*a^2*(a*x)^-5//2*(1+-1*a*x)^(1/2))
@test integrate(x^-4*(a*x)^-1//2*(1+-1*a*x)^-1//2*(1+a*x), x) == :(-208//105*a^3*(a*x)^-1//2*(1+-1*a*x)^(1/2)+-104//105*a^3*(a*x)^-3//2*(1+-1*a*x)^(1/2)+-26//35*a^3*(a*x)^-5//2*(1+-1*a*x)^(1/2)+-2//7*a^3*(a*x)^-7//2*(1+-1*a*x)^(1/2))
@test integrate(x^-5*(a*x)^-1//2*(1+-1*a*x)^-1//2*(1+a*x), x) == :(-544//315*a^4*(a*x)^-1//2*(1+-1*a*x)^(1/2)+-272//315*a^4*(a*x)^-3//2*(1+-1*a*x)^(1/2)+-68//105*a^4*(a*x)^-5//2*(1+-1*a*x)^(1/2)+-34//63*a^4*(a*x)^-7//2*(1+-1*a*x)^(1/2)+-2//9*a^4*(a*x)^-9//2*(1+-1*a*x)^(1/2))
@test integrate(x^-2*(1+x)^-1//2*(-1+x)^-1//2*(-1+2*a*x), x) == :(2*a*arctan((1+x)^(1/2)*(-1+x)^(1/2))+-1*x^-1*(1+x)^(1/2)*(-1+x)^(1/2))
@test integrate(x^-2*(1+x)^-1//2*(-1+x)^-1//2*(-1*(1+-1*a*x)^2+a^2*x^2), x) == :(2*a*arctan((1+x)^(1/2)*(-1+x)^(1/2))+-1*x^-1*(1+x)^(1/2)*(-1+x)^(1/2))
@test integrate((a+b*x)^-1//2*(c+b*x*a^-1*(-1+c))^-1//2*(e+b*x*a^-1*(-1+e))^-1//2*(A+B*x), x) == :(-2*B*a^3//2*b^-2*(1+-1c)^-1//2*(1+-1*e)^-1*Elliptic.E(arcsin(a^-1//2*(1+-1c)^(1/2)*(a+b*x)^(1/2)),(1+-1c)^-1*(1+-1*e))+2*a^(1/2)*b^-2*(1+-1c)^-1//2*(1+-1*e)^-1*(A*(b+-1*b*e)+B*a*e)*Elliptic.F(arcsin(a^-1//2*(1+-1c)^(1/2)*(a+b*x)^(1/2)),(1+-1c)^-1*(1+-1*e)))
@test integrate((a+b*x)^-1//2*(c+d*x)^-1//2*(e+b*x*a^-1*(-1+e))^-1//2*(A+B*x), x) == :(2*a^(1/2)*b^-2*(b*(b*c+-1*a*d)^-1*(c+d*x))^(1/2)*(1+-1*e)^-3//2*(c+d*x)^-1//2*(A*(b+-1*b*e)+B*a*e)*Elliptic.F(arcsin(a^-1//2*(1+-1*e)^(1/2)*(a+b*x)^(1/2)),-1*a*d*(1+-1*e)^-1*(b*c+-1*a*d)^-1)+-2*B*a*b^-2*d^-1//2*(b*(b*c+-1*a*d)^-1*(c+d*x))^(1/2)*(1+-1*e)^-1*(c+d*x)^-1//2*(a*d+-1*b*c)^(1/2)*Elliptic.E(arcsin(d^(1/2)*(a+b*x)^(1/2)*(a*d+-1*b*c)^-1//2),-1*a^-1*d^-1*(1+-1*e)*(b*c+-1*a*d)))
@test integrate((1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^3, x) == :(-1182926269//1603800*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)+-6489123157//699840*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2)+-12243139//356400*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)+-17561//8910*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^2+-427//2970*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^3+2//55*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^4+522167393//139968*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3))
@test integrate((1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^2, x) == :(-5256763//97200*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)+-17746949//29160*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2)+-8141//2700*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)+-61//270*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^2+2//45*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^3+5592499//23328*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3))
@test integrate((1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x), x) == :(-20911//3780*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)+5//28*(1+4x)^3//2*(-5+2x)^3//2*(2+-3x)^(1/2)+136//105*(1+4x)^3//2*(-5+2x)^(1/2)*(2+-3x)^(1/2)+-954811//22680*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2)+72479//4536*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3))
@test integrate((1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2), x) == :(-22//45*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)+1//10*(1+4x)^3//2*(-5+2x)^(1/2)*(2+-3x)^(1/2)+-847//270*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2)+121//108*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3))
@test integrate((1+4x)^(1/2)*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^3, x) == :(46134551//38880*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)+-2161804579//326592*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3)+1//9*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^3+1679//756*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^2+26291//540*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)+2629157597//163296*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2))
@test integrate((1+4x)^(1/2)*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^2, x) == :(73207//1080*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)+-1679161//4536*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3)+1//7*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^2+173//60*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)+8198333//9072*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2))
@test integrate((1+4x)^(1/2)*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x), x) == :(1//4*(1+4x)^3//2*(-5+2x)^(1/2)*(2+-3x)^(1/2)+95//18*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)+-4543//216*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3)+1397//27*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2))
@test integrate((1+4x)^(1/2)*(-5+2x)^-1//2*(2+-3x)^(1/2), x) == :(1//3*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)+-11//9*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3)+55//18*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2))
@test integrate((1+4x)^-1//2*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^3, x) == :(110743//864*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)+-25260049//36288*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3)+5//28*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^2+121//24*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)+15629623//9072*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2))
@test integrate((1+4x)^-1//2*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^2, x) == :(68//9*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)+-17533//432*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3)+1//4*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)+44569//432*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2))
@test integrate((1+4x)^-1//2*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x), x) == :(5//12*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)+-179//72*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3)+241//36*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2))
@test integrate((1+4x)^-1//2*(-5+2x)^-1//2*(2+-3x)^(1/2), x) == :(1//4*22^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.E(arcsin(1//11*11^(1/2)*(1+4x)^(1/2)),3))
@test integrate((1+4x)^-1//2*(-5+2x)^-1//2*(2+-3x)^-1//2*(7+5x)^4, x) == :(-120355//288*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)+-5109835//756*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2)+-305//24*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)+-25//84*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^2+392989907//133056*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3))
@test integrate((1+4x)^-1//2*(-5+2x)^-1//2*(2+-3x)^-1//2*(7+5x)^3, x) == :(-2135//108*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)+-487585//1296*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2)+-5//12*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)+2474201//14256*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3))
@test integrate((1+4x)^-1//2*(-5+2x)^-1//2*(2+-3x)^-1//2*(7+5x)^2, x) == :(-25//36*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)+-2135//108*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2)+24353//2376*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3))
@test integrate((1+4x)^-1//2*(-5+2x)^-1//2*(2+-3x)^-1//2*(7+5x), x) == :(-5//6*11^(1/2)*(-5+2x)^(1/2)*(5+-2x)^-1//2*Elliptic.E(arcsin(2//11*11^(1/2)*(2+-3x)^(1/2)),-1//2)+13//22*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3))
@test integrate((1+4x)^-1//2*(-5+2x)^-1//2*(2+-3x)^-1//2, x) == :(1//33*66^(1/2)*(-5+2x)^-1//2*(5+-2x)^(1/2)*Elliptic.F(arcsin(1//11*33^(1/2)*(1+4x)^(1/2)),1//3))
@test integrate((c+d*x)^-1//2*(e+f*x)^-1//2*(g+h*x)^-1//2*(c*i+d*i*x), x) == :(2*i*f^-1*h^-1//2*(f*(f*g+-1*e*h)^-1*(g+h*x))^(1/2)*(-1*f*(d*e+-1*c*f)^-1*(c+d*x))^-1//2*(c+d*x)^(1/2)*(g+h*x)^-1//2*(e*h+-1*f*g)^(1/2)*Elliptic.E(arcsin(h^(1/2)*(e+f*x)^(1/2)*(e*h+-1*f*g)^-1//2),-1*d*h^-1*(d*e+-1*c*f)^-1*(f*g+-1*e*h)))
@test integrate((c+d*x)^-1//2*(e+f*x)^-1//2*(g+h*x)^-1//2*(a+b*x), x) == :(-2*d^-1*f^-1//2*h^-1*(d*(d*e+-1*c*f)^-1*(e+f*x))^(1/2)*(d*(d*g+-1*c*h)^-1*(g+h*x))^(1/2)*(e+f*x)^-1//2*(g+h*x)^-1//2*(c*f+-1*d*e)^(1/2)*(b*g+-1*a*h)*Elliptic.F(arcsin(f^(1/2)*(c+d*x)^(1/2)*(c*f+-1*d*e)^-1//2),h*f^-1*(d*g+-1*c*h)^-1*(d*e+-1*c*f))+2*b*d^-1*f^-1//2*h^-1*(d*(d*e+-1*c*f)^-1*(e+f*x))^(1/2)*(d*(d*g+-1*c*h)^-1*(g+h*x))^-1//2*(e+f*x)^-1//2*(g+h*x)^(1/2)*(c*f+-1*d*e)^(1/2)*Elliptic.E(arcsin(f^(1/2)*(c+d*x)^(1/2)*(c*f+-1*d*e)^-1//2),h*f^-1*(d*g+-1*c*h)^-1*(d*e+-1*c*f)))
@test integrate((1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-7//2, x) == :(-2852696//11598438735*(1+4x)^(1/2)*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^(1/2)+-2//25*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-5//2+17906//2085525*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-3//2+1426348//2319687747*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-1//2+-48884//220648545*253^(1/2)*((5+-2x)^-1*(7+5x))^-1//2*(-5+2x)^-1//2*(7+5x)^(1/2)*Elliptic.F(arctan((1/2)*2^(1/2)*(1+4x)^(1/2)*(2+-3x)^-1//2),-39//23)+1426348//11598438735*429^(1/2)*((5+-2x)^-1*(2+-3x))^-1//2*((5+-2x)^-1*(7+5x))^(1/2)*(2+-3x)^(1/2)*(7+5x)^-1//2*Elliptic.E(arcsin(1//23*897^(1/2)*(1+4x)^(1/2)*(-5+2x)^-1//2),-23//39))
@test integrate((1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-9//2, x) == :(-65687975672//2257624501329015*(1+4x)^(1/2)*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^(1/2)+-2//35*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-7//2+2558//695175*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-5//2+23758016//57992193675*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-3//2+32843987836//451524900265803*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-1//2+-1212290288//42949018635705*253^(1/2)*((5+-2x)^-1*(7+5x))^-1//2*(-5+2x)^-1//2*(7+5x)^(1/2)*Elliptic.F(arctan((1/2)*2^(1/2)*(1+4x)^(1/2)*(2+-3x)^-1//2),-39//23)+32843987836//2257624501329015*429^(1/2)*((5+-2x)^-1*(2+-3x))^-1//2*((5+-2x)^-1*(7+5x))^(1/2)*(2+-3x)^(1/2)*(7+5x)^-1//2*Elliptic.E(arcsin(1//23*897^(1/2)*(1+4x)^(1/2)*(-5+2x)^-1//2),-23//39))
@test integrate((1+4x)^(1/2)*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^-5//2, x) == :(-9350//3253419*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-1//2+2//117*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-3//2+3740//3253419*(1+4x)^(1/2)*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^(1/2)+44//61893*253^(1/2)*((5+-2x)^-1*(7+5x))^-1//2*(-5+2x)^-1//2*(7+5x)^(1/2)*Elliptic.F(arctan((1/2)*2^(1/2)*(1+4x)^(1/2)*(2+-3x)^-1//2),-39//23)+-1870//3253419*429^(1/2)*((5+-2x)^-1*(2+-3x))^-1//2*((5+-2x)^-1*(7+5x))^(1/2)*(2+-3x)^(1/2)*(7+5x)^-1//2*Elliptic.E(arcsin(1//23*897^(1/2)*(1+4x)^(1/2)*(-5+2x)^-1//2),-23//39))
@test integrate((1+4x)^(1/2)*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^-7//2, x) == :(-20464840//90467822133*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-1//2+-3646//16267095*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-3//2+2//195*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-5//2+8185936//90467822133*(1+4x)^(1/2)*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^(1/2)+111628//1721058651*253^(1/2)*((5+-2x)^-1*(7+5x))^-1//2*(-5+2x)^-1//2*(7+5x)^(1/2)*Elliptic.F(arctan((1/2)*2^(1/2)*(1+4x)^(1/2)*(2+-3x)^-1//2),-39//23)+-4092968//90467822133*429^(1/2)*((5+-2x)^-1*(2+-3x))^-1//2*((5+-2x)^-1*(7+5x))^(1/2)*(2+-3x)^(1/2)*(7+5x)^-1//2*Elliptic.E(arcsin(1//23*897^(1/2)*(1+4x)^(1/2)*(-5+2x)^-1//2),-23//39))
@test integrate((1+4x)^(1/2)*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^-9//2, x) == :(-40944441340//1956607901151813*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-1//2+-3217468//50259901185*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-3//2+2//273*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-7//2+98//1807455*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-5//2+16377776536//1956607901151813*(1+4x)^(1/2)*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^(1/2)+258506776//37222482817611*253^(1/2)*((5+-2x)^-1*(7+5x))^-1//2*(-5+2x)^-1//2*(7+5x)^(1/2)*Elliptic.F(arctan((1/2)*2^(1/2)*(1+4x)^(1/2)*(2+-3x)^-1//2),-39//23)+-8188888268//1956607901151813*429^(1/2)*((5+-2x)^-1*(2+-3x))^-1//2*((5+-2x)^-1*(7+5x))^(1/2)*(2+-3x)^(1/2)*(7+5x)^-1//2*Elliptic.E(arcsin(1//23*897^(1/2)*(1+4x)^(1/2)*(-5+2x)^-1//2),-23//39))
@test integrate((1+4x)^-1//2*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^-5//2, x) == :(-98330//74828637*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-1//2+-10//2691*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-3//2+39332//74828637*(1+4x)^(1/2)*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^(1/2)+716//1423539*253^(1/2)*((5+-2x)^-1*(7+5x))^-1//2*(-5+2x)^-1//2*(7+5x)^(1/2)*Elliptic.F(arctan((1/2)*2^(1/2)*(1+4x)^(1/2)*(2+-3x)^-1//2),-39//23)+-19666//74828637*429^(1/2)*((5+-2x)^-1*(2+-3x))^-1//2*((5+-2x)^-1*(7+5x))^(1/2)*(2+-3x)^(1/2)*(7+5x)^-1//2*Elliptic.E(arcsin(1//23*897^(1/2)*(1+4x)^(1/2)*(-5+2x)^-1//2),-23//39))
@test integrate((1+4x)^-1//2*(-5+2x)^-1//2*(2+-3x)^-1//2*(7+5x)^-1//2, x) == :(2//253*253^(1/2)*((5+-2x)^-1*(7+5x))^-1//2*(-5+2x)^-1//2*(7+5x)^(1/2)*Elliptic.F(arctan((1/2)*2^(1/2)*(1+4x)^(1/2)*(2+-3x)^-1//2),-39//23))
@test integrate((1+4x)^-1//2*(-5+2x)^-1//2*(2+-3x)^-1//2*(7+5x)^-5//2, x) == :(-895300//2319687747*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-1//2+-50//83421*(1+4x)^(1/2)*(-5+2x)^(1/2)*(2+-3x)^(1/2)*(7+5x)^-3//2+358120//2319687747*(1+4x)^(1/2)*(-5+2x)^-1//2*(2+-3x)^(1/2)*(7+5x)^(1/2)+103964//485426799*253^(1/2)*((5+-2x)^-1*(7+5x))^-1//2*(-5+2x)^-1//2*(7+5x)^(1/2)*Elliptic.F(arctan((1/2)*2^(1/2)*(1+4x)^(1/2)*(2+-3x)^-1//2),-39//23)+-179060//2319687747*429^(1/2)*((5+-2x)^-1*(2+-3x))^-1//2*((5+-2x)^-1*(7+5x))^(1/2)*(2+-3x)^(1/2)*(7+5x)^-1//2*Elliptic.E(arcsin(1//23*897^(1/2)*(1+4x)^(1/2)*(-5+2x)^-1//2),-23//39))
@test integrate((a+b*x)^-3//2*(c+d*x)^-1//2*(e+f*x)^-1//2*(g+h*x)^-1//2, x) == :(-2*d*((a+b*x)^-1*(d*e+-1*c*f)^-1*(c+d*x)*(b*e+-1*a*f))^(1/2)*(-1*(a+b*x)^-1*(f*g+-1*e*h)^-1*(g+h*x)*(b*e+-1*a*f))^-1//2*(c+d*x)^-1//2*(g+h*x)^(1/2)*(b*c+-1*a*d)^-1*(b*g+-1*a*h)^-1//2*(f*g+-1*e*h)^-1//2*Elliptic.F(arcsin((a+b*x)^-1//2*(e+f*x)^(1/2)*(b*g+-1*a*h)^(1/2)*(f*g+-1*e*h)^-1//2),-1*(b*g+-1*a*h)^-1*(d*e+-1*c*f)^-1*(b*c+-1*a*d)*(f*g+-1*e*h))+-2*b*((a+b*x)^-1*(d*e+-1*c*f)^-1*(c+d*x)*(b*e+-1*a*f))^-1//2*(-1*(a+b*x)^-1*(f*g+-1*e*h)^-1*(g+h*x)*(b*e+-1*a*f))^(1/2)*(c+d*x)^(1/2)*(g+h*x)^-1//2*(b*c+-1*a*d)^-1*(b*e+-1*a*f)^-1*(b*g+-1*a*h)^-1//2*(f*g+-1*e*h)^(1/2)*Elliptic.E(arcsin((a+b*x)^-1//2*(e+f*x)^(1/2)*(b*g+-1*a*h)^(1/2)*(f*g+-1*e*h)^-1//2),-1*(b*g+-1*a*h)^-1*(d*e+-1*c*f)^-1*(b*c+-1*a*d)*(f*g+-1*e*h)))
@test integrate((a+b*x)^-3//2*(c+d*x)^-3//2*(e+f*x)^-1//2*(g+h*x)^-1//2, x) == :(-2*b^3*(a+b*x)^-1//2*(c+d*x)^(1/2)*(e+f*x)^(1/2)*(g+h*x)^(1/2)*(b*c+-1*a*d)^-2*(b*e+-1*a*f)^-1*(b*g+-1*a*h)^-1+-2*d^3*(a+b*x)^(1/2)*(c+d*x)^-1//2*(e+f*x)^(1/2)*(g+h*x)^(1/2)*(b*c+-1*a*d)^-2*(d*e+-1*c*f)^-1*(d*g+-1*c*h)^-1+-4*b*d*((a+b*x)^-1*(d*e+-1*c*f)^-1*(c+d*x)*(b*e+-1*a*f))^(1/2)*(-1*(a+b*x)^-1*(f*g+-1*e*h)^-1*(g+h*x)*(b*e+-1*a*f))^-1//2*(c+d*x)^-1//2*(g+h*x)^(1/2)*(b*c+-1*a*d)^-2*(b*g+-1*a*h)^-1//2*(f*g+-1*e*h)^-1//2*Elliptic.F(arcsin((a+b*x)^-1//2*(e+f*x)^(1/2)*(b*g+-1*a*h)^(1/2)*(f*g+-1*e*h)^-1//2),-1*(b*g+-1*a*h)^-1*(d*e+-1*c*f)^-1*(b*c+-1*a*d)*(f*g+-1*e*h))+2*b*(a+b*x)^-1//2*(c+d*x)^(1/2)*(e+f*x)^(1/2)*(g+h*x)^(1/2)*(b*c+-1*a*d)^-2*(b*e+-1*a*f)^-1*(b*g+-1*a*h)^-1*(d*e+-1*c*f)^-1*(d*g+-1*c*h)^-1*(b^2*(f*h*c^2+-1*c*d*(e*h+f*g)+2*e*g*d^2)+f*h*a^2*d^2+-1*a*b*d^2*(e*h+f*g))+-2*((a+b*x)^-1*(d*e+-1*c*f)^-1*(c+d*x)*(b*e+-1*a*f))^-1//2*(-1*(a+b*x)^-1*(f*g+-1*e*h)^-1*(g+h*x)*(b*e+-1*a*f))^(1/2)*(c+d*x)^(1/2)*(g+h*x)^-1//2*(b*c+-1*a*d)^-2*(b*e+-1*a*f)^-1*(b*g+-1*a*h)^-1//2*(d*e+-1*c*f)^-1*(d*g+-1*c*h)^-1*(f*g+-1*e*h)^(1/2)*(b^2*(f*h*c^2+-1*c*d*(e*h+f*g)+2*e*g*d^2)+f*h*a^2*d^2+-1*a*b*d^2*(e*h+f*g))*Elliptic.E(arcsin((a+b*x)^-1//2*(e+f*x)^(1/2)*(b*g+-1*a*h)^(1/2)*(f*g+-1*e*h)^-1//2),-1*(b*g+-1*a*h)^-1*(d*e+-1*c*f)^-1*(b*c+-1*a*d)*(f*g+-1*e*h)))
@test integrate((a+b*x)^m*(c+d*x)*(e+f*x)*(g+h*x), x) == :(b^-4*(2+m)^-1*(a+b*x)^(2+m)*(b^2*(c*e*h+c*f*g+d*e*g)+-2*a*b*(c*f*h+d*e*h+d*f*g)+3*d*f*h*a^2)+-1*b^-4*(3+m)^-1*(a+b*x)^(3+m)*(-1*b*(c*f*h+d*e*h+d*f*g)+3*a*d*f*h)+d*f*h*b^-4*(4+m)^-1*(a+b*x)^(4+m)+b^-4*(1+m)^-1*(a+b*x)^(1+m)*(b*c+-1*a*d)*(b*e+-1*a*f)*(b*g+-1*a*h))
@test integrate((a+b*x)^m*(c+d*x)^(-4+-1m)*(e+f*x)*(g+h*x), x) == :(b^-1*d^-2*(3+m)^-1*(a+b*x)^(1+m)*(c+d*x)^(-3+-1m)*(b*c+-1*a*d)^-1*(b*(e*g*d^2+-1*c*d*(e*h+f*g)+-1*f*h*c^2*(2+m))+a*c*d*f*h*(3+m)+-1*d*f*h*x*(3+m)*(b*c+-1*a*d))+b^-1*d^-2*(2+m)^-1*(3+m)^-1*(a+b*x)^(1+m)*(c+d*x)^(-2+-1m)*(b*c+-1*a*d)^-2*(b^2*(2*e*g*d^2+c*d*(1+m)*(e*h+f*g)+f*h*c^2*(2+m^2+3m))+f*h*a^2*d^2*(6+m^2+5m)+-1*a*b*d*(3+m)*(d*(e*h+f*g)+2*c*f*h*(1+m)))+d^-2*(1+m)^-1*(2+m)^-1*(3+m)^-1*(a+b*x)^(1+m)*(c+d*x)^(-1+-1m)*(b*c+-1*a*d)^-3*(b^2*(2*e*g*d^2+c*d*(1+m)*(e*h+f*g)+f*h*c^2*(2+m^2+3m))+f*h*a^2*d^2*(6+m^2+5m)+-1*a*b*d*(3+m)*(d*(e*h+f*g)+2*c*f*h*(1+m))))
@test integrate((a+b*x)^m*(c+d*x)^(-5+-1m)*(e+f*x)*(g+h*x), x) == :((1/2)*b^-1*d^-2*(4+m)^-1*(a+b*x)^(1+m)*(c+d*x)^(-4+-1m)*(b*c+-1*a*d)^-1*(b*(-2*c*d*(e*h+f*g)+2*e*g*d^2+-1*f*h*c^2*(2+m))+a*c*d*f*h*(4+m)+-1*d*f*h*x*(4+m)*(b*c+-1*a*d))+d^-2*(2+m)^-1*(3+m)^-1*(4+m)^-1*(a+b*x)^(1+m)*(c+d*x)^(-2+-1m)*(b*c+-1*a*d)^-3*(b^2*(6*e*g*d^2+f*h*c^2*(2+m^2+3m)+2*c*d*(1+m)*(e*h+f*g))+f*h*a^2*d^2*(12+m^2+7m)+-2*a*b*d*(4+m)*(d*(e*h+f*g)+c*f*h*(1+m)))+(1/2)*b^-1*d^-2*(3+m)^-1*(4+m)^-1*(a+b*x)^(1+m)*(c+d*x)^(-3+-1m)*(b*c+-1*a*d)^-2*(b^2*(6*e*g*d^2+f*h*c^2*(2+m^2+3m)+2*c*d*(1+m)*(e*h+f*g))+f*h*a^2*d^2*(12+m^2+7m)+-2*a*b*d*(4+m)*(d*(e*h+f*g)+c*f*h*(1+m)))+b*d^-2*(1+m)^-1*(2+m)^-1*(3+m)^-1*(4+m)^-1*(a+b*x)^(1+m)*(c+d*x)^(-1+-1m)*(b*c+-1*a*d)^-4*(b^2*(6*e*g*d^2+f*h*c^2*(2+m^2+3m)+2*c*d*(1+m)*(e*h+f*g))+f*h*a^2*d^2*(12+m^2+7m)+-2*a*b*d*(4+m)*(d*(e*h+f*g)+c*f*h*(1+m))))
@test integrate((c+d*x)^(-4+-1m)*(e+f*x)^m*(a+b*x)*(g+h*x), x) == :(-1*d^-2*f^-1*(3+m)^-1*(c+d*x)^(-3+-1m)*(e+f*x)^(1+m)*(d*e+-1*c*f)^-1*(-1*b*c*(d*(f*g+-1*e*h*(3+m))+c*f*h*(2+m))+a*d*f*(d*g+-1*c*h)+b*d*h*x*(3+m)*(d*e+-1*c*f))+d^-2*f^-1*(2+m)^-1*(3+m)^-1*(c+d*x)^(-2+-1m)*(e+f*x)^(1+m)*(d*e+-1*c*f)^-2*(b*(h*c^2*f^2*(2+m^2+3m)+-1*e*d^2*(3+m)*(f*g+-1*e*h*(2+m))+c*d*f*(1+m)*(f*g+-2*e*h*(3+m)))+a*d*f*(d*(2*f*g+-1*e*h*(3+m))+c*f*h*(1+m)))+-1*d^-2*(1+m)^-1*(2+m)^-1*(3+m)^-1*(c+d*x)^(-1+-1m)*(e+f*x)^(1+m)*(d*e+-1*c*f)^-3*(b*(h*c^2*f^2*(2+m^2+3m)+-1*e*d^2*(3+m)*(f*g+-1*e*h*(2+m))+c*d*f*(1+m)*(f*g+-2*e*h*(3+m)))+a*d*f*(d*(2*f*g+-1*e*h*(3+m))+c*f*h*(1+m))))
@test integrate((c+d*x)^(-4+-1m)*(e+f*x)^m*(g+h*x), x) == :(-1*d^-1*(3+m)^-1*(c+d*x)^(-3+-1m)*(e+f*x)^(1+m)*(d*e+-1*c*f)^-1*(d*g+-1*c*h)+d^-1*(2+m)^-1*(3+m)^-1*(c+d*x)^(-2+-1m)*(e+f*x)^(1+m)*(d*e+-1*c*f)^-2*(d*(2*f*g+-1*e*h*(3+m))+c*f*h*(1+m))+-1*f*d^-1*(1+m)^-1*(2+m)^-1*(3+m)^-1*(c+d*x)^(-1+-1m)*(e+f*x)^(1+m)*(d*e+-1*c*f)^-3*(d*(2*f*g+-1*e*h*(3+m))+c*f*h*(1+m)))
@test integrate(x*(1+d*x)^-1//2*(1+-1*d*x)^-1//2*(a+b*x+c*x^2), x) == :((1/2)*b*d^-3*arcsin(d*x)+-1//6*d^-4*(1+-1*d^2*x^2)^(1/2)*(4c+6*a*d^2+3*b*x*d^2)+-1//3*c*d^-2*x^2*(1+-1*d^2*x^2)^(1/2))
@test integrate((1+d*x)^-1//2*(1+-1*d*x)^-1//2*(a+b*x+c*x^2), x) == :((1/2)*d^-3*(c+2*a*d^2)*arcsin(d*x)+-1*b*d^-2*(1+-1*d^2*x^2)^(1/2)+-1//2*c*x*d^-2*(1+-1*d^2*x^2)^(1/2))
@test integrate(x^-1*(1+d*x)^-1//2*(1+-1*d*x)^-1//2*(a+b*x+c*x^2), x) == :(-1*a*arctanh((1+-1*d^2*x^2)^(1/2))+b*d^-1*arcsin(d*x)+-1*c*d^-2*(1+-1*d^2*x^2)^(1/2))
@test integrate(x^-2*(1+d*x)^-1//2*(1+-1*d*x)^-1//2*(a+b*x+c*x^2), x) == :(-1*b*arctanh((1+-1*d^2*x^2)^(1/2))+c*d^-1*arcsin(d*x)+-1*a*x^-1*(1+-1*d^2*x^2)^(1/2))
@test integrate(x^-3*(1+d*x)^-1//2*(1+-1*d*x)^-1//2*(a+b*x+c*x^2), x) == :(-1*(c+(1/2)*a*d^2)*arctanh((1+-1*d^2*x^2)^(1/2))+-1*b*x^-1*(1+-1*d^2*x^2)^(1/2)+-1//2*a*x^-2*(1+-1*d^2*x^2)^(1/2))

max_stars_repo_path: test/1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.4 (a+b x)^m (c+d x)^n (e+f x)^p (g+h x)^q.jl


text:
"A trait for the type of bins a histogram may have."
abstract type BinType end
"Each bin has the same predetermined width."
struct FixedWidth <: BinType end
"Bins have possibly different predetermined widths."
struct VariableWidth <: BinType end
"Bin widths are not known before computing the histogram (i.e. text data). Only 1D histograms are supported."
struct UnboundedWidth <: BinType end

"A trait for the ways the bin search and bin update steps can be parallelized."
abstract type HistogramParallelization end

"No threading nor vectorization."
struct NoParallelization <: HistogramParallelization end

"""
Threads that have private bin data structures that are reduced after their private updates.

Requires these functions to be defined for real-valued histograms:
- `get_subweights(hist)::AbstractArray{Int,N+1}` Returns the weights (i.e. counts) array for an N-dimensional histogram.

Requires these functions to be defined for text histograms:
- `get_subtable(hist)::AbstractVector{AbstractDict{String,Int}}` Returns a vector of independent tables.
"""
struct PrivateThreads <: HistogramParallelization end

"""
SIMD vectorization.

Requires these functions to be defined for real-valued histograms:
- `get_subweights(hist)::AbstractArray{Int,N+1}` Returns the weights (i.e. counts) array for an N-dimensional histogram.

Requires these functions to be defined for text histograms:
- `get_subtable(hist)::AbstractVector{AbstractDict{String,Int}}` Returns a vector of independent tables.
"""
struct SIMD <: HistogramParallelization end

ext: jl


text:
export AbstractLowRankOperator, LowRankOperator

abstract AbstractLowRankOperator{T} <: Operator{T}

immutable LowRankOperator{S<:Space,T} <: AbstractLowRankOperator{T}
    U::Vector{Fun{S,T}}
    V::Vector{Operator{T}}

    function LowRankOperator(U::Vector{Fun{S,T}},V::Vector{Operator{T}})
        @assert all(isafunctional,V)

        @assert length(U) == length(V)
        @assert length(U) > 0
        ds=domainspace(first(V))
        for k=2:length(V)
            @assert domainspace(V[k])==ds
        end
        rs=space(first(U))
        for k=2:length(U)
            @assert space(U[k])==rs
        end
        new(U,V)
    end
end

LowRankOperator{S,T}(U::Vector{Fun{S,T}},V::Vector{Operator{T}})=LowRankOperator{S,T}(U,V)
LowRankOperator{S,T1,T2}(U::Vector{Fun{S,T1}},V::Vector{Operator{T2}})=LowRankOperator(convert(Vector{Fun{S,promote_type(T1,T2)}},U),
                                                                                convert(Vector{Operator{promote_type(T1,T2)}},V))
LowRankOperator{FF<:Fun,FT<:Operator}(U::Vector{FF},V::Vector{FT})=LowRankOperator(U,convert(Vector{Operator{eltype(FT)}},V))

LowRankOperator(B::AbstractVector,S...)=LowRankOperator(convert(Vector{Operator{Float64}},B),S...)

LowRankOperator(A::Fun,B::Operator)=LowRankOperator([A],[B])

datasize(L::LowRankOperator,k) =
    k==1?mapreduce(ncoefficients,max,L.U):mapreduce(bandwidth,max,L.V)
datasize(L::LowRankOperator) = datasize(L,1),datasize(L,2)
bandinds(L::LowRankOperator) = datasize(L,1)-1,datasize(L,2)-1

domainspace(L::LowRankOperator)=domainspace(first(L.V))
rangespace(L::LowRankOperator)=space(first(L.U))
promoterangespace(L::LowRankOperator,sp::Space)=LowRankOperator(map(u->Fun(u,sp),L.U),L.V)
promotedomainspace(L::LowRankOperator,sp::Space)=LowRankOperator(L.U,map(v->promotedomainspace(v,sp),L.V))

function Base.getindex(L::LowRankOperator,k::Integer,j::Integer)
    ret=zero(eltype(L))
    for p=1:length(L.U)
        if k≤ncoefficients(L.U[p])
            ret+=L.U[p].coefficients[k]*L.V[p][j]
        end
    end
    ret
end

Base.rank(L::LowRankOperator) = length(L.U)

-(L::LowRankOperator) = LowRankOperator(-L.U,L.V)

*(L::LowRankOperator,f::Fun) = sum(map((u,v)->u*(v*f),L.U,L.V))

*(A::LowRankOperator,B::LowRankOperator) = LowRankOperator((A.V*B.U.').'*A.U,B.V)
# avoid ambiguituy
for TYP in (:TimesOperator,:PlusOperator,:Conversion,:Operator)
    @eval *(L::LowRankOperator,B::$TYP) = LowRankOperator(L.U,map(v->v*B,L.V))
    @eval *(B::$TYP,L::LowRankOperator) = LowRankOperator(map(u->B*u,L.U),L.V)
end

+(A::LowRankOperator,B::LowRankOperator) = LowRankOperator([A.U;B.U],[A.V;B.V])
-(A::LowRankOperator,B::LowRankOperator) = LowRankOperator([A.U;-B.U],[A.V;B.V])


text:
using Pkg
tempdir = mktempdir()
Pkg.activate(tempdir)
Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "..")))
Pkg.add(["BenchmarkTools", "PkgBenchmark"])
Pkg.resolve()

using Convex: Convex, ProblemDepot
using BenchmarkTools

const SUITE = BenchmarkGroup()

problems =  [
                "constant_fix!_with_complex_numbers",
                "affine_dot_multiply_atom",
                "affine_hcat_atom",
                "affine_trace_atom",
                "exp_entropy_atom",
                "exp_log_perspective_atom",
                "socp_norm_2_atom",
                "socp_quad_form_atom",
                "socp_sum_squares_atom",
                "lp_norm_inf_atom",
                "lp_maximum_atom",
                "sdp_and_exp_log_det_atom",
                "sdp_norm2_atom",
                "sdp_lambda_min_atom",
                "sdp_sum_largest_eigs",
                "mip_integer_variables",
            ]

SUITE["formulation"] = ProblemDepot.benchmark_suite(Convex.conic_problem, problems)


text:
###########################################################################
#
#  Acceptance Rejection Sampling (ARS)
#
#  Parameters:
#    -logCandidate: unscaled candidate
#    -logCandidateScalingFactor: scale factor to ensure the scaled-up
#        logCandidate covers target
#    -scale: for scaling the jumps
#    -tuner: for tuning the scale parameter
#
###########################################################################

###########################################################################
#
# Based on Bolstad: Computational Bayesian Statistics, section 2.1
#
###########################################################################

export ARS

println("Loading ARS(logCandidate, logCandidateScalingFactor, scale, tuner) sampler")

###########################################################################
#                  ARS specific 'tuners'
###########################################################################
abstract ARSTuner <: MCMCTuner

###########################################################################
#                  ARS type
###########################################################################

immutable ARS <: MCMCSampler
  logCandidate::Function
  logCandidateScalingFactor::Float64
  scale::Float64
  tuner::Union(Nothing, ARSTuner)

  function ARS(logCandidate::Function, logCandidateScalingFactor::Float64, x::Float64, t::Union(Nothing, ARSTuner))
    @assert typeof(logCandidate) == Function "Unscaled candidate logCandidate should be a function"
    @assert x>0 "scale should be > 0"
    new(logCandidate, logCandidateScalingFactor, x, t)
  end
end
ARS(logCandidate::Function) = ARS(logCandidate, 1., 1., nothing)
ARS(logCandidate::Function, logCandidateScalingFactor::Float64) = ARS(logCandidate, logCandidateScalingFactor, 1.0, nothing)
ARS(logCandidate::Function, t::ARSTuner) = ARS(logCandidate, 1., 1., t)
ARS(logCandidate::Function; logCandidateScalingFactor::Float64 = 1.0, scale::Float64=1.0, tuner::Union(Nothing, ARSTuner)=nothing) = ARS(logCandidate, logCandidateScalingFactor, scale, tuner)

###########################################################################
#                  ARS task
###########################################################################

# ARS sampling
function SamplerTask(model::MCMCModel, sampler::ARS, runner::MCMCRunner)
	local pars, proposedPars
	local logTarget, proposedLogTarget, proposedLogCandidate
  local scale, weight, mu

	# hook inside Task to allow remote resetting
	task_local_storage(:reset, (resetPars::Vector{Float64}) -> (pars = copy(resetPars); logTarget = model.eval(pars)))

	# initialization
	scale = model.scale .* sampler.scale  # rescale model scale by sampler scale
	pars = copy(model.init)
	logTarget = model.eval(pars)
  @assert isfinite(logTarget) "Initial values out of model support, try other values"

	# main loop
	while true
  	proposedPars = pars + randn(model.size) .* scale
    mu = log(rand())
    proposedLogTarget = model.eval(proposedPars)
		proposedLogCandidate = sampler.logCandidate(proposedPars)
	  weight = proposedLogTarget .- sampler.logCandidateScalingFactor .- proposedLogCandidate
		if weight > mu
			ms = MCMCSample(proposedPars, proposedLogTarget, pars, logTarget,
        {"accept" => true})
	    produce(ms)
	    pars, logTarget = copy(proposedPars), copy(proposedLogTarget)
	  else
			ms = MCMCSample(pars, logTarget, pars, logTarget,
        {"accept" => false})
      produce(ms)
    end
	end
end


text:
# test Utils module
mutable struct TestParams{Tdim} <: AbstractParamType{Tdim}
  time::Timings
end

mutable struct TestData{Tsol, Tres} <: AbstractSolutionData{Tsol, Tres}
  params::TestParams{2}
  M::Array{Float64, 1}
  Minv::Array{Float64, 1}
  comm::MPI.Comm
end

mutable struct TestMesh{Tmsh} <: AbstractMesh{Tmsh}
  jac::AbstractArray{Float64, 2}
  dofs::AbstractArray{Int, 3}
  comm::MPI.Comm
  commsize::Int
  numNodes::Int
  numNodesPerElement::Int
  numEl::Int
  numDofPerNode::Int
  min_node_dist::Float64
  dim::Int
end

mutable struct FakeSBP
end

"""
  Test BufferedIO type
"""
function test_utils_io()
  @testset "---- Testing IO -----" begin
    fname = "iotest.dat"
    f = open(fname, "w")
    f2 = BufferedIO(f)

    println(f2, 1)
    println(f2, 2)

    fstat = stat(fname)
    @test ( fstat.size )== 0
    flush(f2)
  #  close(f)

    data = readdlm(fname)
    @test ( data[1] )== 1
    @test ( data[2] )== 2

    close(f)
  end

  return nothing
end

#test_utils_io()
add_func1!(EulerTests, test_utils_io, [TAG_UTILS, TAG_SHORTTEST])

"""
  Test Timings type
"""
function test_utils_timing()
  @testset "----- Testing Timings -----" begin
    t = Timings()
  #=
    t.t_volume = 1.0
    t.t_face = 2.0
    t.t_source = 3.0
    t.t_sharedface = 4.0
    t.t_bndry = 5.0
    t.t_send = 6.0
    t.t_wait = 7.0
    t.t_allreduce = 8.0
    t.t_jacobian = 9.0
    t.t_solve = 10.0
    t.t_barrier = 11.0
    t.t_barrier2 = 12.0
    t.t_barrier3 = 13.0
    t.t_barriers[1:7] = 14.0:20.0

    write_timings(t, "timing_breakdown")
    vals = readdlm("timing_breakdown.dat")
    for i=1:20
      @test isapprox( vals[i], Float64(i))
    end
  =#

    function mysleep(t)
      sleep(t)
      return 42
    end

    @time_all mysleep(0.0001)
    gc()
    val, telapsed, gc_time, alloced = @time_all mysleep(0.1)
    @test ( val )== 42
    @test isapprox( telapsed, 0.1) atol=0.05
    @test isapprox( gc_time, 0.0)

  end

  return nothing
end

#test_utils_timing()
add_func1!(EulerTests, test_utils_timing, [TAG_UTILS, TAG_SHORTTEST])

"""
  Test various functions: calcNorm, calcMeshH etc.
"""
function test_utils_misc()
  @testset "----- Testing Utility Functions -----" begin
    M = rand(10)
    Minv = 1./M
    params = TestParams{2}(Timings())
    obj = TestData{Float64, Float64}(params, M, Minv, MPI.COMM_WORLD)
    data = rand(10)

    norm1 = calcNorm(obj, data)
    norm2 = sqrt(sum(data.*M.*data))

    @test isapprox( norm1, norm2)

    norm1 = calcNorm(obj, data, strongres=true)
    norm2 = sqrt(sum(data.*Minv.*data))

    data = complex.(rand(10), rand(10))
    norm1 = calcNorm(obj, data)
    norm2 = sqrt(sum(real(data).*M.*real(data)))

    @test isapprox( norm1, norm2)

    norm1 = calcNorm(obj, data, strongres=true)
    norm2 = sqrt(sum(real(data).*Minv.*real(data)))

    @test isapprox( norm1, norm2)

    numEl = 5
    numNodesPerElement = 3
    numNodes = numEl*numNodesPerElement
    numDofPerNode = 1
    jac = zeros(Float64, numNodesPerElement, numEl)
    fill!(jac, 2)
    dofs = zeros(Int, numDofPerNode, numNodesPerElement, numEl)
    for i=1:length(dofs)
      dofs[i] = i
    end
    mesh = TestMesh{Float64}(jac, dofs, MPI.COMM_WORLD, MPI.Comm_size(MPI.COMM_WORLD), numNodes, numNodesPerElement, numEl, numDofPerNode, 0.25, 2)
    opts = Dict{Any, Any}()
    @test isapprox( calcMeshH(mesh, FakeSBP(), obj, opts), mesh.min_node_dist*1./sqrt(2))

  end

  return nothing
end

#test_utils_misc()
add_func1!(EulerTests, test_utils_misc, [TAG_UTILS, TAG_SHORTTEST])

"""
  Test permutation matrices, vectors, functions
"""
function test_utils_perm()
  @testset "----- Testing Permutation Functions -----" begin

    # test permuting rows
    A = [1 2 3 10 ; 4 5 6 11; 7 8 9 12]
    B = zeros(A)
    permvec = [2, 3, 1]
    P = permMatrix(permvec)

    A_orig = copy(A)
    applyPermRowInplace(permvec, A, B)
    A2 = P*A_orig

    for i=1:3
      @test ( A2[i, :] )== A[i, :]
    end

    A = [1 2 3; 4 5 6; 7 8 9]
    A_orig = copy(A)
    permvec = [3, 1, 2]
    P = permMatrix(permvec)

    A2 = P*A
    applyPermRowInplace(permvec, A, B)

    for i=1:3
      @test ( A2[i, :] )== A[i, :]
    end

    #  test inverse permutation
    invperm = zeros(permvec)
    inversePerm(permvec, invperm)
    applyPermRowInplace(invperm, A, B)

    for i=1:length(A)
      @test ( A[i] )== A_orig[i]
    end

    P2 = permMatrix(invperm)

    @test ( P2.' )== P

    # test permuting columns
    A = [1 2 3; 4 5 6; 7 8 9; 10 11 12]
    B = zeros(A)
    A_orig = copy(A)
    P = permMatrix(permvec)

    applyPermColumnInplace(permvec, A, B)
    A2 = A_orig*P

    for i=1:3
      @test ( (A_orig*P)[:, i] )== A[ :, i ]
    end
  end  # end facts block

  return nothing
end

#test_utils_perm()
add_func1!(EulerTests, test_utils_perm, [TAG_UTILS, TAG_SHORTTEST])

function test_utils_projection()
  @testset "----- testing Projection -----" begin
    t = Timings()
    params = TestParams{2}(t)
    params3 = TestParams{3}(t)

    nrm = zeros(2)
    P = zeros(4,4)
    nrm3 = zeros(3)
    P3 = zeros(5,5)

    q = [1.0, 2.0, 3.0, 7.0]
    qprime = zeros(q)
    qprime2 = zeros(q)
    q2 = zeros(q)
    q2a = zeros(q)
    q3 = [1.0, 2.0, 3.0, 4.0, 15.0]
    q3prime = zeros(q3)
    q3prime2 = zeros(q3)
    q32 = zeros(q3)
    q32a = zeros(q3)

    for theta=0.0:0.1:2*pi
      nrm[1] = cos(theta)
      nrm[2] = sin(theta)

      t1, t2 = Utils.getOrthogonalVector(params, nrm)
      dot_val = abs.(t1*nrm[1] + t2*nrm[2])
      # check normal vector calculation
      @test isapprox( sqrt(t1*t1 + t2*t2), 1.0) atol=1e-12
      @test  dot_val  < 1e-13

      # check projection
      getProjectionMatrix(params, nrm, P)

      # check P is unitary
      @test isapprox( norm(inv(P) - P.'), 0.0) atol=1e-12

      # check q1, q4 are same, magnitude of momentum is same
      A_mul_B!(qprime, P, q)
      @test isapprox( q[1], qprime[1]) atol=1e-13
      @test isapprox( q[4], qprime[4]) atol=1e-13
      @test isapprox( q[2]*q[2] + q[3]*q[3], qprime[2]*qprime[2] + qprime[3]*qprime[3]) atol=1e-12

      # check multiplication functions against A_mul_b
      projectToNT(params, P, q, qprime2)
      @test isapprox( qprime, qprime2) atol=1e-12

      # project back
      At_mul_B!(q2, P, qprime)
      @test isapprox( q2, q) atol=1e-12
      projectToXY(params, P, qprime, q2a)
      @test isapprox( q2a, q) atol=1e-12

    end

    # 3D test
    for theta=0:0.1:pi
      for phi=0:0.1:2*pi
        nrm3[1] = sin(theta)*cos(phi)
        nrm3[2] = sin(theta)*sin(phi)
        nrm3[3] = cos(theta)

        t1, t2, t3 = Utils.getOrthogonalVector(params3, nrm3)
        dot_val = abs.(t1*nrm3[1] + t2*nrm3[2] + t3*nrm3[3])
        @test isapprox( t1*t1 + t2*t2 + t3*t3, 1.0) atol=1e-12
        @test isapprox( dot_val, 0.0) atol=1e-13

        getProjectionMatrix(params3, nrm3, P3)

        @test isapprox( norm(inv(P3) - P3.'), 0.0) atol=1e-12

        # check q1, q4 are same, magnitude of momentum is same
        A_mul_B!(q3prime, P3, q3)
        @test isapprox( q3[1], q3prime[1]) atol=1e-13
        @test isapprox( q3[5], q3prime[5]) atol=1e-13
        @test isapprox( q3[2]*q[2] + q3[3]*q3[3] + q3[4]*q3[4], q3prime[2]*q3prime[2] + q3prime[3]*q3prime[3] + q3prime[4]*q3prime[4]) atol=1e-12

        # check multiplication functions against A_mul_b
        projectToNT(params3, P3, q3, q3prime2)
        @test isapprox( q3prime, q3prime2) atol=1e-12

        # project back
        At_mul_B!(q32, P3, q3prime)
        @test isapprox( q32, q3) atol=1e-12
        projectToXY(params3, P3, q3prime, q32a)
        @test isapprox( q32a, q3) atol=1e-12

      end
    end

  end  # end facts block

  return nothing
end

add_func1!(EulerTests, test_utils_projection, [TAG_UTILS, TAG_SHORTTEST])


text:
function capacityio(;
  refreshio = PARAM[:refreshio])
  #provides an easy way to turn off writing the output
  if !PARAM[:refreshanyio]
    map(keys(refreshio) |> collect) do k
      refreshio[k] = false
    end
  end

  refreshio[:comomentumtable] && comomentumtables()
  refreshio[:fundtable] && fundtables()
  refreshio[:lltable] && lltables()
  #refreshio[:comomentumgraphs] && comomentumgraphs()
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
"""
    worldclim(layer::Integer; resolution::Float64=10.0, left=nothing, right=nothing, bottom=nothing, top=nothing)

Download and prepare a WorldClim 2.1 bioclimatic variable, and returns it as an
`SimpleSDMPredictor`. Layers are called by their number, from 1 to 19. The list
of available layers is given in a table below.

The keywords are `resolution`, which must be a floating point value, and either
`2.5`, `5.0`, or `10.0`, as well as optionally `left`, `right`, `bottom`, and
`top`, which will allow to only load parts of the data.

Internally, this function will download the main zip file for the required
resolution from the WordlClim website, extract it, and parse the required
layers.

It is recommended to *keep* the content of the `path` folder, as it will
eliminate the need to download and/or extract the tiff files. For example,
calling `wordlclim(1:19)` will download and extract everything, and future calls
will be much faster.

| Variable | Description                                                |
| ------   | ------                                                     |
| 1        | Annual Mean Temperature                                    |
| 2        | Mean Diurnal Range (Mean of monthly (max temp - min temp)) |
| 3        | Isothermality (BIO2/BIO7) (* 100)                          |
| 4        | Temperature Seasonality (standard deviation *100)          |
| 5        | Max Temperature of Warmest Month                           |
| 6        | Min Temperature of Coldest Month                           |
| 7        | Temperature Annual Range (BIO5-BIO6)                       |
| 8        | Mean Temperature of Wettest Quarter                        |
| 9        | Mean Temperature of Driest Quarter                         |
| 10       | Mean Temperature of Warmest Quarter                        |
| 11       | Mean Temperature of Coldest Quarter                        |
| 12       | Annual Precipitation                                       |
| 13       | Precipitation of Wettest Month                             |
| 14       | Precipitation of Driest Month                              |
| 15       | Precipitation Seasonality (Coefficient of Variation)       |
| 16       | Precipitation of Wettest Quarter                           |
| 17       | Precipitation of Driest Quarter                            |
| 18       | Precipitation of Warmest Quarter                           |
| 19       | Precipitation of Coldest Quarter                           |

Original data: https://www.worldclim.org/data/worldclim21.html
"""
function worldclim(layer::Integer; resolution::Float64=10.0, left=nothing, right=nothing, bottom=nothing, top=nothing)
    return raster(SimpleSDMPredictor, WorldClim(resolution), layer=layer, left=left, right=right, bottom=bottom, top=top)
end

worldclim(layers::Vector{T}; args...) where {T <: Integer} = [worldclim(l; args...) for l in layers]
worldclim(layers::UnitRange{T}; args...) where {T <: Integer} = [worldclim(l; args...) for l in layers]


text:
get_codigo_elemento(::Pais)::String = "81"
get_valor_elemento(p::Pais) = String(p.codigo)

CONST_PAISES = [
	:AD,
	:AE,
	:AF,
	:AG,
	:AI,
	:AL,
	:AM,
	:AO,
	:AQ,
	:AR,
	:AS,
	:AT,
	:AU,
	:AW,
	:AX,
	:AZ,
	:BA,
	:BB,
	:BD,
	:BE,
	:BF,
	:BG,
	:BH,
	:BI,
	:BJ,
	:BL,
	:BM,
	:BN,
	:BO,
	:BQ,
	:BR,
	:BS,
	:BT,
	:BV,
	:BW,
	:BY,
	:BZ,
	:CA,
	:CC,
	:CD,
	:CF,
	:CG,
	:CH,
	:CI,
	:CK,
	:CL,
	:CM,
	:CN,
	:CO,
	:CR,
	:CU,
	:CV,
	:CW,
	:CX,
	:CY,
	:CZ,
	:DE,
	:DJ,
	:DK,
	:DM,
	:DO,
	:DZ,
	:EC,
	:EE,
	:EG,
	:EH,
	:ER,
	:ES,
	:ET,
	:FI,
	:FJ,
	:FK,
	:FM,
	:FO,
	:FR,
	:GA,
	:GB,
	:GD,
	:GE,
	:GF,
	:GG,
	:GH,
	:GI,
	:GL,
	:GM,
	:GN,
	:GP,
	:GQ,
	:GR,
	:GS,
	:GT,
	:GU,
	:GW,
	:GY,
	:HK,
	:HM,
	:HN,
	:HR,
	:HT,
	:HU,
	:ID,
	:IE,
	:IL,
	:IM,
	:IN,
	:IO,
	:IQ,
	:IR,
	:IS,
	:IT,
	:JE,
	:JM,
	:JO,
	:JP,
	:KE,
	:KG,
	:KH,
	:KI,
	:KM,
	:KN,
	:KP,
	:KR,
	:KW,
	:KY,
	:KZ,
	:LA,
	:LB,
	:LC,
	:LI,
	:LK,
	:LR,
	:LS,
	:LT,
	:LU,
	:LV,
	:LY,
	:MA,
	:MC,
	:MD,
	:ME,
	:MF,
	:MG,
	:MH,
	:MK,
	:ML,
	:MM,
	:MN,
	:MO,
	:MP,
	:MQ,
	:MR,
	:MS,
	:MT,
	:MU,
	:MV,
	:MW,
	:MX,
	:MY,
	:MZ,
	:NA,
	:NC,
	:NE,
	:NF,
	:NG,
	:NI,
	:NL,
	:NO,
	:NP,
	:NR,
	:NU,
	:NZ,
	:OM,
	:PA,
	:PE,
	:PF,
	:PG,
	:PH,
	:PK,
	:PL,
	:PM,
	:PN,
	:PR,
	:PS,
	:PT,
	:PW,
	:PY,
	:QA,
	:RE,
	:RO,
	:RS,
	:RU,
	:RW,
	:SA,
	:SB,
	:SC,
	:SD,
	:SE,
	:SG,
	:SH,
	:SI,
	:SJ,
	:SK,
	:SL,
	:SM,
	:SN,
	:SO,
	:SR,
	:SS,
	:ST,
	:SV,
	:SX,
	:SY,
	:SZ,
	:TC,
	:TD,
	:TF,
	:TG,
	:TH,
	:TJ,
	:TK,
	:TL,
	:TM,
	:TN,
	:TO,
	:TR,
	:TT,
	:TV,
	:TW,
	:TZ,
	:UA,
	:UG,
	:UM,
	:US,
	:UY,
	:UZ,
	:VA,
	:VC,
	:VE,
	:VG,
	:VI,
	:VN,
	:VU,
	:WF,
	:WS,
	:YE,
	:YT,
	:ZA,
	:ZM,
	:ZW
]


text:
@testset "BostonHousing" begin
    n_obs = 506
    n_features = 13
    n_targets = 1
    feature_names = ["CRIM", "ZN", "INDUS", "CHAS", "NOX", "RM", "AGE", "DIS", "RAD", "TAX", "PTRATIO", "B", "LSTAT"]
    target_names = ["MEDV"]

    d = BostonHousing()
    test_inmemory_supervised_table_dataset(d;
        n_obs, n_features, n_targets,
        feature_names, target_names)

    d = BostonHousing(as_df=false)
    test_inmemory_supervised_table_dataset(d;
        n_obs, n_features, n_targets,
        feature_names, target_names,
        Tx=Float64, Ty=Float64)
end

@testset "Iris" begin
    n_obs = 150
    n_features = 4
    n_targets = 1

    feature_names = ["sepallength", "sepalwidth", "petallength", "petalwidth"]
    target_names = ["class"]

    d = Iris()
    test_inmemory_supervised_table_dataset(d;
        n_obs, n_features, n_targets,
        feature_names, target_names)

    d = Iris(as_df=false)
    test_inmemory_supervised_table_dataset(d;
        n_obs, n_features, n_targets,
        feature_names, target_names,
        Tx=Float64, Ty=AbstractString)
end

@testset "Mutagenesis" begin
    dtrain = Mutagenesis(split=:train)
    dtest = Mutagenesis(split=:test)
    dval = Mutagenesis(split=:val)
    dall = Mutagenesis(split=:all)

    train_x, train_y = dtrain[]
    test_x, test_y = dtest[]
    val_x, val_y = dval[]
    all_x, all_y = dall[]

    @test length(train_x) == length(train_y) == 100
    @test length(test_x) == length(test_y) == 44
    @test length(val_x) == length(val_y) == 44
    @test length(all_x) == length(all_y) == 188

    # test that label is not contained in features
    @test !any(haskey.(train_x, :mutagenic))
    @test !any(haskey.(test_x, :mutagenic))
    @test !any(haskey.(val_x, :mutagenic))
    # test data is materialized
    @test train_x isa Vector{<:Dict}
    @test test_x isa Vector{<:Dict}
    @test val_x isa Vector{<:Dict}

    x, y = dtrain[1:2]
    @test x isa Vector{Dict{Symbol, Any}}
    @test length(x) == 2
    @test y isa Vector{Int}
    @test length(y) == 2
end

@testset "Titanic" begin
    n_obs = 891
    n_features = 11
    n_targets = 1
    feature_names = ["PassengerId","Pclass", "Name", "Sex", "Age", "SibSp", "Parch", "Ticket", "Fare", "Cabin", "Embarked"]
    target_names = ["Survived"]

    d = Titanic()
    test_inmemory_supervised_table_dataset(d;
        n_obs, n_features, n_targets,
        feature_names, target_names)

    d = Titanic(as_df=false)
    test_inmemory_supervised_table_dataset(d;
        n_obs, n_features, n_targets,
        feature_names, target_names,
        Tx=Any, Ty=Int)

    @test isequal(d[1].features, [1, 3, "Braund, Mr. Owen Harris", "male", 22, 1, 0, "A/5 21171", 7.25, missing, "S"])
end


text:
export check_function, check_method
export StabilityReport, is_stable

"""
    check_function(func, signatures, acceptable_instability=Dict())

Check that the function is stable under each of the given signatures.

Return an array of method signature-`StabilityReport` pairs from
[`check_method`](@ref).
"""
function check_function(func, signatures, acceptable_instability=Dict{Symbol, Type}())
    result = Tuple{Any, StabilityReport}[]
    for params in signatures
        push!(result, (params, check_method(func, params, acceptable_instability)))
    end
    result
end

"""
    check_method(func, signature, acceptable_instability=Dict())

Create a `StabilityReport` object describing the type stability of the method.is

Compute non-concrete types of variables and return value, returning them in
a [`StabilityReport`](@ref) Object

`acceptable_instability`, if present, is a mapping of variables that are
allowed be non-concrete types.  `get` is called with the mapping, the
variable's symbol and `Bool` to get the variable's allowed type.  Additionally,
the return value is checked using `:return` as the symbol.
"""
function check_method(func, signature, acceptable_instability=Dict{Symbol, Type}())
    #Based off julia's code_warntype and show_ir

    unstable_vars_list = Array{Tuple{Symbol, Type}, 1}(undef, 0)

    function var_is_stable(typ, name)
        (isconcretetype(typ) && typ != Core.Box) ||
            begin
                (typ <: get(acceptable_instability, name, Bool))
            end
    end

    if VERSION < v"0.7.0-"
        code = code_typed(func, signature)

        function scan_exprs!(used, exprs)
            for ex in exprs
                if isa(ex, Slot)
                    used[ex.id] = true
                elseif isa(ex, Expr)
                    scan_exprs!(used, ex.args)
                end
            end
        end
    else
        code = code_typed(func, signature; optimize=false)
    end

    if length(code) == 0
        error("No methods found for $func matching $signature")
    elseif length(code) != 1
        warn("Mutliple methods for $func matching $signature")
    end

    for (src, rettyp) in code
        if !var_is_stable(rettyp, :return)
            push!(unstable_vars_list, (:return, rettyp))
        end
        slotnames = Base.sourceinfo_slotnames(src)
        used_slotids = falses(length(slotnames))

        if VERSION < v"0.7.0-"
            scan_exprs!(used_slotids, src.code)
            types = src.slottypes
        else
            stmts = src.code
            ssatypes = src.ssavaluetypes
            types = Vector{Type}(undef, length(slotnames))
            fill!(types, Union{})
            for idx in eachindex(stmts)
                if isa(stmts[idx], Expr) && stmts[idx].head == :(=)
                    if isa(ssatypes[idx], Core.Compiler.Const)
                        typ = typeof(ssatypes[idx].val)
                    else
                        typ = ssatypes[idx]
                    end
                    slotidx = stmts[idx].args[1].id
                    types[slotidx] = Union{types[slotidx], typ}
                    used_slotids[slotidx] = true
                end
            end
        end

        if isa(types, Vector)
            for i = 1:length(slotnames)
                if used_slotids[i]
                    name = Symbol(slotnames[i])
                    typ = types[i]
                    if !var_is_stable(typ, name)
                        push!(unstable_vars_list, (name, typ))
                    end

                    #else not an issue for type stability
                end
            end
        else
            warn("Can't access types of CodeInfo")
        end
    end

    return StabilityReport(unstable_vars_list)
end

"""
    StabilityReport()
    StabilityReport(unstable_variables::Vector{Tuple{Symbol, Type}})

Holds information about the stability of a method.

If `unstable_vars` is present, set the fields.  Otherwise, creates an empty set.

See [`is_stable`](@ref)
"""
struct StabilityReport
    "A list of unstable variables and their values"
    unstable_variables::Dict{Symbol, Type}

    StabilityReport(v::Dict{Symbol, Type}) = new(v)
end

StabilityReport() = StabilityReport(Dict{Symbol, Type}())
StabilityReport(vars) = StabilityReport(Dict{Symbol, Type}(vars))

function Base.:(==)(x::StabilityReport, y::StabilityReport)
    x.unstable_variables == y.unstable_variables
end

"""
    is_stable(report::StabilityReport)::Bool
    is_stable(reports::AbstractArray{StabilityReport})::Bool
    is_stable(reports::AbstractArray{Tuple{<:Any, StabilityReport}})::Bool

Check if the given [`StabilityReport`](@ref)s don't have any unstable types.
"""
is_stable(report::StabilityReport) = length(report.unstable_variables) == 0
is_stable(reports::AbstractArray{StabilityReport}) = all(is_stable.(reports))
is_stable(reports::Set{StabilityReport}) = all(is_stable.(reports))
is_stable(reports::AbstractArray{<:Tuple{<:Any, StabilityReport}}) = all(map(pair->is_stable(pair[2]), reports))

"""
    stability_warn(func_name, report::AbstractArray{Tuple{<:Any,StabilityReport}})

Displays warnings about the function if any of the reports are not stable
"""
function stability_warn(func_name, reports)
    for (args, report) in reports
        if !is_stable(report)
            println(stderr, "$func_name($(join(args, ", "))) is not stable")
            for (var, typ) in report.unstable_variables
                println(stderr, "  $var is of type $typ")
            end
        end
    end
end


text:
using ImageFiltering.Models
using ImageFiltering.Models: preallocate_solve_ROF_PD

@testset "solve_ROF_PD" begin
    # Note: random seed really matters a lot

    @testset "API" begin
        img = rand(MersenneTwister(0), Float64, 64)
        out = solve_ROF_PD(img, 0.2, 30)
        @test eltype(out) == Float32

        out2 = solve_ROF_PD(Float64, img, 0.2, 30)
        @test eltype(out2) == Float64
        @test out ≈ out2

        out3 = similar(img, Float32)
        buffer = preallocate_solve_ROF_PD(img)
        solve_ROF_PD!(out3, buffer, img, 0.2, 30)
        @test out == out3
    end

    @testset "Numerical" begin
        # 2D Gray
        img = restrict(testimage("cameraman"))
        img_noisy = img .+ 0.05randn(MersenneTwister(0), size(img))
        img_smoothed = solve_ROF_PD(img_noisy, 0.05, 20)
        @test ndims(img_smoothed) == 2
        @test eltype(img_smoothed) <: Gray
        @test assess_psnr(img_smoothed, img) > 31.67
        @test assess_ssim(img_smoothed, img) > 0.90

        # 2D RGB
        img = restrict(testimage("lighthouse"))
        img_noisy = img .+ colorview(RGB, ntuple(i->0.05.*randn(MersenneTwister(i), size(img)), 3)...)
        img_smoothed = solve_ROF_PD(img_noisy, 0.03, 20)
        @test ndims(img_smoothed) == 2
        @test eltype(img_smoothed) <: RGB
        @test assess_psnr(img_smoothed, img) > 32.15
        @test assess_ssim(img_smoothed, img) > 0.90

        # 3D Gray
        img = Gray.(restrict(testimage("mri"), (1, 2)))
        img_noisy = img .+ 0.05randn(MersenneTwister(0), size(img))
        img_smoothed = solve_ROF_PD(img_noisy, 0.02, 20)
        @test ndims(img_smoothed) == 3
        @test eltype(img_smoothed) <: Gray
        @test assess_psnr(img_smoothed, img) > 31.78
        @test assess_ssim(img_smoothed, img) > 0.85

        # 3D RGB
        img = RGB.(restrict(testimage("mri"), (1, 2)))
        img_noisy = img .+ colorview(RGB, ntuple(i->0.05.*randn(MersenneTwister(i), size(img)), 3)...)
        img_smoothed = solve_ROF_PD(img_noisy, 0.02, 20)
        @test ndims(img_smoothed) == 3
        @test eltype(img_smoothed) <: RGB
        @test assess_psnr(img_smoothed, img) > 31.17
        @test assess_ssim(img_smoothed, img) > 0.79
    end

    @testset "FixedPointNumbers" begin
        A = rand(N0f8, 20, 20)
        @test solve_ROF_PD(A, 0.01, 5) ≈ solve_ROF_PD(float32.(A), 0.01, 5)
    end

    @testset "OffsetArray" begin
        Ao = OffsetArray(rand(N0f8, 20, 20), -1, -1)
        out = solve_ROF_PD(Ao, 0.01, 5)
        @test axes(out) == axes(Ao)
    end
end


text:
if VERSION < v"0.4.0-"
    hash32(d::Any) = uint32(hash(d))
else
    hash32(d::Any) = hash(d) % Uint32
end

ρ(s::Uint32) = uint32(uint32(leading_zeros(s)) + 0x00000001)

function α(m::Uint32)
    if m == 0x00000010
        return 0.673
    elseif m == 0x00000020
        return 0.697
    elseif m == 0x00000040
        return 0.709
    else # if m >= uint32(128)
        return 0.7213 / (1 + 1.079 / m)
    end
end

type HyperLogLog
    m::Uint32
    M::Vector{Uint32}
    mask::Uint32
    altmask::Uint32
end

function HyperLogLog(b::Integer)
    if !(4 <= b <= 16)
        throw(ArgumentError("b must be an integer between 4 and 16"))
    end

    m = 0x00000001 << b

    M = zeros(Uint32, m)

    mask = 0x00000000
    for i in 1:(b - 1)
        mask |= 0x00000001
        mask <<= 1
    end
    mask |= 0x00000001

    altmask = ~mask

    return HyperLogLog(m, M, mask, altmask)
end

function Base.show(io::IO, counter::HyperLogLog)
    @printf(io, "A HyperLogLog counter w/ %d registers", int(counter.m))
    return
end

function update!(counter::HyperLogLog, v::Any)
    x = hash32(v)
    j = uint32((x & counter.mask) + 0x00000001)
    w = x & counter.altmask
    counter.M[j] = max(counter.M[j], ρ(w))
    return
end

function state(counter::HyperLogLog)
    S = 0.0

    for j in 1:counter.m
        S += 1 / (2^counter.M[j])
    end

    Z = 1 / S

    E = α(counter.m) * uint(counter.m)^2 * Z

    if E <= 5//2 * counter.m
        V = 0
        for j in 1:counter.m
            V += int(counter.M[j] == 0x00000000)
        end
        if V != 0
            E_star = counter.m * log(counter.m / V)
        else
            E_star = E
        end
    elseif E <= 1//30 * 2^32
        E_star = E
    else
        E_star = -2^32 * log(1 - E / (2^32))
    end

    return E_star
end

# TODO: Figure out details here
# function confint(counter::HyperLogLog)
#     e = estimate(counter)
#     delta = e * 1.04 / sqrt(counter.m)
#     return e - delta, e + delta
# end

function Base.copy(counter::HyperLogLog)
    return HyperLogLog(
        counter.m,
        copy(counter.M),
        counter.mask,
        counter.altmask,
    )
end


text:
using DelayEmbeddings: localextrema, exponential_decay_fit
using DynamicalSystemsBase, Test

println("\nTesting delay estimation...")

testval = (val, vmin, vmax) -> @test vmin ≤ val ≤ vmax
diffeq = (atol = 1e-9, rtol = 1e-9, maxiters = typemax(Int))

@testset "Estimate Delay" begin
    # test exponential decay fit
    x = 0:200
    y = @. exp(-x/5)
    @test exponential_decay_fit(x, y, :equal) ≈ 5
    @test exponential_decay_fit(x, y, :small) ≈ 5

    ds = Systems.henon()
    data = trajectory(ds,1000;Δt=1)
    x = data[:,1]
    @test estimate_delay(x,"ac_zero", 0:10) ≤ 2
    @test estimate_delay(x,"ac_min", 0:10)  ≤ 2
    @test estimate_delay(x,"exp_extrema", 0:10)  ≤ 4
    @test estimate_delay(x,"mi_min", 0:10) == 10

    ds = Systems.roessler(ones(3))
    Δt = 0.02
    data = trajectory(ds,500.0;Δt,diffeq...)
    x = data[:,1]
    @test 1.3 ≤ estimate_delay(x,"ac_zero", 1:2:500)*Δt ≤ 1.7
    @test 2.6 ≤ estimate_delay(x,"ac_min", 1:2:500)*Δt  ≤ 3.4
    @test 1.0 ≤ estimate_delay(x,"mi_min", 1:2:500)*Δt ≤ 1.6

    Δt = 0.1
    data = trajectory(ds,2000.0; Δt,diffeq...)
    x = data[:,1]
    @test 1.3 ≤ estimate_delay(x,"ac_zero", 1:1:50)*Δt ≤ 1.7
    @test 2.6 ≤ estimate_delay(x,"ac_min", 1:1:50)*Δt  ≤ 3.4
    @test 1.15 ≤ estimate_delay(x,"mi_min", 1:1:50)*Δt ≤ 1.6

    ds = Systems.lorenz()
    Δt = 0.1
    data = trajectory(ds,5000;Δt,diffeq...)
    x = data[500:end,1]
    @test 0 < estimate_delay(x,"exp_extrema", 0:2:200)  < 200

    # Issue #18
    ds = Systems.gissinger(ones(3)) # 3D continuous chaotic system, also shown in orbit diagrams tutorial
    Δt = 0.05
    data = trajectory(ds, 1000.0, Δt = Δt)
    s = data[:, 1]

    τ = estimate_delay(s, "mi_min", 0:1:400) # this was the erroring line
    @test τ > 0
end


text:
#SeqBounds.jl
struct SeqBoundsResult
    sf::String
    v::Vector
    alpha::Float64
    side::Int
    za::Vector
    zb::Vector
    asfv::Vector
    asfvdiff::Vector
    stdv::Vector
    sdproc::Vector
end
function obf(p, alpha, side)
    2(1 - cdf(Normal(), quantile(Normal(), 1 - alpha/side/2)/sqrt(p)))
end
function pocock(p, alpha, side)
    (alpha/side)*log(1+(exp(1)-1)*p)
end
function power(p, alpha, side; γ = 1)
    alpha/side*p^γ
end

"""
    bounds(v::Vector, alpha::Float64; side = :one, asf = :obf, h::Float64 = 0.05)

* `v` - vector of onformation portion;
* `alpha` - total alpha level;
* `side` - one- or two- sided test (:one or :two);
* `asf` - alpha spending function (:obf, :pocock, :power, :ep);
* `h` - integration grid multiplier.

# Example

bounds([0.2, 0.4, 0.6, 0.8, 1.0], 0.05; side = :one, asf = :obf)

Maximum `v` shoul be less 1, i-th element shod be less i+1.
Alpha shoul be in range > 0.0 and < 1.0.

# asf:

`:obf` - O'Brien-Fleming function

2(1 - cdf(Normal(), quantile(Normal(), 1 - α / side / 2)/sqrt(t)))

`:pocock` - Pocock function.

α / side * log(1 + (exp(1) - 1) * t)

`:power` - power function.

α / side * t^γ

`:ep` - equal parts.

"""
function bounds(v::Vector, alpha::Float64; side = :one, asf = :obf, h::Float64 = 0.05)
    if alpha ≥ 1 || alpha ≤ 0 error("alpha should be in range: 0 < alpha < 1!") end
    if v[end] > 1 error("Last value of v shoul be  ≤ 1!") end
    if length(v) > 1
        for i = 2:length(v)
            if !(v[i] > v[i-1])  error("v[i] shoul be > v[i-1]") end
        end
    end
    if side == :one side = 1 elseif side == :two side  = 2 end
    if asf == :obf
        asfv   = obf.(v, alpha, side) * side
        asfstr = "O'Brien-Fleming"
    elseif asf == :pocock
        asfv   = pocock.(v, alpha, side) * side
        asfstr = "Pocock"
    elseif asf == :power
        asfv   = power.(v, alpha, side) * side
        asfstr = "Power"
    elseif asf == :ep
        asfv   = collect(1:length(v)) * (alpha/length(v))
        asfstr = "Equal parts"
    else
        error("uncknown function")
    end
    #
    zninf       = -8.0
    asfvdiff    = similar(asfv)
    vdiff       = similar(v)
    asfvdiff[1] = asfv[1]
    vdiff[1]    = v[1]
    for i = 2:length(v)
        vdiff[i]    = v[i] - v[i-1]
        asfvdiff[i] = asfv[i] - asfv[i-1]
    end
    stdv     = sqrt.(vdiff)
    sdproc   = sqrt.(v)
    za    = similar(asfv)
    zb    = similar(asfv)
    ya    = similar(asfv)
    yb    = similar(asfv)
    nints = zeros(Int, length(v))
    grids = Vector{StepRangeLen}(undef, length(v))
    #
    zb[1] = quantile(Normal(), 1-asfvdiff[1]/side)
    yb[1] = zb[1]*stdv[1]
    if side == 1 # For one side
        za[1] = zninf
        ya[1] = za[1]*stdv[1]
    else
        za[1] = -zb[1]
        ya[1] = -yb[1]
    end
    nints[1] = ceil((yb[1] - ya[1])/(h*stdv[1]))
    grids[1] = range(ya[1], yb[1], length=nints[1] + 1)
    last     = pdf.(Normal(0, stdv[1]), grids[1])
    for i = 2:length(v)
        if i > 2 # For next step
            nints[i-1] = ceil((yb[i-1]-ya[i-1])/(h*stdv[i-1]))
            grids[i-1] = range(ya[i-1], yb[i-1], length=nints[i-1]+1)
            last     = highordgrid(last, grids[i-1], grids[i-2], stdv[i-1])
        end
        zerof    = x -> asfvdiff[i]/side - integrategrid(x, last, grids[i-1], stdv[i])
        yb[i]    = find_zero(zerof, (0.0, zb[i-1]))
        zb[i]    = yb[i]/sdproc[i]
        if side == 1
            za[i]    = zninf
            ya[i]    = zninf*sdproc[i]
        else
            ya[i]    = -yb[i]
            za[i]    = -zb[i]
        end
    end
    SeqBoundsResult(asfstr, v, alpha, side, za, zb, asfv, asfvdiff, stdv, sdproc)
end

function integrategrid(x, last::Vector{T}, grid, stdv) where T
    hlast = step(grid)
    DIST  = Normal(x, stdv)
    fst   = cdf(DIST, grid[1]) * last[1]
    lst   = cdf(DIST, grid[end]) * last[end]
    sfn   = zero(T)
    for i = 2:length(grid)-1
        sfn += cdf(DIST, grid[i]) * last[i]
    end
    hlast*(sfn + fst/2 + lst/2)
end

function highordgrid(last::Vector{T}, x, grid, stdv) where T
    h    = step(grid)
    rvec = zeros(T, length(x))
    for n = 1:length(x)
        fst = last[1] * pdf(Normal(x[n], stdv), grid[1])
        lst = last[end] * pdf(Normal(x[n], stdv), grid[end])
        for m = 2:length(grid)-1
            rvec[n] += last[m] * pdf(Normal(x[n], stdv), grid[m])
        end
        rvec[n] += fst/2 + lst/2
    end
    rvec * h
end

function Base.show(io::IO, obj::SeqBoundsResult)
    if obj.side == 1
        println(io, """One-sided group sequential design
    Alpha spending function: $(obj.sf),  Alpha = $(obj.alpha)""")
            pretty_table(io, hcat(obj.v, obj.asfv, obj.asfvdiff, obj.zb, ccdf.(Normal(), obj.zb)); header = ["Portion", "Function value", "Spend", "Z", "Nominal p"])
    else
        println(io, """Two-sided group sequential design
    Alpha spending function: $(obj.sf),  Alpha = $(obj.alpha)""")
            pretty_table(io, hcat(obj.v, obj.asfv, obj.asfvdiff, obj.za, obj.zb, ccdf.(Normal(), obj.zb)); header = ["Portion", "Function value", "Spend", "Zl", "Zu", "Nominal p"])
    end
end
################################################################################
#=
Dounds can be found by hcubature with HCubature.jl but it takes more time and memory
=#
#=
f1i(x) = quadgk(ab -> pdf(Normal(0, stdv[1]), ab), ya[1], yb[1], rtol=1e-8)[1]
f1i(yb[1])
f1(x) =  pdf(Normal(0, stdv[1]), x)

f21(x) = quadgk(u -> cdf(Normal(x, stdv[2]), u)*pdf(Normal(0, stdv[1]), u), ya[1], yb[1], rtol=1e-8)[1]
f21(yb[2])

# OR

f22(x) = quadgk(u -> cdf(Normal(x, stdv[2]), u)*f1(u), ya[1], yb[1], rtol=1e-8)[1]
zerof = x ->asfvdiff[2] - f23(x)
fz = find_zero(zerof, (0.0, zb[1]))
zbf = fz/sdproc[2]

f3(x) = hcubature(u -> cdf(Normal(x, stdv[i]), u[1]) * pdf(Normal(u[1], stdv[i-1]), u[2]) * pdf(Normal(0, stdv[1]), u[2]) , [ya[i-1], ya[i-2]], [yb[i-1],yb[i-2]], rtol=1e-8)[1]
f3(yb[3])
zerof = x ->asfvdiff[3] - f3(x)
fz = find_zero(zerof, (0.0, zb[2]))
zbf = fz/sdproc[3]
=#


text:
"Linear two point interpolation"
function interp(x::Real,
        xmin::Real=0, xmax::Real=1,
        ymin::Real=0, ymax::Real=1)
    (x - xmin) * (ymax - ymin) / (xmax - xmin) + ymin
end

"Constant function"
function cst(k::Number)
    return f(args...) = k
end

"Identity function"
id(x::Any) = x


text:
pkg_path = dirname(dirname(pathof(PdfFraud)))
data_path = joinpath(pkg_path, "data")

# first make sure that the proto source files are in place
plugin_path = joinpath(dirname(pathof(ProtoBuf)), "../plugin/protoc-gen-julia")
proto_path = joinpath(data_path, "proto")

"""
	build_proto_definition()

Builds the correct ProtoBuffer Julia definitions.
"""
function build_proto_definition()
	println("Building the protobuffer definitions.")
	cmd = `-I=$(data_path) --plugin=protoc-gen-julia=$(plugin_path) --julia_out=$(proto_path) $(data_path)/sample.proto`
	run(ProtoBuf.protoc(cmd))
	println("Done.")
end

mkpath(proto_path)
if length(readdir(proto_path)) == 0
	build_proto_definition()
end

include(joinpath(proto_path, "google.jl"))
include(joinpath(proto_path, "sample_pb.jl"))

"""
	read_pb(fs::String)

Reads the whole ProtoBuffer data object from file.
"""
function read_pb(fs::String)
	open(fs) do f
		readproto(f, PdfSampleRaw())
	end
end

_hostname = gethostname()
proto_path = ""
extracted_data = ""
if _hostname == "vit-ThinkPad-E470"
	proto_path = "/home/vit/vyzkum/bulletproof/pdf_experiment/data/uloz_to_2019_11_05/"
elseif	_hostname in ["gpu-node", "gpu-titan"]
	proto_path = "/compass/home/skvara/vyzkum/bulletproof/pdf_experiment/data/uloz_to_2019_11_05/"
elseif	_hostname == "tarbik.utia.cas.cz"
	proto_path = "/home/skvara/work/bulletproof/pdf_experiment/data/uloz_to_2019_11_05/"
end
if proto_path != ""
	extracted_data = joinpath(proto_path, "extracted_representations/data.jld2")
else
	@warn "ProtoBuffer data path is not known on this system!"
end

"""
	visual_representation(pdf::PdfSampleRaw[, pageid])

Extracts the raw representation from the PdfSampleRaw object `pdf`.
"""
visual_representation(pdf::PdfSampleRaw, pageid::Int) = pdf.pages[pageid].visual_representation.letter_representation
function visual_representation(pdf::PdfSampleRaw)
	x = []
	for i in 1:length(pdf.pages)
		# we have to do it this way because some pages have empty representations
		y = visual_representation(pdf, i)
		if length(y) != 0
			push!(x, y)
		end
	end
	hcat(x...)
end

"""
	reshape_representation(vector, shape=(48,64))
	reshape_representation(matrix, shape=(48,64))

Given a vector/matrix of floats, transforms this into the correct shape and orientation as a 3D tensor of
pdf representations.
"""
reshape_representation(vec::Vector, shape=(48,64)) = rotl90(reshape(vec, shape...))
reshape_representation(mat::Matrix, shape=(48,64)) =
	cat([reshape_representation(mat[:,i], shape) for i in 1:size(mat,2)]..., dims = 3)

"""
	get_labels([sha256::Vector, ]index_file::String, label_file::String)

Get labels from 2 files. Possibly filter them with a given vector of sha256 strings.
Returns a tuple of (labels, sha256, filenames).
"""
function get_labels(index_file::String, label_file::String)
	# get the index
	f = open(index_file);
	index_entries = map(JSON.parse, readlines(f));
	close(f)
	sha256 = map(x -> x["sha256"], index_entries);
	filenames = map(x -> x["original_name"], index_entries);

	# get the labels
	labeldata = CSV.read(label_file, header=[:label, :file]);

	# connect labels via filenames
	function filename2label(filename)
		label = (filter(row -> row[:file] == filename, labeldata))[!,:label]
		return (length(label) > 0) ? label[1] : nothing
	end

	map(filename2label, filenames), sha256, filenames
end
function get_labels(sha256::Vector, index_file::String, label_file::String)
	# get the connected labelels and shas from the two files
	labels, labeled_sha256, labeled_filenames = get_labels(index_file, label_file)

	# connect labels via sha256
	function sha2label(sha)
		label = labels[labeled_sha256 .== sha]
		return length(label) > 0 ? label[1] : nothing
	end
	function sha2filename(sha)
		filename = labeled_filenames[labeled_sha256 .== sha]
		return length(filename) > 0 ? filename[1] : nothing
	end

	map(sha2label, sha256), sha256, map(sha2filename, sha256)
end


text:
module Elements
import GaussQuadrature

"""
    lglpoints(::Type{T}, N::Integer) where T <: AbstractFloat

returns the points `r` and weights `w` associated with the `N+1`-point
Gauss-Legendre-Lobatto quadrature rule of type `T`

"""
function lglpoints(::Type{T}, N::Integer) where T <: AbstractFloat
  @assert N ≥ 1
  GaussQuadrature.legendre(T, N+1, GaussQuadrature.both)
end

"""
    lgpoints(::Type{T}, N::Integer) where T <: AbstractFloat

returns the points `r` and weights `w` associated with the `N+1`-point
Gauss-Legendre quadrature rule of type `T`
"""
function lgpoints(::Type{T}, N::Integer) where T <: AbstractFloat
  @assert N ≥ 1
  GaussQuadrature.legendre(T, N+1, GaussQuadrature.neither)
end

"""
    baryweights(r)

returns the barycentric weights associated with the array of points `r`

Reference:
  Jean-Paul Berrut & Lloyd N. Trefethen, "Barycentric Lagrange Interpolation",
  SIAM Review 46 (2004), pp. 501-517.
  <https://doi.org/10.1137/S0036144502417715>
"""
function baryweights(r::AbstractVector{T}) where T
  Np = length(r)
  wb = ones(T, Np)

  for j = 1:Np
    for i = 1:Np
      if i != j
        wb[j] = wb[j] * (r[j] - r[i])
      end
    end
    wb[j] = T(1) / wb[j]
  end
  wb
end

"""
    spectralderivative(r::AbstractVector{T},
                       wb=baryweights(r)::AbstractVector{T}) where T

returns the spectral differentiation matrix for a polynomial defined on the
points `r` with associated barycentric weights `wb`

Reference:
  Jean-Paul Berrut & Lloyd N. Trefethen, "Barycentric Lagrange Interpolation",
  SIAM Review 46 (2004), pp. 501-517.
  <https://doi.org/10.1137/S0036144502417715>
"""
function spectralderivative(r::AbstractVector{T},
                            wb=baryweights(r)::AbstractVector{T}) where T
  Np = length(r)
  @assert Np == length(wb)
  D = zeros(T, Np, Np)

  for k = 1:Np
    for j = 1:Np
      if k == j
        for l = 1:Np
          if l != k
            D[j, k] = D[j, k] + T(1) / (r[k] - r[l])
          end
        end
      else
        D[j, k] = (wb[k] / wb[j]) / (r[j] - r[k])
      end
    end
  end
  D
end

"""
    interpolationmatrix(rsrc::AbstractVector{T}, rdst::AbstractVector{T},
                        wbsrc=baryweights(rsrc)::AbstractVector{T}) where T

returns the polynomial interpolation matrix for interpolating between the points
`rsrc` (with associated barycentric weights `wbsrc`) and `rdst`

Reference:
  Jean-Paul Berrut & Lloyd N. Trefethen, "Barycentric Lagrange Interpolation",
  SIAM Review 46 (2004), pp. 501-517.
  <https://doi.org/10.1137/S0036144502417715>
"""
function interpolationmatrix(rsrc::AbstractVector{T}, rdst::AbstractVector{T},
                             wbsrc=baryweights(rsrc)::AbstractVector{T}) where T
  Npdst = length(rdst)
  Npsrc = length(rsrc)
  @assert Npsrc == length(wbsrc)
  I = zeros(T, Npdst, Npsrc)
  for k = 1:Npdst
    for j = 1:Npsrc
      I[k, j] = wbsrc[j] / (rdst[k] - rsrc[j]);
      if !isfinite(I[k,j])
        I[k, :] .= T(0);
        I[k, j] = T(1);
        break
      end
    end
    d = sum(I[k, :]);
    I[k, :] = I[k, :] / d;
  end
  I
end

end # module

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using ITensors
using ITensorMakie
using Graphs
using GLMakie
using LayeredLayouts
using PastaQ: qft

N = 4
gates = qft(N)

s = siteinds("Qubit", N)

U, s̃ = circuit_network(gates, s)
ψ = MPS(s)
ψ̃ = MPS(s̃)
tn = [ψ..., U..., ψ̃...]

edge_labels = (; tags=true, plevs=true)
layout(g) = layered_layout(solve_positions(Zarate(), g))
@visualize fig tn arrow_show = true edge_labels = edge_labels edge_textsize = 20 layout =
  layout
edge_labels = (; plevs=true)
@visualize! fig[1, 2] tn ndims = 3 edge_labels = edge_labels edge_textsize = 20

fig


text:
using Pushover
using Pushover: _sanitize_priority, _crop
using Test
#using Lint  # unsupported Julia 0.7 version as of 2018/08/19

@testset "Pushover" begin
    @testset "_sanitize_priority" begin
        @test _sanitize_priority(-2) == -2
        @test _sanitize_priority(-1) == -1
        @test _sanitize_priority(0) == 0
        @test _sanitize_priority(1) == 1

        @test _sanitize_priority(100) == 0
        @test _sanitize_priority(-100) == 0
        @test _sanitize_priority("abc") == 0
    end

    @testset "_crop" begin
        max_len = 30
        msg = "Lorem ipsum dolor sit amet"
        cropped_msg = _crop(msg, max_len)
        @test cropped_msg == msg
        @test length(cropped_msg) == length(msg)
        @test length(cropped_msg) <= max_len

        msg = "Lorem ipsum dolor sit amet, consectetur adipiscing elit."
        cropped_msg = _crop(msg, max_len)
        @test length(cropped_msg) == max_len
        @test cropped_msg[end-2:end] == "..."
    end

    #@testset "lint" begin
    #    @test isempty(lintpkg("Pushover"))
    #end
end


text:
export t₁, current_epoch, loss_history

# const t₁ = 5.0                      # number of simulation years
const ρ = 900.0                     # Ice density [kg / m^3]
const g = 9.81                    # Gravitational acceleration [m / s^2]
const n = 3.0                       # Glen's flow law exponent
const ϵ = 1e-3                    # small number
const maxA = 8e-17
const minA = 8.5e-20
const maxT = 1.0
const minT = -25.0

# From Cuffey and Paterson
const A_values_sec = [0.0 -2.0 -5.0 -10.0 -15.0 -20.0 -25.0 -30.0 -35.0 -40.0 -45.0 -50.0;
                              2.4e-24 1.7e-24 9.3e-25 3.5e-25 2.1e-25 1.2e-25 6.8e-26 3.7e-26 2.0e-26 1.0e-26 5.2e-27 2.6e-27] # s⁻¹Pa⁻³
const A_values = hcat(A_values_sec[1,:], A_values_sec[2,:].*60*60*24*365.25)'

plots = false                      # Display plots
overwrite_climate = false          # Force re-computing climate data for glaciers
# loss_type = "H"                   # Loss function based on ice thickness
loss_type = "V"                   # Loss function based on ice surface velocities
# loss_type = "HV"                    # Combined loss function based on ice surface velocities and ice thickness
const random_sampling_loss = false  # Use random subset of matrix samples for the loss
const scale_loss = true
const noise = true                  # Add random noise to fake A law
rng_seed() = MersenneTwister(123)   # Random seed

# Machine learning training
current_epoch = 1
loss_history = []


text:
#= This file is being deprecated in favor of using diagnostics.jl =#

initialize_io(self::ForcingBase, Stats) = nothing

initialize_io(self::RadiationBase, Stats::NetCDFIO_Stats) = nothing
io(self::RadiationBase, grid, state, Stats::NetCDFIO_Stats) = nothing

function initialize_io(en::EnvironmentVariables, Stats::NetCDFIO_Stats)
    add_ts(Stats, "env_cloud_base")
    add_ts(Stats, "env_cloud_top")
    add_ts(Stats, "env_cloud_cover")
    add_ts(Stats, "env_lwp")
    add_ts(Stats, "env_iwp")
    return
end

function io(en::EnvironmentVariables, grid, state, Stats::NetCDFIO_Stats)
    # Assuming amximum overlap in environmental clouds
    write_ts(Stats, "env_cloud_cover", en.cloud_cover)
    write_ts(Stats, "env_cloud_base", en.cloud_base)
    write_ts(Stats, "env_cloud_top", en.cloud_top)
    write_ts(Stats, "env_lwp", en.lwp)
    write_ts(Stats, "env_iwp", en.iwp)
    return
end

function initialize_io(precip::PrecipVariables, Stats::NetCDFIO_Stats)
    add_ts(Stats, "rwp_mean")
    add_ts(Stats, "swp_mean")
    add_ts(Stats, "cutoff_precipitation_rate")
    return
end

function io(precip::PrecipVariables, grid, state, Stats::NetCDFIO_Stats)
    write_ts(Stats, "rwp_mean", precip.mean_rwp)
    write_ts(Stats, "swp_mean", precip.mean_rwp)

    #TODO - change to rain rate that depends on rain model choice
    write_ts(Stats, "cutoff_precipitation_rate", precip.cutoff_precipitation_rate)
    return
end

function initialize_io(up::UpdraftVariables, Stats::NetCDFIO_Stats)

    add_ts(Stats, "updraft_cloud_cover")
    add_ts(Stats, "updraft_cloud_base")
    add_ts(Stats, "updraft_cloud_top")
    add_ts(Stats, "updraft_lwp")
    add_ts(Stats, "updraft_iwp")
    return
end

function io(up::UpdraftVariables, grid, state, Stats::NetCDFIO_Stats)
    # Note definition of cloud cover : each updraft is associated with a cloud cover equal to the maximum
    # area fraction of the updraft where ql > 0. Each updraft is assumed to have maximum overlap with respect to
    # itup (i.e. no consideration of tilting due to shear) while the updraft classes are assumed to have no overlap
    # at all. Thus total updraft cover is the sum of each updraft"s cover
    write_ts(Stats, "updraft_cloud_cover", sum(up.cloud_cover))
    write_ts(Stats, "updraft_cloud_base", minimum(abs.(up.cloud_base)))
    write_ts(Stats, "updraft_cloud_top", maximum(abs.(up.cloud_top)))
    write_ts(Stats, "updraft_lwp", up.lwp)
    write_ts(Stats, "updraft_iwp", up.iwp)
    return
end

function initialize_io(gm::GridMeanVariables, Stats::NetCDFIO_Stats)
    add_ts(Stats, "lwp_mean")
    add_ts(Stats, "iwp_mean")
    add_ts(Stats, "cloud_base_mean")
    add_ts(Stats, "cloud_top_mean")
    add_ts(Stats, "cloud_cover_mean")
    return
end

function io(gm::GridMeanVariables, grid, state, Stats::NetCDFIO_Stats)
    write_ts(Stats, "cloud_cover_mean", gm.cloud_cover)

    write_ts(Stats, "lwp_mean", gm.lwp)
    write_ts(Stats, "iwp_mean", gm.iwp)
    write_ts(Stats, "cloud_base_mean", gm.cloud_base)
    write_ts(Stats, "cloud_top_mean", gm.cloud_top)
    return
end

# Initialize the IO pertaining to this class
function initialize_io(edmf::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats)

    initialize_io(edmf.UpdVar, Stats)
    initialize_io(edmf.EnvVar, Stats)
    initialize_io(edmf.Precip, Stats)

    add_ts(Stats, "rd")
    return
end

function io(edmf::EDMF_PrognosticTKE, grid, state, Stats::NetCDFIO_Stats, TS::TimeStepping, param_set)
    io(edmf.UpdVar, grid, state, Stats)
    io(edmf.EnvVar, grid, state, Stats)
    io(edmf.Precip, grid, state, Stats)
    write_ts(Stats, "rd", StatsBase.mean(edmf.pressure_plume_spacing))
    return
end


text:
module Julip

include("jlptypes.jl")
include("util.jl")
include("token.jl")
include("asmtypes.jl")
include("llvmirgen.jl")

function main(args = [])
    user_input = if length(args) != 0 && isfile(args[1])
                     str = open( args[1], "r" ) do fp
                         read(fp, String)
                     end
                 else
                     read(stdin, String)
                 end
    corecode = open(joinpath(dirname(pathof(Julip)), "..", "julipsrc/core.jlp")) do fp
                   read(fp, String)
               end
    tokens = tokenize(corecode * user_input)
    println(setall(tokens))
end

function main(arg::String)
    main([arg])
end

#main(ARGS)

end # module


text:
using CITIfile: rdciti
using Test

s,f = rdciti("seg_list.cit")
@testset "SEG_LIST-style file" begin
    @test size(s) == (1,1,10)
    @test size(f) == (10,)
    @test f == range(1, 4, length=10)
    @test s[1,1,1] ≈ complex(0.86303E-1,-8.98651E-1)
    @test s[1,1,end] ≈ complex(-7.78350E-1,5.72082E-1)
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

name = "SCSBuilder"
version = v"2.1.1"

# Collection of sources required to build SCSBuilder
sources = [
    "https://github.com/cvxgrp/scs/archive/$(version).tar.gz" =>
    "0e20b91e8caf744b84aa985ba4e98cc7235ee33612b2bad2bf31ea5ad4e07d93",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/scs-*
flags="DLONG=1 USE_OPENMP=0"
blasldflags="-L${prefix}/lib"
# see https://github.com/JuliaPackaging/Yggdrasil/blob/0bc1abd56fa176e3d2cc2e48e7bf85a26c948c40/OpenBLAS/build_tarballs.jl#L23
if [[ ${nbits} == 64 ]] && [[ ${target} != aarch64* ]]; then
    flags="${flags} BLAS64=1 BLASSUFFIX=_64_"
    blasldflags+=" -lopenblas64_"
else
    blasldflags+=" -lopenblas"
fi

# Patch to reverse this WIN64 change: https://github.com/cvxgrp/scs/commit/9858d6b562f499de75493b85286276c19ad84c6f#diff-a9dbab3214616022c64ee2656440f544
# Looks like it is known that this change causes trouble witn mingw32 (but not for mingw64?):
# https://github.com/cvxgrp/scs/blob/e6ab81db115bb37502de0a9917041a0bc2ded313/.appveyor.yml#L13-L16
cd include
cp glbopts.h glbopts.h.orig
cat > file.patch <<'END'
--- glbopts.h.orig
+++ glbopts.h
@@ -97,7 +97,7 @@
 #ifdef _WIN64
 /* #include <stdint.h> */
 /* typedef int64_t scs_int; */
-typedef long scs_int;
+typedef __int64 scs_int;
 #else
 typedef long scs_int;
 #endif
END
patch -l glbopts.h.orig file.patch -o glbopts.h
cd ..

make BLASLDFLAGS="${blasldflags}" ${flags} out/libscsdir.${dlext}
make BLASLDFLAGS="${blasldflags}" ${flags} out/libscsindir.${dlext}
mv out/libscs* ${prefix}/lib/

# Copied from https://github.com/JuliaLinearAlgebra/ArpackBuilder/blob/ed55085cd647b39fefb923c3dd95b03591522760/build_tarballs.jl#L72-L77
# Also see https://github.com/JuliaLang/julia/issues/27610
# For now, we'll have to adjust the name of the OpenBLAS library on macOS.
# Eventually, this should be fixed upstream
if [[ ${target} == "x86_64-apple-darwin14" ]]; then
    echo "-- Modifying library name for OpenBLAS"
    install_name_tool -change @rpath/libopenblas64_.0.3.5.dylib @rpath/libopenblas64_.dylib ${prefix}/lib/libscsindir.dylib
    install_name_tool -change @rpath/libopenblas64_.0.3.5.dylib @rpath/libopenblas64_.dylib ${prefix}/lib/libscsdir.dylib
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    FreeBSD(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]
platforms = expand_gcc_versions(platforms)

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libscsindir", :indirect),
    LibraryProduct(prefix, "libscsdir", :direct)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/OpenBLAS-v0.3.5-2/build_OpenBLAS.v0.3.5.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using Polynomials
set_bigfloat_precision(128)
function Baristowa(n,r,s, a,  iter )
	b = zeros(typeof(r), n+1)
	c = zeros(typeof(r), n+1)
	for j = 1:iter
		b[n+1] = a[n+1]
		b[n-1+1] = a[n-1+1]+ r*b[n+1]
		c[n+1] = b[n+1]
		c[n-1+1]= b[n-1+1] + r*c[n+1]
		for i = n-2:-1:0
			b[i+1]= a[i+1] + r*b[i+1+1] + s*b[i+2+1]
			c[i+1] = b[i+1] + r*c[i+1+1] + s*c[i+2+1]
	end
	J = [	c[2+1]	c[3+1];
			c[1+1]	c[2+1]];
	B = [-b[1+1]; -b[0+1]]
	res = \(J,B)
	r += res[1]
	s += res[2]
	end
	return (r,s)
end

(r,s) = Baristowa(4,BigFloat(0.1),BigFloat(0.1),[BigFloat(1.):BigFloat(5.);], 10)

@printf("z^2  %.16f*z - %.16f \n",r,s)
@show ((r+sqrt(complex(r*r+BigFloat(4.)*s)))/BigFloat(2.))

roots(Poly([1.:5.;]))

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
function show(io::IO, x::DoubleFloat{T}) where {T<:IEEEFloat}
    compact = get(io, :compact, true)
    if compact
        str = string(x.hi)
    else
        str = string(x)
    end
    print(io, str)
end

function show(io::IO, x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    compact = get(io, :compact, true)
    if compact
        str = string(x.re.hi, (signbit(x.im.hi) ? " - " : " + "), abs(x.im.hi), "im")
    else
        str = string(x.re, (signbit(x.im.hi) ? " - " : " + "), abs(x.im), "im")
    end
    print(io, str)
end

show(x::DoubleFloat{T}) where {T<:IEEEFloat} = show(Base.stdout, x)

show(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat} = show(Base.stdout, x)

function showtyped(io::IO, x::DoubleFloat{T}) where {T<:IEEEFloat}
    str = stringtyped(x)
    print(io, str)
end

showtyped(x::DoubleFloat{T}) where {T<:IEEEFloat} = showtyped(Base.stdout, x)

showtyped(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat} = showtyped(Base.stdout, x)


text:
function permute_labels!(membership::Vector{Int})
    N = length(membership)
    if maximum(membership) > N || minimum(membership) < 1
        error("Label must between 1 and |V|")
    end
    label_counters = zeros(Int, N)
    j = 1
    for i=1:length(membership)
        k = membership[i]
        if k >= 1
            if label_counters[k] == 0
                # We have seen this label for the first time
                label_counters[k] = j
                k = j
                j += 1
            else
                k = label_counters[k]
            end
        end
        membership[i] = k
    end
    membership
end

function avedegree{V}(g::AbstractGraph{V})
    d = 0.0
    for v in vertices(g)
        d += out_degree(v, g)
    end
    d/num_vertices(g)
end

max_stars_repo_name: afternone/LPA

max_issues_repo_name: afternone/LPA

max_forks_repo_name: afternone/LPA


text:
module PhylogeneticDistance

import StatsBase
using StatsBase
import RCall
using RCall

export
       getLabel,
       getEdgeLength,
       readNewick,
       isRooted,
       GUniFrac,
       KernelMatrix

include("getLabel.jl")
include("getEdgeLength.jl")
include("readNewick.jl")
include("isRooted.jl")
include("GUniFrac.jl")
include("KernelMatrix.jl")
end


text:
"Some utilities"
module Util

import Base: convert

"`BigInt` -> `Vector{UInt8}` conversion"
function convert(::Type{Vector{UInt8}}, x::BigInt)
    sz = div(ndigits(x, base=2) + 8 - 1, 8)
    bytes = UInt8[]
    sizehint!(bytes, sz)
    while x != 0
        push!(bytes, x & 255)
        x >>= 8
    end
    reverse!(bytes)
end

"`Vector{UInt8}` -> `BigInt` conversion"
function convert(::Type{BigInt}, v::Vector{UInt8})
    foldl((acc, b) -> (acc << 8) + b, v; init=big(0))
end

end  # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
"""
armijo(f,fk,dfk,xk,pk)
Backtracked Armijo linesearch
"""
function armijo(f::Function,fk,dfk,xk,pk;t=1.0, maxIter=10, c1=1e-4,b=0.5)
    LS = 1

    while LS<=maxIter
        if f(xk+t*pk)[1] <= fk + t*c1*dot(dfk,pk)
            break
        end
        t *= b
        LS += 1
    end
    if LS>maxIter
    	LS= -1
    	t = 0.0
    end
    return t,LS
end

"""
bfgs(f,df,x)
BFGS method for solving min_x f(x)
"""
function bfgs(f::Function,fdf::Function,x::Vector;H=Matrix(1.0I,length(x),length(x)), maxIter=20,atol=1e-8,out::Int=0,storeInterm::Bool=false,
	lineSearch::Function=(f,fk,dfk,xk,pk,ak)->armijo(f,fk,dfk,xk,pk,maxIter=30,t=ak),cb::Function=()->())

    his = zeros(maxIter,3)
    # I   = speye(length(x))
    X   = (storeInterm) ? zeros(length(x),maxIter) : []
    fk,dfk  = fdf(x)

    i = 1; flag = -1; a0 = 1.0
    while i<=maxIter

        his[i,1:2] = [fk norm(dfk)]
        if storeInterm; X[:,i] = x; end;
        if norm(dfk)<atol
            his  = his[1:i,:]
            flag = 0
            break
        end

        # get search direction
        pk    = - H*dfk
        # line search
        ak,his[i,3] = lineSearch(f,fk,dfk,x,pk,a0)
        if his[i,3]==1
            a0 =5*ak
        else
            a0= 2*ak
        end
        cb(i)
        if out>0
             Printf.@printf( "iter=%4d\t|f|=%1.2e\t|df|=%1.2e\tLS=%d\tmuLS=%1.2e\n", i, his[i,1] ,his[i,2] ,his[i,3],ak)
        end
        if his[i,3]==-1
             flag = -3
             his  = his[1:i,:]
             break;
        end
        x    += ak*pk
        fk,dfnew    = fdf(x)
        sk    = ak*pk
        yk    = dfnew - dfk
        if dot(yk,sk)>0 # ensure that approximate Hessians remain positive definite
        	H     = (I - (sk*yk')/dot(sk,yk)) * H * (I - (yk*sk')/dot(sk,yk)) + (sk*sk')/dot(yk,sk)
        else
            if out>0
                println("bfgs detected negative curvature. Resetting Hessian")
                a0=1.0
            end
            H = Matrix(1.0I,length(x),length(x))
        end
        dfk  = dfnew
        i+=1
    end
    i = min(maxIter,i)

    if out>=0
        if flag==-1
            Printf.@printf("bfgs iterated maxIter (=%d) times but reached only atol of %1.2e instead of tol=%1.2e",i,his[i,2],atol)
        elseif flag==-3
            Printf.@printf("bfgs stopped of a line search fail at iteration %d.",i)
        elseif out>1
            Printf.@printf("bfgs achieved desired atol of %1.2e at iteration %d.",atol,i)
        end
    end

    if storeInterm; X = X[:,1:i]; end
    return x,flag,his,X,H
end

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
#
# `depth_first_constructive` is a very simple recursive implementation
# but it is more difficult to cap the costs so scales very badly
#

function depth_first_constructive(T::Vector{<:ITensor}) where {LabelSetT}
  indsT = [inds(Tₙ) for Tₙ in T]
  return depth_first_constructive(DimT, indsT)
end

function depth_first_constructive(
  ::Type{DimT}, T::Vector{IndexSetT}
) where {IndexSetT<:IndexSet,DimT}
  labels, dims = label_dims(DimT, T)
  nlabels = length(dims)
  if nlabels ≤ 16
    return depth_first_constructive(UInt16, labels, dims)
  elseif nlabels ≤ 32
    return depth_first_constructive(UInt32, labels, dims)
  elseif nlabels ≤ 64
    return depth_first_constructive(UInt64, labels, dims)
  elseif nlabels ≤ 128
    return depth_first_constructive(UInt128, labels, dims)
  else
    return depth_first_constructive(BitSet, labels, dims)
  end
end

function depth_first_constructive(
  ::Type{LabelSetT}, labels::Vector, dims::Vector
) where {LabelSetT}
  return depth_first_constructive(map(label -> bitset(LabelSetT, label), labels), dims)
end

function depth_first_constructive(
  ::Type{LabelSetT}, ::Type{DimT}, T::Vector{<:ITensor}
) where {LabelSetT,DimT}
  indsT = [inds(Tₙ) for Tₙ in T]
  return depth_first_constructive(LabelSetT, DimT, indsT)
end

function depth_first_constructive(
  ::Type{LabelSetT}, ::Type{DimT}, T::Vector{IndexSetT}
) where {IndexSetT<:IndexSet,LabelSetT,DimT}
  labels, dims = label_dims(DimT, T)
  return depth_first_constructive(map(label -> bitset(LabelSetT, label), labels), dims)
end

function depth_first_constructive(T::Vector, ind_dims::Vector)
  optimal_cost = Ref(typemax(eltype(ind_dims)))
  optimal_sequence = Vector{Pair{Int,Int}}(undef, length(T) - 1)
  _depth_first_constructive!(
    optimal_sequence, optimal_cost, Pair{Int,Int}[], T, ind_dims, collect(1:length(T)), 0
  )
  return pair_sequence_to_tree(optimal_sequence, length(T))
end

function _depth_first_constructive!(
  optimal_sequence, optimal_cost, sequence, T, ind_dims, remaining, cost
)
  if length(remaining) == 1
    # Only should get here if the contraction was the best
    # Otherwise it would have hit the `continue` below
    @assert cost ≤ optimal_cost[]
    optimal_cost[] = cost
    optimal_sequence .= sequence
  end
  for aᵢ in 1:(length(remaining) - 1), bᵢ in (aᵢ + 1):length(remaining)
    a = remaining[aᵢ]
    b = remaining[bᵢ]
    current_cost, Tᵈ = contraction_cost(T[a], T[b], ind_dims)
    new_cost = cost + current_cost
    if new_cost ≥ optimal_cost[]
      continue
    end
    new_sequence = push!(copy(sequence), a => b)
    new_T = push!(copy(T), Tᵈ)
    new_remaining = deleteat!(copy(remaining), (aᵢ, bᵢ))
    push!(new_remaining, length(new_T))
    _depth_first_constructive!(
      optimal_sequence, optimal_cost, new_sequence, new_T, ind_dims, new_remaining, new_cost
    )
  end
end


text:
module BoundaryTriangulations

using Gridap

export BoundaryTriangulation

import Gridap: CellPoints
import Gridap: CellBasis
import Gridap: CellGeomap
import Gridap: CellRefFEs

struct BoundaryTriangulation{Z,D} <: Triangulation{Z,D}
  trian::Triangulation{Z,D}
  descriptor::BoundaryDescriptor

  function BoundaryTriangulation(
    trian::Triangulation{Z,D}, descriptor::BoundaryDescriptor) where {Z,D}
    @assert D == Z + 1
    new{Z,D}(trian,descriptor)
  end

end

for op in (:CellPoints,:CellBasis,:CellGeomap,:CellRefFEs)
  @eval begin
    $op(t::BoundaryTriangulation) = $op(t.trian)
  end
end

end # module


text:
import FLOWFarm; const ff = FLOWFarm

turbine_x = [0.0]
turbine_y = [0.0]
turbine_z = [0.0]
turbine_yaw = [0.0]
diam = 80.0
nturbines = 1

# set turbine design parameters
rotor_diameter = zeros(nturbines) .+ diam # m
hub_height = zeros(nturbines) .+ 70.0   # m
cut_in_speed = zeros(nturbines) .+4.  # m/s
cut_out_speed = zeros(nturbines) .+25.  # m/s
rated_speed = zeros(nturbines) .+16.  # m/s
rated_power = zeros(nturbines) .+2.0E6  # W
generator_efficiency = zeros(nturbines) .+0.944

# load power curve
powerdata = readdlm("./inputfiles/niayifar_vestas_v80_power_curve_observed.txt",  ',', skipstart=1)
velpoints = powerdata[:,1]
powerpoints = powerdata[:,2]*1E6

# initialize power model
power_model = ff.PowerModelPowerPoints(velpoints, powerpoints)
power_models = Vector{typeof(power_model)}(undef, nturbines)
for i = 1:nturbines
    power_models[i] = power_model
end

# load thrust curve
ctdata = readdlm("./inputfiles/mfg-ct-vestas-v80-niayifar2016.txt",  ',', skipstart=1)
velpoints = ctdata[:,1]
ctpoints = ctdata[:,2]

# initialize thurst model
ct_model = ff.ThrustModelConstantCt(0.8)
ct_models = Vector{typeof(ct_model)}(undef, nturbines)
for i = 1:nturbines
    ct_models[i] = ct_model
end

sorted_turbine_index = [1]

wind_speed = 9.0
air_density = 1.1716  # kg/m^3
# air_density = 1.2
winddirections = [270.0*pi/180.0]
windspeeds = [wind_speed]
windprobabilities = [1.0]
measurementheights = [hub_height[1]]
wtvelocities = [wind_speed]
ambient_tis = [0.094]
shearexponent = 0.0
wind_shear_model = ff.PowerLawWindShear(shearexponent)
windresource = ff.DiscretizedWindResource(winddirections, windspeeds, windprobabilities, measurementheights, air_density, ambient_tis, wind_shear_model)

# set up wake and related models
wakedeficitmodel = ff.GaussOriginal(0.04)
# wakedeficitmodel = ff.GaussYawVariableSpread()
wakedeflectionmodel = ff.GaussYawVariableSpreadDeflection()
wakecombinationmodel = ff.LinearLocalVelocitySuperposition()
localtimodel = ff.LocalTIModelNoLocalTI()

model_set = ff.WindFarmModelSet(wakedeficitmodel, wakedeflectionmodel, wakecombinationmodel, localtimodel)