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

@test integrate(x^n, x) == :(x^(1+n)*(1+n)^-1)
@test integrate(exp(x), x) == :(exp(x))
@test integrate(x^-1, x) == :(log(x))
@test integrate(a^x, x) == :(a^x*log(a)^-1)
@test integrate(sin(x), x) == :(-1*cos(x))
@test integrate(cos(x), x) == :(sin(x))
@test integrate(sec(x)^2, x) == :(tan(x))
@test integrate(csc(x)^2, x) == :(-1*cot(x))
@test integrate(sec(x)*tan(x), x) == :(sec(x))
@test integrate(cot(x)*csc(x), x) == :(-1*csc(x))
@test integrate(sinh(x), x) == :(cosh(x))
@test integrate(cosh(x), x) == :(sinh(x))
@test integrate(tan(x), x) == :(-1*log(cos(x)))
@test integrate(cot(x), x) == :(log(sin(x)))
@test integrate(x*sin(x), x) == :(-1*x*cos(x)+sin(x))
@test integrate(log(x), x) == :(-1x+x*log(x))
@test integrate(x^2*exp(x), x) == :(2*exp(x)+x^2*exp(x)+-2*x*exp(x))
@test integrate(exp(x)*sin(x), x) == :((1/2)*exp(x)*sin(x)+-1//2*cos(x)*exp(x))
@test integrate(arctan(x), x) == :(-1//2*log(1+x^2)+x*arctan(x))
@test integrate(x*exp(2x), x) == :(-1//4*exp(2x)+(1/2)*x*exp(2x))
@test integrate(x*cos(x), x) == :(x*sin(x)+cos(x))
@test integrate(x*sin(4x), x) == :(1//16*sin(4x)+-1//4*x*cos(4x))
@test integrate(x*log(x), x) == :(-1//4*x^2+(1/2)*x^2*log(x))
@test integrate(x^2*cos(3x), x) == :(-2//27*sin(3x)+1//3*x^2*sin(3x)+2//9*x*cos(3x))
@test integrate(x^2*sin(2x), x) == :(1//4*cos(2x)+(1/2)*x*sin(2x)+-1//2*x^2*cos(2x))
@test integrate(log(x)^2, x) == :(2x+x*log(x)^2+-2*x*log(x))
@test integrate(arcsin(x), x) == :((1+-1*x^2)^(1/2)+x*arcsin(x))
@test integrate(t*cos(t)*sin(t), t) == :(-1//4*t+(1/2)*t*sin(t)^2+1//4*cos(t)*sin(t))
@test integrate(t*sec(t)^2, t) == :(t*tan(t)+log(cos(t)))
@test integrate(t^2*log(t), t) == :(-1//9*t^3+1//3*t^3*log(t))
@test integrate(t^3*exp(t), t) == :(-6*exp(t)+t^3*exp(t)+-3*t^2*exp(t)+6*t*exp(t))
@test integrate(exp(2t)*sin(3t), t) == :(-3//13*cos(3t)*exp(2t)+2//13*exp(2t)*sin(3t))
@test integrate(cos(3t)*exp(-1t), t) == :(-1//10*cos(3t)*exp(-1t)+3//10*exp(-1t)*sin(3t))
@test integrate(y*sinh(y), y) == :(-1*sinh(y)+y*cosh(y))
@test integrate(y*cosh(a*y), y) == :(-1*a^-2*cosh(a*y)+y*a^-1*sinh(a*y))
@test integrate(t*exp(-1t), t) == :(-1*exp(-1t)+-1*t*exp(-1t))
@test integrate(t^(1/2)*log(t), t) == :(-4//9*t^3//2+2//3*t^3//2*log(t))
@test integrate(x*cos(2x), x) == :(1//4*cos(2x)+(1/2)*x*sin(2x))
@test integrate(x^2*exp(-1x), x) == :(-2*exp(-1x)+-1*x^2*exp(-1x)+-2*x*exp(-1x))
@test integrate(arccos(x), x) == :(-1*(1+-1*x^2)^(1/2)+x*arccos(x))
@test integrate(x*csc(x)^2, x) == :(-1*x*cot(x)+log(sin(x)))
@test integrate(cos(5x)*sin(3x), x) == :(-1//16*cos(8x)+1//4*cos(2x))
@test integrate(sin(2x)*sin(4x), x) == :(-1//12*sin(6x)+1//4*sin(2x))
@test integrate(cos(x)*log(sin(x)), x) == :(-1*sin(x)+log(sin(x))*sin(x))
@test integrate(x^3*exp(x^2), x) == :(-1//2*exp(x^2)+(1/2)*x^2*exp(x^2))
@test integrate((3+2x)*exp(x), x) == :(-2*exp(x)+(3+2x)*exp(x))
@test integrate(x*5^x, x) == :(-1*5^x*log(5)^-2+x*5^x*log(5)^-1)
@test integrate(cos(log(x)), x) == :((1/2)*x*cos(log(x))+(1/2)*x*sin(log(x)))
@test integrate(exp(x^(1/2)), x) == :(-2*exp(x^(1/2))+2*x^(1/2)*exp(x^(1/2)))
@test integrate(log(x^(1/2)), x) == :(-1//2*x+x*log(x^(1/2)))
@test integrate(sin(log(x)), x) == :((1/2)*x*sin(log(x))+-1//2*x*cos(log(x)))
@test integrate(sin(x^(1/2)), x) == :(2*sin(x^(1/2))+-2*x^(1/2)*cos(x^(1/2)))
@test integrate(x^5*cos(x^3), x) == :(1//3*cos(x^3)+1//3*x^3*sin(x^3))
@test integrate(x^5*exp(x^2), x) == :((1/2)*x^4*exp(x^2)+-1*x^2*exp(x^2)+exp(x^2))
@test integrate(x*arctan(x), x) == :((1/2)*arctan(x)+-1//2*x+(1/2)*x^2*arctan(x))
@test integrate(x*cos(pi*x), x) == :(pi^-2*cos(pi*x)+x*pi^-1*sin(pi*x))
@test integrate(x^(1/2)*log(x), x) == :(-4//9*x^3//2+2//3*x^3//2*log(x))
@test integrate(sin(3x)^2, x) == :((1/2)*x+-1//6*cos(3x)*sin(3x))
@test integrate(cos(x)^2, x) == :((1/2)*x+(1/2)*cos(x)*sin(x))
@test integrate(cos(x)^4, x) == :(3//8*x+1//4*cos(x)^3*sin(x)+3//8*cos(x)*sin(x))
@test integrate(sin(x)^3, x) == :(-1*cos(x)+1//3*cos(x)^3)
@test integrate(cos(x)^4*sin(x)^3, x) == :(-1//5*cos(x)^5+1//7*cos(x)^7)
@test integrate(cos(x)^3*sin(x)^4, x) == :(-1//7*sin(x)^7+1//5*sin(x)^5)
@test integrate(cos(x)^2*sin(x)^4, x) == :(1//16*x+-1//6*cos(x)^3*sin(x)^3+-1//8*cos(x)^3*sin(x)+1//16*cos(x)*sin(x))
@test integrate(cos(x)^2*sin(x)^2, x) == :(1//8*x+-1//4*cos(x)^3*sin(x)+1//8*cos(x)*sin(x))
@test integrate((1+-1*sin(2x))^2, x) == :(3//2*x+-1//4*cos(2x)*sin(2x)+cos(2x))
@test integrate(cos(x)*sin(x+1//6*pi), x) == :(-1//4*cos(2x+1//6*pi)+1//4*x)
@test integrate(cos(x)^5*sin(x)^5, x) == :(-1//4*sin(x)^8+1//6*sin(x)^6+1//10*sin(x)^10)
@test integrate(sin(x)^6, x) == :(5//16*x+-5//16*cos(x)*sin(x)+-5//24*sin(x)^3*cos(x)+-1//6*sin(x)^5*cos(x))
@test integrate(cos(x)^6, x) == :(5//16*x+1//6*cos(x)^5*sin(x)+5//16*cos(x)*sin(x)+5//24*cos(x)^3*sin(x))
@test integrate(cos(2x)^4*sin(2x)^2, x) == :(1//16*x+-1//12*cos(2x)^5*sin(2x)+1//32*cos(2x)*sin(2x)+1//48*cos(2x)^3*sin(2x))
@test integrate(sin(x)^5, x) == :(-1*cos(x)+-1//5*cos(x)^5+2//3*cos(x)^3)
@test integrate(cos(x)^4*sin(x)^4, x) == :(3//128*x+-1//8*cos(x)^5*sin(x)^3+-1//16*cos(x)^5*sin(x)+1//64*cos(x)^3*sin(x)+3//128*cos(x)*sin(x))
@test integrate(cos(x)^(1/2)*sin(x)^3, x) == :(-2//3*cos(x)^3//2+2//7*cos(x)^7//2)
@test integrate(cos(x)^3*sin(x)^(1/2), x) == :(-2//7*sin(x)^7//2+2//3*sin(x)^3//2)
@test integrate(x^-1//2*cos(x^(1/2))^2, x) == :(x^(1/2)+cos(x^(1/2))*sin(x^(1/2)))
@test integrate(x*sin(x^2)^3, x) == :(-1//2*cos(x^2)+1//6*cos(x^2)^3)
@test integrate(cos(x)^2*tan(x)^3, x) == :((1/2)*cos(x)^2+-1*log(cos(x)))
@test integrate(cot(x)^5*sin(x)^2, x) == :((1/2)*sin(x)^2+-2*log(sin(x))+-1//2*csc(x)^2)
@test integrate(cos(x)^-1*(1+-1*sin(x)), x) == :(log(1+sin(x)))
@test integrate((1+-1*sin(x))^-1, x) == :((1+-1*sin(x))^-1*cos(x))
@test integrate(tan(x)^2, x) == :(-1x+tan(x))
@test integrate(tan(x)^4, x) == :(x+-1*tan(x)+1//3*tan(x)^3)
@test integrate(sec(x)^4, x) == :(1//3*tan(x)^3+tan(x))
@test integrate(sec(x)^6, x) == :(1//5*tan(x)^5+2//3*tan(x)^3+tan(x))
@test integrate(sec(x)^2*tan(x)^4, x) == :(1//5*tan(x)^5)
@test integrate(sec(x)^4*tan(x)^2, x) == :(1//3*tan(x)^3+1//5*tan(x)^5)
@test integrate(sec(x)^3*tan(x), x) == :(1//3*sec(x)^3)
@test integrate(sec(x)^3*tan(x)^3, x) == :(-1//3*sec(x)^3+1//5*sec(x)^5)
@test integrate(tan(x)^5, x) == :(-1*log(cos(x))+-1//2*tan(x)^2+1//4*tan(x)^4)
@test integrate(tan(x)^6, x) == :(-1x+-1//3*tan(x)^3+1//5*tan(x)^5+tan(x))
@test integrate(tan(x)^5*sec(x), x) == :(-2//3*sec(x)^3+1//5*sec(x)^5+sec(x))
@test integrate(sec(x)^3*tan(x)^5, x) == :(-2//5*sec(x)^5+1//3*sec(x)^3+1//7*sec(x)^7)
@test integrate(sec(x)^6*tan(x), x) == :(1//6*sec(x)^6)
@test integrate(sec(x)^6*tan(x)^3, x) == :(-1//6*sec(x)^6+1//8*sec(x)^8)
@test integrate(cot(x)^-1*sec(x)^2, x) == :((1/2)*sec(x)^2)
@test integrate(tan(x)^2*sec(x), x) == :(-1//2*arctanh(sin(x))+(1/2)*sec(x)*tan(x))
@test integrate(cot(x)^2, x) == :(-1x+-1*cot(x))
@test integrate(cot(x)^3, x) == :(-1*log(sin(x))+-1//2*cot(x)^2)
@test integrate(cot(x)^4*csc(x)^4, x) == :(-1//5*cot(x)^5+-1//7*cot(x)^7)
@test integrate(cot(x)^3*csc(x)^4, x) == :(-1//6*csc(x)^6+1//4*csc(x)^4)
@test integrate(csc(x), x) == :(-1*arctanh(cos(x)))
@test integrate(csc(x)^3, x) == :(-1//2*arctanh(cos(x))+-1//2*cot(x)*csc(x))
@test integrate(cos(x)^2*sin(x)^-1, x) == :(-1*arctanh(cos(x))+cos(x))
@test integrate(sin(x)^-4, x) == :(-1*cot(x)+-1//3*cot(x)^3)
@test integrate(sin(2x)*sin(5x), x) == :(-1//14*sin(7x)+1//6*sin(3x))
@test integrate(cos(x)*sin(3x), x) == :(-1//4*cos(2x)+-1//8*cos(4x))
@test integrate(cos(3x)*cos(4x), x) == :((1/2)*sin(x)+1//14*sin(7x))
@test integrate(sin(3x)*sin(6x), x) == :(-1//18*sin(9x)+1//6*sin(3x))
@test integrate(cos(x)^5*sin(x), x) == :(-1//6*cos(x)^6)
@test integrate(cos(x)*cos(2x)*cos(3x), x) == :(1//4*x+1//8*sin(2x)+1//16*sin(4x)+1//24*sin(6x))
@test integrate(sec(x)^-2*(1+-1*tan(x)^2), x) == :(cos(x)*sin(x))
@test integrate(sin(2x)^-1*(cos(x)+sin(x)), x) == :((1/2)*arctanh(sin(x))+-1//2*arctanh(cos(x)))
@test integrate(sin(x)^2*tan(x), x) == :((1/2)*cos(x)^2+-1*log(cos(x)))
@test integrate(cos(x)^2*cot(x)^3, x) == :((1/2)*sin(x)^2+-2*log(sin(x))+-1//2*csc(x)^2)
@test integrate(sec(x)^3*tan(x), x) == :(1//3*sec(x)^3)
@test integrate(sec(x)^3*tan(x)^3, x) == :(-1//3*sec(x)^3+1//5*sec(x)^5)
@test integrate(x^-2*(9+-1*x^2)^(1/2), x) == :(-1*arcsin(1//3*x)+-1*x^-1*(9+-1*x^2)^(1/2))
@test integrate(x^-2*(4+x^2)^-1//2, x) == :(-1//4*x^-1*(4+x^2)^(1/2))
@test integrate(x*(4+x^2)^-1//2, x) == :((4+x^2)^(1/2))
@test integrate((x^2+-1*a^2)^-1//2, x) == :(arctanh(x*(x^2+-1*a^2)^-1//2))
@test integrate(x^3*(9+4*x^2)^-3//2, x) == :(1//16*(9+4*x^2)^(1/2)+9//16*(9+4*x^2)^-1//2)
@test integrate(x*(3+-1*x^2+-2x)^-1//2, x) == :(-1*(3+-1*x^2+-2x)^(1/2)+arcsin(-1//2+-1//2*x))
@test integrate(x^-2*(1+-1*x^2)^-1//2, x) == :(-1*x^-1*(1+-1*x^2)^(1/2))
@test integrate(x^3*(4+-1*x^2)^(1/2), x) == :(-4//3*(4+-1*x^2)^3//2+1//5*(4+-1*x^2)^5//2)
@test integrate(x*(1+-1*x^2)^-1//2, x) == :(-1*(1+-1*x^2)^(1/2))
@test integrate(x*(4+-1*x^2)^(1/2), x) == :(-1//3*(4+-1*x^2)^3//2)
@test integrate((1+-4*x^2)^(1/2), x) == :(1//4*arcsin(2x)+(1/2)*x*(1+-4*x^2)^(1/2))
@test integrate(x^3*(4+x^2)^-1//2, x) == :(-4*(4+x^2)^(1/2)+1//3*(4+x^2)^3//2)
@test integrate((9+x^2)^-1//2, x) == :(arcsinh(1//3*x))
@test integrate((1+x^2)^(1/2), x) == :((1/2)*arcsinh(x)+(1/2)*x*(1+x^2)^(1/2))
@test integrate(x^-3*(-16+x^2)^-1//2, x) == :(1//128*arctan(1//4*(-16+x^2)^(1/2))+1//32*x^-2*(-16+x^2)^(1/2))
@test integrate(x^-4*(x^2+-1*a^2)^(1/2), x) == :(1//3*a^-2*x^-3*(x^2+-1*a^2)^3//2)
@test integrate(x^-1*(-4+9*x^2)^(1/2), x) == :((-4+9*x^2)^(1/2)+-2*arctan((1/2)*(-4+9*x^2)^(1/2)))
@test integrate(x^-2*(-9+16*x^2)^-1//2, x) == :(1//9*x^-1*(-9+16*x^2)^(1/2))
@test integrate(x^2*(a^2+-1*x^2)^-3//2, x) == :(-1*arctan(x*(a^2+-1*x^2)^-1//2)+x*(a^2+-1*x^2)^-1//2)
@test integrate(x^2*(5+-1*x^2)^-1//2, x) == :(5//2*arcsin(1//5*x*5^(1/2))+-1//2*x*(5+-1*x^2)^(1/2))
@test integrate(x^-1*(3+x^2)^-1//2, x) == :(-1//3*3^(1/2)*arctanh(1//3*3^(1/2)*(3+x^2)^(1/2)))
@test integrate(x*(4+x^2)^-5//2, x) == :(-1//3*(4+x^2)^-3//2)
@test integrate(x^3*(4+-9*x^2)^(1/2), x) == :(-4//243*(4+-9*x^2)^3//2+1//405*(4+-9*x^2)^5//2)
@test integrate(x^2*(9+-1*x^2)^(1/2), x) == :(81//8*arcsin(1//3*x)+-9//8*x*(9+-1*x^2)^(1/2)+1//4*x^3*(9+-1*x^2)^(1/2))
@test integrate(5*x*(1+x^2)^(1/2), x) == :(5//3*(1+x^2)^3//2)
@test integrate((-25+4*x^2)^-3//2, x) == :(-1//25*x*(-25+4*x^2)^-1//2)
@test integrate((-1*x^2+2x)^(1/2), x) == :(-1//2*arcsin(1+-1x)+(-1*x^2+2x)^(1/2)*(-1//2+(1/2)*x))
@test integrate((8+x^2+4x)^-1//2, x) == :(arcsinh(1+(1/2)*x))
@test integrate((-8+6x+9*x^2)^-1//2, x) == :(1//3*arctanh((-8+6x+9*x^2)^-1//2*(1+3x)))
@test integrate(x^2*(-1*x^2+4x)^-1//2, x) == :(-6*arcsin(1+-1//2*x)+-3*(-1*x^2+4x)^(1/2)+-1//2*x*(-1*x^2+4x)^(1/2))
@test integrate((2+x^2+2x)^-2, x) == :((1/2)*arctan(1+x)+(4+2*x^2+4x)^-1*(1+x))
@test integrate((5+-1*x^2+-4x)^-5//2, x) == :(1//27*(5+-1*x^2+-4x)^-3//2*(2+x)+1//243*(5+-1*x^2+-4x)^-1//2*(4+2x))
@test integrate((9+-1*exp(2t))^(1/2)*exp(t), t) == :(9//2*arcsin(1//3*exp(t))+(1/2)*(9+-1*exp(2t))^(1/2)*exp(t))
@test integrate((-9+exp(2t))^(1/2), t) == :((-9+exp(2t))^(1/2)+-3*arctan(1//3*(-9+exp(2t))^(1/2)))
@test integrate((a^2+x^2)^-1//2, x) == :(arctanh(x*(a^2+x^2)^-1//2))
@test integrate((-2+x+x^2)^-1*(5+x), x) == :(-1*log(2+x)+2*log(1+-1x))
@test integrate((-1+x)^-1*(x+x^3), x) == :((1/2)*x^2+2x+2*log(1+-1x)+1//3*x^3)
@test integrate((-2x+2*x^3+3*x^2)^-1*(-1+x^2+2x), x) == :((1/2)*log(x)+-1//10*log(2+x)+1//10*log(1+-2x))
@test integrate((1+x^3+-1x+-1*x^2)^-1*(1+x^4+-2*x^2+4x), x) == :(x+(1/2)*x^2+-1*log(1+x)+2*(1+-1x)^-1+log(1+-1x))
@test integrate((x^3+4x)^-1*(4+-1x+2*x^2), x) == :((1/2)*log(4+x^2)+-1//2*arctan((1/2)*x)+log(x))
@test integrate((3+-4x+4*x^2)^-1*(2+-3x+4*x^2), x) == :(x+1//8*log(3+-4x+4*x^2)+1//8*2^(1/2)*arctan((1/2)*2^(1/2)*(1+-2x)))
@test integrate(x^-1*(1+x^2)^-3*(-1+x)^-1*(1+x+x^2)^-1*(1+x^2+x^3), x) == :(-1*log(x)+-1//2*log(1+x+x^2)+1//8*log(1+-1x)+7//16*arctan(x)+15//16*log(1+x^2)+-1*(8+8*x^2)^-1*(3+-3x)+3*x*(16+16*x^2)^-1+-1//3*3^(1/2)*arctan(1//3*3^(1/2)*(1+2x))+1//8*(1+x^2)^-2*(1+x))
@test integrate(x^-1*(1+x^2)^-2*(1+-1*x^3+-3x+2*x^2), x) == :(-2*arctan(x)+-1//2*log(1+x^2)+-1*(2+2*x^2)^-1*(1+2x)+log(x))
@test integrate((1+x^2)^-2, x) == :((1/2)*arctan(x)+x*(2+2*x^2)^-1)
@test integrate((-1+x)^-1*(2+x)^-1, x) == :(-1//3*log(2+x)+1//3*log(1+-1x))
@test integrate(7*(-12+2*x^2+5x)^-1, x) == :(-7//11*log(4+x)+7//11*log(3+-2x))
@test integrate((-1+2x)^-2*(3+2x)^-1*(-4+x^2+3x), x) == :(-9*(32+-64x)^-1+-25//128*log(3+2x)+41//128*log(1+-2x))
@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^4+-1*x^3)^-1, x) == :(x^-1+(1/2)*x^-2+-1*log(x)+log(1+-1x))
@test integrate((x^3+-1x)^-1*(1+x^3+x^4+-1x+-1*x^2), x) == :(x+(1/2)*x^2+(1/2)*log(1+-1*x^2)+-1*log(x))
@test integrate(x^-1*(2+x^2)^-1*(-2+x^2), x) == :(-1*log(x)+log(2+x^2))
@test integrate((1+x^2)^-1*(2+x^2)^-1*(2+x^3+-4*x^2), x) == :(6*arctan(x)+-1//2*log(1+x^2)+-5*2^(1/2)*arctan((1/2)*x*2^(1/2))+log(2+x^2))
@test integrate((1+x^2)^-1*(4+x^2)^-2*(1+x^2+x^4), x) == :(1//9*arctan(x)+25//144*arctan((1/2)*x)+-13*x*(96+24*x^2)^-1)
@test integrate((-3+2x)^-1*(5+x)^-2*(1+x+x^2)^-1*(1+16x), x) == :(-79*(1365+273x)^-1+-481//5586*log(1+x+x^2)+200//3211*log(3+-2x)+2731//24843*log(5+x)+451//8379*3^(1/2)*arctan(1//3*3^(1/2)*(1+2x)))
@test integrate(x^4*(9+x^2)^-3, x) == :(1//8*arctan(1//3*x)+-3*x*(72+8*x^2)^-1+-1//4*x^3*(9+x^2)^-2)
@test integrate(19*x*(-1+x)^-3*(3+4*x^2+5x)^-2, x) == :(-1843*(4416+-4416x)^-1+-399//736*(1+-1x)^-2+-209//4608*log(3+4*x^2+5x)+209//2304*log(1+-1x)+114437//1218816*23^(1/2)*arctan(1//23*23^(1/2)*(5+8x))+1//276*(1+-1x)^-2*(3+4*x^2+5x)^-1*(741+836x))
@test integrate((x^3+x^4+2*x^2)^-1*(1+x^2+x^3), x) == :(-1//2*x^-1+-1//4*log(x)+5//8*log(2+x+x^2)+1//28*7^(1/2)*arctan(1//7*7^(1/2)*(1+2x)))
@test integrate((x^6+-1*x^3)^-1, x) == :((1/2)*x^-2+-1//6*log(1+x+x^2)+1//3*log(1+-1x)+-1//3*3^(1/2)*arctan(1//3*3^(1/2)*(1+2x)))
@test integrate(x^2*(1+x)^-1, x) == :((1/2)*x^2+-1x+log(1+x))
@test integrate(x*(-5+x)^-1, x) == :(x+5*log(5+-1x))
@test integrate((-1+x)^-1*(2+x)^-1*(-1+4x), x) == :(3*log(2+x)+log(1+-1x))
@test integrate((1+x)^-1*(2+x)^-1, x) == :(-1*log(2+x)+log(1+x))
@test integrate((3+2x)^-1*(-5+6x), x) == :(-7*log(3+2x)+3x)
@test integrate((a+x)^-1*(b+x)^-1, x) == :((a+-1b)^-1*log(b+x)+-1*(a+-1b)^-1*log(a+x))
@test integrate((x^2+-1x)^-1*(1+x^2), x) == :(x+-1*log(x)+2*log(1+-1x))
@test integrate((-12+x+x^2)^-1*(1+x^2+x^3+-12x), x) == :((1/2)*x^2+-1//7*log(4+x)+1//7*log(3+-1x))
@test integrate((1+x)^-2*(3+2x), x) == :(-1*(1+x)^-1+2*log(1+x))
@test integrate(x^-1*(1+x)^-1*(3+2x)^-1, x) == :(-1*log(1+x)+1//3*log(x)+2//3*log(3+2x))
@test integrate((x^3+-3x+2*x^2)^-1*(-3+5x+6*x^2), x) == :(2*log(1+-1x)+3*log(3+x)+log(x))
@test integrate(x*(4+x^2+4x)^-1, x) == :(2*(2+x)^-1+log(2+x))
@test integrate((-1+x)^-2*(4+x)^-1, x) == :((5+-5x)^-1+-1//25*log(1+-1x)+1//25*log(4+x))
@test integrate(x^2*(-3+x)^-1*(2+x)^-2, x) == :(4*(10+5x)^-1+9//25*log(3+-1x)+16//25*log(2+x))
@test integrate((x^3+2*x^2)^-1*(-2+3x+5*x^2), x) == :(x^-1+2*log(x)+3*log(2+x))
@test integrate((-6+x+x^3+4*x^2)^-1*(18+-4*x^2+-2x), x) == :(-3*log(3+x)+-2*log(2+x)+log(1+-1x))
@test integrate((4+x^3+3*x^2)^-1*(x^2+2x), x) == :(1//3*log(4+x^3+3*x^2))
@test integrate(x^-2*(-1+x)^-2, x) == :((1+-1x)^-1+-1*x^-1+-2*log(1+-1x)+2*log(x))
@test integrate(x^2*(1+x)^-3, x) == :(2*(1+x)^-1+-1//2*(1+x)^-2+log(1+x))
@test integrate((x^4+-1*x^2)^-1, x) == :(x^-1+-1*arctanh(x))
@test integrate((1+x^4+-1*x^2)^-1*(-1x+2*x^3), x) == :((1/2)*log(1+x^4+-1*x^2))
@test integrate(x^3*(1+x^2)^-1, x) == :((1/2)*x^2+-1//2*log(1+x^2))
@test integrate((2+x^2+2x)^-1*(-1+x), x) == :((1/2)*log(2+x^2+2x)+-2*arctan(1+x))
@test integrate(x*(1+x+x^2)^-1, x) == :((1/2)*log(1+x+x^2)+-1//3*3^(1/2)*arctan(1//3*3^(1/2)*(1+2x)))
@test integrate((5+4x+4*x^2)^-1*(7+4*x^2+5x), x) == :(x+1//8*log(5+4x+4*x^2)+3//8*arctan(1/2+x))
@test integrate((1+x^2)^-1*(-1+x)^-1*(5+-4x+3*x^2), x) == :((1/2)*log(1+x^2)+-3*arctan(x)+2*log(1+-1x))
@test integrate((x^3+3x)^-1*(3+2x), x) == :(-1//2*log(3+x^2)+2//3*3^(1/2)*arctan(1//3*x*3^(1/2))+log(x))
@test integrate((-1+x^3)^-1, x) == :(-1//6*log(1+x+x^2)+1//3*log(1+-1x)+-1//3*3^(1/2)*arctan(1//3*3^(1/2)*(1+2x)))
@test integrate(x^3*(1+x^3)^-1, x) == :(x+-1//3*log(1+x)+1//6*log(1+x^2+-1x)+1//3*3^(1/2)*arctan(1//3*3^(1/2)*(1+-2x)))
@test integrate((1+x^2)^-1*(-1+x)^-2*(-1+x^2+-2x), x) == :((-1+x)^-1+-1//2*log(1+x^2)+arctan(x)+log(1+-1x))
@test integrate(x^4*(-1+x^4)^-1, x) == :(x+-1//2*arctan(x)+-1//2*arctanh(x))
@test integrate((1+x^2)^-1*(2+x^2)^-1*(-4+-1*x^2+3*x^3+6x), x) == :(-3*arctan(x)+3//2*log(1+x^2)+2^(1/2)*arctan((1/2)*x*2^(1/2)))
@test integrate((4+x^4+5*x^2)^-1*(1+x+x^3+-2*x^2), x) == :((1/2)*log(4+x^2)+-3//2*arctan((1/2)*x)+arctan(x))
@test integrate((4+x^2+2x)^-2*(-3+x), x) == :((24+6*x^2+12x)^-1*(-7+-4x)+-2//9*3^(1/2)*arctan(1//3*3^(1/2)*(1+x)))
@test integrate(x^-1*(1+x^2)^-2*(1+x^4), x) == :((1+x^2)^-1+log(x))
@test integrate((2+sin(x)^2+-3*sin(x))^-1*(-3+2*sin(x))*cos(x), x) == :(log(2+sin(x)^2+-3*sin(x)))
@test integrate((5+cos(x)^2)^-1*cos(x)^2*sin(x), x) == :(-1*cos(x)+5^(1/2)*arctan(1//5*5^(1/2)*cos(x)))
@test integrate((-3+x^2+2x)^-1, x) == :(-1//4*log(3+x)+1//4*log(1+-1x))
@test integrate((x^2+-2x)^-1, x) == :((1/2)*log(2+-1x)+-1//2*log(x))
@test integrate((-7+4*x^2+12x)^-1*(1+2x), x) == :(1//8*log(1+-2x)+3//8*log(7+2x))
@test integrate(x*(-1+x+x^2)^-1, x) == :(1//10*(5+5^(1/2))*log(1+5^(1/2)+2x)+1//10*(5+-1*5^(1/2))*log(1+-1*5^(1/2)+2x))
@test integrate((-70+-299x+-286*x^2+-13*x^4+30*x^5+50*x^3)^-1*(-32+-27*x^2+4*x^3+5x), x) == :(-3146//80155*log(7+-3x)+-334//323*log(1+2x)+4822//4879*log(2+5x)+11049//260015*log(5+x+x^2)+3988//260015*19^(1/2)*arctan(1//19*19^(1/2)*(1+2x)))
@test integrate((4+-80*x^3+-80*x^5+-20x+41*x^2+100*x^6+116*x^4)^-1*(8+-13*x^2+-7*x^3+12*x^5), x) == :(5828*(18150+-45375x)^-1+-59096//99825*log(2+-5x)+2843//7986*log(1+2*x^2)+-1*(1452+2904*x^2)^-1*(313+502x)+503//15972*2^(1/2)*arctan(x*2^(1/2)))
@test integrate(x^-1*(4+x)^(1/2), x) == :(-4*arctanh((1/2)*(4+x)^(1/2))+2*(4+x)^(1/2))
@test integrate((x^(1/2)+-1*x^-1//3)^-1, x) == :(2*x^(1/2)+6//5*log(1+-1*x^1//6)+-1*(3//10+-3//10*5^(1/2))*log(2+x^1//6+2*x^1//3+5^(1/2)*x^1//6)+-1*(3//10+3//10*5^(1/2))*log(2+x^1//6+2*x^1//3+-1*5^(1/2)*x^1//6)+-3//5*(10+2*5^(1/2))^(1/2)*arctan((1/2)*(1/2+1//10*5^(1/2))^(1/2)*(1+5^(1/2)+4*x^1//6))+3//5*(10+-2*5^(1/2))^(1/2)*arctan((10+2*5^(1/2))^-1//2*(1+-1*5^(1/2)+4*x^1//6)))
@test integrate((-4*cos(x)+3*sin(x))^-1, x) == :(-1//5*arctanh(3//5*cos(x)+4//5*sin(x)))
@test integrate((1+x^(1/2))^-1, x) == :(-2*log(1+x^(1/2))+2*x^(1/2))
@test integrate((1+x^-1//3)^-1, x) == :(x+-1*log(x)+-3*log(1+x^-1//3)+3*x^1//3+-3//2*x^2//3)
@test integrate(x^(1/2)*(1+x)^-1, x) == :(-2*arctan(x^(1/2))+2*x^(1/2))
@test integrate(x^-1*(1+x)^-1//2, x) == :(-2*arctanh((1+x)^(1/2)))
@test integrate((x+-1*x^1//3)^-1, x) == :(3//2*log(1+-1*x^2//3))
@test integrate((x+-1*(2+x)^(1/2))^-1, x) == :(2//3*log(1+(2+x)^(1/2))+4//3*log(2+-1*(2+x)^(1/2)))
@test integrate(x^2*(-1+x)^-1//2, x) == :(2*(-1+x)^(1/2)+2//5*(-1+x)^5//2+4//3*(-1+x)^3//2)
@test integrate((1+x)^-1*(-1+x)^(1/2), x) == :(2*(-1+x)^(1/2)+-2*2^(1/2)*arctan((1/2)*2^(1/2)*(-1+x)^(1/2)))
@test integrate((1+x^(1/2))^-1//2, x) == :(-4*(1+x^(1/2))^(1/2)+4//3*(1+x^(1/2))^3//2)
@test integrate(x^(1/2)*(x+x^2)^-1, x) == :(2*arctan(x^(1/2)))
@test integrate((-1+x^(1/2))^-1*(1+x^(1/2)), x) == :(x+4*x^(1/2)+4*log(1+-1*x^(1/2)))
@test integrate((-1+x^-1//3)^-1*(1+x^-1//3), x) == :(-1x+-6*x^1//3+-6*log(1+-1*x^1//3)+-3*x^2//3)
@test integrate(x^3*(1+x^2)^-1//3, x) == :(-3//4*(1+x^2)^2//3+3//10*(1+x^2)^5//3)
@test integrate(x^(1/2)*(x^(1/2)+-1*x^-1//3)^-1, x) == :(x+6*x^1//6+6//5*log(1+-1*x^1//6)+-1*(3//10+-3//10*5^(1/2))*log(2+x^1//6+2*x^1//3+-1*5^(1/2)*x^1//6)+-1*(3//10+3//10*5^(1/2))*log(2+x^1//6+2*x^1//3+5^(1/2)*x^1//6)+-3//5*(10+-2*5^(1/2))^(1/2)*arctan((1/2)*(1/2+1//10*5^(1/2))^(1/2)*(1+5^(1/2)+4*x^1//6))+-3//5*(10+2*5^(1/2))^(1/2)*arctan((10+2*5^(1/2))^-1//2*(1+-1*5^(1/2)+4*x^1//6)))
@test integrate((x^(1/2)+x^-1//4)^-1, x) == :(2*x^(1/2)+-2//3*log(1+x^(1/2)+-1*x^1//4)+4//3*log(1+x^1//4)+4//3*3^(1/2)*arctan(1//3*3^(1/2)*(1+-2*x^1//4)))
@test integrate((x^-1//3+x^-1//4)^-1, x) == :(-1x+-12*log(1+x^1//12)+-6*x^1//6+-3*x^1//3+-2*x^(1/2)+4*x^1//4+12*x^1//12+-6//5*x^5//6+-6//7*x^7//6+-3//2*x^2//3+4//3*x^3//4+4//5*x^5//4+12//5*x^5//12+12//7*x^7//12+12//11*x^11//12+12//13*x^13//12)
@test integrate((x^-1*(1+-1x))^(1/2), x) == :(-1*arctan((-1+x^-1)^(1/2))+x*(-1+x^-1)^(1/2))
@test integrate((sin(x)^2+sin(x))^-1*cos(x), x) == :(-1*log(1+sin(x))+log(sin(x)))
@test integrate((2+3*exp(x)+exp(2x))^-1*exp(2x), x) == :(-1*log(1+exp(x))+2*log(2+exp(x)))
@test integrate((1+exp(x))^-1//2, x) == :(-2*arctanh((1+exp(x))^(1/2)))
@test integrate((1+-1*exp(x))^(1/2), x) == :(-2*arctanh((1+-1*exp(x))^(1/2))+2*(1+-1*exp(x))^(1/2))
@test integrate((3+-5*sin(x))^-1, x) == :(-1//4*log(-3*sin((1/2)*x)+cos((1/2)*x))+1//4*log(-1*sin((1/2)*x)+3*cos((1/2)*x)))
@test integrate((cos(x)+sin(x))^-1, x) == :(-1//2*2^(1/2)*arctanh((1/2)*2^(1/2)*(-1*sin(x)+cos(x))))
@test integrate((1+-1*cos(x)+sin(x))^-1, x) == :(-1*log(1+cot((1/2)*x)))
@test integrate((3*sin(x)+4*cos(x))^-1, x) == :(-1//5*arctanh(-4//5*sin(x)+3//5*cos(x)))
@test integrate((sin(x)+tan(x))^-1, x) == :(-1//2*csc(x)^2+-1//2*arctanh(cos(x))+(1/2)*cot(x)*csc(x))
@test integrate((2*sin(x)+sin(2x))^-1, x) == :(1//4*log(tan((1/2)*x))+1//8*tan((1/2)*x)^2)
@test integrate((1+sin(x))^-1*sec(x), x) == :((1/2)*arctanh(sin(x))+-1*(2+2*sin(x))^-1)
@test integrate((a*sin(x)+b*cos(x))^-1, x) == :(-1*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(a*cos(x)+-1*b*sin(x))))
@test integrate((a^2*sin(x)^2+b^2*cos(x)^2)^-1, x) == :(a^-1*b^-1*arctan(a*b^-1*tan(x)))
@test integrate(x*(-1+x^2)^-1, x) == :((1/2)*log(1+-1*x^2))
@test integrate(x^(1/2)*(1+x^(1/2)), x) == :((1/2)*x^2+2//3*x^3//2)
@test integrate((1+-1*cos(x))^-1, x) == :(-1*(1+-1*cos(x))^-1*sin(x))
@test integrate(tan(x)^2*sec(x), x) == :(-1//2*arctanh(sin(x))+(1/2)*sec(x)*tan(x))
@test integrate(sec(x)^3*tan(x)^3, x) == :(-1//3*sec(x)^3+1//5*sec(x)^5)
@test integrate(exp(x^(1/2)), x) == :(-2*exp(x^(1/2))+2*x^(1/2)*exp(x^(1/2)))
@test integrate((x^3+-10x+-3*x^2)^-1*(1+x^5), x) == :(19x+-31//14*log(2+x)+-1//10*log(x)+1//3*x^3+3//2*x^2+3126//35*log(5+-1x))
@test integrate(x^-1*log(x)^-1//2, x) == :(2*log(x)^(1/2))
@test integrate((-3+x)^-1*(5+2x), x) == :(2x+11*log(3+-1x))
@test integrate(exp(x+exp(x)), x) == :(exp(exp(x)))
@test integrate(cos(x)^2*sin(x)^2, x) == :(1//8*x+-1//4*cos(x)^3*sin(x)+1//8*cos(x)*sin(x))
@test integrate((cos(x)+sin(x))^-1*(-1*cos(x)+sin(x)), x) == :(-1*log(cos(x)+sin(x)))
@test integrate(x*(1+-1*x^2)^-1//2, x) == :(-1*(1+-1*x^2)^(1/2))
@test integrate(x^3*log(x), x) == :(-1//16*x^4+1//4*x^4*log(x))
@test integrate((-2+x)^(1/2)*(2+x)^-1, x) == :(-4*arctan((1/2)*(-2+x)^(1/2))+2*(-2+x)^(1/2))
@test integrate(x*(2+x)^-2, x) == :(2*(2+x)^-1+log(2+x))
@test integrate(log(1+x^2), x) == :(-2x+2*arctan(x)+x*log(1+x^2))
@test integrate(x^-1*(1+log(x))^(1/2)*log(x)^-1, x) == :(-2*arctanh((1+log(x))^(1/2))+2*(1+log(x))^(1/2))
@test integrate((1+x^(1/2))^8, x) == :(-2//9*(1+x^(1/2))^9+1//5*(1+x^(1/2))^10)
@test integrate(sec(x)^4*tan(x)^3, x) == :(-1//4*sec(x)^4+1//6*sec(x)^6)
@test integrate(x*(2+x^2+-2x)^-1, x) == :((1/2)*log(2+x^2+-2x)+-1*arctan(1+-1x))
@test integrate(x*arcsin(x), x) == :(-1//4*arcsin(x)+(1/2)*x^2*arcsin(x)+1//4*x*(1+-1*x^2)^(1/2))
@test integrate(x^-1*(9+-1*x^2)^(1/2), x) == :((9+-1*x^2)^(1/2)+-3*arctanh(1//3*(9+-1*x^2)^(1/2)))
@test integrate(x*(2+x^2+3x)^-1, x) == :(-1*log(1+x)+2*log(2+x))
@test integrate(x^2*cosh(x), x) == :(2*sinh(x)+x^2*sinh(x)+-2*x*cosh(x))
@test integrate((x^4+2*x^2+4x)^-1*(1+x+x^3), x) == :(1//4*log(x^4+2*x^2+4x))
@test integrate((1+sin(x)^2)^-1*cos(x), x) == :(arctan(sin(x)))
@test integrate(cos(x^(1/2)), x) == :(2*cos(x^(1/2))+2*x^(1/2)*sin(x^(1/2)))
@test integrate(sin(pi*x), x) == :(-1*pi^-1*cos(pi*x))
@test integrate((1+exp(x))^-1*exp(2x), x) == :(-1*log(1+exp(x))+exp(x))
@test integrate(cos(5x)*exp(3x), x) == :(3//34*cos(5x)*exp(3x)+5//34*exp(3x)*sin(5x))
@test integrate(cos(3x)*cos(5x), x) == :(1//4*sin(2x)+1//16*sin(8x))
@test integrate((1+x+x^2+x^3)^-1, x) == :((1/2)*arctan(x)+(1/2)*log(1+x)+-1//4*log(1+x^2))
@test integrate(x^2*log(1+x), x) == :(-1//3*x+-1//9*x^3+1//3*log(1+x)+1//6*x^2+1//3*x^3*log(1+x))
@test integrate(x^5*exp(-1*x^3), x) == :(-1//3*exp(-1*x^3)+-1//3*x^3*exp(-1*x^3))
@test integrate(tan(4x)^2, x) == :(-1x+1//4*tan(4x))
@test integrate((-5+9*x^2+12x)^-1//2, x) == :(1//3*arctanh((-5+9*x^2+12x)^-1//2*(2+3x)))
@test integrate(x^2*arctan(x), x) == :(-1//6*x^2+1//6*log(1+x^2)+1//3*x^3*arctan(x))
@test integrate(x^-1//3*(1+-1*x^(1/2)), x) == :(-6//7*x^7//6+3//2*x^2//3)
@test integrate((-1*exp(-1x)+exp(x))^-1, x) == :(-1*arctanh(exp(x)))
@test integrate(x*(10+x^4+2*x^2)^-1, x) == :(1//6*arctan(1//3+1//3*x^2))
@test integrate((x+x^-1//3)^-1, x) == :(3//4*log(1+x^4//3))
@test integrate(cos(x)^4*sin(x)^2, x) == :(1//16*x+-1//6*cos(x)^5*sin(x)+1//16*cos(x)*sin(x)+1//24*cos(x)^3*sin(x))
@test integrate((5+-1*x^2+-4x)^-1//2, x) == :(-1*arcsin(-2//3+-1//3*x))
@test integrate(x*(1+(1+-1*x^2)^(1/2)+-1*x^2)^-1, x) == :(-1*log(1+(1+-1*x^2)^(1/2)))
@test integrate((1+cos(x))*csc(x), x) == :(log(1+-1*cos(x)))
@test integrate((-1+exp(2x))^-1*exp(x), x) == :(-1*arctanh(exp(x)))
@test integrate((-8+x^3)^-1, x) == :(-1//24*log(4+x^2+2x)+1//12*log(2+-1x)+-1//12*3^(1/2)*arctan(1//3*3^(1/2)*(1+x)))
@test integrate(x^5*cosh(x), x) == :(-120*cosh(x)+x^5*sinh(x)+-60*x^2*cosh(x)+-5*x^4*cosh(x)+20*x^3*sinh(x)+120*x*sinh(x))
@test integrate(cos(x)^-1*sin(x)^-1*log(tan(x)), x) == :((1/2)*log(tan(x))^2)
@test integrate(x^2+x^3+-2x, x) == :(-1*x^2+1//3*x^3+1//4*x^4)
@test integrate((1+-1*exp(x))^-1*(1+exp(x)), x) == :(x+-2*log(1+-1*exp(x)))
@test integrate(x*(1+x^2)^-1*(4+x^2)^-1, x) == :(-1//6*log(4+x^2)+1//6*log(1+x^2))
@test integrate((4+-5*sin(x))^-1, x) == :(-1//3*log(-2*sin((1/2)*x)+cos((1/2)*x))+1//3*log(-1*sin((1/2)*x)+2*cos((1/2)*x)))
@test integrate(x*(c+x)^1//3, x) == :(3//7*(c+x)^7//3+-3//4*c*(c+x)^4//3)
@test integrate(exp(x^1//3), x) == :(6*exp(x^1//3)+-6*x^1//3*exp(x^1//3)+3*x^2//3*exp(x^1//3))
@test integrate((4+x+(1+x)^(1/2))^-1, x) == :(-2//11*11^(1/2)*arctan(1//11*11^(1/2)*(1+2*(1+x)^(1/2)))+log(4+x+(1+x)^(1/2)))
@test integrate((x^3+-1*x^2)^-1*(1+x^3), x) == :(x+x^-1+-1*log(x)+2*log(1+-1x))
@test integrate((-3+x^2+4x)*sin(2x), x) == :(7//4*cos(2x)+(1/2)*x*sin(2x)+-2*x*cos(2x)+-1//2*x^2*cos(2x)+sin(2x))
@test integrate(cos(cos(x))*sin(x), x) == :(-1*sin(cos(x)))
@test integrate((16+-1*x^2)^-1//2, x) == :(arcsin(1//4*x))
@test integrate(x^3*(1+x)^-10, x) == :(-3//8*(1+x)^-8+-1//6*(1+x)^-6+1//9*(1+x)^-9+3//7*(1+x)^-7)
@test integrate(cot(2x)^3*csc(2x)^3, x) == :(-1//10*csc(2x)^5+1//6*csc(2x)^3)
@test integrate((x+sin(x))^2, x) == :((1/2)*x+2*sin(x)+1//3*x^3+-2*x*cos(x)+-1//2*cos(x)*sin(x))
@test integrate((1+x^2)^-1*exp(arctan(x)), x) == :(exp(arctan(x)))
@test integrate(x^-1*(1+x^4)^-1, x) == :(-1//4*log(1+x^4)+log(x))
@test integrate(t^3*exp(-2t), t) == :(-3//8*exp(-2t)+-3//4*t*exp(-2t)+-3//4*t^2*exp(-2t)+-1//2*t^3*exp(-2t))
@test integrate(t^(1/2)*(1+t^1//3)^-1, t) == :(-6*t^1//6+2*t^(1/2)+6*arctan(t^1//6)+-6//5*t^5//6+6//7*t^7//6)
@test integrate(sin(x)*sin(2x)*sin(3x), x) == :(-1//8*cos(2x)+-1//16*cos(4x)+1//24*cos(6x))
@test integrate(log((1/2)*x), x) == :(-1x+x*log((1/2)*x))
@test integrate(((1+-1x)^-1*(1+x))^(1/2), x) == :(2*arctan(((1+-1x)^-1*(1+x))^(1/2))+-1*((1+-1x)^-1*(1+x))^(1/2)*(1+-1x))
@test integrate(x*(-1+x^2)^-1//2*log(x), x) == :(-1*(-1+x^2)^(1/2)+(-1+x^2)^(1/2)*log(x)+arctan((-1+x^2)^(1/2)))
@test integrate((a^2+x^2)^-1*(a+x), x) == :((1/2)*log(a^2+x^2)+arctan(x*a^-1))
@test integrate((1+x+-1*x^2)^(1/2), x) == :(-5//8*arcsin(1//5*5^(1/2)*(1+-2x))+-1//4*(1+x+-1*x^2)^(1/2)*(1+-2x))
@test integrate(x^4*(16+x^10)^-1, x) == :(1//20*arctan(1//4*x^5))
@test integrate((2+x+x^2)^-1*(2+x), x) == :((1/2)*log(2+x+x^2)+3//7*7^(1/2)*arctan(1//7*7^(1/2)*(1+2x)))
@test integrate(x*sec(x)*tan(x), x) == :(-1*arctanh(sin(x))+x*sec(x))
@test integrate(x*(x^4+-1*a^4)^-1, x) == :(-1//2*a^-2*arctanh(a^-2*x^2))
@test integrate((x^(1/2)+(1+x)^(1/2))^-1, x) == :(-2//3*x^3//2+2//3*(1+x)^3//2)
@test integrate((1+-1*exp(-1x)+2*exp(x))^-1, x) == :(-1//3*log(1+exp(x))+1//3*log(1+-2*exp(x)))
@test integrate(x^-1//2*arctan(x^(1/2)), x) == :(-1*log(1+x)+2*x^(1/2)*arctan(x^(1/2)))
@test integrate(x^-2*log(1+x), x) == :(-1*log(1+x)+-1*x^-1*log(1+x)+log(x))
@test integrate((-1*exp(x)+exp(3x))^-1, x) == :(-1*arctanh(exp(x))+exp(-1x))
@test integrate((1+-1*cos(x)^2)^-1*(1+cos(x)^2), x) == :(-1x+-2*cot(x))
@test integrate(x^-1*(-25+2x)^-1//2, x) == :(2//5*arctan(1//5*(-25+2x)^(1/2)))
@test integrate((9+-1*cos(x)^4)^-1//2*sin(2x), x) == :(-1*arcsin(1//3*cos(x)^2))
@test integrate(x^2*(5+-4*x^2)^-1//2, x) == :(5//16*arcsin(2//5*x*5^(1/2))+-1//8*x*(5+-4*x^2)^(1/2))
@test integrate(x^3*sin(x), x) == :(-6*sin(x)+-1*x^3*cos(x)+3*x^2*sin(x)+6*x*cos(x))
@test integrate(x*(4+x^2+2x)^(1/2), x) == :(-3//2*arcsinh(1//3*3^(1/2)*(1+x))+1//3*(4+x^2+2x)^3//2+-1//2*(4+x^2+2x)^(1/2)*(1+x))
@test integrate(x*(5+x^2)^8, x) == :(1//18*(5+x^2)^9)
@test integrate(cos(x)^2*sin(x)^5, x) == :(-1//3*cos(x)^3+-1//7*cos(x)^7+2//5*cos(x)^5)
@test integrate(cos(4x)*exp(-3x), x) == :(-3//25*cos(4x)*exp(-3x)+4//25*exp(-3x)*sin(4x))
@test integrate(csc((1/2)*x)^3, x) == :(-1*arctanh(cos((1/2)*x))+-1*cot((1/2)*x)*csc((1/2)*x))
@test integrate(x^-2*(-1+9*x^2)^(1/2), x) == :(3*arctanh(3*x*(-1+9*x^2)^-1//2)+-1*x^-1*(-1+9*x^2)^(1/2))
@test integrate(x^-1*(4+-3*x^2)^(1/2), x) == :((4+-3*x^2)^(1/2)+-2*arctanh((1/2)*(4+-3*x^2)^(1/2)))
@test integrate(x^2*exp(3x), x) == :(2//27*exp(3x)+-2//9*x*exp(3x)+1//3*x^2*exp(3x))
@test integrate((1+sin(x))^-1//2*cos(x)*sin(x), x) == :(-2*(1+sin(x))^(1/2)+2//3*(1+sin(x))^3//2)
@test integrate(x*arcsin(x^2), x) == :((1/2)*(1+-1*x^4)^(1/2)+(1/2)*x^2*arcsin(x^2))
@test integrate(x^3*arcsin(x^2), x) == :(-1//8*arcsin(x^2)+1//4*x^4*arcsin(x^2)+1//8*x^2*(1+-1*x^4)^(1/2))
@test integrate(exp(x)*sech(exp(x)), x) == :(arctan(sinh(exp(x))))
@test integrate(x^2*cos(3x), x) == :(-2//27*sin(3x)+1//3*x^2*sin(3x)+2//9*x*cos(3x))
@test integrate((5+-1*x^2+-4x)^(1/2), x) == :(-9//2*arcsin(-2//3+-1//3*x)+(1/2)*(5+-1*x^2+-4x)^(1/2)*(2+x))
@test integrate(x^5*(2^(1/2)+x^2)^-1, x) == :(1//4*x^4+-1//2*2^(1/2)*x^2+log(2^(1/2)+x^2))
@test integrate(sec(x)^5, x) == :(3//8*arctanh(sin(x))+1//4*sec(x)^3*tan(x)+3//8*sec(x)*tan(x))
@test integrate(sin(2x)^6, x) == :(5//16*x+-5//32*cos(2x)*sin(2x)+-5//48*sin(2x)^3*cos(2x)+-1//12*sin(2x)^5*cos(2x))
@test integrate(sin(x)^2*cos(x)*log(sin(x)), x) == :(-1//9*sin(x)^3+1//3*sin(x)^3*log(sin(x)))
@test integrate((1+2*exp(x))^-1*exp(-1x), x) == :(-1*exp(-1x)+-2x+2*log(1+2*exp(x)))
@test integrate((2+3*cos(x))^(1/2)*tan(x), x) == :(-2*(2+3*cos(x))^(1/2)+2*2^(1/2)*arctanh((1/2)*2^(1/2)*(2+3*cos(x))^(1/2)))
@test integrate(x*(x^2+-4x)^-1//2, x) == :((x^2+-4x)^(1/2)+4*arctanh(x*(x^2+-4x)^-1//2))
@test integrate(cos(x)^5, x) == :(-2//3*sin(x)^3+1//5*sin(x)^5+sin(x))
@test integrate(x^4*exp(-1x), x) == :(-24*exp(-1x)+-1*x^4*exp(-1x)+-24*x*exp(-1x)+-12*x^2*exp(-1x)+-4*x^3*exp(-1x))
@test integrate(x^4*(-2+x^10)^-1//2, x) == :(1//5*arctanh(x^5*(-2+x^10)^-1//2))
@test integrate(cos(4+3x)*exp(x), x) == :(1//10*cos(4+3x)*exp(x)+3//10*exp(x)*sin(4+3x))
@test integrate(x^2*arctan(x), x) == :(-1//6*x^2+1//6*log(1+x^2)+1//3*x^3*arctan(x))
@test integrate((-1+exp(2x))^(1/2), x) == :((-1+exp(2x))^(1/2)+-1*arctan((-1+exp(2x))^(1/2)))
@test integrate(exp(sin(x))*sin(2x), x) == :(-2*exp(sin(x))+2*exp(sin(x))*sin(x))
@test integrate(x^2*(5+-1*x^2)^(1/2), x) == :(25//8*arcsin(1//5*x*5^(1/2))+-5//8*x*(5+-1*x^2)^(1/2)+1//4*x^3*(5+-1*x^2)^(1/2))
@test integrate(x^2*(1+x^3)^4, x) == :(1//15*(1+x^3)^5)
@test integrate(cos(x)^3*sin(x)^3, x) == :(-1//6*sin(x)^6+1//4*sin(x)^4)
@test integrate(sec(x)^4*tan(x)^2, x) == :(1//3*tan(x)^3+1//5*tan(x)^5)
@test integrate(x*(1+2x)^(1/2), x) == :(-1//6*(1+2x)^3//2+1//10*(1+2x)^5//2)
@test integrate(sin(x)^4, x) == :(3//8*x+-3//8*cos(x)*sin(x)+-1//4*sin(x)^3*cos(x))
@test integrate(tan(x)^3, x) == :((1/2)*tan(x)^2+log(cos(x)))
@test integrate(x^5*(1+x^2)^(1/2), x) == :(-2//5*(1+x^2)^5//2+1//3*(1+x^2)^3//2+1//7*(1+x^2)^7//2)

max_stars_repo_path: test/0 Independent test suites/Stewart Problems.jl


text:
"""
$SIGNATURES

Specify the folder for the Literate scripts, by default this is `scripts/`.
"""
function literate_folder(rp::String="")
    isempty(rp) && return PATHS[:literate]
    path = joinpath(PATHS[:folder], rp)
    !isdir(path) && error("Specified Literate path not found ($rp -- $path)")
    PATHS[:literate] = path
    return path
end

#
# Convenience functions to work with strings and substrings
#

"""
$(SIGNATURES)

Returns the string corresponding to `s`: `s` itself if it is a string, or the parent string if `s`
is a substring. Do not confuse with `String(s::SubString)` which casts `s` into its own object.

# Example

```julia-repl
julia> a = SubString("hello JuDoc", 3:8);
julia> JuDoc.str(a)
"hello JuDoc"
julia> String(a)
"llo Ju"
```
"""
str(s::String)::String    = s
str(s::SubString)::String = s.string

"""
    subs(s, from, to)
    subs(s, from)
    subs(s, range)
    subs(s)

Convenience functions to take a substring of a string.

# Example
```julia-repl
julia> JuDoc.subs("hello", 2:4)
"ell"
```
"""
subs(s::AS, from::Int, to::Int)::SubString    = SubString(s, from, to)
subs(s::AS, from::Int)::SubString             = subs(s, from, from)
subs(s::AS, range::UnitRange{Int})::SubString = SubString(s, range)
subs(s::AS) = SubString(s)

"""
$(SIGNATURES)

Given a substring `ss`, returns the position in the parent string where the substring starts.

# Example
```julia-repl
julia> ss = SubString("hello", 2:4); JuDoc.from(ss)
2
```
"""
from(ss::SubString)::Int = nextind(str(ss), ss.offset)

"""
$(SIGNATURES)

Given a substring `ss`, returns the position in the parent string where the substring ends.

# Example
```julia-repl
julia> ss = SubString("hello", 2:4); JuDoc.to(ss)
4
```
"""
to(ss::SubString)::Int = ss.offset + ss.ncodeunits

"""
$(SIGNATURES)

Returns the string span of a regex match. Assumes there is no unicode in the match.

# Example
```julia-repl
julia> JuDoc.matchrange(match(r"ell", "hello"))
2:4
```
"""
matchrange(m::RegexMatch)::UnitRange{Int} = m.offset .+ (0:(length(m.match)-1))

"""
$SIGNATURES

In an lxdef or a div, ignore any whitespace at the beginning of lines to avoid ambiguities with indented
code blocks.
"""
function ignore_starting_line_spaces(content::SubString)::AS
    s = strip(content)
    return replace(s, r"\n\s*" => "\n ")
end

# Other convenience functions

"""
$(SIGNATURES)

Convenience function to display a time since `start`.
"""
function time_it_took(start::Float64)
    comp_time = time() - start
    mess = comp_time > 60 ? "$(round(comp_time/60;   digits=1))m" :
           comp_time >  1 ? "$(round(comp_time;      digits=1))s" :
                            "$(round(comp_time*1000; digits=1))ms"
    return "[done $(lpad(mess, 6))]"
end

```
$(SIGNATURES)

Nicer printing of processes.
```
function print_final(startmsg::AS, starttime::Float64)::Nothing
    tit = time_it_took(starttime)
    println("\r$startmsg$tit")
end

"""
$(SIGNATURES)

Convenience function to denote a string as being in a math context in a recursive parsing
situation. These blocks will be processed as math blocks but without adding KaTeX elements to it
given that they are part of a larger context that already has KaTeX elements.
NOTE: this happens when resolving latex commands in a math environment. So for instance if you have
`\$\$ x \\in \\R \$\$` and `\\R` is defined as a command that does `\\mathbb{R}` well that would be
an embedded math environment. These environments are marked as such so that we don't add additional
KaTeX markers around them.
"""
mathenv(s::AS)::String = "_\$>_$(s)_\$<_"

"""
$(SIGNATURES)

Takes a string `s` and replace spaces by underscores so that that we can use it
for hyper-references. So for instance `"aa  bb"` will become `aa_bb`.
It also defensively removes any non-word character so for instance `"aa bb !"` will be `"aa_bb"`
"""
function refstring(s::AS)::String
    # remove html tags
    st = replace(s, r"<[a-z\/]+>"=>"")
    # remove non-word characters
    st = replace(st, r"&#[0-9]+;" => "")
    st = replace(st, r"[^a-zA-Z0-9_\-\s]" => "")
    # replace spaces by dashes
    st = replace(lowercase(strip(st)), r"\s+" => "_")
    # to avoid clashes with numbering of repeated headers, replace
    # double underscores by a single one (see convert_header function)
    st = replace(st, r"__" => "_")
    # in the unlikely event we don't have anything here, return the hash
    # of the original string
    isempty(st) && return string(hash(s))
    return st
end

"""
$(SIGNATURES)

Internal function to take a path and return a unix version of the path (if it isn't already).
Used in [`resolve_assets_rpath`](@ref).
"""
function unixify(rp::String)::String
    cand = Sys.isunix() ? rp : replace(rp, "\\"=>"/")
    isempty(splitext(cand)[2]) || return cand
    endswith(cand, "/") || isempty(cand) || return cand * "/"
    return cand
end

"""
$(SIGNATURES)
Internal function to take a unix path, split it along `/` and re-join it (which will lead to the
same path on unix but not on windows). Only used in [`resolve_assets_rpath`](@ref).
"""
joinrp(rpath::AS) = joinpath(split(rpath, '/')...)

"""
$(SIGNATURES)

Internal function to resolve a relative path. See [`convert_code_block`](@ref) and
[`resolve_lx_input`](@ref).
As an example, let's say you have the path `./blah/blih.txt` somewhere in `src/pages/page1.md`.
The `./` indicates that the path should be reproduced in `/assets/` and so it will lead to
`/assets/pages/blah/blih.txt`.
In the `canonical=false` mode, the returned path is a UNIX path (irrelevant of your platform);
these paths can be given to html anchors (refs, img, ...).
In the `canonical=true` mode, the path is a valid path on the local system. These paths can be
given to Julia to read or write things.
"""
function resolve_assets_rpath(rpath::AS; canonical::Bool=false, code::Bool=false)::String
    @assert length(rpath) > 1 "relative path '$rpath' doesn't seem right"
    if startswith(rpath, "/")
        # this is a full path starting from the website root folder so for instance
        # `/assets/blah/img1.png`; just return that
        canonical || return rpath
        return normpath(joinpath(PATHS[:folder], joinrp(rpath[2:end])))
    elseif startswith(rpath, "./")
        # this is a path relative to the assets folder with the same path as the calling file so
        # for instance if calling from `src/pages/pg1.md` with `./im1.png` it would refer to
        # /assets/pages/pg1/im1.png
        @assert length(rpath) > 2 "relative path '$rpath' doesn't seem right"
        canonical || return "/assets/" * unixify(splitext(CUR_PATH[])[1]) * rpath[3:end]
        return normpath(joinpath(PATHS[:assets], splitext(CUR_PATH[])[1], joinrp(rpath[3:end])))
    end
    if code
        # in the code mode we allow a short, this: `julia:ex` is considered
        # short for `julia:./code/ex`
        return resolve_assets_rpath("./code/" * rpath; canonical=canonical)
    end
    # this is a full path relative to the assets folder for instance `blah/img1.png` would
    # correspond to `/assets/blah/img1.png`
    canonical || return "/assets/" * rpath
    return normpath(joinpath(PATHS[:assets], joinrp(rpath)))
end

ext: jl


text:
using Test
using StaticArrays, JuMP
using Polyhedra

@testset "Default" begin
    solver = MOI.Utilities.MockOptimizer(MOIU.Model{Float64}())
    @testset "Library" begin
        function default_library_test(d)
            @test default_library(d(2), Float64) isa Polyhedra.DefaultLibrary
            @test default_library(d(1), Float64) isa IntervalLibrary
            @test default_library(d(1), Int) isa IntervalLibrary
            @test default_library(d(1), BigInt) isa IntervalLibrary
            @test default_library(d(2), AbstractFloat) isa Polyhedra.DefaultLibrary{Float64}
            @test default_library(d(1), AbstractFloat) isa IntervalLibrary{Float64}
        end
        default_library_test(N -> StaticArrays.Size((N,)))
        default_library_test(identity)
    end
    @testset "LP Solver" begin
        @test Polyhedra.default_solver(hrep([HalfSpace((@SVector [1, 2]), 3)])) === nothing
        @test Polyhedra.linear_objective_solver(hrep([HalfSpace((@SVector [1, 2]), 3)]), solver) === solver
        lib = Polyhedra.DefaultLibrary{Float64}(solver)
        p = polyhedron(convexhull([0.0, 0.0]), lib)
        @test library(p).solver === solver
    end
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
function one_phase_solve(m::JuMP.Model)
    nlp_raw = MathOptNLPModel(m);
    return one_phase_solve(nlp_raw)
end

function one_phase_solve(m::JuMP.Model, pars::Class_parameters)
    nlp_raw = MathOptNLPModel(m);
    return one_phase_solve(nlp_raw, pars::Class_parameters)
end

function one_phase_solve(nlp_raw::NLPModels.AbstractNLPModel)
    pars = Class_parameters()
    return one_phase_solve(nlp_raw, pars)
end

function one_phase_solve(nlp_raw::NLPModels.AbstractNLPModel, pars::Class_parameters)
    nlp = Class_CUTEst(nlp_raw);
	if ncon(nlp) == 0
		throw(ErrorException("Unconstrained minimization problems are unsupported"))
	end
    timer = class_advanced_timer()
    start_advanced_timer(timer)
    start_advanced_timer(timer, "INIT")
    if pars.init.init_style == :gertz
        intial_it = gertz_init(nlp, pars, timer); # Gertz, Michael, Jorge Nocedal, and A. Sartenar. "A starting point strategy for nonlinear interior methods." Applied mathematics letters 17.8 (2004): 945-952.
    elseif pars.init.init_style == :mehrotra
        intial_it = mehrotra_init(nlp, pars, timer);
    elseif pars.init.init_style == :LP
        intial_it = LP_init(nlp, pars, timer);
    else
        error("Init strategy does not exist")
    end
    pause_advanced_timer(timer, "INIT")
    pause_advanced_timer(timer)
    if pars.output_level >= 4
        print_timer_stats(timer)
    end
    start_advanced_timer(timer)
    @assert(is_feasible(intial_it, pars.ls.comp_feas))
    iter, status, hist, t, err = one_phase_IPM(intial_it, pars, timer);

    pause_advanced_timer(timer)

    if pars.output_level >= 3
        print_timer_stats(timer)
    end
    return iter, status, hist, t, err, timer
end

function switching_condition(iter::Class_iterate, last_step_was_superlinear::Bool, pars::Class_parameters)
    # should we take an aggressive step or not?
    is_feas = is_feasible(iter, pars.ls.comp_feas_agg)
    dual_avg = scaled_dual_feas(iter, pars)

    if pars.primal_bounds_dual_feas
      dual_progress = dual_avg < pars.aggressive_dual_threshold * LinearAlgebra.norm(get_primal_res(iter), Inf)
    else
      dual_progress = dual_avg < pars.aggressive_dual_threshold * get_mu(iter)
    end
    delta_small = get_delta(iter) < sqrt(get_mu(iter)) * max(0.1, LinearAlgebra.norm(get_y(iter),Inf))
    lag_grad = LinearAlgebra.norm(eval_grad_lag(iter,get_mu(iter)),1) < sum(iter.point.s .* iter.point.y) + LinearAlgebra.norm(get_grad(iter) + iter.point.mu * eval_grad_r(iter),1) # + LinearAlgebra.norm(get_primal_res(iter), Inf) + 1.0 #+ sqrt(LinearAlgebra.norm(get_y(iter),Inf))

    be_aggressive = is_feas && dual_progress && lag_grad
    be_aggressive |= last_step_was_superlinear && dual_progress && lag_grad

    return be_aggressive
end

function one_phase_IPM(iter::Class_iterate, pars::Class_parameters, timer::class_advanced_timer)
  #####################################################################
  # THE MAIN ALGORITHM
  # input:
  # iter = starting point
  # pars = parameters for running the algorithm
  # timer = code to time the algorithm
  #####################################################################

      # intialize code
      t = 0;
      rpt = 0.0
      progress = Array{alg_history2,1}();
      filter = Array{Class_filter,1}();

      update!(iter, timer, pars) # is this necessary ????

      kkt_solver = pick_KKT_solver(pars);
      initialize!(kkt_solver, iter)

      if pars.output_level >= 1
          head_progress()
      end

      display = pars.output_level >= 1
      record_progress_first_it!(progress, iter, kkt_solver, pars, display)
      if pars.output_level >= 4
          println("")
      end

      init_step_size = 1.0
      status = :success
      agg_next_step = false;
      dir_size_agg = Inf
      dir_size_stable = 0.0
      ls_info = false

      start_time = time()
      last_step_was_superlinear = false
      scale_update = false

      # check termination critieron at starting point
      start_advanced_timer(timer, "misc/terminate")
      status = terminate(iter, pars)
      pause_advanced_timer(timer, "misc/terminate")

      if status != false
        if pars.output_level >= 1
            println("Terminated with ", status)
        end
        return iter, status, progress, t, false
      end

      if time() - start_time > pars.term.max_time
          if pars.output_level >= 1
            println("Terminated due to timeout")
        end
        return iter, :MAX_TIME, progress, t, false
      end

      # run the algorithm
      for t = 1:pars.term.max_it
             @assert(is_feasible(iter, pars.ls.comp_feas))

             for i = 1:pars.max_it_corrections
                         if pars.output_level >= 4
                             println("================================== ITERATION $t, MINOR ITERATION $i ======================================")
                         end

                         if pars.output_level >= 5
                             println("Strict comp = ", maximum(min.(iter.point.s/LinearAlgebra.norm(iter.point.s,Inf),iter.point.y/LinearAlgebra.norm(iter.point.y,Inf))))
                         end

                       tot_num_fac = 0; inertia_num_fac = 0;

                       start_advanced_timer(timer, "misc/checks")
                       be_aggressive = switching_condition(iter, last_step_was_superlinear, pars) # should we take an aggressive step or not?
                       last_step_was_superlinear = false

                       pause_advanced_timer(timer, "misc/checks")

                       if i == 1
                           ##########################################################################################
                           ##### first step of iteration, we need to compute lag hessian and do a factorization #####
                           ##########################################################################################
                           update_H!(iter, timer, pars)
                           @assert(is_updated(iter.cache))

                           # form matrix that we are going to factorize
                           form_system!(kkt_solver, iter, timer)

                           start_advanced_timer(timer, "ipopt_strategy")
                           # Figures out what delta will give
                           fact_succeed, inertia_num_fac, new_delta = ipopt_strategy!(iter, kkt_solver, pars, timer)
                           tot_num_fac = inertia_num_fac
                           old_delta = get_delta(iter)
                           set_delta(iter, new_delta)

                           pause_advanced_timer(timer, "ipopt_strategy")

                           if fact_succeed != :success
                              return iter, :MAX_DELTA, progress, t, false
                           end

                           start_advanced_timer(timer, "STEP")
                           start_advanced_timer(timer, "STEP/first")

                           if pars.output_level >= 5
                             println(pd("**"), pd("status"), pd("delta"), pd("a_P"), pd("dx"), pd("dy"), pd("ds"))
                           end

                           for k = 1:100
                             #status, new_iter, ls_info = take_step!(iter, reduct_factors, kkt_solver, ls_mode, filter, pars, actual_min_step_size, timer)
                             step_status, new_iter, ls_info, reduct_factors = take_step2!(be_aggressive, iter, kkt_solver, filter, pars, timer)

                             if pars.output_level >= 6
                               println("")
                               println(pd("**"), pd(step_status), rd(get_delta(iter)), rd(ls_info.step_size_P), rd(LinearAlgebra.norm(kkt_solver.dir.x,Inf)), rd(LinearAlgebra.norm(kkt_solver.dir.y,Inf)), rd(LinearAlgebra.norm(kkt_solver.dir.s,Inf)))
                             end

                             if step_status == :success
                               break
                             elseif i < 100 && get_delta(iter) < pars.delta.max
                               if pars.test.response_to_failure == :lag_delta_inc
                                 set_delta(iter, max(LinearAlgebra.norm(eval_grad_lag(iter,iter.point.mu),Inf) / LinearAlgebra.norm(kkt_solver.dir.x,Inf),get_delta(iter) * pars.delta.inc, max(pars.delta.start, old_delta * pars.delta.dec)))
                               elseif pars.test.response_to_failure == :default
                                 set_delta(iter, max(get_delta(iter) * pars.delta.inc, max(pars.delta.start, old_delta * pars.delta.dec)))
                               else
                                  error("pars.test.response_to_failure parameter incorrectly set")
                               end
                               inertia = factor!(kkt_solver, get_delta(iter), timer)
                               tot_num_fac += 1
                             elseif LinearAlgebra.norm(comp(iter),Inf) > 1e-14
                                 warn("Error ... large delta causing issues")
                                iter.point.y = iter.point.mu ./ iter.point.s
                                step_status = :success
                                break
                             else
                               pause_advanced_timer(timer, "STEP/first")
                               pause_advanced_timer(timer, "STEP")
                               println("Terminated due to max delta while attempting to take step")
                               println("delta=$(get_delta(iter)), be_aggressive=$be_aggressive, status=$step_status")
                               println("dx = $(LinearAlgebra.norm(kkt_solver.dir.x,2)), dy = $(LinearAlgebra.norm(kkt_solver.dir.y,2)), ds = $(LinearAlgebra.norm(kkt_solver.dir.s,2))")
                               @show reduct_factors #, ls_mode
                               @show ls_info
                               return iter, :MAX_DELTA, progress, t, false
                             end
                           end

                           pause_advanced_timer(timer, "STEP/first")
                           pause_advanced_timer(timer, "STEP")
                       else
                           #######################################
                           ### corrections, reuse factorization ##
                           #######################################
                           start_advanced_timer(timer, "STEP")
                           start_advanced_timer(timer, "STEP/correction")

                           step_status, new_iter, ls_info, reduct_factors = take_step2!(be_aggressive, iter, kkt_solver, filter, pars, timer)

                           if pars.superlinear_theory_mode && be_aggressive
                             if get_mu(new_iter) < get_mu(iter) * 0.1
                                last_step_was_superlinear = true
                             end
                           end

                           pause_advanced_timer(timer, "STEP/correction")
                           pause_advanced_timer(timer, "STEP")
                       end

                       if step_status == :success
                         iter = new_iter
                         if be_aggressive
                           dir_size_agg = LinearAlgebra.norm(kkt_solver.dir.x, 2)
                         end
                       end

                       add!(filter, iter, pars) # update filter

                       start_advanced_timer(timer, "misc/terminate")
                       status = terminate(iter, pars) # check termination criteron
                       pause_advanced_timer(timer, "misc/terminate")

                       start_advanced_timer(timer, "misc/record_progress")
                       # output to the console
                       output_level = pars.output_level
                       display = output_level >= 4 || (output_level >= 3 && i == 1) || (output_level == 2 && t % 10 == 1 && i == 1) || (status != false && output_level >= 1)

                       record_progress!(progress, t, be_aggressive ? "agg" : "stb", iter, kkt_solver, ls_info, reduct_factors, inertia_num_fac, tot_num_fac, pars, display)
                       if pars.output_level >= 4
                           println("")
                       end

                       pause_advanced_timer(timer, "misc/record_progress")
                       @assert(is_updated_correction(iter.cache))
                       check_for_nan(iter.point)

                       # if termination criteron is satisfied stop the algorithm
                       if status != false
                         if pars.output_level >= 1
                             println("Terminated with ", status)
                         end
                         return iter, status, progress, t, false
                       end

                       if time() - start_time > pars.term.max_time
                         if pars.output_level >= 1
                             println("Terminated due to timeout")
                         end
                         return iter, :MAX_TIME, progress, t, false
                       end

                       if step_status != :success
                         break
                       end
             end
      end

  if pars.output_level >= 1
      println("Terminated due to max iterations reached")
  end
  return iter, :MAX_IT, progress, pars.term.max_it, false
end


text:
#=
Template for branch-and-cut method

July 5, 2018
Kibaek Kim
Brian Dandurand
=#

include("../../MaximinOPF/src/MaximinOPF.jl")
include("../../MaximinOPF/src/utils.jl")
using JuMP, MathOptInterface
using Mosek, MosekTools
using CPLEX
using JuMP
using PowerModels
PowerModels.silence()

function solveNodeSP(model_info, ndata; compute_psd_dual=false, io=devnull)
    model = model_info["model"]
    branch_ids = model_info["branch_ids"]

    unfix_vars(model,branch_ids)
    for l in ndata["inactive_branches"]
        fix(variable_by_name(model,"x[$l]_1"), 1; force=true)
    end
    for l in ndata["protected_branches"]
        fix(variable_by_name(model,"x[$l]_1"), 0; force=true)
    end

    JuMP.optimize!(model)
    sg_info = Dict{Tuple{Int64,Int64},Array{Float64,1}}()
    model_info["opt_val"]=JuMP.objective_value(model)
    model_info["solve_status"]=JuMP.termination_status(model)

    model_info["x_soln_str"]=""
    for l in branch_ids
        x_var = variable_by_name(model,"x[$l]_1")
        x_val = JuMP.value(x_var)
        if x_val > 1.0-1.0e-8
	        x_val = 1
            model_info["x_soln_str"] = string(model_info["x_soln_str"]," $l")
        elseif x_val < 1.0e-8
	        x_val = 0
        end
        model_info["x_soln"][l] = x_val
    end
    ndata["x_soln"] = copy(model_info["x_soln"])
    ndata["bound_value"] = model_info["opt_val"]
    return sg_info
end

function applyPrimalHeuristic(model_info)
    model = model_info["model"]
    branch_ids = model_info["branch_ids"]
    unfix_vars(model,branch_ids)
    idx=sortperm(sort(collect(keys(model_info["x_soln"]))), by=kk->model_info["x_soln"][kk]; rev=true)
    n_idx = length(idx)
    n_attack = min(model_info["attacker_budget"],n_idx)
    model_info["x_soln_str"]=""
    for kk in idx
        l=idx[kk]
        if kk <= n_attack && model_info["x_soln"][l] > 0
            model_info["heur_x_soln"][l] = 1
            fix(variable_by_name(model,"x[$l]_1"), 1; force=true)
            model_info["x_soln_str"] = string(model_info["x_soln_str"]," $l")
        else
            model_info["heur_x_soln"][l] = 0
            fix(variable_by_name(model,"x[$l]_1"), 0; force=true)
        end
    end
    JuMP.optimize!(model)
    model_info["heur_opt_val"]=JuMP.objective_value(model)
    model_info["solve_status"]=JuMP.termination_status(model)
    return model_info["heur_opt_val"]
end

function findNextNode(BnBTree::Dict{Int64,Dict{String,Any}})
  nodekey=-1
  weakestUBVal = -1e20
  for kk in keys(BnBTree)
    if BnBTree[kk]["bound_value"] > weakestUBVal
        nodekey = kk
        weakestUBVal = BnBTree[kk]["bound_value"]
    end
  end
  return nodekey,weakestUBVal
end

function findNextIndex(undecided_branches,x_vals)
  return reduce((k1, k2) -> abs(x_vals[k1]-0.5) <= abs(x_vals[k2]-0.5) ? k1 : k2, undecided_branches)
end

function printNode(node::Dict{String,Any}; pretext="",posttext="")
    println(pretext,"Node ",node["node_key"]," has bound ",node["bound_value"],": inactive=",node["inactive_branches"]," protected: ",node["protected_branches"],posttext)
end

function deleteNodesByBound(BnBTree,bdval)
    for bb in keys(BnBTree)
        if BnBTree[bb]["bound_value"] <= bdval
            delete!(BnBTree,bb)
        end
    end
end

function solveMaxminViaBnB(pm_data,pm_form; use_dual_minmax=true)
    #global MAX_TIME
    K=pm_data["attacker_budget"]
    art_bd=10

    start_time = time_ns()
    init_node=Dict("inactive_branches"=>[], "protected_branches"=>[], "bound_value"=>1e20, "attacker_budget"=>K, "node_key"=>0) ### CREATE ROOT NODE

    base_maxmin = MaximinOPF.MaximinOPFModel(pm_data, pm_form; enforce_int=false, rm_rsoc=true, rm_therm_line_lim=false)
    branch_ids=sort(collect(pm_data["undecided_branches"]))
    println("Branch_ids: ",branch_ids)

    psd_optimizer=with_optimizer(Mosek.Optimizer,MSK_IPAR_LOG=0,MSK_IPAR_NUM_THREADS=8)
    JuMP.set_optimizer(base_maxmin,psd_optimizer)

    psd_model_info=Dict{String,Any}("branch_ids"=>branch_ids,"x_soln"=>Dict{Int64,Float64}(),
            "x_soln_str"=>"","heur_x_soln"=>Dict{Int64,Float64}(),"attacker_budget"=>K)
    for l in branch_ids
        psd_model_info["x_soln"][l]=0
    end
    psd_model_info["model"] = base_maxmin

    BnBTree = Dict{Int64,Dict{String,Any}}()
    BnBTree[0] = init_node ### ADD ROOT NODE TO TREE
    nodekey=0

    feasXs = Dict{String,Dict{String,Any}}()
    feasXs[""]=Dict("x_soln"=>Dict{Int64,Float64}(),"x_soln_str"=>"", "bound_value"=>0, "value"=>0) ### CREATE INITIAL INCUMBENT SOLN
    IncX = feasXs[""]
    nXs = 1

    bestUBVal=1e20
    println("Initial processing of root node yields an upper bound of : ",bestUBVal)

    nNodes=1
    maxidx=1
    maxval=-1

    while true
        incumbent_update=""
        currNode = pop!(BnBTree,nodekey)
        if currNode["bound_value"] <= IncX["value"]
            println("\nFathoming node $nodekey due to bound ",currNode["bound_value"]," <= ",IncX["value"])
        else
            printNode(currNode;pretext="\n",posttext="")
            inn_while_cntr = 0
            while true
                solveNodeSP(psd_model_info, currNode; compute_psd_dual=false)

                currNode["bound_value"] = min(psd_model_info["opt_val"],currNode["bound_value"])
                if currNode["bound_value"] <= IncX["value"]
                    println("\tFathoming due to bound ",currNode["bound_value"]," <= ",IncX["value"])
                    break
                end
	            idx=findNextIndex(psd_model_info["branch_ids"],psd_model_info["x_soln"])
                if psd_model_info["x_soln"][idx] == 0 || psd_model_info["x_soln"][idx] == 1
                    x_soln_str = psd_model_info["x_soln_str"]
                    if !haskey(feasXs,x_soln_str)
	                    nXs += 1
                        feasXs[x_soln_str]=Dict("bound_value"=>psd_model_info["opt_val"], "value"=>psd_model_info["opt_val"])
                        feasXs[x_soln_str]["x_soln"] = copy(psd_model_info["x_soln"])
                        feasXs[x_soln_str]["x_soln_str"] = x_soln_str
                        println("New solution has known optimal value: ",feasXs[x_soln_str]["value"])
                    end
	                println("\tFathoming due to optimality. Adding new attack solution: ",currNode["inactive_branches"])
                    feasXs[x_soln_str]["value"] = psd_model_info["opt_val"]
	                if IncX["value"] < feasXs[x_soln_str]["value"]
	                    IncX = feasXs[x_soln_str]
	                    incumbent_update = string("\t***New incumbent***",IncX["x_soln_str"])
                        deleteNodesByBound(BnBTree,IncX["value"])
	                end
                    break
                else
	                # "BRANCH"
                    bd_val=currNode["bound_value"]
	                println("\tAdding two new nodes, the branching index is: ",idx)
                    node0=Dict{String,Any}("node_key"=>nNodes, "bound_value"=>bd_val,"attacker_budget"=>currNode["attacker_budget"],
                        "inactive_branches"=>copy(currNode["inactive_branches"]), "protected_branches"=>copy(currNode["protected_branches"])
                    )
                    push!(node0["protected_branches"],idx)
	                BnBTree[nNodes]=node0
                    nNodes += 1

                    node1=Dict{String,Any}("node_key"=>nNodes, "bound_value"=>bd_val,"attacker_budget"=>currNode["attacker_budget"],
                        "inactive_branches"=>copy(currNode["inactive_branches"]), "protected_branches"=>copy(currNode["protected_branches"])
                    )
                    push!(node1["inactive_branches"],idx)
	                node1["attacker_budget"] -= 1 ### This node will have another branch inactive, debiting from the inherited budget
	                BnBTree[nNodes]=node1
                    nNodes += 1
                    break
                end
            end ### "while loop"
            println("\tNew node bound value: ",currNode["bound_value"])
        end ### "if initial fathoming due to bound"

        if length(BnBTree) > 0
            nodekey,bestUBVal=findNextNode(BnBTree)
        else
            bestUBVal = IncX["value"]
            break
        end
#=
        if (time_ns()-start_time)/1e9 > MAX_TIME
            break
        end
=#
        tree_report=string(" Nodes left: ",length(BnBTree)," out of ", nNodes,)
        println("\tBound gap: [",IncX["value"],", ",bestUBVal,"]",incumbent_update,",\t",tree_report)
        if bestUBVal < IncX["value"]
            println("Terminating since all other nodes will be fathomed by bound")
            break
        end
    end # while

    end_time = time_ns()
    runtime = (end_time-start_time)/1e9

    println("Final best solution: ",IncX["x_soln_str"]," with value ",IncX["value"])
    println("Runtime: ",runtime)
    println("Number of nodes processed: ",nNodes)
    return bestUBVal,nNodes,runtime
end

testcase = Dict(
	"file" => "data/case30.m",
 	"name" => "case30K4",
 	"attack_budget" => 4,
 	"inactive_indices" => [],
 	"protected_indices" => []
	)

pm_data = PowerModels.parse_file(testcase["file"])
pm_data["attacker_budget"] = testcase["attack_budget"] ###Adding another key and entry
pm_data["inactive_branches"] = testcase["inactive_indices"] ###Adding another key and entry
pm_data["protected_branches"] = testcase["protected_indices"] ###Adding another key and entry

solveMaxminViaBnB(pm_data,SparseSDPWRMPowerModel,use_dual_minmax=true)
#solveMaxminViaBnB(pm_data,SDPWRMPowerModel,use_dual_minmax=true)


text:
@testset "Random Symbols" begin
  @test isa(rand(Symbol), Symbol)
  @test length(rand(Symbol) |> string) == 8
  @test length(rand(Symbol, 10) |> string) == 10
end


text:
using Profile, PProf
using Pkg;Pkg.activate(".")
using PairwiseVelocities
using BenchmarkTools
using DelimitedFiles
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--snapshot"
            help = "an option with an argument"
            arg_type = Int
        "--nd"
            help = "an option with an argument"
            arg_type = Int
        "--min_run"
            help = "an option with an argument"
            arg_type = Int
        "--max_run"
            help = "another option with an argument"
            arg_type = Int
    end

    return parse_args(s)
end

parsed_args = parse_commandline()

DATA_DIR = "/cosma7/data/dp004/dc-cues1/DarkQuest/pairwise_velocities/"
snapshot = parsed_args["snapshot"]
boxsize = 2000.
r_max = 80.

number_densities = readdlm("/cosma6/data/dp004/dc-cues1/DarkQuest/xi/log10density_table.dat", ' ', Float32, '\n')
number_density_left = number_densities[parsed_args["nd"],1]
number_density_right = number_densities[parsed_args["nd"],2]
println("number densities = ( ", number_density_left, " , ",number_density_right," ) ")
rbins_c = readdlm("/cosma6/data/dp004/dc-cues1/DarkQuest/xi/separation.dat", '\t', Float64, '\n')
rbins_c = rbins_c[rbins_c .< r_max]

rbins = zeros(Float64, length(rbins_c)+1)

rbins[2:end-1] = rbins_c[1:end-1] + diff(rbins_c)/2.
rbins[1] = 2. *rbins_c[1] - rbins[2]
rbins[end] = 2. *rbins_c[end] - rbins[end-1]

r_max = maximum(rbins)
rbins_c = 0.5*(rbins[2:end] + rbins[1:end-1])

for run in parsed_args["min_run"]:parsed_args["max_run"]
    filename = "run$(run)_nd_$(abs(number_density_left))_$(abs(number_density_right))_snapshot$(snapshot).csv"
    println(filename)
    if number_density_left == number_density_right
        positions, velocities = read_data(
            run, snapshot, 10^number_density_left, boxsize
        )
        positions = convert(Array{Float64}, positions)
        velocities = convert(Array{Float64}, velocities)
        @time moments = PairwiseVelocities.compute_pairwise_velocity_moments(
                            positions,
                            velocities,
                            rbins,
                            boxsize,
        )
    else
        positions_left, velocities_left = read_data(
            run, snapshot, 10^number_density_left, boxsize
        )
        positions_left = convert(Array{Float64}, positions_left)
        velocities_left = convert(Array{Float64}, velocities_left)
        positions_right, velocities_right = read_data(
            run, snapshot, 10^number_density_right, boxsize
        )
        positions_right = convert(Array{Float64}, positions_right)
        velocities_right = convert(Array{Float64}, velocities_right)

        @time moments = PairwiseVelocities.compute_pairwise_velocity_moments(
                            positions_left,
                            velocities_left,
                            positions_right,
                            velocities_right,
                            rbins,
                            boxsize,
        )
    end
    open(DATA_DIR * filename; write=true) do f
        write(f, "# r_c v_r sigma_r sigma_t skewness_r skewness_rt kurtosis_r kurtosis_t kurtosis_rt \n")
        writedlm(f,
                  zip(rbins_c, moments[2][:],
                    moments[3][:], moments[6][:],
                    moments[4][:], moments[7][:],
                    moments[5][:], moments[8][:], moments[9][:]))
    end

    println("Wrote file for $(run)!")
end


text:
using Distributions
export SpfMcmcSampler
export SpfInvolutiveMcmcCentSubseq, SpfInvolutiveMcmcEdit

abstract type SpfMcmcSampler end

struct SpfInvolutiveMcmcEdit <: SpfMcmcSampler
    ν::Int
    desired_samples::Int
    burn_in::Int
    lag::Int
    par_info::Dict
    curr::Vector{Int} # Storage for current val
    prop::Vector{Int} # Storage for proposed val
    ind_del::Vector{Int} # Storage for indexing of deletions
    ind_add::Vector{Int} # Storage for indexing of additions
    vals::Vector{Int} # Storgae for new values to insert
    function SpfInvolutiveMcmcEdit(
        ;ν=4, desired_samples=1000, burn_in=0, lag=1
    )
    # req_samples = burn_in + 1 + (desired_samples - 1) * lag
    curr = Int[]
    prop = Int[]
    ind_del = zeros(Int, ν)
    ind_add = zeros(Int, ν)
    vals = zeros(Int, ν)
    par_info = Dict()
    par_info[:ν] = "(maximum number of edit operations)"
    new(
        ν,
        desired_samples, burn_in, lag, par_info,
        curr, prop, ind_del, ind_add, vals
        )
    end
end

function Base.show(io::IO, sampler::SpfInvolutiveMcmcEdit)
    title = "MCMC Sampler for Spherical Path Family (SPF) via Edit Operations"
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    println(io, "Parameters:")
    num_of_pars = 1
    for par in fieldnames(typeof(sampler))[1:num_of_pars]
        println(io, par, " = $(getfield(sampler, par))  ", sampler.par_info[par])
    end
    println(io, "\nDefault output parameters:")
    for par in fieldnames(typeof(sampler))[(num_of_pars+1):(num_of_pars+3)]
        println(io, par, " = $(getfield(sampler, par))  ")
    end
end

struct SpfInvolutiveMcmcCentSubseq <: SpfMcmcSampler
    p::Real
    ν::Int
    desired_samples::Int
    burn_in::Int
    lag::Int
    par_info::Dict
    function SpfInvolutiveMcmcCentSubseq(
        ;p=0.6, ν=4, desired_samples=1000, burn_in=0, lag=1
    )
    par_info = Dict()
    par_info[:p] = "(preserved subseq size)"
    par_info[:ν] = "(diff in path length)"
    new(p, ν, desired_samples, burn_in, lag, par_info)
    end
end

function Base.show(io::IO, sampler::SpfInvolutiveMcmcCentSubseq)
    title = "MCMC Sampler for Spherical Path Family (SPF) via Subsequence Preservation"
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    println(io, "Parameters:")
    num_of_pars = 2
    for par in fieldnames(typeof(sampler))[1:num_of_pars]
        println(io, par, " = $(getfield(sampler, par))  ", sampler.par_info[par])
    end
    println(io, "\nDefault output parameters:")
    for par in fieldnames(typeof(sampler))[(num_of_pars+1):(num_of_pars+3)]
        println(io, par, " = $(getfield(sampler, par))  ")
    end
end


text:
using SIMDPolynomials
using SIMDPolynomials: divexact, lc, var
using Test

SIMDPolynomials.debugmode() = true

@testset "Fuzz SparsePoly" begin
    x = Uninomial(0, 1)
    for i in 1:10000
        p = sum(rand(-2:2)*x^i for i in 0:10)
        q = sum(rand(-2:2)*x^i for i in 0:5)
        try
            g = gcd(p, q)
            pdg = divexact(p, g)
            qdg = divexact(q, g)
            @test gcd(pdg, qdg) == one(x)
            pq = p * q
            @test divexact(pq, p) == q
            @test divexact(pq, q) == p
        catch
            display(p)
            display(q)
            rethrow()
        end
    end
end

@testset "Fuzz MPoly" begin
    for monos in [[PackedMonomial{4,8}(i) for i in 0:2], [Monomial([i]) for i in 0:2]]
        x, y, z = monos
        for _ in 1:10000
            pterms = [rand(-2:2)*prod(rand([x, y, z])^i for i in 0:3) for j in 0:7]
            p = sum(pterms)
            q = sum(rand(-2:2)*prod(rand([x, y, z])^i for i in 0:2)  for j in 0:5)
            try
                g = gcd(p, q)
                pdg = divexact(p, g)
                qdg = divexact(q, g)
                @test gcd(pdg, qdg) == one(x)
            catch
                display(p)
                display(q)
                rethrow()
            end
        end
    end
end

@testset "pseudorem" begin
    x = Uninomial(0, 1)
    p = 2x^10 + x^7 + 7*x^2 + x + 3x
    q = (p+one(p)) * (p+2one(p)) * (p+3one(p))
    @test pseudorem(q, p) == 12582912*one(p)
    q = x^7 + 20*one(x)
    @test pseudorem(q, p) == q
    @test pseudorem(p, q) == -40*x^3 + 7*x^2 + 4*x - 20*one(x)
    q = x^6 + 23*one(x)
    @test pseudorem(p, q) == -46*x^4 + 7*x^2 - 19*x
end

function test_gcd(x, y)
    g1 = gcd(x, y)
    g2 = gcd(y, x)
    @test sign(lc(g1)) * g1 == sign(lc(g2)) * g2
    return g1
end

@testset "GCD" begin
    x, y, z = [Monomial([i]) for i in 0:2]
    c1 = 10*(x * z + x)
    c2 = 2*(x^2 + z)
    c3 = 2*(2 - z  )
    c4 = 20*(x * z^2)
    e1 = 0
    e2 = 5
    e3 = 7
    e4 = 10
    p = c1 * y^e1 + c2 * y^e2 + c3 * y^e3 + c4 * y^e4
    pp = SIMDPolynomials.SparsePoly(p, y.ids[1])
    @test var(pp) == y.ids[1]
    @test coeffs(pp) == [c4, c3, c2, c1]
    @test pp.exps == [e4, e3, e2, e1]
    q = prod(i->p + i, 0:3)
    @test length(terms(q)) == 262
    for i in 0:3
        @test test_gcd(p + i, q) == p + i
    end

    k = y^2 + 1
    @test test_gcd(x*k, z*k) == k
    @test test_gcd(x*k, (z+1)*k) == k
    @test test_gcd(x*k, p*k) == k

    p = 2*x*y
    q = 2*x*y + x
    @test test_gcd(q, p) == x

    p = x*y + y
    q = -x*y - y
    @test test_gcd(q, p) == p

    @test test_gcd((-x + 1) * (y^2+1), -x+1) == -x + 1
end

@testset "PackedMonomial" begin
    x, y, z, t = [PackedMonomial{4,8}(i) for i in 0:3]
    m = (x^10 * y^3 * z^2 * t^4)^6
    n = (x^9 * y^8 * z^3 * t^3)^6
    g = gcd(m, n)
    @test g == x^54 * y^18 *z^12 * t^18
    @test SIMDPolynomials.degree(g) == 54 + 18 + 12 + 18
    x, y, z, t = [PackedMonomial{4,7}(i) for i in 0:3]
    m = (x^10 * y^3 * z^2 * t^4)^6
    @test_throws Base.OverflowError (x^9 * y^8 * z^3 * t^3)^6
    g = gcd(m, m)
    @test g == m
    @test SIMDPolynomials.degree(g) == 60 + 18 + 12 + 24

    c1 = 10*(x * z + x)
    c2 = 2*(x^2 + z)
    c3 = 2*(2 - z  )
    c4 = 20*(x * z^2)
    e1 = 0
    e2 = 5
    e3 = 7
    e4 = 10
    p = c1 * y^e1 + c2 * y^e2 + c3 * y^e3 + c4 * y^e4
    q = prod(i->p + i, 0:3);
    @test length(terms(q)) == 262
    for i in 0:3
        @test test_gcd(p + i, q) == p + i
    end

    k = y^2 + 1
    @test test_gcd(x*k, z*k) == k
    @test test_gcd(x*k, (z+1)*k) == k
    @test test_gcd(x*k, p*k) == k

    p = 2*x*y
    q = 2*x*y + x
    @test test_gcd(q, p) == x

    p = x*y + y
    q = -x*y - y
    @test test_gcd(q, p) == p

    @test test_gcd((-x + 1) * (y^2+1), -x+1) == -x + 1

    p1 = 1 + x + y + z + t
    n = 15
    p = p1^n;
    @test test_gcd(p, p1) == p1
    k = p;
    for i in 1:n
        k = divexact(k, p1)
    end
    @test k == 1
    q = (p + 1) * p;
    @test gcd(q, p) == p
end


text:
using LinearAlgebra

xrange(::Type{T}, len=11) where {T} = range(-one(T), one(T), length = len)

##################################
# UNNORMALIZED LEGENDRE FUNCTIONS
##################################

# P_0^0 is constant. Verify the output is invariant for several inputs.
@testset "Constant P_0^0 ($T)" for T in NumTypes
    @test @inferred(legendre(LegendreUnitNorm(), 0, 0, T(0.1))) isa T
    @test all(legendre(LegendreUnitNorm(), 0, 0, x) == one(T) for x in xrange(T))
end

# Match P_ℓ^{m=0} terms for ℓ=1...9 where the analytical expressions
# have been taken from a table of equations, assuming x = cos(π/4).
@testset "Analytic checks for P_ℓ^0(cos(π/4)) ($T)" for T in NumTypes
    LMAX = 9
    plm = Matrix{T}(undef, LMAX+1, LMAX+1)
    legendre!(LegendreUnitNorm(), plm, LMAX, LMAX, sqrt(T(2))/2)
    @test plm[2, 1] ≈  sqrt(T(2))/2
    @test plm[3, 1] ≈  T(1)/4
    @test plm[4, 1] ≈ -sqrt(T(2))/8
    @test plm[5, 1] ≈ -T(13)/32
    @test plm[6, 1] ≈ -17sqrt(T(2))/64
    @test plm[7, 1] ≈ -T(19)/128
    @test plm[8, 1] ≈  23sqrt(T(2))/256
    @test plm[9, 1] ≈  T(611)/2048
    @test plm[10,1] ≈  827sqrt(T(2))/4096
end

#################################################
# ORTHONORMAL & SPHERICAL HARMONIC NORMALIZATION
#################################################

# P_0^0 is constant. Verify the output is invariant for several inputs.
@testset "Constant P_0^0 ($T)" for T in NumTypes
    @test @inferred(legendre(LegendreFourPiNorm(), 0, 0, T(0.1))) isa T
    @test @inferred(legendre(LegendreOrthoNorm(),  0, 0, T(0.1))) isa T
    @test @inferred(legendre(LegendreSphereNorm(), 0, 0, T(0.1))) isa T
    @test all(legendre(LegendreFourPiNorm(), 0, 0, x) == one(T)           for x in xrange(T))
    @test all(legendre(LegendreOrthoNorm(),  0, 0, x) == sqrt(T(0.5))     for x in xrange(T))
    @test all(legendre(LegendreSphereNorm(), 0, 0, x) == sqrt(inv(4T(π))) for x in xrange(T))
end

# The initialization condition
#
#   P_ℓ^ℓ(x) = (-1)^ℓ (2ℓ-1)!! (1-x^2)^(ℓ/2)
#
# can be used with extended-precision computations for special values
# of ℓ and m. To account for the spherical harmonic normalization,
# the formula can be simplified to
#
#   λ_ℓ^ℓ(x) = (-1)^ℓ 1/sqrt(4π) sqrt((2ℓ+1)!!/(2ℓ)!!) (1-x^2)^(ℓ/2)
#
# We'll compute the coefficient separately since it's "easy", and then
# we can use convenient angles that also expand to an exact analytic
# value for (1-x^2)^(ℓ/2).
#
# Assuming x = cos(θ), the (1-x^2)^(ℓ/2) term simplifies to an
# expression of \sin(θ):
#
#   (1-x^2)^(ℓ/2) = (sin^2 θ)^(ℓ/2) = (sin θ)^ℓ
#
# For θ = π/4:
#
#   * even ℓ => 1 / 2^(ℓ/2)
#   * odd ℓ => 1 / (2^(k/2) sqrt(2)) where k=floor(ℓ/2).
@testset "Analytic checks for λ_ℓ^m(cos(π/4))" begin
    function SphericalPll_coeff(T, ℓ)
        ratio(x) = sqrt(T(2x + 1) / T(2x))
        factorials = mapreduce(ratio, *, 1:ℓ, init=one(T))
        sign = (ℓ % 2 == 0) ? 1 : -1
        return sign * sqrt(inv(4T(π))) * factorials
    end

    LMAX = 99
    Λ = fill(0.0, LMAX+1, LMAX+1)
    legendre!(LegendreSphereNorm(), Λ, LMAX, LMAX, cosd(45.0))

    @test Λ[3,3]   ≈ SphericalPll_coeff(Float64,2)  / (2^1)
    @test Λ[4,4]   ≈ SphericalPll_coeff(Float64,3)  / (2^1 * sqrt(2))
    @test Λ[21,21] ≈ SphericalPll_coeff(Float64,20) / (2^10)
    @test Λ[22,22] ≈ SphericalPll_coeff(Float64,21) / (2^10 * sqrt(2))
    @test Λ[99,99] ≈ SphericalPll_coeff(Float64,98) / (2^49)
    @test Λ[100,100] ≈ SphericalPll_coeff(Float64,99) / (2^49 * sqrt(2))
end

# The normal P_ℓ^m should be equal to the spherical harmonic normalized
# λ_ℓ^m if we manually normalize them.
@testset "Equality of N_ℓ^m × P_ℓ^m and λ_ℓ^m ($T)" for T in NumTypes
    LMAX = 5
    atol = max(eps(T(100)), eps(100.0)) # maximum(plm_norm) is of order 10²
    plm_norm = zeros(T, LMAX+1, LMAX+1)
    plm_sphr = zeros(T, LMAX+1, LMAX+1)

    lmat = tril(repeat(collect(0:LMAX), 1, LMAX+1))
    mmat = tril(repeat(collect(0:LMAX)', LMAX+1, 1))
    nlm = tril(Nlm.(T, lmat, mmat))
    for x in xrange(T)
        legendre!(LegendreUnitNorm(),   plm_norm, LMAX, LMAX, x)
        legendre!(LegendreSphereNorm(), plm_sphr, LMAX, LMAX, x)
        @test all(isapprox.(nlm.*plm_norm, plm_sphr, atol=atol))
    end
end

# The orthonormal O_ℓ^m should be equal to the spherical harmonic normalized
# λ_ℓ^m if we manually normalize them.
@testset "Equality of (2π)^(-1/2) O_ℓ^m and λ_ℓ^m ($T)" for T in NumTypes
    LMAX = 5
    atol = max(eps(T(100)), eps(100.0)) # maximum(plm_norm) is of order 10²
    plm_orth = zeros(T, LMAX+1, LMAX+1)
    plm_sphr = zeros(T, LMAX+1, LMAX+1)
    for x in xrange(T)
        legendre!(LegendreOrthoNorm(),  plm_orth, LMAX, LMAX, x)
        legendre!(LegendreSphereNorm(), plm_sphr, LMAX, LMAX, x)
        @test all(isapprox.(plm_orth./sqrt(2T(π)), plm_sphr, atol=atol))
    end
end

# The spherical and 4π-normalized functions should be the same up to the
# differing factor of 1/√4π
@testset "Equality of (4π)^(-1/2) F_ℓ^m and λ_ℓ^m ($T)" for T in NumTypes
    LMAX = 5
    atol = max(eps(4T(π)), eps(4.0π))
    plm_4pi  = zeros(T, LMAX+1, LMAX+1)
    plm_sphr = zeros(T, LMAX+1, LMAX+1)
    for x in xrange(T)
        legendre!(LegendreFourPiNorm(), plm_4pi,  LMAX, LMAX, x)
        legendre!(LegendreSphereNorm(), plm_sphr, LMAX, LMAX, x)
        @test all(isapprox.(plm_4pi./sqrt(4T(π)), plm_sphr, atol=atol))
    end
end

################################
# Complex domain
################################

# Trivial case: real-only complex arguments should be identical to the real case
@testset "Complex arguments, real axis ($T)" for T in NumTypes
    LMAX = 10
    x = xrange(T, 100)
    z = complex(x)
    @test Plm.(0:LMAX, 0:LMAX, x) == real.(Plm.(0:LMAX, 0:LMAX, z))
    @test λlm.(0:LMAX, 0:LMAX, x) == real.(λlm.(0:LMAX, 0:LMAX, z))
end

################################
# Other analytic checks / issues
################################

@testset "Numerical stability/accuracy" begin
    # issue CMB.jl#11
    #   Overall accuracy --- fixed by adding FMA y=(1-x^2) and pairwise diagonal iteration
    x = sind(-57.5)
    lmax = 3 * 720
    Λ = λlm.(0:lmax, 0:lmax, x)
    # The amplitude bound of 1.5 is a very rough limit --- simply checking for
    # unbounded growth.
    @test all(abs.(Λ[end,:]) .< 1.5)

    # PR#16
    #   Increase accuracy of μ coefficient for spherical norm.
    #   Also test ν, and both already equivalent for unit norm.
    @testset "coeffs μ, ν ($norm)" for norm in (LegendreSphereNorm(), LegendreUnitNorm())
        lrng = 1:30_000
        μ(T, l) = AssociatedLegendrePolynomials.coeff_μ(norm, T, l)
        ν(T, l) = AssociatedLegendrePolynomials.coeff_ν(norm, T, l)
        @test μ.(Float64, lrng) == Float64.(μ.(BigFloat, lrng))
        @test ν.(Float64, lrng) == Float64.(ν.(BigFloat, lrng))
    end
end


text:
include("sdsos_example.jl")
include("options_pricing.jl")


text:
#!/usr/bin/env julia

using GR

x = [-3.3 + t*0.1 for t in 0:66]
y = [(t^5 - 13*t^3 + 36*t)::Float64 for t in x]
plot(x, y, bgcolor=163)


text:
println([1, 2, 3] .^ 3)
# `.` is a broadcast operator
# a .^ b equal to (^).(a, b)

max_stars_repo_name: dannypsnl/languages-learn

max_issues_repo_name: dannypsnl/languages-learn

max_forks_repo_name: dannypsnl/languages-learn


text:
#        __            ____
#       / /__ _  __   / __/                      __
#      / //_/(_)/ /_ / /  ___   ____ ___  __ __ / /_
#     / ,<  / // __/_\ \ / _ \ / __// _ \/ // // __/
#    /_/|_|/_/ \__//___// .__//_/   \___/\_,_/ \__/
#                      /_/   github.com/KitSprout
#
#   @file    kSerial.py
#   @author  KitSprout
#   @date    Mar-2020
#   @brief   kSerial packet format :
#            byte 1   : header 'K' (75)       [HK]
#            byte 2   : header 'S' (83)       [HS]
#            byte 3   : data bytes (12-bit)   [L ]
#            byte 4   : data type             [T ]
#            byte 5   : parameter 1           [P1]
#            byte 6   : parameter 2           [P2]
#            byte 7   : checksum              [CK]
#             ...
#            byte L-1 : data                  [DN]
#            byte L   : finish '\r' (13)      [ER]

using LibSerialPort
using Printf

typeconv = Dict(
    # string to value
    "UInt8"   => 0, "UInt16"  =>  1, "UInt32"  =>  2, "UInt64" =>  3,
    "Int8"    => 4, "Int16"   =>  5, "Int32"   =>  6, "Int64"  =>  7,
    "Float16" => 9, "Float32" => 10, "Float64" => 11,
    "R0"      => 8, "R1"      => 12, "R2"      => 13, "R3"     => 14, "R4" => 15,
    # value to string
    0 => "UInt8",    1 => "UInt16",   2 => "UInt32",   3 => "UInt64",
    4 => "Int8",     5 => "Int16",    6 => "Int32",    7 => "Int64",
    9 => "Float16", 10 => "Float32", 11 => "Float64",
    8 => "R0",      12 => "R1",      13 => "R2",      14 => "R3",     15 => "R4"
);

# packet = pack();                   ->  R0 (8-bit)                    -> [75, 83, 0, 8, 0, 0, 8, 13]
# packet = pack(param);              ->  only send parameter (8-bit)   -> [75, 83, 0, 8, P1, P2, CK, 13]
# packet = pack(param, type);        ->  user define data type (8-bit) -> [75, 83, L, T, P1, P2, CK, 13]
# packet = pack(param, data);        ->  automatic detect data type    -> [75, 83, L, T, P1, P2, Dn, CK, 13]
# packet = pack(param, data, type);  ->  user define data type         -> [75, 83, L, T, P1, P2, Dn, CK, 13]
function pack(param::Array{UInt8,1}=UInt8[0,0], type::Union{Integer,String}="R0", data::Vector=[])
    # @printf("[DBG] param = %s, type = %s, data = %s\n", typeof(param), typeof(type), typeof(data));
    # @printf("[DBG] param = %d, type = %d, data = %d, size = %d\n", length(param), length(type), length(data), sizeof(data));

    packet_bytes = sizeof(data);
    packet_type  = (packet_bytes >> 8) << 4;
    try
        # check type input
        typeconv[type];
        if typeof(type) == String
            packet_type += typeconv[type];
        else
            packet_type += type;
        end
    catch e
        @printf("can't find dict key-value: %s\n", e);
    end

    if length(data) == 0
        pdata = UInt8[];
    else
        try
            pdata = reinterpret(UInt8, data);
        catch e
            @printf("reinterpret error: %s\n", e);
        end
    end

    packet = UInt8['K', 'S'];               # header 'KS'
    append!(packet, packet_bytes % 256);    # data bytes
    append!(packet, packet_type);           # data type
    append!(packet, param);                 # parameter
    checksum = (packet[3] + packet[4] + packet[5] + packet[6]) % 256;
    append!(packet, checksum);              # checksum
    append!(packet, pdata);                 # data
    append!(packet, '\r');                  # finish '\r'

    return packet;
end

# info, data, last_index = unpack(packet);
# info[:, 1] : data length (bytes)
# info[:, 2] : data type
# info[:, 3] : parameter 1
# info[:, 4] : parameter 2
# info[:, 5] : checksum
function unpack(packet::Array{UInt8,1})
    # @printf("[DBG] packet = %s, %d\n", typeof(packet), length(packet));
    # println(packet);

    info = [];
    data = [];
    last_index = 0;

    packet_bytes = length(packet);

    # check 'KS'
    Ki = findall(x->x == UInt8('K'), packet);
    Ki = Ki[findall(x->x <= packet_bytes, Ki .+ 7)];
    Ki = Ki[findall(x->x == UInt8('S'), packet[Ki .+ 1])];
    if length(Ki) != 0
        # get packet info and checksum
        Ld  = packet[Ki .+ 2];
        Td  = packet[Ki .+ 3];
        P1d = packet[Ki .+ 4];
        P2d = packet[Ki .+ 5];
        CKd = packet[Ki .+ 6];
        # check checksum
        checksum = Ld + Td + P1d + P2d;
        Ki = Ki[findall(x->x == true, CKd .== checksum)];
        Ld = packet[Ki .+ 2];
        Td = packet[Ki .+ 3];
        # get packet data lens
        packet_lens = Ld + (UInt16.(Td) .>> 4) .<< 8;
        # check last packet
        idx = findall(x->x <= packet_bytes, Ki .+ packet_lens .+ 7);
        Er = Ki .+ packet_lens .+ 7;
        # check '\r'
        idx = findall(x->x == UInt8('\r'), packet[Er[idx]]);
        Ki = Ki[idx];
        if length(Ki) != 0
            # get info
            packet_lens  = packet_lens[idx];
            packet_type  = packet[Ki .+ 3] .& 0x0F;
            packet_param = [packet[Ki .+ 4], packet[Ki .+ 5]];
            packet_param = [packet[Ki .+ 4], packet[Ki .+ 5]];
            # info = Int.(hcat([packet_lens, packet_type, packet[Ki .+ 4], packet[Ki .+ 5], packet[Ki .+ 6]]...)'); # 5xN
            info = Int.(hcat([packet_lens, packet_type, packet[Ki .+ 4], packet[Ki .+ 5], packet[Ki .+ 6]]...));    # Nx5
            # get data
            data = [];
            for i = 1 : length(Ki)
                append!(data, [packet[Ki[i] + 6 .+ (1:packet_lens[i])]]);
            end
        end
    end

    return info, data, last_index;
end

function send(s, param, type, data)
    pkdata = pack(param, type, data);
    write(s, pkdata);
    return pkdata, length(pkdata);
end

function recv(s)
    nb = bytesavailable(s);
    packet = readbytes!(s, nb);
    info, data, last_index = unpack(packet);
    return info, data, nb;
end

# pk(UInt8[0x01, 0x00], "R1", [])
function pk(param, type, data)
    println(param, ", ", type, ", ", data);

    packet = pack(param, type, data);
    info, data, last_index = unpack(packet);
    println(UInt8.(info[:][3:4]), ", ", typeconv[info[:][2]], ", ", data[1][:]);
end

function check(delay)
    send(sp, UInt8[0x00, 0x00], "R2", []);
    sleep(delay);
    info, data, nb = recv(sp);
    return data;
end

function scandevice(delay)
    send(sp, UInt8[0xAB, 0x00], "R2", []);
    sleep(delay);
    info, data, nb = recv(sp);
    return data;
end

function scanreg(address::UInt8, delay)
    send(sp, UInt8[0xCB, (address << 1) + 1], "R2", []);
    sleep(delay);
    info, data, nb = recv(sp);
    return data;
end

sp = open("COM3", 115200);
println(check(0.01));
println(scandevice(0.01));
println(scanreg(0x0c, 0.1));
close(sp);


text:
# common car camera tools

using LinearAlgebra
using DataStructures
using RoME, DistributedFactorGraphs

import Rotations as _Rot

## clear old memory

# delete!(vis[:tags])
# delete!(vis[:poses])

##

abstract type SynchronizingBuffer end

struct SynchronizeCarMono <: SynchronizingBuffer
  syncList::Vector{Symbol}
  # sequence #, time nanoseconds, data
  leftFwdCam::CircularBuffer{Tuple{Int, Int64, Any}}
  rightFwdCam::CircularBuffer{Tuple{Int, Int64, Any}}
  camOdo::CircularBuffer{Tuple{Int, Int64, Any}}
  cmdVal::CircularBuffer{Tuple{Int, Int64, Any}}
end

SynchronizeCarMono( len::Int=30;
                    syncList::Vector=Symbol[],
                    leftFwdCam=CircularBuffer{Tuple{Int, Int, Any}}(len),
                    rightFwdCam=CircularBuffer{Tuple{Int, Int, Any}}(len),
                    camOdo=CircularBuffer{Tuple{Int, Int, Any}}(len),
                    cmdVal=CircularBuffer{Tuple{Int, Int, Any}}(len) ) = SynchronizeCarMono(syncList,leftFwdCam,rightFwdCam,camOdo,cmdVal)
#

struct RacecarTools
  tagDetector
end

struct RacecarVisualizations
  gui
  vis
end

struct FrontEndContainer{M,T,A}
  slam::SLAMWrapperLocal
  middleware::M
  synchronizer::T
  tools::A
  datastore
end

FrontEndContainer(s::SLAMWrapperLocal,m::M,t::T,a::A,d) where {M, T, A} = FrontEndContainer{M,T,A}(s,m,t,a,d)

## data management functions

# function getSyncLatestPair(syncHdlrs::Dict; weirdOffset::Dict=Dict())::Tuple
function findSyncLatestIdx( syncz::SynchronizingBuffer;
                            weirdOffset::Dict=Dict(),
                            syncList::Vector{Symbol}=syncz.syncList)
  # get all sequence numbers
  len = length(syncList)
  if len == 1
    # simply return the newest element in the buffer
    sy = getfield(syncz, syncList[1])
    return length(sy)
  elseif len == 0
    @warn "findSyncLatestIdx has nothing to synchronize -- on account of .syncList being empty."
    return 0  # error
  end
  LEN = Vector{Int}(undef, len)
  for idx in 1:len
    LEN[idx] = getfield(syncz, syncList[idx]) |> length
    if LEN[idx] == 0
      return (NTuple{len, Int}(zeros(len)))
    end
  end

  SEQ = Vector{Vector{Int}}(undef, len)
  for idx in 1:len
    SEQ[idx] = (x->getfield(syncz, syncList[idx])[x][1]).(1:LEN[idx])
    if haskey(weirdOffset, syncList[idx])
      SEQ[idx] .+= weirdOffset[syncList[idx]]
    end
  end

  # find greatest common number
  seqLR = intersect(SEQ...)
  if 0 == length(seqLR)
    return NTuple{len,Int}(zeros(len))
  end
  seq = maximum(seqLR)
  # IDX = Vector{Int}(undef, len)

  return NTuple{len, Int}( (y->findfirst(x->x==seq, y)).(SEQ) )
end

function jsonResultsSLAM2D(dfg::AbstractDFG)
  allDicts = []
  allvars = ls(dfg, r"x\d") |> sortDFG

  for vidx in 1:(length(allvars)-1)
    ps = allvars[vidx]
    ns = allvars[vidx+1]
    if 0 == getPPEDict(getVariable(dfg, ns)) |> length || 0 == getPPEDict(getVariable(dfg, ps)) |> length
      continue
    end
    @show "fetchDataElement for $ps"
    if !hasDataEntry(getVariable(dfg,ps), :JOYSTICK_CMD_VALS)
      @error "missing data entry :JOYSTICK_CMD_VALS in $ps of length(allvars)=$(length(allvars))"
      continue
    end
    _,cmdData = getData(dfg, ps, :JOYSTICK_CMD_VALS)
    # cmdData = fetchDataElement(dfg, ps, :JOYSTICK_CMD_VALS)

    # axis: 2:throttle, 4:steering
    cd = Dict{Symbol,Any}(
      :posei => ps,
      :posej => ns,
      :axis_time => (x->x[2]).(cmdData),
      :axis_throttle => (x->x[3].axis[2]).(cmdData),
      :axis_steer => (x->x[3].axis[4]).(cmdData),
      :world_posei => getPPE(dfg,ps).suggested,
      :world_posej => getPPE(dfg,ns).suggested,
      :ti => getTimestamp(getVariable(dfg, ps)),
      :tj => getTimestamp(getVariable(dfg, ns)),
    )
    DT = (cd[:tj]-cd[:ti]).value*1e-3
    cd[:deltaT] = DT
    deltax = (cd[:world_posej][1] - cd[:world_posei][1])
    deltay = (cd[:world_posej][2] - cd[:world_posei][2])
    velx = deltax./DT
    vely = deltay./DT
    wVelx = isnan(velx) ? 0.0 : velx
    wVely = isnan(vely) ? 0.0 : vely
    # rotate velocity from world to body frame
    cd[:bodyi_delta_pose] = se2vee(SE2(cd[:world_posei])\SE2(cd[:world_posej]))
    biRw = TU.R(-cd[:world_posei][3])
    bVel = biRw*[wVelx;wVely]
    velx = bVel[1]
    vely = bVel[2]
    if vidx in [1;2]
      velx = 0.0
      vely = 0.0
    end
    cd[:body_veli] = [velx;vely]
    push!(allDicts, cd)
  end
  return allDicts
end

## Gadfly to video

# reproject a bearing range onto (assumed level) image.

## sequence of localProduct plots to monitor progression of a solve

# kys = filter(x->!(x in [:graphinit;:default]), listSolveKeys(fg) |> collect |>sortDFG)

# PL = kys .|> x->plotSLAM2D(fg, solveKey=x, xmin=-8,xmax=6,ymin=-4,ymax=8, contour=false);

# PL = kys .|> x->plotLocalProduct(fg, :x7, solveKey=x, dims=[1;2]);
# coo = Coord.Cartesian(xmin=-8,xmax=6,ymin=-4,ymax=6)

# PL .|> x->(x.coord=coo);

# PLi = convert.(Matrix{RGB}, PL);

# writevideo("/tmp/slam2d__.avi", PLi, fps=5)

#


text:
"""
  Profiling for all matpower instances
    To be adapted...
"""
include(joinpath(ROOT,"..","src_PowSysMod", "PowSysMod_body.jl"))

struct Measurement
  nb_buses
  time
  bytes
  gctime
  nbpoolallocs
  nbmalloc
  summarysize
end

## Input
files_m = [split(name, ".")[1] for name in readdir(joinpath("instances", "matpower"))]
files_dat = [split(name, ".")[1] for name in readdir(joinpath("instances", "matpower_QCQP"))]

instances = intersect(files_m, files_dat)

nb_repeat = 3

instancename = sort(collect(instances))[1]

data_OPFpbs = Dict{String, Measurement}()
data_Problems = Dict{String, Measurement}()
data_errors = Dict{String, Float64}()

instnb = 1
for instancename in sort(collect(instances))
  println("\n-- Working on $instancename ($instnb/$(length(instances)))")
  nb_buses = parse(matchall(r"\d+", instancename)[1])

  ## Building OPF_problems
  t = bytes = gctime = nbpoolallocs = nbmalloc = sumsize = 0
  OPFpbs = build_OPFproblems(MatpowerSimpleInput, [joinpath("instances", "matpower", "$instancename.m")])
  for i=1:nb_repeat
    print("$i ")
    OPFpbs, t_, bytes_, gctime_, gc_ = @timed build_OPFproblems(MatpowerSimpleInput, [joinpath("instances", "matpower", "$instancename.m")])
    t += t_/nb_repeat
    bytes += bytes_/nb_repeat
    gctime += gctime_/nb_repeat
    nbpoolallocs += gc_.poolalloc/nb_repeat
    nbmalloc += gc_.malloc/nb_repeat
    sumsize += Base.summarysize(OPFpbs)/nb_repeat
  end
  data_OPFpbs[instancename] = Measurement(nb_buses, t, bytes, gctime, nbpoolallocs, nbmalloc, sumsize)

  for (busname, busdata) in OPFpbs["BaseCase"].ds.bus
    if haskey(busdata, "Gen_reference")
      OPFpbs["BaseCase"].mc.node_formulations[busname]["Gen_reference"] = :None
    end
  end

  ## Building OPF problem
  t = bytes = gctime = nbpoolallocs = nbmalloc = sumsize = 0
  pb = build_Problem!(OPFpbs, "BaseCase")
  for i=1:nb_repeat
    print("$i ")
    pb, t_, bytes_, gctime_, gc_ = @timed build_Problem!(OPFpbs, "BaseCase")
    t += t_/nb_repeat
    bytes += bytes_/nb_repeat
    gctime += gctime_/nb_repeat
    nbpoolallocs += gc_.poolalloc/nb_repeat
    nbmalloc += gc_.malloc/nb_repeat
    sumsize += Base.summarysize(pb)/nb_repeat
  end
  data_Problems[instancename] = Measurement(nb_buses, t, bytes, gctime, nbpoolallocs, nbmalloc, sumsize)

  export_to_dat(pb, "my$instancename.dat")
  filename_dat = joinpath("instances", "matpower_QCQP", "$instancename.dat")
  error, _ = compare_dat("my$instancename.dat", filename_dat)
  rm("my$instancename.dat")
  data_errors[instancename] = error

  instnb+=1
end

# output
filename = "buildOPFpbs_$(Dates.format(now(), "yy_u_dd_HH_MM")).csv"
touch(filename)

open(filename, "w") do outfile
  print(outfile, "instancename;nb_buses;time;bytes;gctime;nbpoolallocs;nbmalloc;summarysize;dat_error\n")

  insancenames_ord = sort(collect(keys(data_OPFpbs)), by=x->data_OPFpbs[x].nb_buses)

  for instancename in insancenames_ord
    measure = data_OPFpbs[instancename]
    print(outfile, "$instancename;$(measure.nb_buses);$(measure.time);$(measure.bytes);$(measure.gctime);$(measure.nbpoolallocs);$(measure.nbmalloc);$(measure.summarysize);$(data_errors[instancename])\n")
  end
end

filename = "build_Problem_$(Dates.format(now(), "yy_u_dd_HH_MM")).csv"
touch(filename)

open(filename, "w") do outfile
  print(outfile, "instancename;nb_buses;time;bytes;gctime;nbpoolallocs;nbmalloc;summarysize;dat_error\n")

  insancenames_ord = sort(collect(keys(data_Problems)), by=x->data_Problems[x].nb_buses)

  for instancename in insancenames_ord
    measure = data_Problems[instancename]
    print(outfile, "$instancename;$(measure.nb_buses);$(measure.time);$(measure.bytes);$(measure.gctime);$(measure.nbpoolallocs);$(measure.nbmalloc);$(measure.summarysize);$(data_errors[instancename])\n")
  end
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
# load required modules
using DataFrames
using Gadfly
using Random
using RDatasets

# set seed for reproducible sequence of numbers
Random.seed!(42);

# load iris dataset
iris = dataset("datasets", "iris")

# scratter plot of sepal width vs sepal length
# with species as color, label and legend
plot(iris, x=:SepalLength, y=:SepalWidth, color=:Species, Geom.point,
    Guide.xlabel("Sepal Length"), Guide.ylabel("Sepal Width"),
    Guide.title("Sepal Width vs Sepal Length"))

# stack horizontally two scatter plots
figa = plot(iris, x=:SepalLength, y=:SepalWidth, color=:Species)
figb = plot(iris, x=:PetalLength, y=:PetalWidth, color=:Species)
hstack(figa, figb)

# boxplot of sepal length
plot(iris, y=:SepalLength, Geom.boxplot, Guide.ylabel("Sepal Length"))

# boxplot of sepal width vs sepal length
plot(iris, x=:SepalLength, y=:SepalWidth, Geom.boxplot)

# densityplot of sepal length to show
# density of the sepal length over the particular range of each species
plot(iris, x=:SepalLength, color=:Species, Geom.density, Guide.xlabel("Sepal Length"))

# densityplot of sepal width to show
# density of the sepal width over the particular range of each species
plot(iris, x=:SepalWidth, color=:Species, Geom.density, Guide.xlabel("Sepal Width"))

# histograms of sepal length distribution of each species
plot(iris, x=:SepalLength, color=:Species, Geom.histogram, Guide.xlabel("Sepal Length"))

# histograms of sepal width distribution of each species
plot(iris, x=:SepalWidth, color=:Species, Geom.histogram, Guide.xlabel("Sepal Width"))

# visualize sepal length density of each species
plot(
    iris,
    x=:Species, y=:SepalLength, color=:Species,
    Geom.violin,
    Guide.ylabel("Sepal Length"), Guide.xlabel("Species"),
    Guide.title("Species Sepal Length Density")
)

# visualize sepal width density of each species
plot(
    iris,
    x=:Species, y=:SepalWidth, color=:Species,
    Geom.violin,
    Guide.ylabel("Sepal Width"), Guide.xlabel("Species"),
    Guide.title("Species Sepal Width Density")
)

# visualize sepal length density of each species
plot(
    iris,
    x=:Species, y=:SepalLength, color=:Species,
    Geom.beeswarm,
    Guide.ylabel("Sepal Length"), Guide.xlabel("Species"),
    Guide.title("Species Sepal Length Density")
)

# visualize sepal width density of each species
plot(
    iris,
    x=:Species, y=:SepalWidth, color=:Species,
    Geom.beeswarm,
    Guide.ylabel("Sepal Width"), Guide.xlabel("Species"),
    Guide.title("Species Sepal Width Density")
)

max_stars_repo_name: lykmapipo/data-science-learning

max_issues_repo_name: lykmapipo/data-science-learning

max_forks_repo_name: lykmapipo/data-science-learning