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

@test integrate(sech(a+b*x), x) == :(b^-1*arctan(sinh(a+b*x)))
@test integrate(sech(a+b*x)^2, x) == :(b^-1*tanh(a+b*x))
@test integrate(sech(a+b*x)^3, x) == :((1/2)*b^-1*arctan(sinh(a+b*x))+(1/2)*b^-1*sech(a+b*x)*tanh(a+b*x))
@test integrate(sech(a+b*x)^4, x) == :(b^-1*tanh(a+b*x)+-1//3*b^-1*tanh(a+b*x)^3)
@test integrate(sech(a+b*x)^5, x) == :(3//8*b^-1*arctan(sinh(a+b*x))+1//4*b^-1*sech(a+b*x)^3*tanh(a+b*x)+3//8*b^-1*sech(a+b*x)*tanh(a+b*x))
@test integrate(sech(a+b*x)^6, x) == :(b^-1*tanh(a+b*x)+-2//3*b^-1*tanh(a+b*x)^3+1//5*b^-1*tanh(a+b*x)^5)
@test integrate(sech(7x)^4, x) == :(-1//21*tanh(7x)^3+1//7*tanh(7x))
@test integrate(sech(pi*x)^6, x) == :(pi^-1*tanh(pi*x)+-2//3*pi^-1*tanh(pi*x)^3+1//5*pi^-1*tanh(pi*x)^5)
@test integrate(sech(a+b*x)^5//2, x) == :(2//3*b^-1*sech(a+b*x)^3//2*sinh(a+b*x)+-2//3*im*b^-1*cosh(a+b*x)^(1/2)*sech(a+b*x)^(1/2)*Elliptic.F((1/2)*im*(a+b*x),2))
@test integrate(sech(a+b*x)^3//2, x) == :(2*b^-1*sech(a+b*x)^(1/2)*sinh(a+b*x)+2*im*b^-1*cosh(a+b*x)^(1/2)*sech(a+b*x)^(1/2)*Elliptic.E((1/2)*im*(a+b*x),2))
@test integrate(sech(a+b*x)^(1/2), x) == :(-2*im*b^-1*cosh(a+b*x)^(1/2)*sech(a+b*x)^(1/2)*Elliptic.F((1/2)*im*(a+b*x),2))
@test integrate(sech(a+b*x)^-1//2, x) == :(-2*im*b^-1*cosh(a+b*x)^(1/2)*sech(a+b*x)^(1/2)*Elliptic.E((1/2)*im*(a+b*x),2))
@test integrate(sech(a+b*x)^-3//2, x) == :(2//3*b^-1*sech(a+b*x)^-1//2*sinh(a+b*x)+-2//3*im*b^-1*cosh(a+b*x)^(1/2)*sech(a+b*x)^(1/2)*Elliptic.F((1/2)*im*(a+b*x),2))
@test integrate(sech(a+b*x)^-5//2, x) == :(2//5*b^-1*sech(a+b*x)^-3//2*sinh(a+b*x)+-6//5*im*b^-1*cosh(a+b*x)^(1/2)*sech(a+b*x)^(1/2)*Elliptic.E((1/2)*im*(a+b*x),2))
@test integrate((b*sech(c+d*x))^7//2, x) == :(2//5*b*d^-1*(b*sech(c+d*x))^5//2*sinh(c+d*x)+6//5*b^3*d^-1*(b*sech(c+d*x))^(1/2)*sinh(c+d*x)+6//5*im*b^4*d^-1*(b*sech(c+d*x))^-1//2*cosh(c+d*x)^-1//2*Elliptic.E((1/2)*im*(c+d*x),2))
@test integrate((b*sech(c+d*x))^5//2, x) == :(2//3*b*d^-1*(b*sech(c+d*x))^3//2*sinh(c+d*x)+-2//3*im*b^2*d^-1*(b*sech(c+d*x))^(1/2)*cosh(c+d*x)^(1/2)*Elliptic.F((1/2)*im*(c+d*x),2))
@test integrate((b*sech(c+d*x))^3//2, x) == :(2*b*d^-1*(b*sech(c+d*x))^(1/2)*sinh(c+d*x)+2*im*b^2*d^-1*(b*sech(c+d*x))^-1//2*cosh(c+d*x)^-1//2*Elliptic.E((1/2)*im*(c+d*x),2))
@test integrate((b*sech(c+d*x))^(1/2), x) == :(-2*im*d^-1*(b*sech(c+d*x))^(1/2)*cosh(c+d*x)^(1/2)*Elliptic.F((1/2)*im*(c+d*x),2))
@test integrate((b*sech(c+d*x))^-1//2, x) == :(-2*im*d^-1*(b*sech(c+d*x))^-1//2*cosh(c+d*x)^-1//2*Elliptic.E((1/2)*im*(c+d*x),2))
@test integrate((b*sech(c+d*x))^-3//2, x) == :(2//3*b^-1*d^-1*(b*sech(c+d*x))^-1//2*sinh(c+d*x)+-2//3*im*b^-2*d^-1*(b*sech(c+d*x))^(1/2)*cosh(c+d*x)^(1/2)*Elliptic.F((1/2)*im*(c+d*x),2))
@test integrate((b*sech(c+d*x))^-5//2, x) == :(2//5*b^-1*d^-1*(b*sech(c+d*x))^-3//2*sinh(c+d*x)+-6//5*im*b^-2*d^-1*(b*sech(c+d*x))^-1//2*cosh(c+d*x)^-1//2*Elliptic.E((1/2)*im*(c+d*x),2))
@test integrate((b*sech(c+d*x))^-7//2, x) == :(2//7*b^-1*d^-1*(b*sech(c+d*x))^-5//2*sinh(c+d*x)+10//21*b^-3*d^-1*(b*sech(c+d*x))^-1//2*sinh(c+d*x)+-10//21*im*b^-4*d^-1*(b*sech(c+d*x))^(1/2)*cosh(c+d*x)^(1/2)*Elliptic.F((1/2)*im*(c+d*x),2))
@test integrate((sech(a+b*x)^2)^7//2, x) == :(5//16*b^-1*arcsin(tanh(a+b*x))+1//6*b^-1*(sech(a+b*x)^2)^5//2*tanh(a+b*x)+5//16*b^-1*(sech(a+b*x)^2)^(1/2)*tanh(a+b*x)+5//24*b^-1*(sech(a+b*x)^2)^3//2*tanh(a+b*x))
@test integrate((sech(a+b*x)^2)^5//2, x) == :(3//8*b^-1*arcsin(tanh(a+b*x))+1//4*b^-1*(sech(a+b*x)^2)^3//2*tanh(a+b*x)+3//8*b^-1*(sech(a+b*x)^2)^(1/2)*tanh(a+b*x))
@test integrate((sech(a+b*x)^2)^3//2, x) == :((1/2)*b^-1*arcsin(tanh(a+b*x))+(1/2)*b^-1*(sech(a+b*x)^2)^(1/2)*tanh(a+b*x))
@test integrate((sech(a+b*x)^2)^(1/2), x) == :(b^-1*arcsin(tanh(a+b*x)))
@test integrate((sech(a+b*x)^2)^-1//2, x) == :(b^-1*(sech(a+b*x)^2)^-1//2*tanh(a+b*x))
@test integrate((sech(a+b*x)^2)^-3//2, x) == :(1//3*b^-1*(sech(a+b*x)^2)^-3//2*tanh(a+b*x)+2//3*b^-1*(sech(a+b*x)^2)^-1//2*tanh(a+b*x))
@test integrate((sech(a+b*x)^2)^-5//2, x) == :(1//5*b^-1*(sech(a+b*x)^2)^-5//2*tanh(a+b*x)+4//15*b^-1*(sech(a+b*x)^2)^-3//2*tanh(a+b*x)+8//15*b^-1*(sech(a+b*x)^2)^-1//2*tanh(a+b*x))
@test integrate((sech(a+b*x)^2)^-7//2, x) == :(1//7*b^-1*(sech(a+b*x)^2)^-7//2*tanh(a+b*x)+6//35*b^-1*(sech(a+b*x)^2)^-5//2*tanh(a+b*x)+8//35*b^-1*(sech(a+b*x)^2)^-3//2*tanh(a+b*x)+16//35*b^-1*(sech(a+b*x)^2)^-1//2*tanh(a+b*x))
@test integrate((a*sech(x)^2)^5//2, x) == :(3//8*a^5//2*arctan(a^(1/2)*(a*sech(x)^2)^-1//2*tanh(x))+1//4*a*(a*sech(x)^2)^3//2*tanh(x)+3//8*a^2*(a*sech(x)^2)^(1/2)*tanh(x))
@test integrate((a*sech(x)^2)^3//2, x) == :((1/2)*a^3//2*arctan(a^(1/2)*(a*sech(x)^2)^-1//2*tanh(x))+(1/2)*a*(a*sech(x)^2)^(1/2)*tanh(x))
@test integrate((a*sech(x)^2)^(1/2), x) == :(a^(1/2)*arctan(a^(1/2)*(a*sech(x)^2)^-1//2*tanh(x)))
@test integrate((a*sech(x)^2)^-1//2, x) == :((a*sech(x)^2)^-1//2*tanh(x))
@test integrate((a*sech(x)^2)^-3//2, x) == :(1//3*(a*sech(x)^2)^-3//2*tanh(x)+2//3*a^-1*(a*sech(x)^2)^-1//2*tanh(x))
@test integrate((a*sech(x)^2)^-5//2, x) == :(1//5*(a*sech(x)^2)^-5//2*tanh(x)+4//15*a^-1*(a*sech(x)^2)^-3//2*tanh(x)+8//15*a^-2*(a*sech(x)^2)^-1//2*tanh(x))
@test integrate((a*sech(x)^2)^-7//2, x) == :(1//7*(a*sech(x)^2)^-7//2*tanh(x)+6//35*a^-1*(a*sech(x)^2)^-5//2*tanh(x)+8//35*a^-2*(a*sech(x)^2)^-3//2*tanh(x)+16//35*a^-3*(a*sech(x)^2)^-1//2*tanh(x))
@test integrate((a*sech(x)^3)^5//2, x) == :(154//585*a^2*(a*sech(x)^3)^(1/2)*tanh(x)+2//13*a^2*(a*sech(x)^3)^(1/2)*sech(x)^4*tanh(x)+22//117*a^2*(a*sech(x)^3)^(1/2)*sech(x)^2*tanh(x)+154//195*a^2*(a*sech(x)^3)^(1/2)*cosh(x)*sinh(x)+154//195*im*a^2*(a*sech(x)^3)^(1/2)*cosh(x)^3//2*Elliptic.E((1/2)*im*x,2))
@test integrate((a*sech(x)^3)^3//2, x) == :(10//21*a*(a*sech(x)^3)^(1/2)*sinh(x)+2//7*a*(a*sech(x)^3)^(1/2)*sech(x)*tanh(x)+-10//21*im*a*(a*sech(x)^3)^(1/2)*cosh(x)^3//2*Elliptic.F((1/2)*im*x,2))
@test integrate((a*sech(x)^3)^(1/2), x) == :(2*(a*sech(x)^3)^(1/2)*cosh(x)*sinh(x)+2*im*(a*sech(x)^3)^(1/2)*cosh(x)^3//2*Elliptic.E((1/2)*im*x,2))
@test integrate((a*sech(x)^3)^-1//2, x) == :(2//3*(a*sech(x)^3)^-1//2*tanh(x)+-2//3*im*(a*sech(x)^3)^-1//2*cosh(x)^-3//2*Elliptic.F((1/2)*im*x,2))
@test integrate((a*sech(x)^3)^-3//2, x) == :(14//45*a^-1*(a*sech(x)^3)^-1//2*sinh(x)+2//9*a^-1*(a*sech(x)^3)^-1//2*cosh(x)^2*sinh(x)+-14//15*im*a^-1*(a*sech(x)^3)^-1//2*cosh(x)^-3//2*Elliptic.E((1/2)*im*x,2))
@test integrate((a*sech(x)^3)^-5//2, x) == :(26//77*a^-2*(a*sech(x)^3)^-1//2*tanh(x)+2//15*a^-2*(a*sech(x)^3)^-1//2*cosh(x)^5*sinh(x)+26//165*a^-2*(a*sech(x)^3)^-1//2*cosh(x)^3*sinh(x)+78//385*a^-2*(a*sech(x)^3)^-1//2*cosh(x)*sinh(x)+-26//77*im*a^-2*(a*sech(x)^3)^-1//2*cosh(x)^-3//2*Elliptic.F((1/2)*im*x,2))
@test integrate((a*sech(x)^4)^7//2, x) == :(a^3*(a*sech(x)^4)^(1/2)*cosh(x)*sinh(x)+-2*a^3*(a*sech(x)^4)^(1/2)*sinh(x)^2*tanh(x)+3*a^3*(a*sech(x)^4)^(1/2)*sinh(x)^2*tanh(x)^3+-20//7*a^3*(a*sech(x)^4)^(1/2)*sinh(x)^2*tanh(x)^5+-6//11*a^3*(a*sech(x)^4)^(1/2)*sinh(x)^2*tanh(x)^9+1//13*a^3*(a*sech(x)^4)^(1/2)*sinh(x)^2*tanh(x)^11+5//3*a^3*(a*sech(x)^4)^(1/2)*sinh(x)^2*tanh(x)^7)
@test integrate((a*sech(x)^4)^5//2, x) == :(a^2*(a*sech(x)^4)^(1/2)*cosh(x)*sinh(x)+-4//3*a^2*(a*sech(x)^4)^(1/2)*sinh(x)^2*tanh(x)+-4//7*a^2*(a*sech(x)^4)^(1/2)*sinh(x)^2*tanh(x)^5+1//9*a^2*(a*sech(x)^4)^(1/2)*sinh(x)^2*tanh(x)^7+6//5*a^2*(a*sech(x)^4)^(1/2)*sinh(x)^2*tanh(x)^3)
@test integrate((a*sech(x)^4)^3//2, x) == :(a*(a*sech(x)^4)^(1/2)*cosh(x)*sinh(x)+-2//3*a*(a*sech(x)^4)^(1/2)*sinh(x)^2*tanh(x)+1//5*a*(a*sech(x)^4)^(1/2)*sinh(x)^2*tanh(x)^3)
@test integrate((a*sech(x)^4)^(1/2), x) == :((a*sech(x)^4)^(1/2)*cosh(x)*sinh(x))
@test integrate((a*sech(x)^4)^-1//2, x) == :((1/2)*(a*sech(x)^4)^-1//2*tanh(x)+(1/2)*x*(a*sech(x)^4)^-1//2*sech(x)^2)
@test integrate((a*sech(x)^4)^-3//2, x) == :(5//16*a^-1*(a*sech(x)^4)^-1//2*tanh(x)+1//6*a^-1*(a*sech(x)^4)^-1//2*cosh(x)^3*sinh(x)+5//16*x*a^-1*(a*sech(x)^4)^-1//2*sech(x)^2+5//24*a^-1*(a*sech(x)^4)^-1//2*cosh(x)*sinh(x))
@test integrate((a*sech(x)^4)^-5//2, x) == :(63//256*a^-2*(a*sech(x)^4)^-1//2*tanh(x)+1//10*a^-2*(a*sech(x)^4)^-1//2*cosh(x)^7*sinh(x)+9//80*a^-2*(a*sech(x)^4)^-1//2*cosh(x)^5*sinh(x)+21//128*a^-2*(a*sech(x)^4)^-1//2*cosh(x)*sinh(x)+21//160*a^-2*(a*sech(x)^4)^-1//2*cosh(x)^3*sinh(x)+63//256*x*a^-2*(a*sech(x)^4)^-1//2*sech(x)^2)
@test integrate((a+a*sech(x))^-1*sinh(x)^4, x) == :(-1//3*a^-1*sinh(x)^3+-1//8*x*a^-1+-1//8*a^-1*cosh(x)*sinh(x)+1//4*a^-1*cosh(x)^3*sinh(x))
@test integrate((a+a*sech(x))^-1*sinh(x)^3, x) == :(-1//2*a^-1*sinh(x)^2+1//3*a^-1*cosh(x)^3)
@test integrate((a+a*sech(x))^-1*sinh(x)^2, x) == :((1/2)*x*a^-1+-1*a^-1*sinh(x)+(1/2)*a^-1*cosh(x)*sinh(x))
@test integrate((a+a*sech(x))^-1*sinh(x), x) == :(a^-1*cosh(x)+-1*a^-1*log(1+cosh(x)))
@test integrate((a+a*sech(x))^-1*csch(x), x) == :((1/2)*a^-1*csch(x)^2+-1//2*a^-1*arctanh(cosh(x))+-1//2*a^-1*coth(x)*csch(x))
@test integrate((a+a*sech(x))^-1*csch(x)^2, x) == :(-1//3*a^-1*coth(x)^3+1//3*a^-1*csch(x)^3)
@test integrate((a+a*sech(x))^-1*csch(x)^3, x) == :(1//4*a^-1*csch(x)^4+1//8*a^-1*arctanh(cosh(x))+-1//4*a^-1*csch(x)^3*coth(x)+-1//8*a^-1*coth(x)*csch(x))
@test integrate((a+a*sech(x))^-1*csch(x)^4, x) == :(-1//5*a^-1*coth(x)^5+1//3*a^-1*coth(x)^3+1//5*a^-1*csch(x)^5)
@test integrate((a+b*sech(x))^-1*sinh(x)^4, x) == :(-1//12*a^-2*sinh(x)^3*(4b+-3*a*cosh(x))+1//8*x*a^-5*(3*a^4+8*b^4+-12*a^2*b^2)+1//8*a^-4*(8*b*(a^2+-1*b^2)+-1*a*(-4*b^2+3*a^2)*cosh(x))*sinh(x)+-2*b*a^-5*(a+b)^3//2*(a+-1b)^3//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*sinh(x)^3, x) == :(1//3*a^-1*cosh(x)^3+-1*a^-3*(a^2+-1*b^2)*cosh(x)+-1//2*b*a^-2*cosh(x)^2+b*a^-4*(a^2+-1*b^2)*log(b+a*cosh(x)))
@test integrate((a+b*sech(x))^-1*sinh(x)^2, x) == :(-1//2*x*a^-3*(a^2+-2*b^2)+-1//2*a^-2*(2b+-1*a*cosh(x))*sinh(x)+2*b*a^-3*(a+b)^(1/2)*(a+-1b)^(1/2)*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*sinh(x), x) == :(a^-1*cosh(x)+-1*b*a^-2*log(b+a*cosh(x)))
@test integrate((a+b*sech(x))^-1*csch(x), x) == :((2a+2b)^-1*log(1+-1*cosh(x))+-1*(-2b+2a)^-1*log(1+cosh(x))+b*(a^2+-1*b^2)^-1*log(b+a*cosh(x)))
@test integrate((a+b*sech(x))^-1*csch(x)^2, x) == :((a^2+-1*b^2)^-1*(b+-1*a*cosh(x))*csch(x)+2*a*b*(a+b)^-3//2*(a+-1b)^-3//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*csch(x)^3, x) == :((-2*b^2+2*a^2)^-1*csch(x)^2*(b+-1*a*cosh(x))+-1//4*a*(a+b)^-2*log(1+-1*cosh(x))+1//4*a*(a+-1b)^-2*log(1+cosh(x))+-1*b*a^2*(a^2+-1*b^2)^-2*log(b+a*cosh(x)))
@test integrate((a+b*sech(x))^-1*csch(x)^4, x) == :((-3*b^2+3*a^2)^-1*csch(x)^3*(b+-1*a*cosh(x))+-1//3*(a^2+-1*b^2)^-2*(3*b*a^2+-1*a*(b^2+2*a^2)*cosh(x))*csch(x)+-2*b*a^3*(a+b)^-5//2*(a+-1b)^-5//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+a*sech(x))^-1*cosh(x)^4, x) == :(-4*a^-1*sinh(x)+-4//3*a^-1*sinh(x)^3+15//8*x*a^-1+-1*(a+a*sech(x))^-1*cosh(x)^3*sinh(x)+5//4*a^-1*cosh(x)^3*sinh(x)+15//8*a^-1*cosh(x)*sinh(x))
@test integrate((a+a*sech(x))^-1*cosh(x)^3, x) == :(4*a^-1*sinh(x)+-3//2*x*a^-1+4//3*a^-1*sinh(x)^3+-1*(a+a*sech(x))^-1*cosh(x)^2*sinh(x)+-3//2*a^-1*cosh(x)*sinh(x))
@test integrate((a+a*sech(x))^-1*cosh(x)^2, x) == :(-2*a^-1*sinh(x)+3//2*x*a^-1+-1*(a+a*sech(x))^-1*cosh(x)*sinh(x)+3//2*a^-1*cosh(x)*sinh(x))
@test integrate((a+a*sech(x))^-1*cosh(x), x) == :(-1*x*a^-1+-1*(a+a*sech(x))^-1*sinh(x)+2*a^-1*sinh(x))
@test integrate((a+a*sech(x))^-1*sech(x), x) == :((a+a*sech(x))^-1*tanh(x))
@test integrate((a+a*sech(x))^-1*sech(x)^2, x) == :(a^-1*arctan(sinh(x))+-1*(a+a*sech(x))^-1*tanh(x))
@test integrate((a+a*sech(x))^-1*sech(x)^3, x) == :(a^-1*tanh(x)+(a+a*sech(x))^-1*tanh(x)+-1*a^-1*arctan(sinh(x)))
@test integrate((a+a*sech(x))^-1*sech(x)^4, x) == :(-2*a^-1*tanh(x)+3//2*a^-1*arctan(sinh(x))+-1*(a+a*sech(x))^-1*sech(x)^2*tanh(x)+3//2*a^-1*sech(x)*tanh(x))
@test integrate((a+a*sech(c+d*x))^-1, x) == :(x*a^-1+-1*d^-1*(a+a*sech(c+d*x))^-1*tanh(c+d*x))
@test integrate((a+-1*a*sech(c+d*x))^-1, x) == :(x*a^-1+-1*d^-1*(a+-1*a*sech(c+d*x))^-1*tanh(c+d*x))
@test integrate((a+a*sech(c+d*x))^5//2, x) == :(2*a^5//2*d^-1*arctanh(a^(1/2)*(a+a*sech(c+d*x))^-1//2*tanh(c+d*x))+2//3*a^2*d^-1*(a+a*sech(c+d*x))^(1/2)*tanh(c+d*x)+14//3*a^3*d^-1*(a+a*sech(c+d*x))^-1//2*tanh(c+d*x))
@test integrate((a+a*sech(c+d*x))^3//2, x) == :(2*a^3//2*d^-1*arctanh(a^(1/2)*(a+a*sech(c+d*x))^-1//2*tanh(c+d*x))+2*a^2*d^-1*(a+a*sech(c+d*x))^-1//2*tanh(c+d*x))
@test integrate((a+a*sech(c+d*x))^(1/2), x) == :(2*a^(1/2)*d^-1*arctanh(a^(1/2)*(a+a*sech(c+d*x))^-1//2*tanh(c+d*x)))
@test integrate((a+a*sech(c+d*x))^-1//2, x) == :(2*a^-1//2*d^-1*arctanh(a^(1/2)*(a+a*sech(c+d*x))^-1//2*tanh(c+d*x))+-1*2^(1/2)*a^-1//2*d^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+a*sech(c+d*x))^-1//2*tanh(c+d*x)))
@test integrate((a+a*sech(c+d*x))^-3//2, x) == :(2*a^-3//2*d^-1*arctanh(a^(1/2)*(a+a*sech(c+d*x))^-1//2*tanh(c+d*x))+-1//2*d^-1*(a+a*sech(c+d*x))^-3//2*tanh(c+d*x)+-5//4*2^(1/2)*a^-3//2*d^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+a*sech(c+d*x))^-1//2*tanh(c+d*x)))
@test integrate((a+-1*a*sech(c+d*x))^(1/2), x) == :(2*a^(1/2)*d^-1*arctanh(a^(1/2)*(a+-1*a*sech(c+d*x))^-1//2*tanh(c+d*x)))
@test integrate((a+-1*a*sech(c+d*x))^-1//2, x) == :(2*a^-1//2*d^-1*arctanh(a^(1/2)*(a+-1*a*sech(c+d*x))^-1//2*tanh(c+d*x))+-1*2^(1/2)*a^-1//2*d^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+-1*a*sech(c+d*x))^-1//2*tanh(c+d*x)))
@test integrate((3+3*sech(x))^(1/2), x) == :(2*3^(1/2)*arctanh((1+sech(x))^-1//2*tanh(x)))
@test integrate((3+-3*sech(x))^(1/2), x) == :(2*3^(1/2)*arctanh((1+-1*sech(x))^-1//2*tanh(x)))
@test integrate((a+b*sech(c+d*x))^4, x) == :(x*a^4+1//3*b^2*d^-1*(a+b*sech(c+d*x))^2*tanh(c+d*x)+1//3*b^2*d^-1*(2*b^2+17*a^2)*tanh(c+d*x)+2*a*b*d^-1*(b^2+2*a^2)*arctan(sinh(c+d*x))+4//3*a*b^3*d^-1*sech(c+d*x)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x))^3, x) == :(x*a^3+(1/2)*b*d^-1*(b^2+6*a^2)*arctan(sinh(c+d*x))+(1/2)*b^2*d^-1*(a+b*sech(c+d*x))*tanh(c+d*x)+5//2*a*b^2*d^-1*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x))^2, x) == :(x*a^2+b^2*d^-1*tanh(c+d*x)+2*a*b*d^-1*arctan(sinh(c+d*x)))
@test integrate(a+b*sech(c+d*x), x) == :(a*x+b*d^-1*arctan(sinh(c+d*x)))
@test integrate((a+b*sech(c+d*x))^-1, x) == :(x*a^-1+-2*b*a^-1*d^-1*(a+b)^-1//2*(a+-1b)^-1//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*c+(1/2)*d*x)))
@test integrate((a+b*sech(c+d*x))^-2, x) == :(x*a^-2+a^-1*b^2*d^-1*(a+b*sech(c+d*x))^-1*(a^2+-1*b^2)^-1*tanh(c+d*x)+-2*b*a^-2*d^-1*(a+b)^-3//2*(a+-1b)^-3//2*(-1*b^2+2*a^2)*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*c+(1/2)*d*x)))
@test integrate((a+b*sech(c+d*x))^-3, x) == :(x*a^-3+(1/2)*a^-1*b^2*d^-1*(a+b*sech(c+d*x))^-2*(a^2+-1*b^2)^-1*tanh(c+d*x)+(1/2)*a^-2*b^2*d^-1*(a+b*sech(c+d*x))^-1*(a^2+-1*b^2)^-2*(-2*b^2+5*a^2)*tanh(c+d*x)+-1*b*a^-3*d^-1*(a+b)^-5//2*(a+-1b)^-5//2*(2*b^4+6*a^4+-5*a^2*b^2)*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*c+(1/2)*d*x)))
@test integrate((a+b*sech(x))^-1*cosh(x)^4, x) == :(1//4*a^-1*cosh(x)^3*sinh(x)+1//8*x*a^-5*(3*a^4+8*b^4+4*a^2*b^2)+-1//3*b*a^-4*(2*a^2+3*b^2)*sinh(x)+-1//3*b*a^-2*cosh(x)^2*sinh(x)+1//8*a^-3*(3*a^2+4*b^2)*cosh(x)*sinh(x)+-2*a^-5*b^5*(a+b)^-1//2*(a+-1b)^-1//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*cosh(x)^3, x) == :(1//3*a^-1*cosh(x)^2*sinh(x)+1//3*a^-3*(2*a^2+3*b^2)*sinh(x)+-1//2*b*x*a^-4*(a^2+2*b^2)+-1//2*b*a^-2*cosh(x)*sinh(x)+2*a^-4*b^4*(a+b)^-1//2*(a+-1b)^-1//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*cosh(x)^2, x) == :((1/2)*x*a^-3*(a^2+2*b^2)+(1/2)*a^-1*cosh(x)*sinh(x)+-1*b*a^-2*sinh(x)+-2*a^-3*b^3*(a+b)^-1//2*(a+-1b)^-1//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*cosh(x), x) == :(a^-1*sinh(x)+-1*b*x*a^-2+2*a^-2*b^2*(a+b)^-1//2*(a+-1b)^-1//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*sech(x), x) == :(2*(a+b)^-1//2*(a+-1b)^-1//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*sech(x)^2, x) == :(b^-1*arctan(sinh(x))+-2*a*b^-1*(a+b)^-1//2*(a+-1b)^-1//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*sech(x)^3, x) == :(b^-1*tanh(x)+-1*a*b^-2*arctan(sinh(x))+2*a^2*b^-2*(a+b)^-1//2*(a+-1b)^-1//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*sech(x)^4, x) == :((1/2)*b^-1*sech(x)*tanh(x)+(1/2)*b^-3*(b^2+2*a^2)*arctan(sinh(x))+-1*a*b^-2*tanh(x)+-2*a^3*b^-3*(a+b)^-1//2*(a+-1b)^-1//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+a*sech(x))^-1*tanh(x)^6, x) == :(x*a^-1+-3//8*a^-1*arctan(sinh(x))+-1//8*a^-1*(8+-3*sech(x))*tanh(x)+-1//12*a^-1*tanh(x)^3*(4+-3*sech(x)))
@test integrate((a+a*sech(x))^-1*tanh(x)^5, x) == :(a^-1*log(cosh(x))+a^-1*sech(x)+(1/2)*a^-1*sech(x)^2+-1//3*a^-1*sech(x)^3)
@test integrate((a+a*sech(x))^-1*tanh(x)^4, x) == :(x*a^-1+-1//2*a^-1*arctan(sinh(x))+-1//2*a^-1*(2+-1*sech(x))*tanh(x))
@test integrate((a+a*sech(x))^-1*tanh(x)^3, x) == :(a^-1*log(cosh(x))+a^-1*sech(x))
@test integrate((a+a*sech(x))^-1*tanh(x)^2, x) == :(x*a^-1+-1*a^-1*arctan(sinh(x)))
@test integrate((a+a*sech(x))^-1*tanh(x), x) == :(a^-1*log(1+cosh(x)))
@test integrate((a+a*sech(x))^-1*coth(x), x) == :((1/2)*a^-1*(1+cosh(x))^-1+1//4*a^-1*log(1+-1*cosh(x))+3//4*a^-1*log(1+cosh(x)))
@test integrate((a+a*sech(x))^-1*coth(x)^2, x) == :(x*a^-1+-1//3*a^-1*coth(x)^3*(1+-1*sech(x))+-1//3*a^-1*(3+-2*sech(x))*coth(x))
@test integrate((a+a*sech(x))^-1*coth(x)^3, x) == :(-1//8*a^-1*(1+cosh(x))^-2+1//8*a^-1*(1+-1*cosh(x))^-1+3//4*a^-1*(1+cosh(x))^-1+5//16*a^-1*log(1+-1*cosh(x))+11//16*a^-1*log(1+cosh(x)))
@test integrate((a+a*sech(x))^-1*coth(x)^4, x) == :(x*a^-1+-1//5*a^-1*coth(x)^5*(1+-1*sech(x))+-1//15*a^-1*coth(x)^3*(5+-4*sech(x))+-1//15*a^-1*(15+-8*sech(x))*coth(x))
@test integrate((a+b*sech(x))^-1*tanh(x)^7, x) == :(a^-1*log(cosh(x))+1//5*b^-1*sech(x)^5+b^-5*(a^4+3*b^4+-3*a^2*b^2)*sech(x)+-1//4*a*b^-2*sech(x)^4+1//3*b^-3*sech(x)^3*(a^2+-3*b^2)+-1*a^-1*b^-6*(a^2+-1*b^2)^3*log(a+b*sech(x))+-1//2*a*b^-4*sech(x)^2*(a^2+-3*b^2))
@test integrate((a+b*sech(x))^-1*tanh(x)^6, x) == :(x*a^-1+-3//8*b^-1*arctan(sinh(x))+a*b^-2*tanh(x)+-1*b^-5*(a^4+3*b^4+-3*a^2*b^2)*arctan(sinh(x))+-3//8*b^-1*sech(x)*tanh(x)+-1//2*b^-3*(a^2+-3*b^2)*arctan(sinh(x))+-1//3*a*b^-2*tanh(x)^3+-1//4*b^-1*sech(x)^3*tanh(x)+a*b^-4*(a^2+-3*b^2)*tanh(x)+-1//2*b^-3*(a^2+-3*b^2)*sech(x)*tanh(x)+2*a^-1*b^-5*(a+b)^5//2*(a+-1b)^5//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*tanh(x)^5, x) == :(a^-1*log(cosh(x))+-1//3*b^-1*sech(x)^3+(1/2)*a*b^-2*sech(x)^2+-1*b^-3*(a^2+-2*b^2)*sech(x)+a^-1*b^-4*(a^2+-1*b^2)^2*log(a+b*sech(x)))
@test integrate((a+b*sech(x))^-1*tanh(x)^4, x) == :(x*a^-1+(1/2)*b^-1*sech(x)*tanh(x)+(1/2)*b^-3*(-3*b^2+2*a^2)*arctan(sinh(x))+-1*a*b^-2*tanh(x)+-2*a^-1*b^-3*(a+b)^3//2*(a+-1b)^3//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*tanh(x)^3, x) == :(a^-1*log(cosh(x))+b^-1*sech(x)+a^-1*(1+-1*a^2*b^-2)*log(a+b*sech(x)))
@test integrate((a+b*sech(x))^-1*tanh(x)^2, x) == :(x*a^-1+-1*b^-1*arctan(sinh(x))+2*a^-1*b^-1*(a+b)^(1/2)*(a+-1b)^(1/2)*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*tanh(x), x) == :(a^-1*log(a+b*sech(x))+a^-1*log(cosh(x)))
@test integrate((a+b*sech(x))^-1*coth(x), x) == :(a^-1*log(cosh(x))+(-2b+2a)^-1*log(1+sech(x))+(2a+2b)^-1*log(1+-1*sech(x))+-1*a^-1*b^2*(a^2+-1*b^2)^-1*log(a+b*sech(x)))
@test integrate((a+b*sech(x))^-1*coth(x)^2, x) == :(a*x*(a^2+-1*b^2)^-1+b*(a^2+-1*b^2)^-1*csch(x)+-1*a*(a^2+-1*b^2)^-1*coth(x)+-1*x*a^-1*b^2*(a^2+-1*b^2)^-1+2*a^-1*b^3*(a+b)^-3//2*(a+-1b)^-3//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*coth(x)^3, x) == :(a^-1*log(cosh(x))+-1*(1+-1*sech(x))^-1*(4a+4b)^-1+-1*(1+sech(x))^-1*(-4b+4a)^-1+1//4*(a+b)^-2*(2a+3b)*log(1+-1*sech(x))+1//4*(a+-1b)^-2*(-3b+2a)*log(1+sech(x))+a^-1*b^4*(a^2+-1*b^2)^-2*log(a+b*sech(x)))
@test integrate((a+b*sech(x))^-1*coth(x)^4, x) == :(a*x*(a^2+-1*b^2)^-1+b*(a^2+-1*b^2)^-1*csch(x)+b*(-3*b^2+3*a^2)^-1*csch(x)^3+-1*a*(a^2+-1*b^2)^-1*coth(x)+-1*a*(-3*b^2+3*a^2)^-1*coth(x)^3+-1*b^3*(a^2+-1*b^2)^-2*csch(x)+a*b^2*(a^2+-1*b^2)^-2*coth(x)+x*a^-1*b^4*(a^2+-1*b^2)^-2+-1*a*x*b^2*(a^2+-1*b^2)^-2+-2*a^-1*b^5*(a+b)^-5//2*(a+-1b)^-5//2*arctan((a+b)^-1//2*(a+-1b)^(1/2)*tanh((1/2)*x)))
@test integrate((a+b*sech(x))^-1*coth(x)^5, x) == :(a^-1*log(cosh(x))+-1*(1+-1*sech(x))^-2*(16a+16b)^-1+-1*(1+sech(x))^-2*(-16b+16a)^-1+-1//16*(1+-1*sech(x))^-1*(a+b)^-2*(5a+7b)+-1//16*(1+sech(x))^-1*(a+-1b)^-2*(-7b+5a)+1//16*(a+b)^-3*(8*a^2+15*b^2+21*a*b)*log(1+-1*sech(x))+1//16*(a+-1b)^-3*(8*a^2+15*b^2+-21*a*b)*log(1+sech(x))+-1*a^-1*b^6*(a^2+-1*b^2)^-3*log(a+b*sech(x)))
@test integrate((a+b*sech(c+d*x))^(1/2)*tanh(c+d*x)^5, x) == :(-2*d^-1*(a+b*sech(c+d*x))^(1/2)+2*a^(1/2)*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2))+-2//9*b^-4*d^-1*(a+b*sech(c+d*x))^9//2+-1//5*b^-4*d^-1*(a+b*sech(c+d*x))^5//2*(-4*b^2+6*a^2)+6//7*a*b^-4*d^-1*(a+b*sech(c+d*x))^7//2+2//3*a*b^-4*d^-1*(a+b*sech(c+d*x))^3//2*(a^2+-2*b^2))
@test integrate((a+b*sech(c+d*x))^(1/2)*tanh(c+d*x)^3, x) == :(-2*d^-1*(a+b*sech(c+d*x))^(1/2)+2*a^(1/2)*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2))+2//5*b^-2*d^-1*(a+b*sech(c+d*x))^5//2+-2//3*a*b^-2*d^-1*(a+b*sech(c+d*x))^3//2)
@test integrate((a+b*sech(c+d*x))^(1/2)*tanh(c+d*x), x) == :(-2*d^-1*(a+b*sech(c+d*x))^(1/2)+2*a^(1/2)*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2)))
@test integrate((a+b*sech(c+d*x))^(1/2)*coth(c+d*x), x) == :(-1*d^-1*(a+b)^(1/2)*arctanh((a+b)^-1//2*(a+b*sech(c+d*x))^(1/2))+-1*d^-1*(a+-1b)^(1/2)*arctanh((a+-1b)^-1//2*(a+b*sech(c+d*x))^(1/2))+2*a^(1/2)*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2)))
@test integrate((a+b*sech(c+d*x))^(1/2)*coth(c+d*x)^3, x) == :(2*a^(1/2)*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2))+-1//2*d^-1*(a+b*sech(c+d*x))^(1/2)*coth(c+d*x)^2+-1*a*d^-1*(a+b)^-1//2*arctanh((a+b)^-1//2*(a+b*sech(c+d*x))^(1/2))+-1*a*d^-1*(a+-1b)^-1//2*arctanh((a+-1b)^-1//2*(a+b*sech(c+d*x))^(1/2))+-3//4*b*d^-1*(a+b)^-1//2*arctanh((a+b)^-1//2*(a+b*sech(c+d*x))^(1/2))+3//4*b*d^-1*(a+-1b)^-1//2*arctanh((a+-1b)^-1//2*(a+b*sech(c+d*x))^(1/2)))
@test integrate((a+b*sech(c+d*x))^-1//2*tanh(c+d*x)^5, x) == :(2*a^-1//2*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2))+-2//7*b^-4*d^-1*(a+b*sech(c+d*x))^7//2+-1//3*b^-4*d^-1*(a+b*sech(c+d*x))^3//2*(-4*b^2+6*a^2)+6//5*a*b^-4*d^-1*(a+b*sech(c+d*x))^5//2+2*a*b^-4*d^-1*(a+b*sech(c+d*x))^(1/2)*(a^2+-2*b^2))
@test integrate((a+b*sech(c+d*x))^-1//2*tanh(c+d*x)^3, x) == :(2*a^-1//2*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2))+2//3*b^-2*d^-1*(a+b*sech(c+d*x))^3//2+-2*a*b^-2*d^-1*(a+b*sech(c+d*x))^(1/2))
@test integrate((a+b*sech(c+d*x))^-1//2*tanh(c+d*x), x) == :(2*a^-1//2*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2)))
@test integrate((a+b*sech(c+d*x))^-1//2*coth(c+d*x), x) == :(-1*d^-1*(a+b)^-1//2*arctanh((a+b)^-1//2*(a+b*sech(c+d*x))^(1/2))+-1*d^-1*(a+-1b)^-1//2*arctanh((a+-1b)^-1//2*(a+b*sech(c+d*x))^(1/2))+2*a^-1//2*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2)))
@test integrate((a+b*sech(c+d*x))^-1//2*coth(c+d*x)^3, x) == :(-1*d^-1*(a+b)^-1//2*arctanh((a+b)^-1//2*(a+b*sech(c+d*x))^(1/2))+-1*d^-1*(a+-1b)^-1//2*arctanh((a+-1b)^-1//2*(a+b*sech(c+d*x))^(1/2))+2*a^-1//2*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2))+-1*d^-1*(1+-1*sech(c+d*x))^-1*(a+b*sech(c+d*x))^(1/2)*(4a+4b)^-1+-1*d^-1*(1+sech(c+d*x))^-1*(a+b*sech(c+d*x))^(1/2)*(-4b+4a)^-1+-1//4*b*d^-1*(a+b)^-3//2*arctanh((a+b)^-1//2*(a+b*sech(c+d*x))^(1/2))+1//4*b*d^-1*(a+-1b)^-3//2*arctanh((a+-1b)^-1//2*(a+b*sech(c+d*x))^(1/2)))
@test integrate((a+b*sech(c+d*x))^-3//2*tanh(c+d*x)^5, x) == :(2*a^-3//2*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2))+-2//5*b^-4*d^-1*(a+b*sech(c+d*x))^5//2+-1*b^-4*d^-1*(a+b*sech(c+d*x))^(1/2)*(-4*b^2+6*a^2)+2*a*b^-4*d^-1*(a+b*sech(c+d*x))^3//2+-2*a^-1*b^-4*d^-1*(a+b*sech(c+d*x))^-1//2*(a^2+-1*b^2)^2)
@test integrate((a+b*sech(c+d*x))^-3//2*tanh(c+d*x)^3, x) == :(2*a^-3//2*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2))+2*b^-2*d^-1*(a+b*sech(c+d*x))^(1/2)+a^-1*b^-2*d^-1*(a+b*sech(c+d*x))^-1//2*(-2*b^2+2*a^2))
@test integrate((a+b*sech(c+d*x))^-3//2*tanh(c+d*x), x) == :(-2*a^-1*d^-1*(a+b*sech(c+d*x))^-1//2+2*a^-3//2*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2)))
@test integrate((a+b*sech(c+d*x))^-3//2*coth(c+d*x), x) == :(-1*d^-1*(a+b)^-3//2*arctanh((a+b)^-1//2*(a+b*sech(c+d*x))^(1/2))+-1*d^-1*(a+-1b)^-3//2*arctanh((a+-1b)^-1//2*(a+b*sech(c+d*x))^(1/2))+2*a^-3//2*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2))+2*a^-1*b^2*d^-1*(a+b*sech(c+d*x))^-1//2*(a^2+-1*b^2)^-1)
@test integrate((a+b*sech(c+d*x))^-3//2*coth(c+d*x)^3, x) == :(2*a^-3//2*d^-1*arctanh(a^-1//2*(a+b*sech(c+d*x))^(1/2))+-1//2*d^-1*(a+b)^-5//2*(2a+3b)*arctanh((a+b)^-1//2*(a+b*sech(c+d*x))^(1/2))+-1//2*d^-1*(a+-1b)^-5//2*(-3b+2a)*arctanh((a+-1b)^-1//2*(a+b*sech(c+d*x))^(1/2))+-1//4*b*d^-1*(a+b)^-5//2*arctanh((a+b)^-1//2*(a+b*sech(c+d*x))^(1/2))+-1//4*d^-1*(1+-1*sech(c+d*x))^-1*(a+b)^-2*(a+b*sech(c+d*x))^(1/2)+-1//4*d^-1*(1+sech(c+d*x))^-1*(a+-1b)^-2*(a+b*sech(c+d*x))^(1/2)+1//4*b*d^-1*(a+-1b)^-5//2*arctanh((a+-1b)^-1//2*(a+b*sech(c+d*x))^(1/2))+-2*a^-1*b^4*d^-1*(a+b*sech(c+d*x))^-1//2*(a^2+-1*b^2)^-2)
@test integrate((sech(a*c+b*c*x)^2)^7//2*exp(c*(a+b*x)), x) == :(48*b^-1*c^-1*(sech(a*c+b*c*x)^2)^(1/2)*(1+exp(2*c*(a+b*x)))^-4*cosh(a*c+b*c*x)+-192//5*b^-1*c^-1*(sech(a*c+b*c*x)^2)^(1/2)*(1+exp(2*c*(a+b*x)))^-5*cosh(a*c+b*c*x)+-64//3*b^-1*c^-1*(sech(a*c+b*c*x)^2)^(1/2)*(1+exp(2*c*(a+b*x)))^-3*cosh(a*c+b*c*x)+32//3*b^-1*c^-1*(sech(a*c+b*c*x)^2)^(1/2)*(1+exp(2*c*(a+b*x)))^-6*cosh(a*c+b*c*x))
@test integrate((sech(a*c+b*c*x)^2)^5//2*exp(c*(a+b*x)), x) == :(-8*b^-1*c^-1*(sech(a*c+b*c*x)^2)^(1/2)*(1+exp(2*c*(a+b*x)))^-2*cosh(a*c+b*c*x)+-4*b^-1*c^-1*(sech(a*c+b*c*x)^2)^(1/2)*(1+exp(2*c*(a+b*x)))^-4*cosh(a*c+b*c*x)+32//3*b^-1*c^-1*(sech(a*c+b*c*x)^2)^(1/2)*(1+exp(2*c*(a+b*x)))^-3*cosh(a*c+b*c*x))
@test integrate((sech(a*c+b*c*x)^2)^3//2*exp(c*(a+b*x)), x) == :(2*b^-1*c^-1*(sech(a*c+b*c*x)^2)^(1/2)*(1+exp(2*c*(a+b*x)))^-2*cosh(a*c+b*c*x)*exp(4*c*(a+b*x)))
@test integrate((sech(a*c+b*c*x)^2)^(1/2)*exp(c*(a+b*x)), x) == :(b^-1*c^-1*(sech(a*c+b*c*x)^2)^(1/2)*cosh(a*c+b*c*x)*log(1+exp(2*c*(a+b*x))))
@test integrate((sech(a*c+b*c*x)^2)^-1//2*exp(c*(a+b*x)), x) == :((1/2)*x*(sech(a*c+b*c*x)^2)^-1//2*sech(a*c+b*c*x)+1//4*b^-1*c^-1*(sech(a*c+b*c*x)^2)^-1//2*exp(2*c*(a+b*x))*sech(a*c+b*c*x))
@test integrate((sech(a*c+b*c*x)^2)^-3//2*exp(c*(a+b*x)), x) == :(3//8*x*(sech(a*c+b*c*x)^2)^-1//2*sech(a*c+b*c*x)+-1//16*b^-1*c^-1*(sech(a*c+b*c*x)^2)^-1//2*exp(-2*c*(a+b*x))*sech(a*c+b*c*x)+1//32*b^-1*c^-1*(sech(a*c+b*c*x)^2)^-1//2*exp(4*c*(a+b*x))*sech(a*c+b*c*x)+3//16*b^-1*c^-1*(sech(a*c+b*c*x)^2)^-1//2*exp(2*c*(a+b*x))*sech(a*c+b*c*x))
@test integrate((sech(a*c+b*c*x)^2)^-5//2*exp(c*(a+b*x)), x) == :(5//16*x*(sech(a*c+b*c*x)^2)^-1//2*sech(a*c+b*c*x)+-5//64*b^-1*c^-1*(sech(a*c+b*c*x)^2)^-1//2*exp(-2*c*(a+b*x))*sech(a*c+b*c*x)+-1//128*b^-1*c^-1*(sech(a*c+b*c*x)^2)^-1//2*exp(-4*c*(a+b*x))*sech(a*c+b*c*x)+1//192*b^-1*c^-1*(sech(a*c+b*c*x)^2)^-1//2*exp(6*c*(a+b*x))*sech(a*c+b*c*x)+5//32*b^-1*c^-1*(sech(a*c+b*c*x)^2)^-1//2*exp(2*c*(a+b*x))*sech(a*c+b*c*x)+5//128*b^-1*c^-1*(sech(a*c+b*c*x)^2)^-1//2*exp(4*c*(a+b*x))*sech(a*c+b*c*x))
@test integrate(x^5*sech(2*log(c*x))^-1//2, x) == :(1//7*x^6*sech(2*log(c*x))^-1//2+2//21*c^-4*x^2*sech(2*log(c*x))^-1//2+1//21*c^-5*x^-1*((c^2+x^-2)^-2*(c^4+x^-4))^(1/2)*(c^4+x^-4)^-1*sech(2*log(c*x))^-1//2*(c^2+x^-2)*Elliptic.F(2*arccot(c*x),1/2))
@test integrate(x^4*sech(2*log(c*x))^-1//2, x) == :(1//6*c^-4*x^5*sech(2*log(c*x))^-1//2*(c^4+x^-4))
@test integrate(x^3*sech(2*log(c*x))^-1//2, x) == :(1//5*x^4*sech(2*log(c*x))^-1//2+2//5*c^-4*sech(2*log(c*x))^-1//2+-2//5*c^-4*x^-2*(c^2+x^-2)^-1*sech(2*log(c*x))^-1//2+-1//5*c^-3*x^-1*((c^2+x^-2)^-2*(c^4+x^-4))^(1/2)*(c^4+x^-4)^-1*sech(2*log(c*x))^-1//2*(c^2+x^-2)*Elliptic.F(2*arccot(c*x),1/2)+2//5*c^-3*x^-1*((c^2+x^-2)^-2*(c^4+x^-4))^(1/2)*(c^4+x^-4)^-1*sech(2*log(c*x))^-1//2*(c^2+x^-2)*Elliptic.E(2*arccot(c*x),1/2))
@test integrate(x^2*sech(2*log(c*x))^-1//2, x) == :(1//4*x^3*sech(2*log(c*x))^-1//2+1//4*c^-4*x^-1*(1+c^-4*x^-4)^-1//2*sech(2*log(c*x))^-1//2*arctanh((1+c^-4*x^-4)^(1/2)))
@test integrate(x*sech(2*log(c*x))^-1//2, x) == :(1//3*x^2*sech(2*log(c*x))^-1//2+-1//3*c^-1*x^-1*((c^2+x^-2)^-2*(c^4+x^-4))^(1/2)*(c^4+x^-4)^-1*sech(2*log(c*x))^-1//2*(c^2+x^-2)*Elliptic.F(2*arccot(c*x),1/2))
@test integrate(sech(2*log(c*x))^-1//2, x) == :((1/2)*x*sech(2*log(c*x))^-1//2+-1//2*c^-2*x^-1*(1+c^-4*x^-4)^-1//2*sech(2*log(c*x))^-1//2*arccsch(c^2*x^2))
@test integrate(x^-1*sech(2*log(c*x))^(1/2), x) == :(-1*im*cosh(2*log(c*x))^(1/2)*sech(2*log(c*x))^(1/2)*Elliptic.F(im*log(c*x),2))
@test integrate(x^-2*sech(2*log(c*x))^(1/2), x) == :(-1//2*x*c^2*(1+c^-4*x^-4)^(1/2)*sech(2*log(c*x))^(1/2)*arccsch(c^2*x^2))
@test integrate(x^-3*sech(2*log(c*x))^(1/2), x) == :(-1*(c^2+x^-2)^-1*sech(2*log(c*x))^(1/2)*(c^4+x^-4)+c*x*((c^2+x^-2)^-2*(c^4+x^-4))^(1/2)*sech(2*log(c*x))^(1/2)*(c^2+x^-2)*Elliptic.E(2*arccot(c*x),1/2)+-1//2*c*x*((c^2+x^-2)^-2*(c^4+x^-4))^(1/2)*sech(2*log(c*x))^(1/2)*(c^2+x^-2)*Elliptic.F(2*arccot(c*x),1/2))
@test integrate(x^-4*sech(2*log(c*x))^(1/2), x) == :(x*sech(2*log(c*x))^(1/2)*(-1//2*c^4+-1//2*x^-4))
@test integrate(x^-5*sech(2*log(c*x))^(1/2), x) == :(sech(2*log(c*x))^(1/2)*(-1//3*c^4+-1//3*x^-4)+1//6*x*c^3*((c^2+x^-2)^-2*(c^4+x^-4))^(1/2)*sech(2*log(c*x))^(1/2)*(c^2+x^-2)*Elliptic.F(2*arccot(c*x),1/2))
@test integrate(x^8*sech(2*log(c*x))^-3//2, x) == :(1//12*x^9*sech(2*log(c*x))^-3//2+x^5*(16*c^4+16*x^-4)^-1*sech(2*log(c*x))^-3//2+1//32*x*c^-4*(c^4+x^-4)^-1*sech(2*log(c*x))^-3//2+-1//32*c^-12*x^-3*(1+c^-4*x^-4)^-3//2*sech(2*log(c*x))^-3//2*arctanh((1+c^-4*x^-4)^(1/2)))
@test integrate(x^7*sech(2*log(c*x))^-3//2, x) == :(1//11*x^8*sech(2*log(c*x))^-3//2+6*x^4*(77*c^4+77*x^-4)^-1*sech(2*log(c*x))^-3//2+4//77*c^-4*(c^4+x^-4)^-1*sech(2*log(c*x))^-3//2+2//77*c^-5*x^-3*((c^2+x^-2)^-2*(c^4+x^-4))^(1/2)*(c^4+x^-4)^-2*sech(2*log(c*x))^-3//2*(c^2+x^-2)*Elliptic.F(2*arccot(c*x),1/2))
@test integrate(x^6*sech(2*log(c*x))^-3//2, x) == :(1//10*c^-4*x^7*sech(2*log(c*x))^-3//2*(c^4+x^-4))
@test integrate(x^5*sech(2*log(c*x))^-3//2, x) == :(1//9*x^6*sech(2*log(c*x))^-3//2+2*x^2*(15*c^4+15*x^-4)^-1*sech(2*log(c*x))^-3//2+4//15*c^-4*x^-2*(c^4+x^-4)^-1*sech(2*log(c*x))^-3//2+-4//15*c^-4*x^-4*(c^2+x^-2)^-1*(c^4+x^-4)^-1*sech(2*log(c*x))^-3//2+-2//15*c^-3*x^-3*((c^2+x^-2)^-2*(c^4+x^-4))^(1/2)*(c^4+x^-4)^-2*sech(2*log(c*x))^-3//2*(c^2+x^-2)*Elliptic.F(2*arccot(c*x),1/2)+4//15*c^-3*x^-3*((c^2+x^-2)^-2*(c^4+x^-4))^(1/2)*(c^4+x^-4)^-2*sech(2*log(c*x))^-3//2*(c^2+x^-2)*Elliptic.E(2*arccot(c*x),1/2))
@test integrate(x^4*sech(2*log(c*x))^-3//2, x) == :(1//8*x^5*sech(2*log(c*x))^-3//2+3*x*(16*c^4+16*x^-4)^-1*sech(2*log(c*x))^-3//2+3//16*c^-8*x^-3*(1+c^-4*x^-4)^-3//2*sech(2*log(c*x))^-3//2*arctanh((1+c^-4*x^-4)^(1/2)))
@test integrate(x^3*sech(2*log(c*x))^-3//2, x) == :(2*(7*c^4+7*x^-4)^-1*sech(2*log(c*x))^-3//2+1//7*x^4*sech(2*log(c*x))^-3//2+-2//7*c^-1*x^-3*((c^2+x^-2)^-2*(c^4+x^-4))^(1/2)*(c^4+x^-4)^-2*sech(2*log(c*x))^-3//2*(c^2+x^-2)*Elliptic.F(2*arccot(c*x),1/2))
@test integrate(x^2*sech(2*log(c*x))^-3//2, x) == :(1//6*x^3*sech(2*log(c*x))^-3//2+x^-1*(2*c^4+2*x^-4)^-1*sech(2*log(c*x))^-3//2+-1//2*c^-6*x^-3*(1+c^-4*x^-4)^-3//2*sech(2*log(c*x))^-3//2*arccsch(c^2*x^2))
@test integrate(x*sech(2*log(c*x))^-3//2, x) == :(1//5*x^2*sech(2*log(c*x))^-3//2+6*x^-2*(5*c^4+5*x^-4)^-1*sech(2*log(c*x))^-3//2+-12*x^-4*(c^2+x^-2)^-1*(5*c^4+5*x^-4)^-1*sech(2*log(c*x))^-3//2+-6//5*c*x^-3*((c^2+x^-2)^-2*(c^4+x^-4))^(1/2)*(c^4+x^-4)^-2*sech(2*log(c*x))^-3//2*(c^2+x^-2)*Elliptic.F(2*arccot(c*x),1/2)+12//5*c*x^-3*((c^2+x^-2)^-2*(c^4+x^-4))^(1/2)*(c^4+x^-4)^-2*sech(2*log(c*x))^-3//2*(c^2+x^-2)*Elliptic.E(2*arccot(c*x),1/2))
@test integrate(sech(2*log(c*x))^-3//2, x) == :(1//4*x*sech(2*log(c*x))^-3//2+-3*x^-3*(4*c^4+4*x^-4)^-1*sech(2*log(c*x))^-3//2+3//4*c^-4*x^-3*(1+c^-4*x^-4)^-3//2*sech(2*log(c*x))^-3//2*arctanh((1+c^-4*x^-4)^(1/2)))
@test integrate(x^-1*sech(2*log(c*x))^3//2, x) == :(sech(2*log(c*x))^(1/2)*sinh(2*log(c*x))+im*cosh(2*log(c*x))^(1/2)*sech(2*log(c*x))^(1/2)*Elliptic.E(im*log(c*x),2))
@test integrate(x^-2*sech(2*log(c*x))^3//2, x) == :(x^3*sech(2*log(c*x))^3//2*((1/2)*c^4+(1/2)*x^-4))
@test integrate(x^-3*sech(2*log(c*x))^3//2, x) == :(x^2*sech(2*log(c*x))^3//2*((1/2)*c^4+(1/2)*x^-4)+-1//4*c^-1*x^3*((c^2+x^-2)^-2*(c^4+x^-4))^(1/2)*sech(2*log(c*x))^3//2*(c^2+x^-2)*(c^4+x^-4)*Elliptic.F(2*arccot(c*x),1/2))
@test integrate(x^-4*sech(2*log(c*x))^3//2, x) == :(x*sech(2*log(c*x))^3//2*((1/2)*c^4+(1/2)*x^-4)+-1//2*c^6*x^3*(1+c^-4*x^-4)^3//2*sech(2*log(c*x))^3//2*arccsch(c^2*x^2))
@test integrate((1+-1*b^2*n^2)*sech(a+b*log(c*x^n))+2*b^2*n^2*sech(a+b*log(c*x^n))^3, x) == :(x*sech(a+b*log(c*x^n))+b*n*x*sech(a+b*log(c*x^n))*tanh(a+b*log(c*x^n)))
@test integrate(sech(a+2*log(c*x^(1/2)))^3, x) == :(2*c^6*(c^4+x^-2*exp(-2a))^-2*exp(-1a))
@test integrate(sech(a+2*log(c*x^-1//2))^3, x) == :(2*c^2*(c^4*x^-2+exp(-2a))^-2*exp(-3a))
@test integrate(sech(a+n^-1*(-2+p)^-1*log(c*x^n))^p, x) == :(x*(c*x^n)^(-2*n^-1*(2+-1p)^-1)*(2+-2p)^-1*sech(a+-1*n^-1*(2+-1p)^-1*log(c*x^n))^p*(1+(c*x^n)^(2*n^-1*(2+-1p)^-1)*exp(-2a))*(2+-1p)*exp(2a))
@test integrate(sech(a+-1*n^-1*(-2+p)^-1*log(c*x^n))^p, x) == :(x*(2+-2p)^-1*sech(a+n^-1*(2+-1p)^-1*log(c*x^n))^p*(1+(c*x^n)^(-2*n^-1*(2+-1p)^-1)*exp(-2a))*(2+-1p))
@test integrate(x^-1*sech(a+b*log(c*x^n)), x) == :(b^-1*n^-1*arctan(sinh(a+b*log(c*x^n))))
@test integrate(x^-1*sech(a+b*log(c*x^n))^2, x) == :(b^-1*n^-1*tanh(a+b*log(c*x^n)))
@test integrate(x^-1*sech(a+b*log(c*x^n))^3, x) == :((1/2)*b^-1*n^-1*arctan(sinh(a+b*log(c*x^n)))+(1/2)*b^-1*n^-1*sech(a+b*log(c*x^n))*tanh(a+b*log(c*x^n)))
@test integrate(x^-1*sech(a+b*log(c*x^n))^4, x) == :(b^-1*n^-1*tanh(a+b*log(c*x^n))+-1//3*b^-1*n^-1*tanh(a+b*log(c*x^n))^3)
@test integrate(x^-1*sech(a+b*log(c*x^n))^5, x) == :(3//8*b^-1*n^-1*arctan(sinh(a+b*log(c*x^n)))+1//4*b^-1*n^-1*sech(a+b*log(c*x^n))^3*tanh(a+b*log(c*x^n))+3//8*b^-1*n^-1*sech(a+b*log(c*x^n))*tanh(a+b*log(c*x^n)))
@test integrate(x^-1*sech(a+b*log(c*x^n))^5//2, x) == :(2//3*b^-1*n^-1*sech(a+b*log(c*x^n))^3//2*sinh(a+b*log(c*x^n))+-2//3*im*b^-1*n^-1*cosh(a+b*log(c*x^n))^(1/2)*sech(a+b*log(c*x^n))^(1/2)*Elliptic.F((1/2)*im*(a+b*log(c*x^n)),2))
@test integrate(x^-1*sech(a+b*log(c*x^n))^3//2, x) == :(2*b^-1*n^-1*sech(a+b*log(c*x^n))^(1/2)*sinh(a+b*log(c*x^n))+2*im*b^-1*n^-1*cosh(a+b*log(c*x^n))^(1/2)*sech(a+b*log(c*x^n))^(1/2)*Elliptic.E((1/2)*im*(a+b*log(c*x^n)),2))
@test integrate(x^-1*sech(a+b*log(c*x^n))^(1/2), x) == :(-2*im*b^-1*n^-1*cosh(a+b*log(c*x^n))^(1/2)*sech(a+b*log(c*x^n))^(1/2)*Elliptic.F((1/2)*im*(a+b*log(c*x^n)),2))
@test integrate(x^-1*sech(a+b*log(c*x^n))^-1//2, x) == :(-2*im*b^-1*n^-1*cosh(a+b*log(c*x^n))^(1/2)*sech(a+b*log(c*x^n))^(1/2)*Elliptic.E((1/2)*im*(a+b*log(c*x^n)),2))
@test integrate(x^-1*sech(a+b*log(c*x^n))^-3//2, x) == :(2//3*b^-1*n^-1*sech(a+b*log(c*x^n))^-1//2*sinh(a+b*log(c*x^n))+-2//3*im*b^-1*n^-1*cosh(a+b*log(c*x^n))^(1/2)*sech(a+b*log(c*x^n))^(1/2)*Elliptic.F((1/2)*im*(a+b*log(c*x^n)),2))
@test integrate(x^-1*sech(a+b*log(c*x^n))^-5//2, x) == :(2//5*b^-1*n^-1*sech(a+b*log(c*x^n))^-3//2*sinh(a+b*log(c*x^n))+-6//5*im*b^-1*n^-1*cosh(a+b*log(c*x^n))^(1/2)*sech(a+b*log(c*x^n))^(1/2)*Elliptic.E((1/2)*im*(a+b*log(c*x^n)),2))

max_stars_repo_path: test/6 Hyperbolic functions/6.5 Hyperbolic secant/6.5.3 Hyperbolic secant functions.jl


text:
function sl_sip(h, b, K, CV, c, μ::Vector, start_inventory, LT::Int = 0; lostsales = false, horizon = length(μ), policy = sSPolicy(), periods = horizon)
    item = EndProduct(
        Market(b, CVNormal{CV}, horizon, 0, μ, lostsales = lostsales),
        Inventory(h, start_inventory),
        Supplier(K,c, leadtime = LeadTime(LT, 0)),
        policy = policy
    )
    InventorySystem(periods, [item])
end

function sl_sip(h, b, K, CV, c, μ::Distribution, start_inventory, LT::Int = 0; lostsales = false, policy = sSPolicy(), horizon::Int, periods::Int)
    item = EndProduct(
        Market(b, CVNormal{CV}, horizon, 0, μ, lostsales = lostsales),
        Inventory(h, start_inventory),
        Supplier(K,c, leadtime = LeadTime(LT, 0)),
        policy = policy
    )
    InventorySystem(periods, [item])
end

"""
InventoryModels.assembly10(holding_costs = ones(10), fixed_costs = 100*ones(10), variable_costs = zeros(10), leadtimes=ones(10), backorder_cost = 10, demand_forecast = [20,40,60,40,20,40,60], CV = 0.25)

order capacities are decentralized, this does not implements the "shared ressources" generalization of Thevenin et al.
"""
#=
function assembly10(;holding_costs, fixed_costs, variable_costs, leadtimes, backorder_cost, demand_forecast, CV, order_capacities = fill(Inf, 10), inventory_capacities = fill(Inf, 10), initlevels = zeros(10))
    suppliers = []
    lts = []
    items = []
    assembs = []
    for i in 5:10
        sup = Supplier(fixed_costs[i], variable_costs[i], capacity = order_capacities[i])
        lt = LeadTime(leadtimes[i], 0., sup)
        item = Item(holding_costs[i], sSPolicy(), initlevels[i], lt, capacity = inventory_capacities[i])
        push!(suppliers, sup)
        push!(lts, lt)
        push!(items, item)
    end
    sfitems = Iterators.Stateful(items)
    for i in 2:4
        ass = Assembly(fixed_costs[i], variable_costs[i], popfirst!(sfitems) => 1, popfirst!(sfitems) => 1, capacity = order_capacities[i])
        lt = LeadTime(leadtimes[i], 0., ass)
        item = Item(holding_costs[i], sSPolicy(), initlevels[i], lt, capacity = inventory_capacities[i])
        push!(assembs, ass)
        push!(lts, lt)
        push!(items, item)
    end
    i = 1
    ass = Assembly(fixed_costs[i], variable_costs[i], popfirst!(sfitems) => 1, popfirst!(sfitems) => 1, popfirst!(sfitems) => 1, capacity = order_capacities[i])
    lt = LeadTime(leadtimes[i], 0., ass)
    item = Item(holding_costs[i], sSPolicy(), initlevels[i], lt, capacity = inventory_capacities[i])
    push!(assembs, ass)
    push!(lts, lt)
    push!(items, item)

    market = Market(backorder_cost, CVNormal{CV}, length(demand_forecast), items[1], 0., [demand_forecast; demand_forecast])
    InventorySystem(length(demand_forecast), suppliers..., assembs..., lts..., items..., market)
end

function general10(;holding_costs, fixed_costs, variable_costs, leadtimes, backorder_costs, demand_forecasts, CVs, order_capacities = fill(Inf, 10), inventory_capacities = fill(Inf, 10), initlevels = zeros(10))
    suppliers = []
    lts = []
    items = []
    assembs = []
    markets = []
    for i in 8:10
        sup = Supplier(fixed_costs[i], variable_costs[i], capacity = order_capacities[i])
        lt = LeadTime(leadtimes[i], 0., sup)
        item = Item(holding_costs[i], sSPolicy(), initlevels[i], lt, capacity = inventory_capacities[i])
        push!(suppliers, sup)
        push!(lts, lt)
        push!(items, item)
    end
    for (j, i) in enumerate(5:6)
        ass = Assembly(fixed_costs[i], variable_costs[i], items[j] => 1, items[j+1] => 1, capacity = order_capacities[i])
        lt = LeadTime(leadtimes[i], 0., ass)
        item = Item(holding_costs[i], sSPolicy(), initlevels[i], lt, capacity = inventory_capacities[i])
        push!(assembs, ass)
        push!(lts, lt)
        push!(items, item)
    end
    i = 7
    ass = Assembly(fixed_costs[i], variable_costs[i], items[3] => 1, capacity = order_capacities[i])
    lt = LeadTime(leadtimes[i], 0., ass)
    item = Item(holding_costs[i], sSPolicy(), initlevels[i], lt, capacity = inventory_capacities[i])
    push!(assembs, ass)
    push!(lts, lt)
    push!(items, item)
    i = 1
    ass = Assembly(fixed_costs[i], variable_costs[i], items[4] => 1, capacity = order_capacities[i])
    lt = LeadTime(leadtimes[i], 0., ass)
    item = Item(holding_costs[i], sSPolicy(), initlevels[i], lt, capacity = inventory_capacities[i])
    push!(assembs, ass)
    push!(lts, lt)
    push!(items, item)
    for i in 2:3
        ass = Assembly(fixed_costs[i], variable_costs[i], items[i+2] => 1, items[i+3] => 1, capacity = order_capacities[i])
        lt = LeadTime(leadtimes[i], 0., ass)
        item = Item(holding_costs[i], sSPolicy(), initlevels[i], lt, capacity = inventory_capacities[i])
        push!(assembs, ass)
        push!(lts, lt)
        push!(items, item)
    end
    i = 4
    ass = Assembly(fixed_costs[i], variable_costs[i], items[6] => 1, capacity = order_capacities[i])
    lt = LeadTime(leadtimes[i], 0., ass)
    item = Item(holding_costs[i], sSPolicy(), initlevels[i], lt, capacity = inventory_capacities[i])
    push!(assembs, ass)
    push!(lts, lt)
    push!(items, item)
    for i in 1:4
        df = demand_forecasts[i]
        push!(markets, Market(backorder_costs[i], CVNormal{CVs[i]}, length(df), items[i+6], 0., [df; df]))
    end
    InventorySystem(minimum(length.(demand_forecasts)), suppliers..., lts..., assembs..., items..., markets...)
end
=#

ext: jl


text:
#=
    Contents: Constructors for ConnectionPool instances.
=#

type ConnectionPool
    connection_prototype       # A disconnected instance of the connection
    target_pool_size::Int64    # Target number of connections in the pool
    peak_size::Int64           # Peak number of connections in the pool
    unoccupied::Set            # The set of connections in the pool that are not being used
    occupied::Set              # The set of connections in the pool that are being used
    ms_wait::Int64             # If all connections are busy, how long to wait (ms) before trying to connect again
    n_tries::Int64             # How many times to retry acquiring a connection

    """
    On creation, target_lb connections will be created
    """
    function ConnectionPool(connection, target_pool_size::Int64, peak_size::Int64, ms_wait::Int64, n_tries::Int64)
	assert(target_pool_size <= peak_size)
	assert(isfinite(peak_size))
	define_new_connection(connection)    # Includes the new_connection method for the connection type
	unoccupied = Set{typeof(connection)}()
	occupied   = Set{typeof(connection)}()
	disconnect(connection)
	new(connection, target_pool_size, peak_size, unoccupied, occupied, ms_wait, n_tries)
    end
end

function define_new_connection(connection)    # Includes the new_connection method for the connection type
    tp = string(typeof(connection))
    if tp == "Redis.RedisConnection"
	include(joinpath(dirname(@__FILE__), "supported_databases/redis_cp.jl"))
    else
	error("Unrecognized database type $tp.")
    end
end

# EOF


text:
"""
    test_statistic(ψ_1, B_11, A_21, B_21)

Computes the test statistic for the score test
"""
function test_statistic!(
    st::Union{WSVarScoreTest{T},WSVarScoreTestInvariant{T}},
    ψ_1::AbstractVector{T},
    B_11::AbstractMatrix{T},
    A_21::AbstractMatrix{T},
    B_21::AbstractMatrix{T},
    nm::WSVarLmmModel{T},
    m::Integer,
    tmp_rr::AbstractMatrix{T},
    tmp_sr::AbstractMatrix{T},
    tmp_r::AbstractVector{T}
) where {T <: BlasReal}
    # Vψ_1 = B_11 - transpose(A_21) * nm.Ainv * B_21 -
    #    transpose(B_21) * nm.Ainv * A_21 +
    #    transpose(A_21) * st.AinvBAinv * A_21
    tmp_rr .= B_11
    # println("B11: ", tmp_rr)
    # println("Ainv: ", nm.Ainv)
    mul!(tmp_sr, nm.Ainv, B_21)
    mul!(tmp_rr, transpose(A_21), tmp_sr, -one(T), one(T))
    # println("B11 - A21^T Ainv B21: ", tmp_rr)
    mul!(tmp_sr, nm.Ainv, A_21)
    mul!(tmp_rr, transpose(B_21), tmp_sr, -one(T), one(T))
    # println("B11 - A21^T Ainv B21 - B21^T Ainv A21: ", tmp_rr)
    mul!(tmp_sr, st.AinvBAinv, A_21)
    mul!(tmp_rr, transpose(A_21), tmp_sr, one(T), one(T))
    # println("Vψ_1: ", tmp_rr)

    eigfact = eigen!(Symmetric(tmp_rr))
    mul!(tmp_r, transpose(eigfact.vectors), ψ_1, one(T) / sqrt(m), zero(T))
    # println("ψ1: ", ψ_1)
    # println("ev^T ψ1: ", tmp_r)
    # v1 = one(T) / sqrt(m) * transpose(eigfact.vectors) * ψ_1
    atol = 1e-8 # tolerance for determining rank
    s = convert(Int, sign(tmp_r[1]))
    rk = 0 # rank
    ts = zero(T)
    for j in 1:length(ψ_1)
        if eigfact.values[j] > atol
            ts += abs2(tmp_r[j]) / eigfact.values[j]
            # println("ts added: ", abs2(tmp_r[j]) / eigfact.values[j])
            rk += 1
        end
    end
    if rk == 1
        return ts, s == 0 ? 1 : s
    else
        return ts, max(rk, 1)
    end
end

function test_statistic_single!(
    st::Union{WSVarScoreTest{T},WSVarScoreTestInvariant{T}},
    ψ_1::AbstractVector{T},
    B_11::AbstractMatrix{T},
    A_21::AbstractMatrix{T},
    B_21::AbstractMatrix{T},
    nm::WSVarLmmModel{T},
    m::Integer,
    tmp_rr::AbstractMatrix{T},
    tmp_sr::AbstractMatrix{T},
    tmp_r::AbstractVector{T}
) where {T <: BlasReal}
    @assert size(tmp_rr, 1) == 1
    # Vψ_1 = B_11 - transpose(A_21) * nm.Ainv * B_21 -
    #    transpose(B_21) * nm.Ainv * A_21 +
    #    transpose(A_21) * st.AinvBAinv * A_21
    tmp_rr .= B_11
    # println("B11: ", tmp_rr)
    # println("Ainv: ", nm.Ainv)
    mul!(tmp_sr, nm.Ainv, B_21)
    mul!(tmp_rr, transpose(A_21), tmp_sr, -one(T), one(T))
    # println("B11 - A21^T Ainv B21: ", tmp_rr)
    mul!(tmp_sr, nm.Ainv, A_21)
    mul!(tmp_rr, transpose(B_21), tmp_sr, -one(T), one(T))
    # println("B11 - A21^T Ainv B21 - B21^T Ainv A21: ", tmp_rr)
    mul!(tmp_sr, st.AinvBAinv, A_21)
    mul!(tmp_rr, transpose(A_21), tmp_sr, one(T), one(T))
    # println("Vψ_1: ", tmp_rr)

    ts = ψ_1[1] / (sqrt(m))# * sqrt(tmp_rr[1]))
    return ts
end

function zs!(st::Union{WSVarScoreTest{T},WSVarScoreTestInvariant{T}}
    ) where {T <: BlasReal}
    @assert st.r_X1 == 1
    @assert st.r_W1 == 1
    nm = st.nullmodel
    r_X1 = 1

    # for X1
    ψ_1p  = @view st.ψ_1[1:r_X1]
    B_11p = @view st.B_11[1:r_X1, 1:r_X1]
    A_21p = @view st.A_21[:, 1:r_X1]
    B_21p = @view st.B_21[:, 1:r_X1]
    z1 = test_statistic_single!(st, ψ_1p, B_11p, A_21p, B_21p, nm, st.m,
        st.tmp_rx1rx1, st.tmp_srx1, st.tmp_rx1)

    # for W1
    ψ_1p  = @view st.ψ_1[(r_X1 + 1):end]
    B_11p = @view st.B_11[(r_X1 + 1):end, (r_X1 + 1):end]
    A_21p = @view st.A_21[:, (r_X1 + 1):end]
    B_21p = @view st.B_21[:, (r_X1 + 1):end]
    z2 = test_statistic_single!(st, ψ_1p, B_11p, A_21p, B_21p, nm, st.m,
        st.tmp_rw1rw1, st.tmp_srw1, st.tmp_rw1)
    return z1, z2
end

"""
    pvalues!(st::WSVarScoreTest)

Returns three p-values for the score test:
for mean (`p1`), within-sample variability (`p2`), and both (`p3`).
-1 is returned for any skipped tests.
"""
function pvalues!(st::Union{WSVarScoreTest{T},WSVarScoreTestInvariant{T}}
) where {T <: BlasReal}
    nm = st.nullmodel
    r_X1, r_W1, r = st.r_X1, st.r_W1, st.r
    if st.r_X1 > 0
        ψ_1p  = @view st.ψ_1[1:r_X1]
        B_11p = @view st.B_11[1:r_X1, 1:r_X1]
        A_21p = @view st.A_21[:, 1:r_X1]
        B_21p = @view st.B_21[:, 1:r_X1]
        v1, r1 = test_statistic!(st, ψ_1p, B_11p, A_21p, B_21p, nm, st.m,
            st.tmp_rx1rx1, st.tmp_srx1, st.tmp_rx1)
        p1 = v1 ≤ 0 ? 1.0 : ccdf(Chisq(abs(r1)), v1)
    else
        p1 = -one(T)
        r1 = 0
    end

    if st.r_W1 > 0
        ψ_1p  = @view st.ψ_1[(r_X1 + 1):end]
        B_11p = @view st.B_11[(r_X1 + 1):end, (r_X1 + 1):end]
        A_21p = @view st.A_21[:, (r_X1 + 1):end]
        B_21p = @view st.B_21[:, (r_X1 + 1):end]
        v2, r2 = test_statistic!(st, ψ_1p, B_11p, A_21p, B_21p, nm, st.m,
            st.tmp_rw1rw1, st.tmp_srw1, st.tmp_rw1)
        p2 = v2 ≤ 0 ? 1.0 : ccdf(Chisq(abs(r2)), v2)
    else
        p2 = -one(T)
        r2 = 0
    end
    r2 = abs(r2) == 1 ? -r2 : r2
    p1, p2, r1, r2
end


text:
# function that handles prepares the queues and ressource flow for the
# next minute of the simulationn
function prepare_next_minute!(sim_data::DataFrame,
                              incidents::DataFrame,
                              ressource_flow::Array{Int64,3},
                              incident_queue::Array{Union{Missing,Int64},3},
                              queue_used::Array{Int64,2},
                              queue_change::Array{Bool,2},
                              mnt::Int64,
                              drop_incident::Int64)
    # if a location has unused ressources and a backlog, the ressources are assigned
    # to work on the backlog as long as there is a backlog left
    for i = 1:size(ressource_flow,2)
        if ressource_flow[mnt,i,1] > 0 && ressource_flow[mnt,i,7] > 0
            ressource_flow[mnt,i,5] += min(ressource_flow[mnt,i,1],ressource_flow[mnt,i,7])
            ressource_flow[mnt,i,1] -= min(ressource_flow[mnt,i,1],ressource_flow[mnt,i,7])
        end
        if ressource_flow[mnt,i,7] > 0
            ressource_flow[mnt,i,7] = max(ressource_flow[mnt,i,7] - ressource_flow[mnt,i,5], 0)
        else
            ressource_flow[mnt,i,1] += ressource_flow[mnt,i,5]
            ressource_flow[mnt,i,5]  = 0
        end
        if ressource_flow[mnt,i,6] / sum(ressource_flow[mnt,i,1:6]) < patrol_ratio && ressource_flow[mnt,i,1] > 1
            ressource_flow[mnt,i,6] += 1
            ressource_flow[mnt,i,1] -= 1
        elseif (ressource_flow[mnt,i,6] + 1) / sum(ressource_flow[mnt,i,1:6]) > patrol_ratio && ressource_flow[mnt,i,6] > 0
            ressource_flow[mnt,i,6] -= 1
            ressource_flow[mnt,i,1] += 1
        end
    end
    # drop all fulfilled incidents and the incidents that pass the "drop_incident" threshold from
    # the incident queues
    for i = 1:size(incident_queue,2)
        for p = 1:size(incident_queue,3)
            for j = 1:size(incident_queue,1)
                if ismissing(incident_queue[j,i,p]) == false
                    current_case = incident_queue[j,i,p]
                    if  sim_data[current_case,:cars_missing] == 0 ||
                        (ismissing(sim_data[current_case,:dispatch_minute_first]) == false &&
                        mnt - sim_data[current_case,:dispatch_minute_first] >= incidents[current_case,:length]) ||
                        mnt - sim_data[current_case,:incident_minute] >= drop_incident
                            incident_queue[j,i,p] = missing
                            queue_change[i,p] = true
                            queue_used[i,p] -= 1
                    end
                end
            end
        end
    end
    # move the current ressource status to the next minute
    if mnt + 1 <= size(ressource_flow,1)
        for i = 1:size(ressource_flow,2)
            for j = 1:size(ressource_flow,3)
                ressource_flow[mnt+1,i,j] += ressource_flow[mnt,i,j]
            end
        end
    end
end


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

# Real
function realeigtest(n, iter)
    A = rand(n,n)
    d = Vector{eltype(A)}
    v = similar(A)
    for i = 1:iter
        d,v = eig(A)
    end
    d,v
end

# Symmetric
function symeigtest(n, iter)
    A = rand(n,n)
    A = A + adjoint(A)
    d = Vector{eltype(A)}
    v = similar(A)
    for i = 1:iter
        d,v = eig(A)
    end
    d,v
end

# Hermitian
function hermitianeigtest(n, iter)
    A = rand(n,n) + im*rand(n,n)
    A = A + adjoint(A)
    d = Vector{eltype(A)}
    v = similar(A)
    for i = 1:iter
        d,v = eig(A)
    end
    d,v
end

problemsizes = [(2, 10^4, "tiny"), (2^4, 10^3, "small"), (2^6, 10^2, "medium"), (2^8, 5, "large"), (2^10, 1, "huge")]
testdata = [(realeigtest, "realeig", "Real matrix eigenfactorization", problemsizes),
            (symeigtest, "symeig", "Symmetric matrix eigenfactorization", problemsizes),
            (hermitianeigtest, "hermitianeig", "Hermitian matrix eigenfactorization", problemsizes)]
include("../perfgeneric.jl")

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
# complex dense, projected w/ transient
# Compare to figures in S&P, chap. 22.
# ∃ small discrepancies for large ϵ,
# but we are consistent w/ EigTool.
using Pseudospectra, Test

@testset "Projection/Transients" begin
    # for a sensible projection, must provide axes
    opts = Dict{Symbol,Any}(:ax => [-1,0.3,-1.2,0.2],:npts=>100,:proj_lev=>1.5)
    A = Pseudospectra.orrsommerfeld(200,10000,1.02)

    ps_data = new_matrix(A,opts)
    driver!(ps_data,opts,gs)
    @test iscomputed(ps_data)

    mtxexpsplot(gs,ps_data,1.0,100,lbmethod=:old)

    # additional plots to compare w/ book Fig. 22.6:

    # modeplot(ps_data,gs,0,0.0-0.2im)
    # modeplot(ps_data,gs,1,-0.1-0.8im)
end


text:
"""
    ForeverVector{T}

Vector of values of type `T` that reuses the value at the last index for
subsequent indices.
"""
struct ForeverVector{T}
    "Elements of the vector."
    xs::Vector{T}
end

"""
    ForeverVector(x)

Create a `ForeverVector` with value(s) `x`.
"""
ForeverVector(x) = ForeverVector([x])

Base.getindex(fv::ForeverVector, i::Int) =
    i <= length(fv.xs) ? fv.xs[i] : fv.xs[end]

"""
    SweepSettings

Parameters for a single DMRG sweep.
"""
struct SweepSettings
    "Number of diagonalization iterations in a sweep."
    num_iters::Int

    "Minimum bond dimension to keep when truncating."
    bond_min::Int
    "Maximum bond dimension to keep when truncating."
    bond_max::Int
    "Maximum eigenvalue weight to drop when truncating."
    cutoff_max::Float64
end

"""
    SweepSchedule

Schedule of parameters for multiple DMRG sweeps.
"""
struct SweepSchedule
    "Maximum number of sweeps."
    max_sweeps::Int
    "Wavefunction tolerance."
    tolerance::Float64

    # SweepSettings parameters.
    num_iters::ForeverVector{Int}
    bond_min::ForeverVector{Int}
    bond_max::ForeverVector{Int}
    cutoff_max::ForeverVector{Float64}
end

"""
    SweepSchedule(max_sweeps::Int; tolerance=nothing, kwargs...)

Create a `SweepSchedule` for `max_sweeps` sweeps, with optional parameters
specified in `kwargs`.
"""
function SweepSchedule(max_sweeps::Int; tolerance=nothing, kwargs...)
    max_sweeps >= 1 || throw(DomainError(max_sweeps, "At least 1 sweep."))

    if isnothing(tolerance)
        # Default.
        tolerance = 1e-5
    elseif tolerance < 0.0
        throw(DomainError(tolerance, "Non-negative tolerance."))
    end

    kwargs = Dict(kwargs)

    num_iters = ForeverVector(get(kwargs, :num_iters, 2))
    bond_min = ForeverVector(get(kwargs, :bond_min, 5))
    bond_max = ForeverVector(get(kwargs, :bond_max, 100))
    cutoff_max = ForeverVector(get(kwargs, :cutoff_max, 1e-10))

    SweepSchedule(max_sweeps, tolerance,
                  num_iters, bond_min, bond_max, cutoff_max)
end

Base.getindex(sch::SweepSchedule, i::Int) =
    SweepSettings(sch.num_iters[i],
                  sch.bond_min[i],
                  sch.bond_max[i],
                  sch.cutoff_max[i])


text:
glbvsc    = UInt(1) #dataObject.defaultSpaceCnt
glbvTcNum = UInt(2) #dataObject.defaultTimeCnt

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
module TestUtils
    using Test
    using StaticArrays

    using iLQGames:
        GeneralGame,
        QuadCache,
        n_states,
        n_controls,
        quadraticize!,
        _quadraticize_ad,
        player_costs

    export
        quad_sanity_check

    # TODO: move to general game test utils
    function quad_sanity_check(g::GeneralGame)
        nx = n_states(g)
        nu = n_controls(g)
        qcache = zero(QuadCache{n_states(g), n_controls(g)})
        @testset "Quadratization Sanity Check" begin
            for i in 1:100
                x = SVector{nx, Float64}(randn(nx))
                u = SVector{nu, Float64}(randn(nu))
                t = 10.0

                for pc in player_costs(g)
                    qc_manual = quadraticize!(qcache, pc, g, x, u, t)
                    qc_ad = _quadraticize_ad(pc, g, x, u, t)

                    @test isapprox(qc_manual.Q, qc_ad.Q)
                    @test isapprox(qc_manual.l, qc_ad.l)
                    @test isapprox(qc_manual.R, qc_ad.R)
                end
            end
        end;
    end
end


text:
### GAMC state subtypes

abstract type GAMCState{F<:VariateForm} <: LMCSamplerState{F} end

## MuvGAMCState holds the internal state ("local variables") of the GAMC sampler for multivariate parameters

mutable struct MuvGAMCState <: GAMCState{Multivariate}
  proposal::Union{MultivariateGMM, AbstractMvNormal, Void}
  pstate::ParameterState{Continuous, Multivariate}
  tune::MCTunerState
  sqrttunestep::Real
  ratio::Real
  μ::RealVector
  lastmean::RealVector
  secondlastmean::RealVector
  newinvtensor::RealMatrix
  oldinvtensor::RealMatrix
  cholinvtensor::RealLowerTriangular
  newfirstterm::RealVector
  oldfirstterm::RealVector
  sqrtminorscale::Real
  w::RealVector
  presentupdatetensor::Bool
  pastupdatetensor::Bool
  count::Integer
  updatetensorcount::Integer
  diagnosticindices::Dict{Symbol, Integer}

  function MuvGAMCState(
    proposal::Union{MultivariateGMM, AbstractMvNormal, Void},
    pstate::ParameterState{Continuous, Multivariate},
    tune::MCTunerState,
    sqrttunestep::Real,
    ratio::Real,
    μ::RealVector,
    lastmean::RealVector,
    secondlastmean::RealVector,
    newinvtensor::RealMatrix,
    oldinvtensor::RealMatrix,
    cholinvtensor::RealLowerTriangular,
    newfirstterm::RealVector,
    oldfirstterm::RealVector,
    sqrtminorscale::Real,
    w::RealVector,
    presentupdatetensor::Bool,
    pastupdatetensor::Bool,
    count::Integer,
    updatetensorcount::Integer,
    diagnosticindices::Dict{Symbol, Integer}
  )
    if !isnan(ratio)
      @assert 0 < ratio < 1 "Acceptance ratio should be between 0 and 1"
      @assert sqrttunestep > 0 "Square root of tuned drift step is not positive"
    end

    if !isnan(sqrtminorscale)
      @assert sqrtminorscale >= 0 "Scaling of stabilizing covariance should be non-negative"
    end

    if !isnan(w[1])
      @assert w[1] > 0 "Weight of core mixture component must be positive"
    end
    if !isnan(w[2])
      @assert w[2] >= 0 "Weight of minor mixture component must be non-negative"
    end

    @assert count >= 0 "Number of iterations (count) should be non-negative"

    new(
      proposal,
      pstate,
      tune,
      sqrttunestep,
      ratio,
      μ,
      lastmean,
      secondlastmean,
      newinvtensor,
      oldinvtensor,
      cholinvtensor,
      newfirstterm,
      oldfirstterm,
      sqrtminorscale,
      w,
      presentupdatetensor,
      pastupdatetensor,
      count,
      updatetensorcount,
      diagnosticindices
    )
  end
end

MuvGAMCState(
  pstate::ParameterState{Continuous, Multivariate},
  sqrtminorscale::Real,
  w::RealVector,
  tune::MCTunerState=GAMCTune(),
  lastmean::RealVector=Array{eltype(pstate)}(pstate.size)
) =
  MuvGAMCState(
  nothing,
  pstate,
  tune,
  NaN,
  NaN,
  Array{eltype(pstate)}(pstate.size),
  Array{eltype(pstate)}(pstate.size),
  Array{eltype(pstate)}(pstate.size),
  Array{eltype(pstate)}(pstate.size, pstate.size),
  Array{eltype(pstate)}(pstate.size, pstate.size),
  RealLowerTriangular(Array{eltype(pstate)}(pstate.size, pstate.size)),
  Array{eltype(pstate)}(pstate.size),
  Array{eltype(pstate)}(pstate.size),
  sqrtminorscale,
  w,
  true,
  false,
  0,
  0,
  Dict{Symbol, Integer}()
)

### Geometric adaptive Monte Carlo (GAMC)

struct GAMC <: LMCSampler
  update!::Function
  transform::Union{Function, Void}
  driftstep::Real
  minorscale::Real # Scaling factor of covariance matrix of the stabilizing mixture component
  c::Real # Non-negative constant with relative small value that determines the mixture weight of the stabilizing component
  t0::Integer

  function GAMC(update!::Function, transform::Union{Function, Void}, driftstep::Real, minorscale::Real, c::Real, t0::Integer)
    @assert driftstep > 0 "Drift step is not positive"
    @assert 0 < minorscale "Constant minorscale must be positive, got $minorscale"
    @assert 0 <= c "Constant c must be non-negative"
    @assert t0 > 0 "t0 is not positive"
    new(update!, transform, driftstep, minorscale, c, t0)
  end
end

GAMC(;
  update::Function=rand_update!,
  transform::Union{Function, Void}=nothing,
  driftstep::Real=1.,
  minorscale::Real=1.,
  c::Real=0.05,
  t0::Integer=3
) =
  GAMC(update, transform, driftstep, minorscale, c, t0)

### Initialize GAMC sampler

## Initialize parameter state

function initialize!(
  pstate::ParameterState{Continuous, Multivariate},
  parameter::Parameter{Continuous, Multivariate},
  sampler::GAMC,
  outopts::Dict
)
  parameter.uptotensorlogtarget!(pstate)
  if sampler.transform != nothing
    pstate.tensorlogtarget = sampler.transform(pstate.tensorlogtarget)
  end
  pstate.diagnosticvalues[1] = false
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"
  @assert all(isfinite.(pstate.gradlogtarget)) "Gradient of log-target not finite: initial values out of support"
  @assert all(isfinite.(pstate.tensorlogtarget)) "Tensor of log-target not finite: initial values out of support"

  if !isempty(outopts[:diagnostics])
    pstate.diagnostickeys = copy(outopts[:diagnostics])
    pstate.diagnosticvalues = Array{Any}(length(pstate.diagnostickeys))
  end
end

## Initialize GAMC state

set_gmm(
  sampler::GAMC,
  pstate::ParameterState{Continuous, Multivariate},
  C::RealMatrix,
  corescale::Real,
  sqrtminorscale::Real,
  w::RealVector
) =
  MixtureModel([MvNormal(pstate.value, corescale*C), MvNormal(pstate.value, sqrtminorscale)], w)

set_gmm!(sstate::MuvGAMCState, sampler::GAMC, pstate::ParameterState{Continuous, Multivariate}) =
  sstate.proposal = set_gmm(
    sampler, pstate, sstate.oldinvtensor, sstate.tune.totaltune.step, sstate.sqrtminorscale, sstate.w
  )

tuner_state(sampler::GAMC, tuner::GAMCTuner) =
  GAMCTune(
    BasicMCTune(NaN, 0, 0, tuner.smmalatuner.period),
    BasicMCTune(NaN, 0, 0, tuner.amtuner.period),
    BasicMCTune(sampler.driftstep, 0, 0, tuner.totaltuner.period)
  )

function sampler_state(
  parameter::Parameter{Continuous, Multivariate},
  sampler::GAMC,
  tuner::MCTuner,
  pstate::ParameterState{Continuous, Multivariate},
  vstate::VariableStateVector,
  diagnostickeys::Vector{Symbol}
)
  oldinvtensor = inv(pstate.tensorlogtarget)
  w = [1-sampler.c, sampler.c]

  sstate = MuvGAMCState(
    generate_empty(pstate, parameter.diffmethods, parameter.diffopts),
    sqrt(sampler.minorscale),
    w,
    tuner_state(sampler, tuner),
    copy(pstate.value)
  )

  sstate.sqrttunestep = sqrt(sampler.driftstep)
  sstate.oldinvtensor[:, :] = inv(pstate.tensorlogtarget)
  sstate.cholinvtensor[:, :] = ctranspose(chol(Hermitian(sstate.oldinvtensor)))
  sstate.oldfirstterm[:] = sstate.oldinvtensor*pstate.gradlogtarget
  set_diagnosticindices!(sstate, [:accept], diagnostickeys)

  sstate
end

### Reset GAMC sampler

## Reset parameter state

function reset!(
  pstate::ParameterState{Continuous, Multivariate},
  x::RealVector,
  parameter::Parameter{Continuous, Multivariate},
  sampler::GAMC
)
  pstate.value = copy(x)
  parameter.uptotensorlogtarget!(pstate)
  pstate.diagnosticvalues[1] = false
end

## Reset sampler state

function reset!(
  sstate::MuvGAMCState,
  pstate::ParameterState{Continuous, Multivariate},
  parameter::Parameter{Continuous, Multivariate},
  sampler::GAMC,
  tuner::MCTuner
)
  reset!(sstate.tune, sampler, tuner)
  sstate.sqrttunestep = sqrt(sampler.driftstep)
  sstate.lastmean = copy(pstate.value)
  sstate.presentupdatetensor = true
  sstate.pastupdatetensor = false
  sstate.count = 0
end

### GAMC schedules

am_only_update!(sstate::MuvGAMCState, pstate::ParameterState{Continuous, Multivariate}, i::Integer, tot::Integer) =
  sstate.presentupdatetensor = false

cos_update!(
  sstate::MuvGAMCState,
  pstate::ParameterState{Continuous, Multivariate},
  i::Integer,
  tot::Integer,
  a::Real=10.,
  b::Real=0.2,
  c::Real=0.2
) =
  sstate.presentupdatetensor = rand(Bernoulli(b*cos(a*i*pi/tot)+c))

mod_update!(
  sstate::MuvGAMCState,
  pstate::ParameterState{Continuous, Multivariate},
  i::Integer,
  tot::Integer,
  n::Integer=10
) =
  sstate.presentupdatetensor = mod(i, n) == 0 ? true : false

rand_exp_decay_update!(
  sstate::MuvGAMCState,
  pstate::ParameterState{Continuous, Multivariate},
  i::Integer,
  tot::Integer,
  a::Real=10.,
  b::Real=0.
) =
  sstate.presentupdatetensor = rand(Bernoulli(exp_decay(i, tot, a, b)))

rand_lin_decay_update!(
  sstate::MuvGAMCState,
  pstate::ParameterState{Continuous, Multivariate},
  i::Integer,
  tot::Integer,
  a::Real=1e-2,
  b::Real=0.
) =
  sstate.presentupdatetensor = rand(Bernoulli(lin_decay(i, tot, a, b)))

rand_pow_decay_update!(
  sstate::MuvGAMCState,
  pstate::ParameterState{Continuous, Multivariate},
  i::Integer,
  tot::Integer,
  a::Real=1e-3,
  b::Real=0.
) =
  sstate.presentupdatetensor = rand(Bernoulli(pow_decay(i, tot, a, b)))

rand_quad_decay_update!(
  sstate::MuvGAMCState,
  pstate::ParameterState{Continuous, Multivariate},
  i::Integer,
  tot::Integer,
  a::Real=1e-5,
  b::Real=0.
) =
  sstate.presentupdatetensor = rand(Bernoulli(quad_decay(i, tot, a, b)))

rand_update!(
  sstate::MuvGAMCState,
  pstate::ParameterState{Continuous, Multivariate},
  i::Integer,
  tot::Integer,
  p::Real=0.5
) =
  sstate.presentupdatetensor = rand(Bernoulli(p))

smmala_only_update!(sstate::MuvGAMCState, pstate::ParameterState{Continuous, Multivariate}, i::Integer, tot::Integer) =
  sstate.presentupdatetensor = true

### GAMC show methods

show(io::IO, sampler::GAMC) = print(
  io,
  "GAMC sampler: drift step = ",
  sampler.driftstep,
  ", scaling of stabilizing covariance = ",
  sampler.minorscale,
  ", weight of stabilizing mixture component = ",
  sampler.c
)


text:
@export_const Float begin
    metre = 1e2
    meter = metre
    cm = metre / 1e2
    mm = metre / 1e3
    um = metre / 1e6
    nm = metre / 1e9

    second = 1e3
    ms = second / 1e3
    Hz = 1 / second

    voltage = 1e3
    mV = voltage / 1e3

    ampere = 1e6
    mA = ampere / 1e3
    uA = ampere / 1e6
    nA = ampere / 1e9

    farad=1e6
    uF = farad/1e6
    ufarad = uF

    siemens = 1e3
    mS = siemens / 1e3
    msiemens = mS
    nS = siemens / 1e9
    nsiemens = nS
end


text:
################################################################################
#
#       DedekindCriterion.jl : Dedekinds Criterion for maximality
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

function dedekind_test(O::NfOrd, p::fmpz, compute_order::Type{Val{S}} = Val{true}) where S
  !isequation_order(O) && error("Order must be an equation order")

  if rem(discriminant(O), p^2) != 0
    if compute_order == Val{true}
      return true, O
    else
      return true
    end
  end

  Zy, y = PolynomialRing(FlintZZ, "y")
  Kx, x = PolynomialRing(GF(p, cached=false), "x", cached=false)

  f = nf(O).pol

  Zyf = Zy(f)

  fmodp = Kx(Zyf)

  fac = factor_squarefree(fmodp)

  g = prod(x for x in keys(fac.fac))
  h = divexact(fmodp,g)

  # first build 1/p ( f - g*h)
  gZ = lift(Zy,g)
  hZ = lift(Zy, h)

  g1 = Zyf - gZ*hZ

  for i in 0:degree(g1)
    q,r = divrem(coeff(g1,i),p)
    @assert r == 0
    setcoeff!(g1,i,q)
  end

  g1modp = Kx(g1)
  U = gcd(gcd(g, h), g1modp)

  if compute_order == Val{false}
    if isone(U)
      return true
    else
      return false
    end
  else

    pmaximal = (isone(U))

    if pmaximal
      return true, O
    end

    @hassert :NfOrd 1 rem(fmodp, U) == zero(Kx)
    U = divexact(fmodp, U)

    @hassert :NfOrd 1 rem(O.disc, p^2) == 0

    alpha = nf(O)(parent(f)(lift(Zy,U)))

    # build the new basis matrix
    # we have to take the representation matrix of alpha!
    # concatenating the coefficient vector won't help
    Malpha, d = representation_matrix_q(alpha)
    @assert isone(d)
    n = _hnf_modular_eldiv(Malpha, p, :lowerleft)
    b = FakeFmpqMat(n, p)
    @hassert :NfOrd 1 defines_order(nf(O), b)[1]
    OO = Order(nf(O), b, check = false)

    OO.isequation_order = false

    OO.disc = divexact(O.disc, p^(2*(degree(O)-degree(U))))

    return false, OO
  end
end

dedekind_test(O::NfOrd, p::Integer) = dedekind_test(O, FlintZZ(p))

dedekind_ispmaximal(O::NfOrd, p::fmpz) = dedekind_test(O, p, Val{false})

dedekind_ispmaximal(O::NfOrd, p::Integer) = dedekind_ispmaximal(O, FlintZZ(p))

dedekind_poverorder(O::NfOrd, p::fmpz) = dedekind_test(O, p)[2]

dedekind_poverorder(O::NfOrd, p::Integer) = dedekind_poverorder(O, FlintZZ(p))

###############################################################################
#
#  Non-prime case
#
###############################################################################

function dedekind_test_composite(O::NfOrd, p::fmpz)
  @assert isequation_order(O)

  Zy = PolynomialRing(FlintZZ, "y")[1]
  R = ResidueRing(FlintZZ, p, cached = false)
  Rx = PolynomialRing(R, "x", cached=false)[1]

  f = Zy(nf(O).pol)

  fmodp = Rx(f)

  # Now, I would like to have the squarefree factorization of fmodp
  # I go for the f/gcd(f,f')

  divs, gcdfderf = _gcd_with_failure(fmodp, derivative(fmodp))

  if !isone(divs)
    return gcd(divs, p), O
  end

  sqff = divexact(fmodp, gcdfderf)

  # first build 1/p ( f - g*h)
  gZ = lift(Zy,sqff)
  hZ = lift(Zy,gcdfderf)

  g1 = f - gZ*hZ
  g1 = divexact(g1, p)

  g1modp = Rx(g1)

  divs, par1 = _gcd_with_failure(gcdfderf, sqff)
  if !isone(divs)
    return gcd(divs, p), O
  end
  divs, U = _gcd_with_failure(par1, g1modp)
  if !isone(divs)
    return gcd(divs, p), O
  end

  if isone(U)
    return fmpz(1), O
  end

  U = divexact(fmodp, U)
  alpha = nf(O)(parent(f)(lift(Zy,U)))

  Malpha, d = representation_matrix_q(alpha)
  hnf_modular_eldiv!(Malpha, p, :lowerleft)
  b = FakeFmpqMat(Malpha, p)

  @hassert :NfOrd 1 defines_order(nf(O), b)[1]
  OO = NfOrd(nf(O), b)
  temp = divexact(b.den^degree(O), prod_diagonal(b.num))
  fl, qq = divides(discriminant(O), temp^2)
  @assert fl
  OO.disc = qq
  @hassert :NfOrd 1 discriminant(basis(OO)) == OO.disc
  OO.index = temp

  return fmpz(1), OO
end


text:
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.
using ISA, RecipesBase

# Load Colormap
include(raw"colorMaps/cubeYF.jl")
include(raw"colorMaps/viridis.jl")

"""
    plot(z::AMFMmodel; timeaxis = 0.0:0.005:1.0)
    plot(ψ::AMFMcomp;  timeaxis = 0.0:0.005:1.0)

Create a 3D Argand Digram ...

# Examples
```@example
using ISA, Plots
ψ₀ = AMFMcomp(t->10t,t->25cos(t),0.0)
plot(ψ₀)
```
"""

# 3D Argand Digram
function viewAngle3DArgand(view)
   normalized = uppercase(view)
   if normalized == "DEFAULT"
      return "time(s)","imag","real",(45,45),true,(true,true,true)
   elseif normalized == "TI" || normalized == "IT"
      return "time(s)","imag","",(0,90),true,(true,true,false)
   elseif normalized == "TR" || normalized == "RT"
      return "time(s)","","real",(0,0),false,(true,false,true)
   elseif normalized == "IR" || normalized == "RI"
      return "","Imag","real",(90,0),false,(false,true,true)
   else return error("Invalid Argand Graph. Choose from \"TR\", \"TI\", or \"IR\"")
   end
end

function colorSelect(CM)
   if isa(CM,Vector{RGB{Float64}})
      if length(CM)==256
         return CM
      else
         return error("Invalid color map. Colormap must be of lengh 256.")
      end
   else
      normalized = uppercase(CM)
      if normalized == "CUBEYF"
         return cubeYF()
      elseif normalized == "VIRIDIS"
         return viridis()
      else return error("Invalid color map. Choose from \"cubeYF\" or \"viridis\"")
      end
   end
end

@recipe function temp(ψ::AMFMcomp; timeaxis = 0.0:0.005:1.0,
                      FreqUnits = "rad/s", view="default",
                      colorMap="VIRIDIS")
   xguide --> viewAngle3DArgand(view)[1]
   yguide --> viewAngle3DArgand(view)[2]
   zguide --> viewAngle3DArgand(view)[3]
   camera --> viewAngle3DArgand(view)[4]
   ymirror --> viewAngle3DArgand(view)[5]
   xticks --> viewAngle3DArgand(view)[6][1]
   yticks --> viewAngle3DArgand(view)[6][2]
   zticks --> viewAngle3DArgand(view)[6][3]
   background_color --> colorSelect(colorMap)[1]
   foreground_color --> :white
   legend --> false
   framestyle --> :origin
   t = timeaxis
   linealpha --> sign.(abs.( ψ(timeaxis) ))
   a_max = maximum(abs.(ψ.C.a.(t)))
   clim = (0,1)
   seriescolor := colorSelect(colorMap)[ max.(min.(round.(Int, abs.(ψ.(t)) .* 256/a_max ),256),50) ]
   t,imag(ψ.(t)),real(ψ.(t))
end

@recipe function temp(z::AMFMmodel; timeaxis = 0.0:0.005:1.0,
                      FreqUnits = "rad/s", view="default",
                      colorMap="VIRIDIS")
   xguide --> viewAngle3DArgand(view)[1]
   yguide --> viewAngle3DArgand(view)[2]
   zguide --> viewAngle3DArgand(view)[3]
   camera --> viewAngle3DArgand(view)[4]
   ymirror --> viewAngle3DArgand(view)[5]
   xticks --> viewAngle3DArgand(view)[6][1]
   yticks --> viewAngle3DArgand(view)[6][2]
   zticks --> viewAngle3DArgand(view)[6][3]
   background_color --> colorSelect(colorMap)[1]
   foreground_color --> :white
   legend --> false
   framestyle --> :origin
   t = timeaxis
   linealpha --> sign.(abs.(z(timeaxis)))
   a_max = maximum(abs.(z(t)))
   clim = (0,1)
   seriescolor := colorSelect(colorMap)[ max.(min.(round.(Int, abs.(z(t)) .* 256/a_max ),256),50) ]
   timeaxis, imag(z(t)), real(z(t))
end

"""
    plot(S::compSet;     timeaxis = 0.0:0.005:1.0, FreqUnits = "rad/s")
    plot(𝐶::AMFMtriplet; timeaxis = 0.0:0.005:1.0, FreqUnits = "rad/s")

Create a 3D Instantaneous Spectrum...

# Examples
```@example
using ISA, Plots

```
"""

# 3D IS Plot
function getFnorm(FreqUnits)
   if FreqUnits == "rad/s"
      Fnorm = 1
   elseif FreqUnits == "Hz"
      Fnorm = 1/2π
   else
      error("invalid FreqUnits")
   end
end
function viewAngleIS(view,FreqUnits)
   normalized = uppercase(view)
   if normalized == "DEFAULT"
      return "time(s)","freq("*FreqUnits*")","real",(20,80),true,0,(true,true,true)
   elseif normalized == "TF" || normalized == "FT"
      return "time(s)","freq("*FreqUnits*")"," ",(0,90),true,-90,(true,true,false)
   elseif normalized == "TR" || normalized == "RT"
      return "time(s)"," ","real",(0,0),false,0,(true,false,true)
   elseif normalized == "FR" || normalized == "RF"
      return "","freq("*FreqUnits*")","real",(90,0),false,0,(false,true,true)
   else return error("Invalid IS Graph. Choose from \"TR\", \"TF\", or \"FR\"")
   end
end

@recipe function temp(𝐶::AMFMtriplet; timeaxis = 0.0:0.005:1.0,
                      FreqUnits = "rad/s",view="default",
                      colorMap="VIRIDIS")
   xguide --> viewAngleIS(view,FreqUnits)[1]
   yguide --> viewAngleIS(view,FreqUnits)[2]
   zguide --> viewAngleIS(view,FreqUnits)[3]
   camera --> viewAngleIS(view,FreqUnits)[4]
   ymirror --> viewAngleIS(view,FreqUnits)[5]
   yrotation --> viewAngleIS(view,FreqUnits)[6]
   xticks --> viewAngleIS(view,FreqUnits)[7][1]
   yticks --> viewAngleIS(view,FreqUnits)[7][2]
   zticks --> viewAngleIS(view,FreqUnits)[7][3]
   background_color --> colorSelect(colorMap)[1]
   foreground_color --> :white
   legend --> false
   framestyle --> :origin
   Fnorm = getFnorm(FreqUnits)
   t = timeaxis
   linealpha --> sign.(abs.( 𝐶.a.(timeaxis) ))
   a_max = maximum(abs.(𝐶.a.(t)))
   clim = (0,1)
   seriescolor := colorSelect(colorMap)[ max.(min.(round.(Int, abs.(AMFMcomp(𝐶).(t)) .* 256/a_max ),256),50) ]
   t,𝐶.ω.(t),real.(AMFMcomp(𝐶).(t))
end

@recipe function temp(S::compSet; timeaxis = 0.0:0.005:1.0,
                      FreqUnits = "rad/s",view="default",
                      colorMap="VIRIDIS",realProj=false)
   xguide --> viewAngleIS(view,FreqUnits)[1]
   yguide --> viewAngleIS(view,FreqUnits)[2]
   zguide --> viewAngleIS(view,FreqUnits)[3]
   camera --> viewAngleIS(view,FreqUnits)[4]
   ymirror --> viewAngleIS(view,FreqUnits)[5]
   yrotation --> viewAngleIS(view,FreqUnits)[6]
   xticks --> viewAngleIS(view,FreqUnits)[7][1]
   yticks --> viewAngleIS(view,FreqUnits)[7][2]
   zticks --> viewAngleIS(view,FreqUnits)[7][3]
   background_color --> colorSelect(colorMap)[1]
   foreground_color --> :white
   legend --> false
   framestyle --> :origin
   Fnorm = getFnorm(FreqUnits)
   t = timeaxis
   a_max = maximum([maximum(abs.(S.S[k].a.(t))) for k in 1:length(S.S)])
   clim = (0,1)
   realProj ? projection = zeros(length(t)) : nothing
   realProj ? minVector = zeros(length(S.S)) : nothing
   for k in 1:length(S.S)
      seriescolor := colorSelect(colorMap)[ max.(min.(round.(Int, abs.(S.S[k].a.(t)) .* 256/a_max ),256),50) ]
      #linealpha --> max.(min.( abs.(z.S.S[1].a.(t)).^(1/2) .* 1/a_max ,1),0)

      @series begin
         linealpha --> sign.(abs.( S.S[k].a.(timeaxis) ))
         timeaxis, Fnorm.*S.S[k].ω.(t), real(AMFMcomp(S.S[k]).(t))
      end
      realProj ? minVector[k] = minimum(Fnorm.*S.S[k].ω.(t)) : nothing
      realProj ? projection += real(AMFMcomp(S.S[k]).(t)) : nothing
   end
   if realProj
      append!(minVector,0)
      @series begin
         seriescolor := :white
         timeaxis, zeros(length(t)).+minimum(minVector), projection
      end
   end
end

# Models
# 3D Argand Digram
@recipe function temp(𝚿::numComp;FreqUnits = "rad/s",
                      view="default", colorMap = "VIRIDIS")
   xguide --> viewAngle3DArgand(view)[1]
   yguide --> viewAngle3DArgand(view)[2]
   zguide --> viewAngle3DArgand(view)[3]
   camera --> viewAngle3DArgand(view)[4]
   ymirror --> viewAngle3DArgand(view)[5]
   xticks --> viewAngle3DArgand(view)[6][1]
   yticks --> viewAngle3DArgand(view)[6][2]
   zticks --> viewAngle3DArgand(view)[6][3]
   background_color --> colorSelect(colorMap)[1]
   linealpha --> sign.(abs.(𝚿(𝚿.t)))
   foreground_color --> :white
   legend --> false
   framestyle --> :origin
   a_max = maximum(abs.(𝚿.t))
   clim = (0,1)
   seriescolor := colorSelect(colorMap)[ max.(min.(round.(Int, abs.(reverse(𝚿.t)) .* 256/a_max ),256),50) ]
   𝚿.t,imag(𝚿.Ψ),real(𝚿.Ψ)
end

@recipe function temp(𝐳::numModel; FreqUnits = "rad/s",
                      view="default", colorMap = "VIRIDIS")
   xguide --> viewAngle3DArgand(view)[1]
   yguide --> viewAngle3DArgand(view)[2]
   zguide --> viewAngle3DArgand(view)[3]
   camera --> viewAngle3DArgand(view)[4]
   ymirror --> viewAngle3DArgand(view)[5]
   xticks --> viewAngle3DArgand(view)[6][1]
   yticks --> viewAngle3DArgand(view)[6][2]
   zticks --> viewAngle3DArgand(view)[6][3]
   background_color --> colorSelect(colorMap)[1]
   tt = 𝐳.𝚿ₖ[1].t
   linealpha --> sign.(abs.( 𝐳(tt)))
   foreground_color --> :white
   legend --> false
   framestyle --> :origin
   a_max = maximum(abs.(𝐳.(tt)))
   clim = (0,1)
   seriescolor := colorSelect(colorMap)[ max.(min.(round.(Int, abs.(𝐳(tt)) .* 256/a_max ),256),50) ]
   tt, imag(𝐳(tt)), real(𝐳(tt))
end

# 3D IS Plot
@recipe function temp(𝐂::numTriplet;FreqUnits = "rad/s",
                      view="default",colorMap = "VIRIDIS")
   xguide --> viewAngleIS(view,FreqUnits)[1]
   yguide --> viewAngleIS(view,FreqUnits)[2]
   zguide --> viewAngleIS(view,FreqUnits)[3]
   camera --> viewAngleIS(view,FreqUnits)[4]
   ymirror --> viewAngleIS(view,FreqUnits)[5]
   yrotation --> viewAngleIS(view,FreqUnits)[6]
   xticks --> viewAngleIS(view,FreqUnits)[7][1]
   yticks --> viewAngleIS(view,FreqUnits)[7][2]
   zticks --> viewAngleIS(view,FreqUnits)[7][3]
   background_color --> colorSelect(colorMap)[1]
   foreground_color --> :white
   legend --> false
   framestyle --> :origin
   Fnorm = getFnorm(FreqUnits)
   a_max = maximum(abs.(𝐂.a))
   clim = (0,1)
   rm = (!).(isnan.(𝐂.ω))# remove "NaN" values in 𝐂.ω
   seriescolor := colorSelect(colorMap)[ max.(min.(round.(Int, abs.((𝐂.a)) .* 256/a_max ),256),50) ][rm]
   𝐂.t[rm],𝐂.ω[rm],real.(𝐂.Ψ)[rm]
end

@recipe function temp(𝐒::numSet;FreqUnits = "rad/s",
                      view="default",colorMap = "VIRIDIS",realProj=false)
   xguide --> viewAngleIS(view,FreqUnits)[1]
   yguide --> viewAngleIS(view,FreqUnits)[2]
   zguide --> viewAngleIS(view,FreqUnits)[3]
   camera --> viewAngleIS(view,FreqUnits)[4]
   ymirror --> viewAngleIS(view,FreqUnits)[5]
   yrotation --> viewAngleIS(view,FreqUnits)[6]
   xticks --> viewAngleIS(view,FreqUnits)[7][1]
   yticks --> viewAngleIS(view,FreqUnits)[7][2]
   zticks --> viewAngleIS(view,FreqUnits)[7][3]
   background_color --> colorSelect(colorMap)[1]
   foreground_color --> :white
   legend --> false
   framestyle --> :origin
   Fnorm = getFnorm(FreqUnits)
   a_max = maximum([maximum(abs.(𝐒.S[k].a)) for k in 1:length(𝐒.S)])
   timeaxis = 𝐒.S[1].t
   clim = (0,1)
   realProj ? projection = zeros(length(timeaxis)) : nothing
   realProj ? minVector = zeros(length(𝐒.S)) : nothing
   for k in 1:length(𝐒.S)
      rm = (!).(isnan.(𝐒.S[k].ω))# remove "NaN" values in 𝐂.ω
      seriescolor := colorSelect(colorMap)[ max.(min.(round.(Int, abs.((𝐒.S[k].a)) .* 256/a_max ),256),50) ][rm]

      @series begin
         𝐒.S[k].t[rm],𝐒.S[k].ω[rm],real.(𝐒.S[k].Ψ)[rm]
      end
      realProj ? minVector[k] = sort(unique(𝐒.S[k].ω))[1] : nothing
      realProj ? projection += real.(𝐒.S[k].Ψ) : nothing
   end
   if realProj
      append!(minVector,0)
      @series begin
         seriescolor := :white
         timeaxis, zeros(length(timeaxis)).+minimum(minVector), projection
      end
   end
end


text:
"""
    Kmeans(K; max=300, th=1e-4)
Kmeans method.

# Parameters:
- `K`: number of class
- `max`: maximum number of repitition
- `th`: converge threshold

# Example
```jldoctest
julia> model = Kmenas(3)
Kmeans{3}(Matrix{Float64}(undef, 0, 0), 100000000, 0.0001)

julia> using HorseML.Clustering: fit!

julia> fit!(model, x)

julia> model.labels |> size
(100, 3)
```
"""
mutable struct Kmeans{K}
    μ::Matrix{Float64}
    max::Int64
    th::Float64
    labels::Vector{Int}
    Kmeans(K; max = 300, th = 1e-4) = new{K}(Array{Float64}(undef, 0, 0), max, th, [])
end
n_clusters(model::Kmeans{K}) where K = K

function clustering(μ, x, K)
    r = zeros(size(x, 1), K)
    for n in 1 : size(x, 1)
        diff = Matrix(undef, size(μ)...)
        for k in 1 : K
            diff[k, :] = x[n, :] - μ[k, :]
        end
        distance = @. abs(diff)
        r[n, argmin(vec(sum(distance, dims=2)))] = 1
    end
    return r
end

function update(μ, x, K)
    r = clustering(μ, x, K)
    μ = similar(μ)
    for k in 1 : K
        μ[k, :] = sum(x .* r[:, k], dims=1) / sum(r[:, k])
    end
    return μ
end

"""
    fit!(model, x)
fit Kmeans or GMM model with data.

size of `x` is (numer of data, number of features).
"""
function fit!(model::Kmeans{K}, x) where {K}
    coverge(x, th) = @. abs(x) < th
    μ = kmeanspp(x, K)
    for _ in 1 : model.max
        μ_new = update(μ, x, K)
        if coverge.(μ_new - μ, model.th) == trues(size(μ)...)
            model.μ = μ_new
            labels = Array{Int}(undef, size(x, 1))
            r = clustering(model.μ, x, K)
            for i in 1 : size(r, 2)
                labels[findall(r[:, i].==1)] .= i
            end
            model.labels = labels
            return
        end
        μ = μ_new
    end
    @warn "Not Converged!"
    model.μ = μ
    labels = Array{Int}(undef, size(x, 1))
    r = clustering(model.μ, x, K)
    for i in 1 : size(r, 2)
        labels[findall(r[:, i].==1)] .= i
    end
    model.labels = labels
end

#initialize μ(centeroid vectors) with keams++
function kmeanspp(x, K)
    μ = Array{Float64, 2}(undef, K, size(x, 2))
    #pick first centeroid randomly
    μ[1, :] = x[rand(1:size(x, 1)), :]

    for i in 2 : K
        #compute the distance between each point and each centeroid
        distances = compute_distance(x, μ[1:i-1, :], size(x, 1), i-1)
        #compute the distance to the nearest centeroid
        closest_dist_sq = minimum(distances.^2, dims=2)
        #add the squares of each distances
        weights = sum(closest_dist_sq)

        rand_vals = rand() * weights
        #compute the culative sum of the squares of distances
        #Get the index of the data point where the cumulative sum and the value of rand_val are closest
        candidate_ids = searchsortedfirst(vec(cumsum(closest_dist_sq, dims=1)), rand_vals)
        #add centeroidx
        μ[i, :] = x[candidate_ids, :]
    end
    return μ
end

function compute_distance(x, μ, n, nc)
    dists = Array{Float64}(undef, n, nc)
    for i in 1 : nc
        dists[:, i] = sum((x.-μ[i, :]').^2, dims=2)
    end
    return dists
end


text:
"""
    Chi(ν)

The *Chi distribution* `ν` degrees of freedom has probability density function

```math
f(x; \\nu) = \\frac{1}{\\Gamma(\\nu/2)} 2^{1 - \\nu/2} x^{\\nu-1} e^{-x^2/2}, \\quad x > 0
```

It is the distribution of the square-root of a [`Chisq`](@ref) variate.

```julia
Chi(ν)       # Chi distribution with ν degrees of freedom

params(d)    # Get the parameters, i.e. (ν,)
dof(d)       # Get the degrees of freedom, i.e. ν
```

External links

* [Chi distribution on Wikipedia](http://en.wikipedia.org/wiki/Chi_distribution)

"""
struct Chi{T<:Real} <: ContinuousUnivariateDistribution
    ν::T
    Chi{T}(ν::T) where {T} = new{T}(ν)
end

function Chi(ν::Real; check_args::Bool=true)
    @check_args Chi (ν, ν > zero(ν))
    return Chi{typeof(ν)}(ν)
end

Chi(ν::Integer; check_args::Bool=true) = Chi(float(ν); check_args=check_args)

@distr_support Chi 0.0 Inf

### Conversions
convert(::Type{Chi{T}}, ν::Real) where {T<:Real} = Chi(T(ν))
convert(::Type{Chi{T}}, d::Chi{S}) where {T <: Real, S <: Real} = Chi(T(d.ν), check_args=false)

#### Parameters

dof(d::Chi) = d.ν
params(d::Chi) = (d.ν,)
@inline partype(d::Chi{T}) where {T<:Real} = T

#### Statistics

mean(d::Chi) = (h = d.ν/2; sqrt2 * exp(loggamma(h + 1//2) - loggamma(h)))

var(d::Chi) = d.ν - mean(d)^2
_chi_skewness(μ::Real, σ::Real) = (σ2 = σ^2; σ3 = σ2 * σ; (μ / σ3) * (1 - 2σ2))

function skewness(d::Chi)
    μ = mean(d)
    σ = sqrt(d.ν - μ^2)
    _chi_skewness(μ, σ)
end

function kurtosis(d::Chi)
    μ = mean(d)
    σ = sqrt(d.ν - μ^2)
    γ = _chi_skewness(μ, σ)
    (2/σ^2) * (1 - μ * σ * γ - σ^2)
end

entropy(d::Chi{T}) where {T<:Real} = (ν = d.ν;
    loggamma(ν/2) - T(logtwo)/2 - ((ν - 1)/2) * digamma(ν/2) + ν/2)

function mode(d::Chi; check_args::Bool=true)
    ν = d.ν
    @check_args(
        Chi,
        (ν, ν >= 1, "Chi distribution has no mode when ν < 1"),
    )
    sqrt(ν - 1)
end

#### Evaluation

function logpdf(d::Chi, x::Real)
    ν, _x = promote(d.ν, x)
    xsq = _x^2
    val = (xlogy(ν - 1, xsq / 2) - xsq + logtwo) / 2 - loggamma(ν / 2)
    return x < zero(x) ? oftype(val, -Inf) : val
end

gradlogpdf(d::Chi{T}, x::Real) where {T<:Real} = x >= 0 ? (d.ν - 1) / x - x : zero(T)

for f in (:cdf, :ccdf, :logcdf, :logccdf)
    @eval $f(d::Chi, x::Real) = $f(Chisq(d.ν; check_args=false), max(x, 0)^2)
end

for f in (:quantile, :cquantile, :invlogcdf, :invlogccdf)
    @eval $f(d::Chi, p::Real) = sqrt($f(Chisq(d.ν; check_args=false), p))
end

#### Sampling

rand(rng::AbstractRNG, d::Chi) =
    (ν = d.ν; sqrt(rand(rng, Gamma(ν / 2.0, 2.0one(ν)))))

struct ChiSampler{S <: Sampleable{Univariate,Continuous}} <:
    Sampleable{Univariate,Continuous}
    s::S
end

rand(rng::AbstractRNG, s::ChiSampler) = sqrt(rand(rng, s.s))

sampler(d::Chi) = ChiSampler(sampler(Chisq(d.ν)))


text:
include("parserMomkpZL.jl")
include("parserMomkpPG.jl")
include("vOptMomkp.jl")
include("displayGraphic.jl")
include("VNS.jl")
include("MOKP.jl")
include("tools.jl")
using CPLEX
name = "knapsack.500.2"
fname = "instancesZL/" * name

momkpZL = readInstanceMOMKPformatZL(false, fname)
dat1 = _bi01IP(momkpZL.P[1:2, :], momkpZL.W, momkpZL.ω)
solverSelected = CPLEX.Optimizer
YN, XE = vSolveBi01IP(solverSelected, dat1.C, dat1.A, dat1.b)
Yn_int = Vector{Vector{Int}}(undef, length(YN))
for j = 1:length(YN)
    Yn_int[j] = [0, 0]
    Yn_int[j][1] = round(Int, YN[j][1])
    Yn_int[j][2] = round(Int, YN[j][2])
end
TOOL_write_YN(Yn_int, name)


text:
function main_carlo(args)
    println("="^70)
    println("Début de l'action carlo")

    # Affiche l'ensemble des options du programme (cf fichier args.jl)
    # Args.show_args()   # Déjà fait dans Args.jl

    inst = Instance(args[:infile])

    # Si :itermax vaut 0 (i.e automatique) on choisit une valeur "pertinante"
    itermax_default = 500*inst.nb_planes
    itermax = args[:itermax] == 0 ? itermax_default : args[:itermax]

    # Le timingSolver sous-jacent est choisi en fonction de l'option
    # --timing-algo-solver (alias -t )
    # (Accessible par : Args.get(:timing_algo_solver))
    cursol = Solution(inst, update=false)
    # Attention : bestsol ne doit pas recopier cursol car cursol ne sera pas la
    # premiere solution explorée !
    bestsol = Solution(inst, update=false)

    @show cursol
    @show bestsol
    @show itermax
    @show Args.get(:level)

    ms_start = ms() # seconde depuis le démarrage avec précision à la ms
    for i in 1:itermax
        # On contruit un ordre "intelligent"

        # Cas particulier : si iterman==1 : on fait le meilleur coup possibe
        # en triant sur le target
        if itermax == 1
            initial_sort!(cursol, presort=:target)
        else
            shuffle!(cursol, do_update=false)
        end

        # Puis on lance l'algo de timing (presélectionné selon l'option --algo)
        # en fonction de l'ordre des avions ainsi contruit.
        solve!(cursol)

        lg(i," ")
        if  cursol.cost < bestsol.cost
            # On mémorise ce record
            copy!(bestsol, cursol)
            write(bestsol)
            # On enregistre cette nouvelle solution dans un fichier

            if lg1()
                # Ce code n'est exécuter que si --level vaut au moins 1
                print(i, ":" , bestsol.cost)
                if lg2()
                    # Ce code n'est exécuter que si --level vaut au moins 2
                    print(to_s(bestsol))
                end
                println()
            end
        end
    end
    ms_stop = ms()

    # le print_sol final n'est exécuté que si lg1() retourne true
    lg1() && print_sol(bestsol)

    nb_calls = bestsol.solver.nb_calls
    nb_infeasable = bestsol.solver.nb_infeasable
    nb_sec = round(ms_stop - ms_start, digits=3)
    nb_call_per_sec = round(nb_calls/nb_sec, digits=3)
    println("Performance: ")
    println("  nb_calls=$nb_calls")
    println("  nb_infeasable=$nb_infeasable")
    println("  nb_sec=$nb_sec")
    println("  => nb_call_per_sec = $nb_call_per_sec call/sec")

    println("Fin de l'action carlo")
end

main_carlo(Args.args)


text:
using Test
using DataFrames
using Deldir
using Gadfly
using GeometricalPredicates

include("../src/scaler.jl")
include("../src/edges.jl")

bbox = [0.0 0.0; 1.0 1.0]

function plotexample()
    x = [0.2, 0.3, 0.4, 0.4, 0.5, 0.6]
    y = [0.0, 0.2, 0.0, 0.4, 0.2, 0.0]
    del, vor, sum = deldir(x, y)
    x, y = edges(del)
    println(del)
    println("-----------------------------")
    println(sum)
    Gadfly.plot(x = x, y = y, Geom.path, Coord.cartesian(fixed = true))
end

@testset "Compute centroid from scaled triangle" begin
    a = Point(0.0, 0.0)
    b = Point(0.5, 0.0)
    c = Point(0.0, 0.5)
    triangle = [a, b, c]
    scaler = Scaler(bbox)
    expectedcenterpoint = Point(0.16666666666666666, 0.16666666666666666)
    @test unscaledcenterpoint(triangle, scaler) == expectedcenterpoint
end

@testset "Build vectoriezed edges" begin
    a = Point(0.0, 0.0)
    b = Point(0.5, 0.0)
    c = Point(0.0, 0.5)
    triangle = [a, b, c]
    expectedvectorizededges = [
        [0.0 0.0; 0.5 0.0],
        [0.5 0.0; 0.0 0.5],
        [0.0 0.5; 0.0 0.0],
        [0.0 0.0; 0.0 0.5],
        [0.0 0.5; 0.5 0.0],
        [0.5 0.0; 0.0 0.0]
    ]
    vedges = vectorizededges(triangle)
    @test size(vedges, 1) == size(findall(in(expectedvectorizededges), vedges), 1)
end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
module Derivatives

using Sylvia

export derivative,
    a, b, c, d,
    ctx

const ctx = Sylvia.Context(Sylvia.@__context__)

function derivative end
Sylvia.@register ctx derivative 2

@sym [Number] a b c d

@scope ctx let
    @sym [Wild{Function}] f
    @sym [Wild{Number}] x y z

    @! set -(-x) --> x

    @! set derivative(f(x), y) --> derivative(f(x), x) * derivative(x, y)
    @! set derivative(x, x) --> 1
    @! set derivative(-x, y) --> -derivative(x, y)
    @! set derivative(x + y, z) --> derivative(x, z) + derivative(y, z)
    @! set derivative(x - y, z) --> derivative(x, z) - derivative(y, z)
    @! set derivative(x*y, z) --> derivative(x, z)*y + x*derivative(y, z)
    @! set derivative(x^y, x) --> y*x^(1 - y)
    @! set derivative(cos(x), x) --> -sin(x)
    @! set derivative(sin(x), x) --> cos(x)
end

end