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

@test integrate(sinh(c+d*x)^4*(a+b*sech(c+d*x)^2), x) == :(x*(-3//2*b+3//8*a)+b*d^-1*tanh(c+d*x)+-1//8*d^-1*(-4b+5a)*cosh(c+d*x)*sinh(c+d*x)+1//4*a*d^-1*cosh(c+d*x)^3*sinh(c+d*x))
@test integrate(sinh(c+d*x)^3*(a+b*sech(c+d*x)^2), x) == :(b*d^-1*sech(c+d*x)+-1*d^-1*(a+-1b)*cosh(c+d*x)+1//3*a*d^-1*cosh(c+d*x)^3)
@test integrate(sinh(c+d*x)^2*(a+b*sech(c+d*x)^2), x) == :(x*(b+-1//2*a)+-1*b*d^-1*tanh(c+d*x)+(1/2)*a*d^-1*cosh(c+d*x)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)*sinh(c+d*x), x) == :(a*d^-1*cosh(c+d*x)+-1*b*d^-1*sech(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)*csch(c+d*x), x) == :(b*d^-1*sech(c+d*x)+-1*d^-1*(a+b)*arctanh(cosh(c+d*x)))
@test integrate(csch(c+d*x)^2*(a+b*sech(c+d*x)^2), x) == :(-1*b*d^-1*tanh(c+d*x)+-1*d^-1*(a+b)*coth(c+d*x))
@test integrate(csch(c+d*x)^3*(a+b*sech(c+d*x)^2), x) == :((1/2)*d^-1*(a+3b)*arctanh(cosh(c+d*x))+-1*b*d^-1*sech(c+d*x)+-1//2*d^-1*(a+b)*coth(c+d*x)*csch(c+d*x))
@test integrate(csch(c+d*x)^4*(a+b*sech(c+d*x)^2), x) == :(b*d^-1*tanh(c+d*x)+d^-1*(a+2b)*coth(c+d*x)+-1//3*d^-1*coth(c+d*x)^3*(a+b))
@test integrate((a+b*sech(c+d*x)^2)^2*sinh(c+d*x)^4, x) == :(x*(b^2+3//8*a^2+-3*a*b)+-1//3*b^2*d^-1*tanh(c+d*x)^3+-1//4*d^-1*(a^2+4*b^2+-8*a*b)*tanh(c+d*x)+1//4*a^2*d^-1*sinh(c+d*x)^4*tanh(c+d*x)+-1//8*a*d^-1*(a+-8b)*cosh(c+d*x)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*sinh(c+d*x)^3, x) == :(1//3*a^2*d^-1*cosh(c+d*x)^3+1//3*b^2*d^-1*sech(c+d*x)^3+b*d^-1*(-1b+2a)*sech(c+d*x)+-1*a*d^-1*(a+-2b)*cosh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*sinh(c+d*x)^2, x) == :(-1//2*a*x*(a+-4b)+1//3*b^2*d^-1*tanh(c+d*x)^3+(1/2)*a*d^-1*(a+-4b)*tanh(c+d*x)+(1/2)*a^2*d^-1*sinh(c+d*x)^2*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*sinh(c+d*x), x) == :(a^2*d^-1*cosh(c+d*x)+-1//3*b^2*d^-1*sech(c+d*x)^3+-2*a*b*d^-1*sech(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*csch(c+d*x), x) == :(-1*d^-1*(a+b)^2*arctanh(cosh(c+d*x))+1//3*b^2*d^-1*sech(c+d*x)^3+b*d^-1*(b+2a)*sech(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*csch(c+d*x)^2, x) == :(-1*d^-1*(a+b)^2*coth(c+d*x)+1//3*b^2*d^-1*tanh(c+d*x)^3+-2*b*d^-1*(a+b)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*csch(c+d*x)^3, x) == :((1/2)*d^-1*(a+b)*(a+5b)*arctanh(cosh(c+d*x))+-1//3*b*d^-1*(5b+6a)*sech(c+d*x)+-1//6*d^-1*(3*a^2+5*b^2+6*a*b)*coth(c+d*x)*csch(c+d*x)+1//3*b^2*d^-1*csch(c+d*x)^2*sech(c+d*x)^3)
@test integrate((a+b*sech(c+d*x)^2)^2*csch(c+d*x)^4, x) == :(-1//3*b^2*d^-1*tanh(c+d*x)^3+-1//3*d^-1*(a+b)^2*coth(c+d*x)^3+b*d^-1*(2a+3b)*tanh(c+d*x)+d^-1*(a+b)*(a+3b)*coth(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^3*sinh(c+d*x)^4, x) == :(3//8*a*x*(a^2+8*b^2+-12*a*b)+-3//8*a*d^-1*(a^2+8*b^2+-12*a*b)*tanh(c+d*x)+-1//40*b^2*d^-1*tanh(c+d*x)^5*(-48b+15a)+1//4*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^3*sinh(c+d*x)^3*cosh(c+d*x)+1//8*b*d^-1*tanh(c+d*x)^3*(-8*b^2+6*a^2+-23*a*b)+-1//8*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^2*sinh(c+d*x)^2*(-6b+3a)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^3*sinh(c+d*x)^3, x) == :(1//3*a^3*d^-1*cosh(c+d*x)^3+1//5*b^3*d^-1*sech(c+d*x)^5+-1*a^2*d^-1*(a+-3b)*cosh(c+d*x)+1//3*b^2*d^-1*sech(c+d*x)^3*(-1b+3a)+3*a*b*d^-1*(a+-1b)*sech(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^3*sinh(c+d*x), x) == :(a^3*d^-1*cosh(c+d*x)+-1//5*b^3*d^-1*sech(c+d*x)^5+-1*a*b^2*d^-1*sech(c+d*x)^3+-3*b*a^2*d^-1*sech(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^3*csch(c+d*x), x) == :(-1*d^-1*(a+b)^3*arctanh(cosh(c+d*x))+1//5*b^3*d^-1*sech(c+d*x)^5+b*d^-1*(b^2+3*a^2+3*a*b)*sech(c+d*x)+1//3*b^2*d^-1*sech(c+d*x)^3*(b+3a))
@test integrate((a+b*sech(c+d*x)^2)^3*csch(c+d*x)^2, x) == :(-1*d^-1*(a+b)^3*coth(c+d*x)+-1//5*b^3*d^-1*tanh(c+d*x)^5+b^2*d^-1*tanh(c+d*x)^3*(a+b)+-3*b*d^-1*(a+b)^2*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^3*csch(c+d*x)^3, x) == :((1/2)*d^-1*(a+b)^2*(a+7b)*arctanh(cosh(c+d*x))+-1//2*d^-1*(a+b)^2*(a+7b)*sech(c+d*x)+-1//6*b*d^-1*sech(c+d*x)^3*(6*a^2+7*b^2+15*a*b)+-1//10*b^2*d^-1*sech(c+d*x)^5*(5a+7b)+-1//2*d^-1*(b+a*cosh(c+d*x)^2)^2*csch(c+d*x)^2*sech(c+d*x)^5*(a+b))
@test integrate((a+b*sech(c+d*x)^2)^3*csch(c+d*x)^4, x) == :(-1//3*d^-1*(a+b)^3*coth(c+d*x)^3+1//5*b^3*d^-1*tanh(c+d*x)^5+d^-1*(a+b)^2*(a+4b)*coth(c+d*x)+-1//3*b^2*d^-1*tanh(c+d*x)^3*(3a+4b)+3*b*d^-1*(a+b)*(a+2b)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-1*sinh(c+d*x)^4, x) == :(1//8*x*a^-3*(3*a^2+8*b^2+12*a*b)+1//4*a^-1*d^-1*cosh(c+d*x)^3*sinh(c+d*x)+-1*a^-3*b^(1/2)*d^-1*(a+b)^3//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+-1//8*a^-2*d^-1*(4b+5a)*cosh(c+d*x)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-1*sinh(c+d*x)^3, x) == :(1//3*a^-1*d^-1*cosh(c+d*x)^3+-1*a^-2*d^-1*(a+b)*cosh(c+d*x)+a^-5//2*b^(1/2)*d^-1*(a+b)*arctan(a^(1/2)*b^-1//2*cosh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*sinh(c+d*x)^2, x) == :(-1//2*x*a^-2*(a+2b)+(1/2)*a^-1*d^-1*cosh(c+d*x)*sinh(c+d*x)+a^-2*b^(1/2)*d^-1*(a+b)^(1/2)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*sinh(c+d*x), x) == :(a^-1*d^-1*cosh(c+d*x)+-1*a^-3//2*b^(1/2)*d^-1*arctan(a^(1/2)*b^-1//2*cosh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*csch(c+d*x), x) == :(-1*d^-1*(a+b)^-1*arctanh(cosh(c+d*x))+a^-1//2*b^(1/2)*d^-1*(a+b)^-1*arctan(a^(1/2)*b^-1//2*cosh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*csch(c+d*x)^2, x) == :(-1*d^-1*(a+b)^-1*coth(c+d*x)+b^(1/2)*d^-1*(a+b)^-3//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*csch(c+d*x)^3, x) == :((1/2)*d^-1*(a+b)^-2*(a+-1b)*arctanh(cosh(c+d*x))+-1*d^-1*(2a+2b)^-1*coth(c+d*x)*csch(c+d*x)+-1*a^(1/2)*b^(1/2)*d^-1*(a+b)^-2*arctan(a^(1/2)*b^-1//2*cosh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*csch(c+d*x)^4, x) == :(-1*d^-1*(3a+3b)^-1*coth(c+d*x)^3+a*d^-1*(a+b)^-2*coth(c+d*x)+-1*a*b^(1/2)*d^-1*(a+b)^-5//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-2*sinh(c+d*x)^4, x) == :(1//8*x*a^-4*(3*a^2+24*b^2+24*a*b)+1//4*a^-1*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*cosh(c+d*x)^3*sinh(c+d*x)+-3//2*a^-4*b^(1/2)*d^-1*(a+b)^(1/2)*(a+2b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+-3//8*b*a^-3*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*(3a+4b)*tanh(c+d*x)+-1//8*a^-2*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*(5a+6b)*cosh(c+d*x)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-2*sinh(c+d*x)^3, x) == :(1//3*a^-2*d^-1*cosh(c+d*x)^3+-1*a^-3*d^-1*(a+2b)*cosh(c+d*x)+(1/2)*a^-7//2*b^(1/2)*d^-1*(3a+5b)*arctan(a^(1/2)*b^-1//2*cosh(c+d*x))+-1//2*b*a^-3*d^-1*(b+a*cosh(c+d*x)^2)^-1*(a+b)*cosh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-2*sinh(c+d*x)^2, x) == :(-1//2*x*a^-3*(a+4b)+b*a^-2*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*tanh(c+d*x)+(1/2)*a^-1*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*cosh(c+d*x)*sinh(c+d*x)+(1/2)*a^-3*b^(1/2)*d^-1*(a+b)^-1//2*(3a+4b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-2*sinh(c+d*x), x) == :(3//2*a^-2*d^-1*cosh(c+d*x)+-3//2*a^-5//2*b^(1/2)*d^-1*arctan(a^(1/2)*b^-1//2*cosh(c+d*x))+-1//2*a^-1*d^-1*(b+a*cosh(c+d*x)^2)^-1*cosh(c+d*x)^3)
@test integrate((a+b*sech(c+d*x)^2)^-2*csch(c+d*x), x) == :(-1*d^-1*(a+b)^-2*arctanh(cosh(c+d*x))+(1/2)*a^-3//2*b^(1/2)*d^-1*(a+b)^-2*(b+3a)*arctan(a^(1/2)*b^-1//2*cosh(c+d*x))+-1//2*b*a^-1*d^-1*(a+b)^-1*(b+a*cosh(c+d*x)^2)^-1*cosh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-2*csch(c+d*x)^2, x) == :(-3//2*d^-1*(a+b)^-2*coth(c+d*x)+d^-1*(2a+2b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*coth(c+d*x)+3//2*b^(1/2)*d^-1*(a+b)^-5//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-2*csch(c+d*x)^3, x) == :((1/2)*d^-1*(a+b)^-3*(a+-3b)*arctanh(cosh(c+d*x))+-1*d^-1*(b+a*cosh(c+d*x)^2)^-1*(2a+2b)^-1*coth(c+d*x)*csch(c+d*x)+-1//2*d^-1*(a+b)^-2*(b+a*cosh(c+d*x)^2)^-1*(a+-1b)*cosh(c+d*x)+-1//2*a^-1//2*b^(1/2)*d^-1*(a+b)^-3*(-1b+3a)*arctan(a^(1/2)*b^-1//2*cosh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-2*csch(c+d*x)^4, x) == :(-1//3*d^-1*(a+b)^-2*coth(c+d*x)^3+d^-1*(a+b)^-3*(a+-1b)*coth(c+d*x)+-1//2*b^(1/2)*d^-1*(a+b)^-7//2*(-2b+3a)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+-1//2*a*b*d^-1*(a+b)^-3*(a+b+-1*b*tanh(c+d*x)^2)^-1*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*sinh(c+d*x)^4, x) == :(1//8*x*a^-5*(3*a^2+48*b^2+36*a*b)+1//4*a^-1*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*cosh(c+d*x)^3*sinh(c+d*x)+-3//2*b*a^-4*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*(a+2b)*tanh(c+d*x)+-3//8*a^-5*b^(1/2)*d^-1*(a+b)^-1//2*(5*a^2+16*b^2+20*a*b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+-1//8*b*a^-3*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*(7a+12b)*tanh(c+d*x)+-1//8*a^-2*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*(5a+8b)*cosh(c+d*x)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*sinh(c+d*x)^3, x) == :(1//3*a^-3*d^-1*cosh(c+d*x)^3+-1*a^-4*d^-1*(a+3b)*cosh(c+d*x)+5//8*a^-9//2*b^(1/2)*d^-1*(3a+7b)*arctan(a^(1/2)*b^-1//2*cosh(c+d*x))+-1//8*b*a^-4*d^-1*(b+a*cosh(c+d*x)^2)^-1*(9a+13b)*cosh(c+d*x)+1//4*a^-4*b^2*d^-1*(b+a*cosh(c+d*x)^2)^-2*(a+b)*cosh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*sinh(c+d*x)^2, x) == :(-1//2*x*a^-4*(a+6b)+(1/2)*a^-1*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*cosh(c+d*x)*sinh(c+d*x)+3//4*b*a^-2*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*tanh(c+d*x)+1//8*a^-4*b^(1/2)*d^-1*(a+b)^-3//2*(15*a^2+24*b^2+40*a*b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+1//8*b*a^-3*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*(11a+12b)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*sinh(c+d*x), x) == :(15//8*a^-3*d^-1*cosh(c+d*x)+-15//8*a^-7//2*b^(1/2)*d^-1*arctan(a^(1/2)*b^-1//2*cosh(c+d*x))+-5//8*a^-2*d^-1*(b+a*cosh(c+d*x)^2)^-1*cosh(c+d*x)^3+-1//4*a^-1*d^-1*(b+a*cosh(c+d*x)^2)^-2*cosh(c+d*x)^5)
@test integrate((a+b*sech(c+d*x)^2)^-3*csch(c+d*x), x) == :(-1*d^-1*(a+b)^-3*arctanh(cosh(c+d*x))+-1//4*b*a^-1*d^-1*(a+b)^-1*(b+a*cosh(c+d*x)^2)^-2*cosh(c+d*x)^3+1//8*a^-5//2*b^(1/2)*d^-1*(a+b)^-3*(3*b^2+15*a^2+10*a*b)*arctan(a^(1/2)*b^-1//2*cosh(c+d*x))+-1//8*b*a^-2*d^-1*(a+b)^-2*(b+a*cosh(c+d*x)^2)^-1*(3b+7a)*cosh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*csch(c+d*x)^2, x) == :(-15//8*d^-1*(a+b)^-3*coth(c+d*x)+d^-1*(4a+4b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*coth(c+d*x)+5//8*d^-1*(a+b)^-2*(a+b+-1*b*tanh(c+d*x)^2)^-1*coth(c+d*x)+15//8*b^(1/2)*d^-1*(a+b)^-7//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-3*csch(c+d*x)^3, x) == :((1/2)*d^-1*(a+b)^-4*(a+-5b)*arctanh(cosh(c+d*x))+-1*d^-1*(b+a*cosh(c+d*x)^2)^-2*(2a+2b)^-1*coth(c+d*x)^2*cosh(c+d*x)+-1//8*a^-1*d^-1*(a+b)^-3*(b+a*cosh(c+d*x)^2)^-1*(-1*b^2+4*a^2+-9*a*b)*cosh(c+d*x)+-1//8*a^-3//2*b^(1/2)*d^-1*(a+b)^-4*(-1*b^2+15*a^2+-10*a*b)*arctan(a^(1/2)*b^-1//2*cosh(c+d*x))+1//4*b*a^-1*d^-1*(a+b)^-2*(b+a*cosh(c+d*x)^2)^-2*(-1b+2a)*cosh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*csch(c+d*x)^4, x) == :(-1//3*d^-1*(a+b)^-3*coth(c+d*x)^3+d^-1*(a+b)^-4*(a+-2b)*coth(c+d*x)+-1//8*b^(1/2)*d^-1*(a+b)^-9//2*(-20b+15a)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+-1//4*a*b*d^-1*(a+b)^-3*(a+b+-1*b*tanh(c+d*x)^2)^-2*tanh(c+d*x)+-1//8*b*d^-1*(a+b)^-4*(a+b+-1*b*tanh(c+d*x)^2)^-1*(-4b+7a)*tanh(c+d*x))
@test integrate(cosh(c+d*x)^4*(a+b*sech(c+d*x)^2), x) == :(1//8*x*(3a+4b)+1//4*a*d^-1*cosh(c+d*x)^3*sinh(c+d*x)+1//8*d^-1*(3a+4b)*cosh(c+d*x)*sinh(c+d*x))
@test integrate(cosh(c+d*x)^3*(a+b*sech(c+d*x)^2), x) == :(d^-1*(a+b)*sinh(c+d*x)+1//3*a*d^-1*sinh(c+d*x)^3)
@test integrate(cosh(c+d*x)^2*(a+b*sech(c+d*x)^2), x) == :((1/2)*x*(a+2b)+(1/2)*a*d^-1*cosh(c+d*x)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)*cosh(c+d*x), x) == :(a*d^-1*sinh(c+d*x)+b*d^-1*arctan(sinh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)*sech(c+d*x), x) == :((1/2)*d^-1*(b+2a)*arctan(sinh(c+d*x))+(1/2)*b*d^-1*sech(c+d*x)*tanh(c+d*x))
@test integrate(sech(c+d*x)^3*(a+b*sech(c+d*x)^2), x) == :(1//8*d^-1*(3b+4a)*arctan(sinh(c+d*x))+1//4*b*d^-1*sech(c+d*x)^3*tanh(c+d*x)+1//8*d^-1*(3b+4a)*sech(c+d*x)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*cosh(c+d*x)^4, x) == :(x*(b^2+3//8*a^2+a*b)+1//4*a*d^-1*cosh(c+d*x)^3*(a+b+-1*b*tanh(c+d*x)^2)*sinh(c+d*x)+3//8*a*d^-1*(a+2b)*cosh(c+d*x)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*cosh(c+d*x)^3, x) == :(b^2*d^-1*arctan(sinh(c+d*x))+1//3*a^2*d^-1*sinh(c+d*x)^3+a*d^-1*(a+2b)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*cosh(c+d*x)^2, x) == :(b^2*d^-1*tanh(c+d*x)+(1/2)*a*x*(a+4b)+(1/2)*a^2*d^-1*cosh(c+d*x)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*cosh(c+d*x), x) == :(a^2*d^-1*sinh(c+d*x)+(1/2)*b*d^-1*(b+4a)*arctan(sinh(c+d*x))+(1/2)*b^2*d^-1*sech(c+d*x)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*sech(c+d*x), x) == :(1//8*d^-1*(3*b^2+8*a^2+8*a*b)*arctan(sinh(c+d*x))+1//4*b*d^-1*sech(c+d*x)^3*(a+b+a*sinh(c+d*x)^2)*tanh(c+d*x)+3//8*b*d^-1*(b+2a)*sech(c+d*x)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*sech(c+d*x)^2, x) == :(d^-1*(a+b)^2*tanh(c+d*x)+1//5*b^2*d^-1*tanh(c+d*x)^5+-2//3*b*d^-1*tanh(c+d*x)^3*(a+b))
@test integrate((a+b*sech(c+d*x)^2)^2*sech(c+d*x)^3, x) == :(1//16*d^-1*(5*b^2+8*a^2+12*a*b)*arctan(sinh(c+d*x))+1//16*d^-1*(5*b^2+8*a^2+12*a*b)*sech(c+d*x)*tanh(c+d*x)+1//6*b*d^-1*sech(c+d*x)^5*(a+b+a*sinh(c+d*x)^2)*tanh(c+d*x)+1//24*b*d^-1*sech(c+d*x)^3*(5b+8a)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*sech(c+d*x)^4, x) == :(d^-1*(a+b)^2*tanh(c+d*x)+-1//7*b^2*d^-1*tanh(c+d*x)^7+-1//3*d^-1*tanh(c+d*x)^3*(a+b)*(a+3b)+1//5*b*d^-1*tanh(c+d*x)^5*(2a+3b))
@test integrate((a+b*sech(c+d*x)^2)^3*cosh(c+d*x)^4, x) == :(b^3*d^-1*tanh(c+d*x)+3//8*a*x*(a^2+8*b^2+4*a*b)+1//4*a^3*d^-1*cosh(c+d*x)^3*sinh(c+d*x)+3//8*a^2*d^-1*(a+4b)*cosh(c+d*x)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^3*cosh(c+d*x)^3, x) == :(1//3*a^3*d^-1*sinh(c+d*x)^3+a^2*d^-1*(a+3b)*sinh(c+d*x)+(1/2)*b^2*d^-1*(b+6a)*arctan(sinh(c+d*x))+(1/2)*b^3*d^-1*sech(c+d*x)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^3*cosh(c+d*x)^2, x) == :((1/2)*x*a^2*(a+6b)+-1//3*b^3*d^-1*tanh(c+d*x)^3+b^2*d^-1*(b+3a)*tanh(c+d*x)+(1/2)*a^3*d^-1*cosh(c+d*x)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^3*cosh(c+d*x), x) == :(a^3*d^-1*sinh(c+d*x)+1//4*b^3*d^-1*sech(c+d*x)^3*tanh(c+d*x)+3//8*b*d^-1*(b^2+8*a^2+4*a*b)*arctan(sinh(c+d*x))+3//8*b^2*d^-1*(b+4a)*sech(c+d*x)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^3*sech(c+d*x), x) == :(1//16*d^-1*(b+2a)*(5*b^2+8*a^2+8*a*b)*arctan(sinh(c+d*x))+1//6*b*d^-1*(a+b+a*sinh(c+d*x)^2)^2*sech(c+d*x)^5*tanh(c+d*x)+1//48*b*d^-1*(15*b^2+44*a^2+44*a*b)*sech(c+d*x)*tanh(c+d*x)+5//24*b*d^-1*sech(c+d*x)^3*(b+2a)*(a+b+a*sinh(c+d*x)^2)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^3*sech(c+d*x)^2, x) == :(d^-1*(a+b)^3*tanh(c+d*x)+-1//7*b^3*d^-1*tanh(c+d*x)^7+-1*b*d^-1*(a+b)^2*tanh(c+d*x)^3+3//5*b^2*d^-1*tanh(c+d*x)^5*(a+b))
@test integrate((a+b*sech(c+d*x)^2)^3*sech(c+d*x)^3, x) == :(1//128*d^-1*(35*b^3+64*a^3+120*a*b^2+144*b*a^2)*arctan(sinh(c+d*x))+1//128*d^-1*(35*b^3+64*a^3+120*a*b^2+144*b*a^2)*sech(c+d*x)*tanh(c+d*x)+1//8*b*d^-1*(a+b+a*sinh(c+d*x)^2)^2*sech(c+d*x)^7*tanh(c+d*x)+1//192*b*d^-1*sech(c+d*x)^3*(35*b^2+72*a^2+92*a*b)*tanh(c+d*x)+1//48*b*d^-1*sech(c+d*x)^5*(7b+12a)*(a+b+a*sinh(c+d*x)^2)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^3*sech(c+d*x)^4, x) == :(d^-1*(a+b)^3*tanh(c+d*x)+1//9*b^3*d^-1*tanh(c+d*x)^9+-1//3*d^-1*(a+b)^2*tanh(c+d*x)^3*(a+4b)+-1//7*b^2*d^-1*tanh(c+d*x)^7*(3a+4b)+3//5*b*d^-1*tanh(c+d*x)^5*(a+b)*(a+2b))
@test integrate((a+b*sech(c+d*x)^2)^-1*cosh(c+d*x)^4, x) == :(1//8*x*a^-3*(3*a^2+8*b^2+-4*a*b)+1//4*a^-1*d^-1*cosh(c+d*x)^3*sinh(c+d*x)+-1*a^-3*b^5//2*d^-1*(a+b)^-1//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+1//8*a^-2*d^-1*(-4b+3a)*cosh(c+d*x)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-1*cosh(c+d*x)^3, x) == :(1//3*a^-1*d^-1*sinh(c+d*x)^3+a^-2*d^-1*(a+-1b)*sinh(c+d*x)+a^-5//2*b^2*d^-1*(a+b)^-1//2*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*cosh(c+d*x)^2, x) == :((1/2)*x*a^-2*(a+-2b)+(1/2)*a^-1*d^-1*cosh(c+d*x)*sinh(c+d*x)+a^-2*b^3//2*d^-1*(a+b)^-1//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*cosh(c+d*x), x) == :(a^-1*d^-1*sinh(c+d*x)+-1*b*a^-3//2*d^-1*(a+b)^-1//2*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*sech(c+d*x), x) == :(a^-1//2*d^-1*(a+b)^-1//2*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*sech(c+d*x)^2, x) == :(b^-1//2*d^-1*(a+b)^-1//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*sech(c+d*x)^3, x) == :(b^-1*d^-1*arctan(sinh(c+d*x))+-1*a^(1/2)*b^-1*d^-1*(a+b)^-1//2*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*sech(c+d*x)^4, x) == :(b^-1*d^-1*tanh(c+d*x)+-1*a*b^-3//2*d^-1*(a+b)^-1//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*sech(c+d*x)^5, x) == :((1/2)*b^-1*d^-1*sech(c+d*x)*tanh(c+d*x)+-1//2*b^-2*d^-1*(-1b+2a)*arctan(sinh(c+d*x))+a^3//2*b^-2*d^-1*(a+b)^-1//2*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*sech(c+d*x)^6, x) == :(-1//3*b^-1*d^-1*tanh(c+d*x)^3+-1*b^-2*d^-1*(a+-1b)*tanh(c+d*x)+a^2*b^-5//2*d^-1*(a+b)^-1//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-2*cosh(c+d*x)^3, x) == :(1//3*a^-2*d^-1*sinh(c+d*x)^3+a^-3*d^-1*(a+-2b)*sinh(c+d*x)+(1/2)*a^-7//2*b^2*d^-1*(a+b)^-3//2*(5b+6a)*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x))+-1//2*a^-3*b^3*d^-1*(a+b)^-1*(a+b+a*sinh(c+d*x)^2)^-1*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-2*cosh(c+d*x)^2, x) == :((1/2)*x*a^-3*(a+-4b)+(1/2)*a^-1*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*cosh(c+d*x)*sinh(c+d*x)+(1/2)*a^-3*b^3//2*d^-1*(a+b)^-3//2*(4b+5a)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+(1/2)*b*a^-2*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*(a+2b)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-2*cosh(c+d*x), x) == :(a^-2*d^-1*sinh(c+d*x)+(1/2)*a^-2*b^2*d^-1*(a+b)^-1*(a+b+a*sinh(c+d*x)^2)^-1*sinh(c+d*x)+-1//2*b*a^-5//2*d^-1*(a+b)^-3//2*(3b+4a)*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-2*sech(c+d*x), x) == :((1/2)*a^-3//2*d^-1*(a+b)^-3//2*(b+2a)*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x))+-1//2*b*a^-1*d^-1*(a+b)^-1*(a+b+a*sinh(c+d*x)^2)^-1*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-2*sech(c+d*x)^2, x) == :(d^-1*(2a+2b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*tanh(c+d*x)+(1/2)*b^-1//2*d^-1*(a+b)^-3//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-2*sech(c+d*x)^3, x) == :(d^-1*(2a+2b)^-1*(a+b+a*sinh(c+d*x)^2)^-1*sinh(c+d*x)+(1/2)*a^-1//2*d^-1*(a+b)^-3//2*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-2*sech(c+d*x)^4, x) == :((1/2)*b^-3//2*d^-1*(a+b)^-3//2*(a+2b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+-1//2*a*b^-1*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-2*sech(c+d*x)^5, x) == :(b^-2*d^-1*arctan(sinh(c+d*x))+-1//2*a*b^-1*d^-1*(a+b)^-1*(a+b+a*sinh(c+d*x)^2)^-1*sinh(c+d*x)+-1//2*a^(1/2)*b^-2*d^-1*(a+b)^-3//2*(2a+3b)*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-2*sech(c+d*x)^6, x) == :(b^-2*d^-1*tanh(c+d*x)+(1/2)*a^2*b^-2*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*tanh(c+d*x)+-1//2*a*b^-5//2*d^-1*(a+b)^-3//2*(3a+4b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-2*sech(c+d*x)^7, x) == :(-1//2*b^-3*d^-1*(-1b+4a)*arctan(sinh(c+d*x))+(1/2)*b^-1*d^-1*(a+b+a*sinh(c+d*x)^2)^-1*sech(c+d*x)*tanh(c+d*x)+(1/2)*a^3//2*b^-3*d^-1*(a+b)^-3//2*(4a+5b)*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x))+(1/2)*a*b^-2*d^-1*(a+b)^-1*(a+b+a*sinh(c+d*x)^2)^-1*(b+2a)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*cosh(c+d*x)^2, x) == :((1/2)*x*a^-4*(a+-6b)+(1/2)*a^-1*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*cosh(c+d*x)*sinh(c+d*x)+1//8*a^-4*b^3//2*d^-1*(a+b)^-5//2*(24*b^2+35*a^2+56*a*b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+1//4*b*a^-2*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*(2a+3b)*tanh(c+d*x)+1//8*b*a^-3*d^-1*(a+b)^-2*(a+b+-1*b*tanh(c+d*x)^2)^-1*(a+4b)*(3b+4a)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*cosh(c+d*x), x) == :(a^-3*d^-1*sinh(c+d*x)+-3//8*b*a^-7//2*d^-1*(a+b)^-5//2*((b+2a)^2+4*(a+b)^2)*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x))+-1//4*a^-3*b^3*d^-1*(a+b)^-1*(a+b+a*sinh(c+d*x)^2)^-2*sinh(c+d*x)+3//8*a^-3*b^2*d^-1*(a+b)^-2*(a+b+a*sinh(c+d*x)^2)^-1*(3b+4a)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*sech(c+d*x), x) == :(1//8*a^-5//2*d^-1*(a+b)^-5//2*(3*b^2+8*a^2+8*a*b)*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x))+-3//8*b*a^-2*d^-1*(a+b)^-2*(a+b+a*sinh(c+d*x)^2)^-1*(b+2a)*sinh(c+d*x)+-1//4*b*a^-1*d^-1*(a+b)^-1*(a+b+a*sinh(c+d*x)^2)^-2*cosh(c+d*x)^2*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*sech(c+d*x)^2, x) == :(d^-1*(4a+4b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*tanh(c+d*x)+3//8*b^-1//2*d^-1*(a+b)^-5//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+3//8*d^-1*(a+b)^-2*(a+b+-1*b*tanh(c+d*x)^2)^-1*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*sech(c+d*x)^3, x) == :(1//8*a^-3//2*d^-1*(a+b)^-5//2*(b+4a)*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x))+-1//4*b*a^-1*d^-1*(a+b)^-1*(a+b+a*sinh(c+d*x)^2)^-2*sinh(c+d*x)+1//8*a^-1*d^-1*(a+b)^-2*(a+b+a*sinh(c+d*x)^2)^-1*(b+4a)*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*sech(c+d*x)^4, x) == :(1//8*b^-3//2*d^-1*(a+b)^-5//2*(a+4b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+-1//4*a*b^-1*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*tanh(c+d*x)+1//8*b^-1*d^-1*(a+b)^-2*(a+b+-1*b*tanh(c+d*x)^2)^-1*(a+4b)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*sech(c+d*x)^5, x) == :(d^-1*(4a+4b)^-1*(a+b+a*sinh(c+d*x)^2)^-2*sinh(c+d*x)+3//8*a^-1//2*d^-1*(a+b)^-5//2*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x))+3//8*d^-1*(a+b)^-2*(a+b+a*sinh(c+d*x)^2)^-1*sinh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*sech(c+d*x)^6, x) == :(1//8*b^-5//2*d^-1*(a+b)^-5//2*(3*a^2+8*b^2+8*a*b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+-3//8*a*b^-2*d^-1*(a+b)^-2*(a+b+-1*b*tanh(c+d*x)^2)^-1*(a+2b)*tanh(c+d*x)+-1//4*a*b^-1*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*sech(c+d*x)^2*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*sech(c+d*x)^7, x) == :(b^-3*d^-1*arctan(sinh(c+d*x))+-1//4*a*b^-1*d^-1*(a+b)^-1*(a+b+a*sinh(c+d*x)^2)^-2*sinh(c+d*x)+-1//8*a^(1/2)*b^-3*d^-1*(a+b)^-5//2*(8*a^2+15*b^2+20*a*b)*arctan(a^(1/2)*(a+b)^-1//2*sinh(c+d*x))+-1//8*a*b^-2*d^-1*(a+b)^-2*(a+b+a*sinh(c+d*x)^2)^-1*(4a+7b)*sinh(c+d*x))
@test integrate(tanh(c+d*x)^4*(a+b*sech(c+d*x)^2), x) == :(a*x+-1*a*d^-1*tanh(c+d*x)+-1//3*a*d^-1*tanh(c+d*x)^3+1//5*b*d^-1*tanh(c+d*x)^5)
@test integrate(tanh(c+d*x)^3*(a+b*sech(c+d*x)^2), x) == :(a*d^-1*log(cosh(c+d*x))+(1/2)*d^-1*sech(c+d*x)^2*(a+-1b)+1//4*b*d^-1*sech(c+d*x)^4)
@test integrate(tanh(c+d*x)^2*(a+b*sech(c+d*x)^2), x) == :(a*x+-1*a*d^-1*tanh(c+d*x)+1//3*b*d^-1*tanh(c+d*x)^3)
@test integrate((a+b*sech(c+d*x)^2)*tanh(c+d*x), x) == :(a*d^-1*log(cosh(c+d*x))+-1//2*b*d^-1*sech(c+d*x)^2)
@test integrate(a+b*sech(c+d*x)^2, x) == :(a*x+b*d^-1*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)*coth(c+d*x), x) == :(d^-1*(a+b)*log(sinh(c+d*x))+-1*b*d^-1*log(cosh(c+d*x)))
@test integrate(coth(c+d*x)^2*(a+b*sech(c+d*x)^2), x) == :(a*x+-1*d^-1*(a+b)*coth(c+d*x))
@test integrate(coth(c+d*x)^3*(a+b*sech(c+d*x)^2), x) == :(a*d^-1*log(sinh(c+d*x))+-1//2*d^-1*csch(c+d*x)^2*(a+b))
@test integrate(coth(c+d*x)^4*(a+b*sech(c+d*x)^2), x) == :(a*x+-1*a*d^-1*coth(c+d*x)+-1//3*d^-1*coth(c+d*x)^3*(a+b))
@test integrate(coth(c+d*x)^5*(a+b*sech(c+d*x)^2), x) == :(a*d^-1*log(sinh(c+d*x))+-1//2*d^-1*csch(c+d*x)^2*(b+2a)+-1//4*d^-1*csch(c+d*x)^4*(a+b))
@test integrate((a+b*sech(c+d*x)^2)^2*tanh(c+d*x)^4, x) == :(x*a^2+-1*a^2*d^-1*tanh(c+d*x)+-1//3*a^2*d^-1*tanh(c+d*x)^3+-1//7*b^2*d^-1*tanh(c+d*x)^7+1//5*b*d^-1*tanh(c+d*x)^5*(b+2a))
@test integrate((a+b*sech(c+d*x)^2)^2*tanh(c+d*x)^3, x) == :(a^2*d^-1*log(cosh(c+d*x))+1//6*b^2*d^-1*sech(c+d*x)^6+(1/2)*a*d^-1*sech(c+d*x)^2*(a+-2b)+1//4*b*d^-1*sech(c+d*x)^4*(-1b+2a))
@test integrate((a+b*sech(c+d*x)^2)^2*tanh(c+d*x)^2, x) == :(x*a^2+-1*a^2*d^-1*tanh(c+d*x)+-1//5*b^2*d^-1*tanh(c+d*x)^5+1//3*b*d^-1*tanh(c+d*x)^3*(b+2a))
@test integrate((a+b*sech(c+d*x)^2)^2*tanh(c+d*x), x) == :(a^2*d^-1*log(cosh(c+d*x))+-1//4*b^2*d^-1*sech(c+d*x)^4+-1*a*b*d^-1*sech(c+d*x)^2)
@test integrate((a+b*sech(c+d*x)^2)^2, x) == :(x*a^2+-1//3*b^2*d^-1*tanh(c+d*x)^3+b*d^-1*(b+2a)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*coth(c+d*x), x) == :(d^-1*(a+b)^2*log(sinh(c+d*x))+(1/2)*b^2*d^-1*sech(c+d*x)^2+-1*b*d^-1*(b+2a)*log(cosh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^2*coth(c+d*x)^2, x) == :(x*a^2+-1*b^2*d^-1*tanh(c+d*x)+-1*d^-1*(a+b)^2*coth(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^2*coth(c+d*x)^3, x) == :(b^2*d^-1*log(cosh(c+d*x))+d^-1*(a^2+-1*b^2)*log(sinh(c+d*x))+-1//2*d^-1*(a+b)^2*csch(c+d*x)^2)
@test integrate((a+b*sech(c+d*x)^2)^2*coth(c+d*x)^4, x) == :(x*a^2+-1*d^-1*(a^2+-1*b^2)*coth(c+d*x)+-1//3*d^-1*(a+b)^2*coth(c+d*x)^3)
@test integrate((a+b*sech(c+d*x)^2)^2*coth(c+d*x)^5, x) == :(a^2*d^-1*log(sinh(c+d*x))+-1//4*d^-1*(a+b)^2*csch(c+d*x)^4+-1*a*d^-1*csch(c+d*x)^2*(a+b))
@test integrate((a+b*sech(c+d*x)^2)^2*coth(c+d*x)^6, x) == :(x*a^2+-1*a^2*d^-1*coth(c+d*x)+-1//3*d^-1*coth(c+d*x)^3*(a^2+-1*b^2)+-1//5*d^-1*(a+b)^2*coth(c+d*x)^5)
@test integrate((a+b*sech(c+d*x)^2)^2*coth(c+d*x)^7, x) == :(a^2*d^-1*log(sinh(c+d*x))+-1//4*d^-1*(a+b)^2*csch(c+d*x)^4+-1*a*d^-1*csch(c+d*x)^2*(a+b)+-1*d^-1*(b+a*cosh(c+d*x)^2)^3*(6a+6b)^-1*csch(c+d*x)^6)
@test integrate((a+b*sech(c+d*x)^2)^3*tanh(c+d*x)^4, x) == :(x*a^3+-1*a^3*d^-1*tanh(c+d*x)+-1//3*a^3*d^-1*tanh(c+d*x)^3+1//9*b^3*d^-1*tanh(c+d*x)^9+-1//7*b^2*d^-1*tanh(c+d*x)^7*(2b+3a)+1//5*b*d^-1*tanh(c+d*x)^5*(b^2+3*a^2+3*a*b))
@test integrate((a+b*sech(c+d*x)^2)^3*tanh(c+d*x)^3, x) == :(a^3*d^-1*log(cosh(c+d*x))+-1//6*b^3*d^-1*sech(c+d*x)^6+-3//2*b*a^2*d^-1*sech(c+d*x)^2+-3//4*a*b^2*d^-1*sech(c+d*x)^4+1//8*b^-1*d^-1*(b+a*cosh(c+d*x)^2)^4*sech(c+d*x)^8)
@test integrate((a+b*sech(c+d*x)^2)^3*tanh(c+d*x)^2, x) == :(x*a^3+-1*a^3*d^-1*tanh(c+d*x)+1//7*b^3*d^-1*tanh(c+d*x)^7+-1//5*b^2*d^-1*tanh(c+d*x)^5*(2b+3a)+1//3*b*d^-1*tanh(c+d*x)^3*(b^2+3*a^2+3*a*b))
@test integrate((a+b*sech(c+d*x)^2)^3*tanh(c+d*x), x) == :(a^3*d^-1*log(cosh(c+d*x))+-1//6*b^3*d^-1*sech(c+d*x)^6+-3//2*b*a^2*d^-1*sech(c+d*x)^2+-3//4*a*b^2*d^-1*sech(c+d*x)^4)
@test integrate((a+b*sech(c+d*x)^2)^3, x) == :(x*a^3+1//5*b^3*d^-1*tanh(c+d*x)^5+b*d^-1*(b^2+3*a^2+3*a*b)*tanh(c+d*x)+-1//3*b^2*d^-1*tanh(c+d*x)^3*(2b+3a))
@test integrate((a+b*sech(c+d*x)^2)^3*coth(c+d*x), x) == :(d^-1*(a+b)^3*log(sinh(c+d*x))+1//4*b^3*d^-1*sech(c+d*x)^4+(1/2)*b^2*d^-1*sech(c+d*x)^2*(b+3a)+-1*b*d^-1*(b^2+3*a^2+3*a*b)*log(cosh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^3*coth(c+d*x)^2, x) == :(x*a^3+-1*d^-1*(a+b)^3*coth(c+d*x)+1//3*b^3*d^-1*tanh(c+d*x)^3+-1*b^2*d^-1*(2b+3a)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^3*coth(c+d*x)^3, x) == :(-1//2*b^3*d^-1*sech(c+d*x)^2+-1//2*d^-1*(a+b)^3*csch(c+d*x)^2+b^2*d^-1*(2b+3a)*log(cosh(c+d*x))+d^-1*(a+b)^2*(a+-2b)*log(sinh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^3*coth(c+d*x)^4, x) == :(x*a^3+b^3*d^-1*tanh(c+d*x)+-1//3*d^-1*(a+b)^3*coth(c+d*x)^3+-1*d^-1*(a+b)^2*(a+-2b)*coth(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^3*coth(c+d*x)^5, x) == :(d^-1*(a^3+b^3)*log(sinh(c+d*x))+-1*b^3*d^-1*log(cosh(c+d*x))+-1//4*d^-1*(a+b)^3*csch(c+d*x)^4+-1//2*d^-1*(a+b)^2*csch(c+d*x)^2*(-1b+2a))
@test integrate((a+b*sech(c+d*x)^2)^3*coth(c+d*x)^6, x) == :(x*a^3+-1*d^-1*(a^3+b^3)*coth(c+d*x)+-1//5*d^-1*(a+b)^3*coth(c+d*x)^5+-1//3*d^-1*(a+b)^2*coth(c+d*x)^3*(a+-2b))
@test integrate((a+b*sech(c+d*x)^2)^3*coth(c+d*x)^7, x) == :(a^3*d^-1*log(sinh(c+d*x))+-1//6*d^-1*(a+b)^3*csch(c+d*x)^6+-3//2*a^2*d^-1*csch(c+d*x)^2*(a+b)+-3//4*a*d^-1*(a+b)^2*csch(c+d*x)^4)
@test integrate((a+b*sech(c+d*x)^2)^4, x) == :(x*a^4+-1//7*b^4*d^-1*tanh(c+d*x)^7+-1//3*b^2*d^-1*tanh(c+d*x)^3*(3*b^2+6*a^2+8*a*b)+1//5*b^3*d^-1*tanh(c+d*x)^5*(3b+4a)+b*d^-1*(b+2a)*(b^2+2*a^2+2*a*b)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^5, x) == :(x*a^5+1//9*b^5*d^-1*tanh(c+d*x)^9+b*d^-1*(b^4+5*a^4+5*a*b^3+10*b*a^3+10*a^2*b^2)*tanh(c+d*x)+-1//3*b^2*d^-1*tanh(c+d*x)^3*(4*b^3+10*a^3+15*a*b^2+20*b*a^2)+-1//7*b^4*d^-1*tanh(c+d*x)^7*(4b+5a)+1//5*b^3*d^-1*tanh(c+d*x)^5*(6*b^2+10*a^2+15*a*b))
@test integrate((a+b*sech(c+d*x)^2)^-1*tanh(c+d*x)^5, x) == :(-1//2*b^-1*d^-1*sech(c+d*x)^2+-1*b^-2*d^-1*(a+2b)*log(cosh(c+d*x))+(1/2)*a^-1*b^-2*d^-1*(a+b)^2*log(b+a*cosh(c+d*x)^2))
@test integrate((a+b*sech(c+d*x)^2)^-1*tanh(c+d*x)^4, x) == :(x*a^-1+b^-1*d^-1*tanh(c+d*x)+-1*a^-1*b^-3//2*d^-1*(a+b)^3//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*tanh(c+d*x)^3, x) == :(-1*b^-1*d^-1*log(cosh(c+d*x))+(1/2)*a^-1*b^-1*d^-1*(a+b)*log(b+a*cosh(c+d*x)^2))
@test integrate((a+b*sech(c+d*x)^2)^-1*tanh(c+d*x)^2, x) == :(x*a^-1+-1*a^-1*b^-1//2*d^-1*(a+b)^(1/2)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*tanh(c+d*x), x) == :((1/2)*a^-1*d^-1*log(b+a*cosh(c+d*x)^2))
@test integrate((a+b*sech(c+d*x)^2)^-1, x) == :(x*a^-1+-1*a^-1*b^(1/2)*d^-1*(a+b)^-1//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*coth(c+d*x), x) == :(d^-1*(a+b)^-1*log(sinh(c+d*x))+(1/2)*b*a^-1*d^-1*(a+b)^-1*log(b+a*cosh(c+d*x)^2))
@test integrate((a+b*sech(c+d*x)^2)^-1*coth(c+d*x)^2, x) == :(x*a^-1+-1*d^-1*(a+b)^-1*coth(c+d*x)+-1*a^-1*b^3//2*d^-1*(a+b)^-3//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-1*coth(c+d*x)^3, x) == :(-1*d^-1*(2a+2b)^-1*csch(c+d*x)^2+d^-1*(a+b)^-2*(a+2b)*log(sinh(c+d*x))+(1/2)*a^-1*b^2*d^-1*(a+b)^-2*log(b+a*cosh(c+d*x)^2))
@test integrate((a+b*sech(c+d*x)^2)^-1*coth(c+d*x)^4, x) == :(x*a^-1+-1*d^-1*(3a+3b)^-1*coth(c+d*x)^3+-1*d^-1*(a+b)^-2*(a+2b)*coth(c+d*x)+-1*a^-1*b^5//2*d^-1*(a+b)^-5//2*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-2*tanh(c+d*x)^5, x) == :(b^-2*d^-1*log(cosh(c+d*x))+(1/2)*d^-1*(a^-2+-1*b^-2)*log(b+a*cosh(c+d*x)^2)+(1/2)*a^-2*b^-1*d^-1*(a+b)^2*(b+a*cosh(c+d*x)^2)^-1)
@test integrate((a+b*sech(c+d*x)^2)^-2*tanh(c+d*x)^4, x) == :(x*a^-2+(1/2)*a^-2*b^-3//2*d^-1*(a+b)^(1/2)*(a+-2b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+-1//2*a^-1*b^-1*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*(a+b)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-2*tanh(c+d*x)^3, x) == :((1/2)*a^-2*d^-1*log(b+a*cosh(c+d*x)^2)+(1/2)*a^-2*d^-1*(b+a*cosh(c+d*x)^2)^-1*(a+b))
@test integrate((a+b*sech(c+d*x)^2)^-2*tanh(c+d*x)^2, x) == :(x*a^-2+-1//2*a^-1*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*tanh(c+d*x)+-1//2*a^-2*b^-1//2*d^-1*(a+b)^-1//2*(a+2b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-2*tanh(c+d*x), x) == :((1/2)*a^-2*d^-1*log(b+a*cosh(c+d*x)^2)+(1/2)*b*a^-2*d^-1*(b+a*cosh(c+d*x)^2)^-1)
@test integrate((a+b*sech(c+d*x)^2)^-2, x) == :(x*a^-2+-1//2*b*a^-1*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*tanh(c+d*x)+-1//2*a^-2*b^(1/2)*d^-1*(a+b)^-3//2*(2b+3a)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-2*coth(c+d*x), x) == :(d^-1*(a+b)^-2*log(sinh(c+d*x))+(1/2)*a^-2*b^2*d^-1*(a+b)^-1*(b+a*cosh(c+d*x)^2)^-1+(1/2)*b*a^-2*d^-1*(a+b)^-2*(b+2a)*log(b+a*cosh(c+d*x)^2))
@test integrate((a+b*sech(c+d*x)^2)^-2*coth(c+d*x)^2, x) == :(x*a^-2+-1//2*a^-1*d^-1*(a+b)^-2*(-1b+2a)*coth(c+d*x)+-1//2*b*a^-1*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*coth(c+d*x)+-1//2*a^-2*b^3//2*d^-1*(a+b)^-5//2*(2b+5a)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-2*coth(c+d*x)^3, x) == :(-1//2*d^-1*(a+b)^-2*csch(c+d*x)^2+d^-1*(a+b)^-3*(a+3b)*log(sinh(c+d*x))+(1/2)*a^-2*b^3*d^-1*(a+b)^-2*(b+a*cosh(c+d*x)^2)^-1+(1/2)*a^-2*b^2*d^-1*(a+b)^-3*(b+3a)*log(b+a*cosh(c+d*x)^2))
@test integrate((a+b*sech(c+d*x)^2)^-2*coth(c+d*x)^4, x) == :(x*a^-2+-1//2*a^-1*d^-1*(a+b)^-3*(-1*b^2+2*a^2+6*a*b)*coth(c+d*x)+-1//6*a^-1*d^-1*(a+b)^-2*coth(c+d*x)^3*(-3b+2a)+-1//2*b*a^-1*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*coth(c+d*x)^3+-1//2*a^-2*b^5//2*d^-1*(a+b)^-7//2*(2b+7a)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x)))
@test integrate((a+b*sech(c+d*x)^2)^-3*tanh(c+d*x)^6, x) == :(x*a^-3+-1//4*a^-1*b^-1*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*tanh(c+d*x)^3*(a+b)+-1//8*a^-3*b^-5//2*d^-1*(a+b)^(1/2)*(3*a^2+8*b^2+-4*a*b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+1//8*a^-2*b^-2*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*(a+b)*(-4b+3a)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*tanh(c+d*x)^5, x) == :((1/2)*a^-3*d^-1*log(b+a*cosh(c+d*x)^2)+a^-3*d^-1*(b+a*cosh(c+d*x)^2)^-1*(a+b)+-1//4*a^-3*d^-1*(a+b)^2*(b+a*cosh(c+d*x)^2)^-2)
@test integrate((a+b*sech(c+d*x)^2)^-3*tanh(c+d*x)^4, x) == :(x*a^-3+-1//4*a^-1*b^-1*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*(a+b)*tanh(c+d*x)+1//8*a^-3*b^-3//2*d^-1*(a+b)^-1//2*(a^2+-8*b^2+-4*a*b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+1//8*a^-2*b^-1*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*(a+-4b)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*tanh(c+d*x)^3, x) == :((1/2)*a^-3*d^-1*log(b+a*cosh(c+d*x)^2)+(1/2)*a^-3*d^-1*(b+a*cosh(c+d*x)^2)^-1*(a+2b)+-1//4*b*a^-3*d^-1*(b+a*cosh(c+d*x)^2)^-2*(a+b))
@test integrate((a+b*sech(c+d*x)^2)^-3*tanh(c+d*x)^2, x) == :(x*a^-3+-1//4*a^-1*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*tanh(c+d*x)+-1//8*a^-3*b^-1//2*d^-1*(a+b)^-3//2*(3*a^2+8*b^2+12*a*b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+-1//8*a^-2*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-1*(3a+4b)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*tanh(c+d*x), x) == :((1/2)*a^-3*d^-1*log(b+a*cosh(c+d*x)^2)+b*a^-3*d^-1*(b+a*cosh(c+d*x)^2)^-1+-1//4*a^-3*b^2*d^-1*(b+a*cosh(c+d*x)^2)^-2)
@test integrate((a+b*sech(c+d*x)^2)^-3, x) == :(x*a^-3+-1//4*b*a^-1*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*tanh(c+d*x)+-1//8*a^-3*b^(1/2)*d^-1*(a+b)^-5//2*(8*b^2+15*a^2+20*a*b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+-1//8*b*a^-2*d^-1*(a+b)^-2*(a+b+-1*b*tanh(c+d*x)^2)^-1*(4b+7a)*tanh(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*coth(c+d*x), x) == :(d^-1*(a+b)^-3*log(sinh(c+d*x))+-1//4*a^-3*b^3*d^-1*(a+b)^-1*(b+a*cosh(c+d*x)^2)^-2+(1/2)*b*a^-3*d^-1*(a+b)^-3*(b^2+3*a^2+3*a*b)*log(b+a*cosh(c+d*x)^2)+(1/2)*a^-3*b^2*d^-1*(a+b)^-2*(b+a*cosh(c+d*x)^2)^-1*(2b+3a))
@test integrate((a+b*sech(c+d*x)^2)^-3*coth(c+d*x)^2, x) == :(x*a^-3+-1//8*a^-2*d^-1*(a+b)^-3*(-4*b^2+8*a^2+-11*a*b)*coth(c+d*x)+-1//4*b*a^-1*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*coth(c+d*x)+-1//8*a^-3*b^3//2*d^-1*(a+b)^-7//2*(8*b^2+35*a^2+28*a*b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+-1//8*b*a^-2*d^-1*(a+b)^-2*(a+b+-1*b*tanh(c+d*x)^2)^-1*(4b+9a)*coth(c+d*x))
@test integrate((a+b*sech(c+d*x)^2)^-3*coth(c+d*x)^3, x) == :(-1//2*d^-1*(a+b)^-3*csch(c+d*x)^2+d^-1*(a+b)^-4*(a+4b)*log(sinh(c+d*x))+-1//4*a^-3*b^4*d^-1*(a+b)^-2*(b+a*cosh(c+d*x)^2)^-2+a^-3*b^3*d^-1*(a+b)^-3*(b+a*cosh(c+d*x)^2)^-1*(b+2a)+(1/2)*a^-3*b^2*d^-1*(a+b)^-4*(b^2+6*a^2+4*a*b)*log(b+a*cosh(c+d*x)^2))
@test integrate((a+b*sech(c+d*x)^2)^-3*coth(c+d*x)^4, x) == :(x*a^-3+-1//8*a^-2*d^-1*(a+b)^-4*(-4*b^3+8*a^3+-15*a*b^2+32*b*a^2)*coth(c+d*x)+-1//24*a^-2*d^-1*(a+b)^-3*coth(c+d*x)^3*(-12*b^2+8*a^2+-39*a*b)+-1//4*b*a^-1*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-2*coth(c+d*x)^3+-1//8*a^-3*b^5//2*d^-1*(a+b)^-9//2*(8*b^2+63*a^2+36*a*b)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+-1//8*b*a^-2*d^-1*(a+b)^-2*(a+b+-1*b*tanh(c+d*x)^2)^-1*coth(c+d*x)^3*(4b+11a))
@test integrate((a+b*sech(c+d*x)^2)^-4, x) == :(x*a^-4+-1//6*b*a^-1*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-3*tanh(c+d*x)+-1//16*a^-4*b^(1/2)*d^-1*(a+b)^-7//2*(16*b^3+35*a^3+56*a*b^2+70*b*a^2)*arctanh(b^(1/2)*(a+b)^-1//2*tanh(c+d*x))+-1//16*b*a^-3*d^-1*(a+b)^-3*(a+b+-1*b*tanh(c+d*x)^2)^-1*(8*b^2+19*a^2+22*a*b)*tanh(c+d*x)+-1//24*b*a^-2*d^-1*(a+b)^-2*(a+b+-1*b*tanh(c+d*x)^2)^-2*(6b+11a)*tanh(c+d*x))
@test integrate((1+-1*sech(x)^2)^3//2, x) == :(-1//2*(tanh(x)^2)^3//2*coth(x)+(tanh(x)^2)^(1/2)*coth(x)*log(cosh(x)))
@test integrate((1+-1*sech(x)^2)^(1/2), x) == :((tanh(x)^2)^(1/2)*coth(x)*log(cosh(x)))
@test integrate((1+-1*sech(x)^2)^-1//2, x) == :((tanh(x)^2)^-1//2*log(sinh(x))*tanh(x))
@test integrate((-1+sech(x)^2)^3//2, x) == :((1/2)*(-1*tanh(x)^2)^(1/2)*tanh(x)+-1*(-1*tanh(x)^2)^(1/2)*coth(x)*log(cosh(x)))
@test integrate((-1+sech(x)^2)^(1/2), x) == :((-1*tanh(x)^2)^(1/2)*coth(x)*log(cosh(x)))
@test integrate((-1+sech(x)^2)^-1//2, x) == :((-1*tanh(x)^2)^-1//2*log(sinh(x))*tanh(x))
@test integrate((a+b*sech(x)^2)^(1/2)*tanh(x)^5, x) == :(-1*(a+b*sech(x)^2)^(1/2)+a^(1/2)*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1//5*b^-2*(a+b*sech(x)^2)^5//2+1//3*b^-2*(a+b*sech(x)^2)^3//2*(a+2b))
@test integrate((a+b*sech(x)^2)^(1/2)*tanh(x)^4, x) == :(a^(1/2)*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1//4*(a+b+-1*b*tanh(x)^2)^(1/2)*tanh(x)^3+-1//8*b^-3//2*(a^2+-3*b^2+6*a*b)*arctan(b^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+1//8*b^-1*(a+b+-1*b*tanh(x)^2)^(1/2)*(a+-3b)*tanh(x))
@test integrate((a+b*sech(x)^2)^(1/2)*tanh(x)^3, x) == :(-1*(a+b*sech(x)^2)^(1/2)+a^(1/2)*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+1//3*b^-1*(a+b*sech(x)^2)^3//2)
@test integrate((a+b*sech(x)^2)^(1/2)*tanh(x)^2, x) == :(a^(1/2)*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1//2*(a+b+-1*b*tanh(x)^2)^(1/2)*tanh(x)+-1//2*b^-1//2*(a+-1b)*arctan(b^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x)))
@test integrate((a+b*sech(x)^2)^(1/2)*tanh(x), x) == :(-1*(a+b*sech(x)^2)^(1/2)+a^(1/2)*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2)))
@test integrate((a+b*sech(x)^2)^(1/2), x) == :(a^(1/2)*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+b^(1/2)*arctan(b^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x)))
@test integrate((a+b*sech(x)^2)^(1/2)*coth(x), x) == :(a^(1/2)*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1*(a+b)^(1/2)*arctanh((a+b)^-1//2*(a+b*sech(x)^2)^(1/2)))
@test integrate((a+b*sech(x)^2)^(1/2)*coth(x)^2, x) == :(a^(1/2)*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1*(a+b+-1*b*tanh(x)^2)^(1/2)*coth(x))
@test integrate((a+b*sech(x)^2)^(1/2)*coth(x)^3, x) == :(a^(1/2)*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1//2*(a+b*sech(x)^2)^(1/2)*coth(x)^2+-1//2*(a+b)^-1//2*(b+2a)*arctanh((a+b)^-1//2*(a+b*sech(x)^2)^(1/2)))
@test integrate((a+b*sech(x)^2)^(1/2)*coth(x)^4, x) == :(a^(1/2)*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1//3*(a+b+-1*b*tanh(x)^2)^(1/2)*coth(x)^3+-1*(3a+3b)^-1*(a+b+-1*b*tanh(x)^2)^(1/2)*(2b+3a)*coth(x))
@test integrate((a+b*sech(x)^2)^(1/2)*coth(x)^5, x) == :(a^(1/2)*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1//4*(a+b*sech(x)^2)^(1/2)*coth(x)^4+-1//8*(a+b)^-3//2*(3*b^2+8*a^2+12*a*b)*arctanh((a+b)^-1//2*(a+b*sech(x)^2)^(1/2))+-1*(a+b*sech(x)^2)^(1/2)*(8a+8b)^-1*coth(x)^2*(3b+4a))
@test integrate((a+b*sech(x)^2)^3//2*tanh(x)^3, x) == :(-1//3*(a+b*sech(x)^2)^3//2+a^3//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1*a*(a+b*sech(x)^2)^(1/2)+1//5*b^-1*(a+b*sech(x)^2)^5//2)
@test integrate((a+b*sech(x)^2)^3//2*tanh(x)^2, x) == :(a^3//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1//8*b^-1//2*(-1*b^2+3*a^2+-6*a*b)*arctan(b^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1//8*(a+b+-1*b*tanh(x)^2)^(1/2)*(b+5a)*tanh(x)+1//4*b*(a+b+-1*b*tanh(x)^2)^(1/2)*tanh(x)^3)
@test integrate((a+b*sech(x)^2)^3//2*tanh(x), x) == :(-1//3*(a+b*sech(x)^2)^3//2+a^3//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1*a*(a+b*sech(x)^2)^(1/2))
@test integrate((a+b*sech(x)^2)^3//2, x) == :(a^3//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+(1/2)*b*(a+b+-1*b*tanh(x)^2)^(1/2)*tanh(x)+(1/2)*b^(1/2)*(b+3a)*arctan(b^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x)))
@test integrate((a+b*sech(x)^2)^3//2*coth(x), x) == :(b*(a+b*sech(x)^2)^(1/2)+a^3//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1*(a+b)^3//2*arctanh((a+b)^-1//2*(a+b*sech(x)^2)^(1/2)))
@test integrate((a+b*sech(x)^2)^3//2*coth(x)^2, x) == :(a^3//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1*b^3//2*arctan(b^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1*(a+b+-1*b*tanh(x)^2)^(1/2)*(a+b)*coth(x))
@test integrate((a+b*sech(c+d*x)^2)^5//2, x) == :(a^5//2*d^-1*arctanh(a^(1/2)*(a+b+-1*b*tanh(c+d*x)^2)^-1//2*tanh(c+d*x))+1//4*b*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^3//2*tanh(c+d*x)+1//8*b^(1/2)*d^-1*(3*b^2+15*a^2+10*a*b)*arctan(b^(1/2)*(a+b+-1*b*tanh(c+d*x)^2)^-1//2*tanh(c+d*x))+1//8*b*d^-1*(a+b+-1*b*tanh(c+d*x)^2)^(1/2)*(3b+7a)*tanh(c+d*x))
@test integrate((a+b*sech(x)^2)^-1//2*tanh(x)^5, x) == :(a^-1//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1//3*b^-2*(a+b*sech(x)^2)^3//2+b^-2*(a+b*sech(x)^2)^(1/2)*(a+2b))
@test integrate((a+b*sech(x)^2)^-1//2*tanh(x)^4, x) == :(a^-1//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+(1/2)*b^-1*(a+b+-1*b*tanh(x)^2)^(1/2)*tanh(x)+-1//2*b^-3//2*(a+3b)*arctan(b^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x)))
@test integrate((a+b*sech(x)^2)^-1//2*tanh(x)^3, x) == :(a^-1//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+b^-1*(a+b*sech(x)^2)^(1/2))
@test integrate((a+b*sech(x)^2)^-1//2*tanh(x)^2, x) == :(a^-1//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1*b^-1//2*arctan(b^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x)))
@test integrate((a+b*sech(x)^2)^-1//2*tanh(x), x) == :(a^-1//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2)))
@test integrate((a+b*sech(x)^2)^-1//2, x) == :(a^-1//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x)))
@test integrate((a+b*sech(x)^2)^-1//2*coth(x), x) == :(a^-1//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1*(a+b)^-1//2*arctanh((a+b)^-1//2*(a+b*sech(x)^2)^(1/2)))
@test integrate((a+b*sech(x)^2)^-1//2*coth(x)^2, x) == :(a^-1//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1*(a+b)^-1*(a+b+-1*b*tanh(x)^2)^(1/2)*coth(x))
@test integrate((a+b*sech(x)^2)^-1//2*coth(x)^3, x) == :(a^-1//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1*(a+b*sech(x)^2)^(1/2)*(2a+2b)^-1*coth(x)^2+-1//2*(a+b)^-3//2*(2a+3b)*arctanh((a+b)^-1//2*(a+b*sech(x)^2)^(1/2)))
@test integrate((a+b*sech(x)^2)^-3//2*tanh(x)^5, x) == :(a^-3//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1*b^-2*(a+b*sech(x)^2)^(1/2)+-1*a^-1*b^-2*(a+b)^2*(a+b*sech(x)^2)^-1//2)
@test integrate((a+b*sech(x)^2)^-3//2*tanh(x)^4, x) == :(a^-3//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+b^-3//2*arctan(b^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1*a^-1*b^-1*(a+b+-1*b*tanh(x)^2)^-1//2*(a+b)*tanh(x))
@test integrate((a+b*sech(x)^2)^-3//2*tanh(x)^3, x) == :(a^-3//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1*a^-1*b^-1*(a+b*sech(x)^2)^-1//2*(a+b))
@test integrate((a+b*sech(x)^2)^-3//2*tanh(x)^2, x) == :(a^-3//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1*a^-1*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))
@test integrate((a+b*sech(x)^2)^-3//2*tanh(x), x) == :(a^-3//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1*a^-1*(a+b*sech(x)^2)^-1//2)
@test integrate((a+b*sech(x)^2)^-3//2, x) == :(a^-3//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1*b*a^-1*(a+b)^-1*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))
@test integrate((a+b*sech(x)^2)^-3//2*coth(x), x) == :(a^-3//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1*(a+b)^-3//2*arctanh((a+b)^-1//2*(a+b*sech(x)^2)^(1/2))+-1*b*a^-1*(a+b)^-1*(a+b*sech(x)^2)^-1//2)
@test integrate((a+b*sech(x)^2)^-3//2*coth(x)^2, x) == :(a^-3//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1*b*a^-1*(a+b)^-1*(a+b+-1*b*tanh(x)^2)^-1//2*coth(x)+-1*a^-1*(a+b)^-2*(a+b+-1*b*tanh(x)^2)^(1/2)*(a+-1b)*coth(x))
@test integrate((a+b*sech(x)^2)^-5//2*tanh(x)^6, x) == :(a^-5//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1*b^-5//2*arctan(b^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1*(a+b+-1*b*tanh(x)^2)^-1//2*(a^-2+-1*b^-2)*tanh(x)+-1//3*a^-1*b^-1*(a+b+-1*b*tanh(x)^2)^-3//2*tanh(x)^3*(a+b))
@test integrate((a+b*sech(x)^2)^-5//2*tanh(x)^5, x) == :(a^-5//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1*(a+b*sech(x)^2)^-1//2*(a^-2+-1*b^-2)+-1//3*a^-1*b^-2*(a+b)^2*(a+b*sech(x)^2)^-3//2)
@test integrate((a+b*sech(x)^2)^-5//2*tanh(x)^4, x) == :(a^-5//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1//3*a^-1*b^-1*(a+b+-1*b*tanh(x)^2)^-3//2*(a+b)*tanh(x)+1//3*a^-2*b^-1*(a+b+-1*b*tanh(x)^2)^-1//2*(a+-3b)*tanh(x))
@test integrate((a+b*sech(x)^2)^-5//2*tanh(x)^3, x) == :(a^-5//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1*a^-2*(a+b*sech(x)^2)^-1//2+-1//3*a^-1*b^-1*(a+b*sech(x)^2)^-3//2*(a+b))
@test integrate((a+b*sech(x)^2)^-5//2*tanh(x)^2, x) == :(a^-5//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1//3*a^-1*(a+b+-1*b*tanh(x)^2)^-3//2*tanh(x)+-1//3*a^-2*(a+b)^-1*(a+b+-1*b*tanh(x)^2)^-1//2*(2a+3b)*tanh(x))
@test integrate((a+b*sech(x)^2)^-5//2*tanh(x), x) == :(a^-5//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1*a^-2*(a+b*sech(x)^2)^-1//2+-1//3*a^-1*(a+b*sech(x)^2)^-3//2)
@test integrate((a+b*sech(x)^2)^-5//2, x) == :(a^-5//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1//3*b*a^-1*(a+b)^-1*(a+b+-1*b*tanh(x)^2)^-3//2*tanh(x)+-1//3*b*a^-2*(a+b)^-2*(a+b+-1*b*tanh(x)^2)^-1//2*(3b+5a)*tanh(x))
@test integrate((a+b*sech(x)^2)^-5//2*coth(x), x) == :(a^-5//2*arctanh(a^-1//2*(a+b*sech(x)^2)^(1/2))+-1*(a+b)^-5//2*arctanh((a+b)^-1//2*(a+b*sech(x)^2)^(1/2))+-1//3*b*a^-1*(a+b)^-1*(a+b*sech(x)^2)^-3//2+-1*b*a^-2*(a+b)^-2*(a+b*sech(x)^2)^-1//2*(b+2a))
@test integrate((a+b*sech(x)^2)^-5//2*coth(x)^2, x) == :(a^-5//2*arctanh(a^(1/2)*(a+b+-1*b*tanh(x)^2)^-1//2*tanh(x))+-1//3*b*a^-1*(a+b)^-1*(a+b+-1*b*tanh(x)^2)^-3//2*coth(x)+-1//3*b*a^-2*(a+b)^-2*(a+b+-1*b*tanh(x)^2)^-1//2*(3b+7a)*coth(x)+-1//3*a^-2*(a+b)^-3*(a+b+-1*b*tanh(x)^2)^(1/2)*(a+-3b)*(b+3a)*coth(x))
@test integrate((a+b*sech(c+d*x)^2)^-7//2, x) == :(a^-7//2*d^-1*arctanh(a^(1/2)*(a+b+-1*b*tanh(c+d*x)^2)^-1//2*tanh(c+d*x))+-1//5*b*a^-1*d^-1*(a+b)^-1*(a+b+-1*b*tanh(c+d*x)^2)^-5//2*tanh(c+d*x)+-1//15*b*a^-3*d^-1*(a+b)^-3*(a+b+-1*b*tanh(c+d*x)^2)^-1//2*(15*b^2+33*a^2+40*a*b)*tanh(c+d*x)+-1//15*b*a^-2*d^-1*(a+b)^-2*(a+b+-1*b*tanh(c+d*x)^2)^-3//2*(5b+9a)*tanh(c+d*x))

max_stars_repo_path: test/6 Hyperbolic functions/6.5 Hyperbolic secant/6.5.7 (d hyper)^m (a+b (c sech)^n)^p.jl


text:
using JuliaExampleTest
using Test
@testset "JuliaExampleTest.jl" begin
    # Write your tests here.
    @test subfunction(2,1) == 8
    @test subfunction(2,3) == 12
end

@testset "Derivative Test" begin
    @test derivative_of_my_f(2,1) == 3
end

ext: jl


text:
PS-2

using Optim,DataFrames, CSV, HTTP,GLM,FreqTables

#Question 1

f(x) = -x[1]^4-10x[1]^3-2x[1]^2-3x[1]-2
	minusf(x) = x[1]^4+10x[1]^3+2x[1]^2+3x[1]+2
	startval = rand(1)   # random starting value
	result = optimize(minusf, startval, BFGS())
	println(result)

    #Question 2

    url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS1-julia-intro/nlsw88.csv"
	df = CSV.read(HTTP.get(url).body, DataFrame)
	X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
	y = df.married.==1

	function ols(beta, X, y)
		ssr = (y.-X*beta)'*(y.-X*beta)
		return ssr
	end

	beta_hat_ols = optimize(b -> ols(b, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
	println(beta_hat_ols.minimizer)

	bols = inv(X'*X)*X'*y
	df.white = df.race.==1
	bols_lm = lm(@formula(married ~ age + white + collgrad), df)

    #Question 3

    function logit(alpha, X, d)
		log_mat=[d[i]*log((1+exp(-(X*alpha)[i]))^-1)+(1-d[i])*log(1-(1+exp(-(X*alpha)[i]))^-1) for i in 1:size(d,1)]
	#     p=(1+exp(-X*alpha))^-1
		loglike=-sum(log_mat)

		return loglike
	end
	beta_loglike = optimize(b -> logit(b, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
	println(beta_loglike.minimizer)

    #Question 4

    df.white = df.race.==1
	α̂ = glm(@formula(married ~ age + white + collgrad), df, Binomial(), LogitLink())
	println(α̂)

    #Question 5

    	freqtable(df, :occupation) # note small number of obs in some occupations
	df = dropmissing(df, :occupation)
	df[df.occupation.>7,:occupation] .= 7
	freqtable(df, :occupation) # problem solved

	X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
	y = df.occupation

	function mnlogit(alpha, x, d)
		c=[i==j ? 1 : 0 for i in d , j in 1:6]
		log_mat=[log(exp.(x[i,:]'*alpha[:,j])/(1+sum(exp.(x[i,:]'*alpha)))) for i in 1:size(x,1) ,j in 1:6]
		loglike=-sum(c.*log_mat)

		return loglike
	end

	beta_mnlog = optimize(b -> mnlogit(b, X, y), rand(size(X,2),6), LBFGS(), Optim.Options(g_tol=1e-7, iterations=100_000, show_trace=true))
	println(beta_mnlog.minimizer)


text:
"""
    AllParticles{APC} <: AbstractParticles where {APC<:AbstractParticleCollection}

Combines the particles of different types from a particle collection lazily.

The property arrays of the collection's particles are concatenated lazily with
[LazyArrays.jl](https://github.com/JuliaArrays/LazyArrays.jl).
Scalar properties (e.g. mass for equal-mass particles) are repeated according to the particle number of the given
particle type with [FillArrays.jl](https://github.com/JuliaArrays/FillArrays.jl).
Properties that only exist for one of the particle types are repeated as `missing` values in the same way as
scalar properties.

# Usage
The individual properties of `AllParticles` can be converted as regular arrays with `Array(ap.pos)`
and all or selected particle properties can be converted to a `Particles` object with `Particles(ap)`,
where all lazy arrays are materialized as `Array`s.
`AllParticles` cannot be sorted, filtered, or copied.
In-place modifications are possible, such as [`rotate!`](@ref), [`translate!`](@ref), or [`to_comoving`](@ref).

!!! warning "In-place modifications"
    Modifying an `AllParticles` object in-place results in the modification of the underlying particle collection!
    To prevent bugs or unexpected behavior, it is advised to modify the particle collection instead and then
    obtain the object through `pc.all` or `AllParticles(pc)` (these are equivalent).

The lazy concatenation occurs only when the property is needed and is only performed once.
Note that this means that calling `pc.all.prop` multiple times results in the given property being lazily
concatenated every time. As long as the property array pointers of the particle collection are not modified
(e.g. by `sort!` or `filter!`), it is possible to reuse an object `ap = pc.all` without the need of regenerating
the concatenated arrays.

!!! warning "Modifying the underlying particle collection"
    After modifying the property array pointers of a particle collection (e.g. by `sort!` or `filter!`), any
    already created `AllParticles` objects from that collection may have lazy arrays pointing to the old,
    non-modified arrays. When this happens, create a new object through `pc.all` or `AllParticles(pc)`.

# Example
```julia
julia> pc::ParticleCollection
ParticleCollection
dm: 100 Particles
 id mass pos
gas: 50 Particles
 id mass pos temp

julia> ap = AllParticles(pc) # equivalent to ap = pc.all
all: 150 Particles
 id mass pos temp

julia> p = Particles(ap) # p is detached from pc and therefore modifiable
all: 150 Particles
 id mass pos temp

julia> sort!(p, :id)
```
"""
struct AllParticles{APC} <: AbstractParticles where {APC<:AbstractParticleCollection}
    particle_collection::APC
    props::Dict{Symbol,Any}

    AllParticles(pc::APC) where {APC} = new{APC}(pc, Dict{Symbol,Any}())
end

"""
    particle_collection(p::AllParticles)

Returns the underlying particle collection.

This is not exported.
"""
particle_collection(p::AllParticles) = p.particle_collection

function Base.getindex(ap::AllParticles, sym::Symbol)
    # get property index if property already saved or if key does not exist at all to raise Dict error
    if haskey(ap.props, sym) || !haskey(ap, sym)
        return ap.props[sym]
    end

    pc = ap.particle_collection

    # get array dimensions
    dims = 0 # 0 for vectors; for arrays: size(arr, 1)
    for key in keys(pc)
        p = pc[key]
        if haskey(p, sym)
            p[sym] isa AbstractVector && break
            if p[sym] isa AbstractMatrix
                dims = size(p[sym], 1)
                break
            end
        end
    end

    # create array of particle property arrays
    arrs = []

    # iterate through particles
    for p in values(pc)
        # fill array with missing if property does not exist
        if haskey(p, sym)
            # fill array with scalar value if property is not an array
            if p[sym] isa AbstractArray
                arr = p[sym]
            else
                if dims == 0
                    arr = Fill(p[sym], particle_number(p))
                else
                    arr = Fill(p[sym], dims, particle_number(p))
                end
            end
            push!(arrs, arr)
        else
            if dims == 0
                arr = Fill(missing, particle_number(p))
            else
                arr = Fill(missing, dims, particle_number(p))
            end
            push!(arrs, arr)
        end
    end

    cat_func = dims == 0 ? vcat : hcat
    lazy = ApplyArray(cat_func, arrs...)

    ap.props[sym] = lazy
    return lazy
end

function Base.setindex!(::AllParticles, _, ::Symbol)
    error(
        "Setting a property of an AllParticles object is not allowed. Try converting it via Particles(allparticles) first.",
    )
end

function Base.keys(ap::AllParticles)
    keys_arr = [keys(p) for p in values(ap.particle_collection)]
    length(keys_arr) == 0 && return Symbol[]
    length(keys_arr) == 1 && return keys_arr[1]
    return union(keys_arr...) |> collect
end

Base.haskey(p::AllParticles, key) = haskey(p.props, key) || key in keys(p)

Base.values(p::AllParticles) = Iterators.map(key -> getindex(p, key), keys(p))

function Base.copy(::AllParticles)
    error(
        "Copying an AllParticles object is not allowed. Try converting it via Particles(allparticles) instead.",
    )
    ## NOTE: potential way of making copy work
    # pnew = Particles(p.particle_collection)
    # copy!(pnew.props, p.props)
    # return pnew
end

function Base.copy!(::AllParticles, ::AllParticles, props=nothing)
    error(
        "Copying particles to an AllParticles object is not allowed. Try converting it via Particles(allparticles) instead.",
    )
end

function Base.copy!(::AllParticles, ::AbstractParticles, props=nothing)
    error(
        "Copying particles to an AllParticles object is not allowed. Try converting it via Particles(allparticles) instead.",
    )
end

"""
    Base.copy!(dst::AbstractParticles, src::AllParticles[, props])

Copies the materialized particle properties to another particle object.

To only copy certain properties, `props` can be passed as a vector of `Symbol`s.
The `dst` cannot be of type `AllParticles`.
"""
function Base.copy!(dst::AbstractParticles, src::AllParticles, props=nothing)
    empty!(dst)
    for key in keys(src)
        if isnothing(props) || key in props
            val = src[key]
            dst[key] = Array(val)
        end
    end
    return dst
end

function Base.empty(::AllParticles)
    error(
        "Creating an empty AllParticles object is not supported. Try emptying the underlying particle collection instead.",
    )
end

function Base.empty!(::AllParticles)
    error(
        "Emptying an AllParticles object is not allowed. Try emptying the underlying particle collection instead.",
    )
end

function Base.isempty(p::AllParticles)
    for p in values(p.particle_collection)
        isempty(p) || return false
    end
    return true
end

Base.:(==)(p1::AllParticles, p2::AllParticles) = p1.particle_collection == p2.particle_collection

particle_name(::AllParticles) = "Particles"

function particle_number(ap::AllParticles)
    n = 0
    for p in values(ap.particle_collection)
        n += particle_number(p)
    end
    return n
end

function Base.propertynames(p::AllParticles; private=false)
    if private
        return [fieldnames(AllParticles) |> collect; keys(p) |> collect]
    else
        return keys(p) |> collect
    end
end

function Base.show(io::IO, mime::MIME"text/plain", p::AllParticles)
    printstyled(io, "all"; bold=true)
    print(io, ": ")
    show_properties(io, mime, p)
end

"""
    Particles(p::AllParticles, props=nothing)

Materializes the particle properties in a new particle object of type `:all`.

To only copy certain properties, `props` can be passed as a tuple of `Symbol`s.
"""
Particles(p::AllParticles, props=nothing) = copy!(Particles(:all), p, props)

function applyind!(::AllParticles, ::AbstractVector)
    error(
        "Indexing into an AllParticles object is not allowed. Try converting it via Particles(allparticles) first.",
    )
end

applyind(p::AllParticles, ind::AbstractVector) = applyind!(p, ind)


text:
function GlobalMappingFunctionModel(arg0::jdouble, arg1::jdouble)
    return GlobalMappingFunctionModel((jdouble, jdouble), arg0, arg1)
end

function GlobalMappingFunctionModel(arg0::jdouble, arg1::jdouble, arg2::TimeScale)
    return GlobalMappingFunctionModel((jdouble, jdouble, TimeScale), arg0, arg1, arg2)
end

function get_parameters_drivers(obj::GlobalMappingFunctionModel)
    return jcall(obj, "getParametersDrivers", List, ())
end

function mapping_factors(obj::GlobalMappingFunctionModel, arg0::RealFieldElement, arg1::RealFieldElement, arg2::Vector{RealFieldElement}, arg3::FieldAbsoluteDate)
    return jcall(obj, "mappingFactors", Vector{RealFieldElement}, (RealFieldElement, RealFieldElement, Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1, arg2, arg3)
end

function mapping_factors(obj::GlobalMappingFunctionModel, arg0::jdouble, arg1::jdouble, arg2::Vector{jdouble}, arg3::AbsoluteDate)
    return jcall(obj, "mappingFactors", Vector{jdouble}, (jdouble, jdouble, Vector{jdouble}, AbsoluteDate), arg0, arg1, arg2, arg3)
end


text:
@testset "Constant variables: " for solver in solvers

    @testset "Issue #166" begin
        # Issue #166
        α = Variable(5)
        fix!(α, ones(5,1))

        # has const vexity, but not at the head
        c = (rand(5,5) * α) * ones(1,5)

        β = Variable(5)
        β.value = ones(5)

        problem = Problem{BigFloat}(:minimize, sum(c * β), [β >= 0])
        solve!(problem, solver)
        @test problem.optval ≈ evaluate(sum(c * β)) atol=TOL
        @test problem.optval ≈ 0.0 atol=TOL
        @test β.value ≈ zeros(5) atol=TOL
    end

    @testset "Issue #228" begin
        x = Variable(2)
        y = Variable(2)
        fix!(x, [1 1]')
        prob = Problem{BigFloat}(:minimize, y'*(x+[2 2]'), [y>=0])
        solve!(prob, solver)
        @test prob.optval ≈ 0.0 atol = TOL

        prob = Problem{BigFloat}(:minimize, x'*y, [y>=0])
        solve!(prob, solver)
        @test prob.optval ≈ 0.0 atol = TOL
    end

end


text:
module VMLinuxBindings
    using CBinding

    @static if Sys.islinux()
        try
            vmlinux_iob = IOBuffer()
            run(pipeline(`bpftool btf dump file /sys/kernel/btf/vmlinux format c`; stdout=vmlinux_iob))
            vmlinux_str = String(take!(vmlinux_iob))

            vmlinux_path = joinpath(@__DIR__, "vmlinux.h")
            open(vmlinux_path, "w") do io
                write(io, vmlinux_str)
                flush(io)
            end

            @eval c`-std=c99 -fparse-all-comments -I$(@__DIR__)`
            @eval c"#include \"vmlinux.h\""
        catch err
            @warn "Failed to load VMLinux bindings!"
            @debug "VMLinux error: $err"
        end
    end
end


text:
module AlgoptTests

import Algopt

# Make sure to not start with an outdated registry
# https://github.com/JuliaLang/Pkg.jl/blob/master/test/runtests.jl
rm(joinpath(@__DIR__, "registries"); force = true, recursive = true)

tests = [
    "uni_derivative.jl",
    "bracketing.jl",
    "local_descent.jl",
    "first_order.jl",
    "test_functions.jl",
]

for test in tests
    include(test)
end

# clean up locally cached registry
# https://github.com/JuliaLang/Pkg.jl/blob/master/test/runtests.jl
rm(joinpath(@__DIR__, "registries"); force = true, recursive = true)

end # module


text:
using Documenter, SDFReader

makedocs(;
    modules=[SDFReader],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/SebastianM-C/SDFReader.jl/blob/{commit}{path}#L{line}",
    sitename="SDFReader.jl",
    authors="Sebastian Micluța-Câmpeanu <m.c.sebastian95@gmail.com>",
    assets=String[],
)


text:
#if !@isdefined  PDMats
#   using PDMats
#end

function calc_mean_stddev(pop::abc_population_type) #; param_active = nothing)
  theta_mean =  sum(pop.theta.*pop.weights',dims=2) # weighted mean for parameters
  theta_stddev = sqrt(var_weighted(pop.theta'.-theta_mean',pop.weights))  # scaled, weighted covar for parameters
  return (theta_mean,theta_stddev)
end

function make_proposal_dist_gaussian_full_covar(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false, # param_active = nothing,
         max_maha_distsq_per_dim::Real = 4.0)
  theta_mean =  sum(pop.theta.*pop.weights',dims=2) # weighted mean for parameters
  rawtau = cov_weighted(pop.theta'.-theta_mean',pop.weights)  # scaled, weighted covar for parameters
  tau = tau_factor*rawtau
  #tau = tau_factor*make_matrix_pd(rawtau)
  if verbose
    println("theta_mean = ", theta_mean)
    println("pop.theta = ", pop.theta)
    println("pop.weights = ", pop.weights)
    println("tau = ", tau)
  end
  covar = PDMat(tau)
  max_maha_distsq = 4.0*size(theta_mean,1)
  if max_maha_distsq_per_dim > 0
     max_maha_distsq = max_maha_distsq_per_dim*size(theta_mean,1)
     sampler = GaussianMixtureModelCommonCovarTruncated(pop.theta,pop.weights,covar,max_maha_distsq)
  else
     sampler = GaussianMixtureModelCommonCovar(pop.theta,pop.weights,covar)
  end
  return sampler
end

function make_proposal_dist_gaussian_diag_covar(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false, # param_active = nothing,
         max_maha_distsq_per_dim::Real = 4.0)
  theta_mean =  sum(pop.theta.*pop.weights',dims=2) # weighted mean for parameters
  tau = tau_factor*var_weighted(pop.theta'.-theta_mean',pop.weights)  # scaled, weighted covar for parameters
  if verbose
    println("theta_mean = ", theta_mean)
    println("pop.theta = ", pop.theta)
    println("pop.weights = ", pop.weights)
    println("tau = ", tau)
  end
  covar = PDiagMat(tau)
  max_maha_distsq = 4.0*size(theta_mean,1)
  if max_maha_distsq_per_dim > 0
      max_maha_distsq = max_maha_distsq_per_dim*size(theta_mean,1)
      sampler = GaussianMixtureModelCommonCovarTruncated(pop.theta,pop.weights,covar,max_maha_distsq)
  else
     sampler = GaussianMixtureModelCommonCovar(pop.theta,pop.weights,covar)
  end
  return sampler
end

function make_proposal_dist_gaussian_subset_full_covar(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false, param_active::Vector{Int64} = collect(1:size(pop.covar,1)) )
  weights = ones(pop.weights)./length(pop.weights) # pop.weights
  #weights = pop.weights
  theta_mean =  sum(pop.theta.*weights',dims=2) # weighted mean for parameters
  rawtau = cov_weighted(pop.theta'.-theta_mean',weights)  # scaled, weighted covar for parameters
  tau = tau_factor*rawtau
  #tau = tau_factor*make_matrix_pd(rawtau)
  if verbose
    println("theta_mean = ", theta_mean)
    println("pop.theta = ", pop.theta)
    println("pop.weights = ", pop.weights)
    println("weights = ", weights)
    println("tau = ", tau)
  end
  covar = tau # PDMat(tau[param_active])
  sampler = GaussianMixtureModelCommonCovarSubset(pop.theta,weights,covar,param_active)
end

function make_proposal_dist_gaussian_subset_diag_covar(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false, param_active::Vector{Int64} = collect(1:size(pop.covar,1)) )
  weights = ones(pop.weights)./length(pop.weights) # pop.weights
  theta_mean =  sum(pop.theta.*weights',dims=2) # weighted mean for parameters
  tau = tau_factor*var_weighted(pop.theta'.-theta_mean',weights)  # scaled, weighted covar for parameters
  if verbose
    println("theta_mean = ", theta_mean)
    println("pop.theta = ", pop.theta)
    println("pop.weights = ", pop.weights)
    println("weights = ", weights)
    println("tau = ", tau)
  end
  covar = tau
  sampler = GaussianMixtureModelCommonCovarSubset(pop.theta,weights,covar,param_active)
end

function make_proposal_dist_gaussian_rand_subset_full_covar(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false, num_param_active::Integer = 2)
  nparam = size(pop.theta,1)
  if num_param_active == nparam
     param_active = 1:nparam
  else
     param_active = union(sample(1:size(pop.theta,1),num_param_active,replace=false))
     println("# Proposals to perturb parameters: ", param_active)
  end
  make_proposal_dist_gaussian_subset_full_covar(pop,tau_factor, verbose=verbose, param_active=param_active )
end

function make_proposal_dist_gaussian_rand_subset_diag_covar(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false, num_param_active::Integer = 2)
  nparam = size(pop.theta,1)
  if num_param_active == nparam
     param_active = 1:nparam
  else
     param_active = union(sample(1:size(pop.theta,1),num_param_active,replace=false))
     println("# Proposals to perturb parameters: ", param_active)
  end
  make_proposal_dist_gaussian_subset_diag_covar(pop,tau_factor, verbose=verbose, param_active=param_active )
end

function make_proposal_dist_gaussian_rand_subset_neighbors_diag_covar(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false, num_param_active::Integer = 2)
  nparam = size(pop.theta,1)
  if num_param_active == nparam
     param_active = 1:nparam
  else
     param_active_first = rand(1:nparam)
     param_active = collect(param_active_first:(param_active_first+num_param_active-1))
     for i in 1:length(param_active)
         if param_active[i] > nparam
            param_active[i] -= nparam
         end
     end
     println("# Proposals to perturb parameters: ", param_active)
  end
  make_proposal_dist_gaussian_subset_diag_covar(pop,tau_factor, verbose=verbose, param_active=param_active )
end

function make_proposal_dist_gaussian_rand_subset_neighbors_full_covar(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false, num_param_active::Integer = 2)
  nparam = size(pop.theta,1)
  if num_param_active == nparam
     param_active = 1:nparam
  else
     param_active_first = rand(1:nparam)
     param_active = collect(param_active_first:(param_active_first+num_param_active-1))
     for i in 1:length(param_active)
         if param_active[i] > nparam
            param_active[i] -= nparam
         end
     end
     println("# Proposals to perturb parameters: ", param_active)
  end
  make_proposal_dist_gaussian_subset_full_covar(pop,tau_factor, verbose=verbose, param_active=param_active )
end

make_proposal_dist_gaussian_cycle_subset_neighbors_cycle_start_idx_next = 1
make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next = make_proposal_dist_gaussian_cycle_subset_neighbors_cycle_start_idx_next

function make_proposal_dist_gaussian_cycle_subset_neighbors_diag_covar(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false, num_param_active::Integer = 2)
  global make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next, make_proposal_dist_gaussian_cycle_subset_neighbors_cycle_start_idx_next
  nparam = size(pop.theta,1)
  if num_param_active >= nparam
     param_active = 1:nparam
  else
     param_active_first = ( (make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next-1) % nparam ) + 1
     param_active = collect(param_active_first:(param_active_first+num_param_active-1))
     for i in 1:length(param_active)
         if param_active[i] > nparam
            param_active[i] -= nparam
         end
     end
     println("# Proposals to perturb parameters: ", param_active)
     make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next += num_param_active
     if make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next > nparam
        make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next = make_proposal_dist_gaussian_cycle_subset_neighbors_cycle_start_idx_next
        make_proposal_dist_gaussian_cycle_subset_neighbors_cycle_start_idx_next = ( (make_proposal_dist_gaussian_cycle_subset_neighbors_cycle_start_idx_next) % nparam ) + 1
     end
  end
  make_proposal_dist_gaussian_subset_diag_covar(pop,tau_factor, verbose=verbose, param_active=param_active )
end

function make_proposal_dist_gaussian_cycle_subset_neighbors_full_covar(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false, num_param_active::Integer = 2)
  global make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next, make_proposal_dist_gaussian_cycle_subset_neighbors_cycle_start_idx_next
  nparam = size(pop.theta,1)
  if num_param_active == nparam
     param_active = 1:nparam
  else
     param_active_first = ( (make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next-1) % nparam ) + 1
     param_active = collect(param_active_first:(param_active_first+num_param_active-1))
     for i in 1:length(param_active)
         if param_active[i] > nparam
            param_active[i] -= nparam
         end
     end
     println("# Proposals to perturb parameters: ", param_active)
     make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next += num_param_active
     if make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next > nparam
        make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next = make_proposal_dist_gaussian_cycle_subset_neighbors_cycle_start_idx_next
        make_proposal_dist_gaussian_cycle_subset_neighbors_cycle_start_idx_next = ( (make_proposal_dist_gaussian_cycle_subset_neighbors_cycle_start_idx_next) % nparam ) + 1
     end
  end
  make_proposal_dist_gaussian_subset_full_covar(pop,tau_factor, verbose=verbose, param_active=param_active )
end

function make_proposal_dist_gaussian_subset_zoomed_diag_covar(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false, param_active::Vector{Int64} = collect(1:size(pop.covar,1)), inactive_scale_factor::Float64 = 0.25, max_maha_distsq_per_dim::Real = 4.0 )
  weights = pop.weights
  # weights = ones(pop.weights)./length(pop.weights) # pop.weights
  theta_mean =  sum(pop.theta.*weights',dims=2) # weighted mean for parameters
  tau = tau_factor*var_weighted(pop.theta'.-theta_mean',weights)  # scaled, weighted covar for parameters
  if verbose
    println("theta_mean = ", theta_mean)
    println("pop.theta = ", pop.theta)
    println("pop.weights = ", pop.weights)
    #println("weights = ", weights)
    println("tau = ", tau)
  end
  for i in 1:length(param_active)
    @assert 1<=param_active[i]<=length(tau)
    tau[param_active[i]] *= inactive_scale_factor
  end
  covar = PDiagMat(tau)
  max_maha_distsq = 4.0
  if max_maha_distsq_per_dim > 0
      max_maha_distsq = max_maha_distsq_per_dim*size(theta_mean,1)
      sampler = GaussianMixtureModelCommonCovarTruncated(pop.theta,pop.weights,covar,max_maha_distsq)
  else
     sampler = GaussianMixtureModelCommonCovar(pop.theta,pop.weights,covar)
  end
  return sampler
end

function make_proposal_dist_gaussian_cycle_zoomed_subset_neighbors_diag_covar(pop::abc_population_type, tau_factor::Float64; verbose::Bool = false, num_param_active::Integer = 2, inactive_scale_factor::Float64 = 0.25,  max_maha_distsq_per_dim::Real = 4.0 )
  global make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next, make_proposal_dist_gaussian_cycle_subset_neighbors_cycle_start_idx_next
  nparam = size(pop.theta,1)
  if num_param_active == nparam
     param_active = 1:nparam
  else
     param_active_first = ( (make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next-1) % nparam ) + 1
     param_active = collect(param_active_first:(param_active_first+num_param_active-1))
     for i in 1:length(param_active)
         if param_active[i] > nparam
            param_active[i] -= nparam
         end
     end
     println("# Proposals to perturb parameters: ", param_active)
     make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next += num_param_active
     if make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next > nparam
        make_proposal_dist_gaussian_cycle_subset_neighbors_param_idx_next = make_proposal_dist_gaussian_cycle_subset_neighbors_cycle_start_idx_next
        make_proposal_dist_gaussian_cycle_subset_neighbors_cycle_start_idx_next = ( (make_proposal_dist_gaussian_cycle_subset_neighbors_cycle_start_idx_next) % nparam ) + 1
     end
  end
  make_proposal_dist_gaussian_subset_zoomed_diag_covar(pop,tau_factor, verbose=verbose, param_active=param_active, inactive_scale_factor=inactive_scale_factor, max_maha_distsq_per_dim=max_maha_distsq_per_dim )
end


text:
using Test
using TestSetExtensions
using Qaintessent

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "view basic" begin

    N = 5
    cgc = [
        circuit_gate(3, HadamardGate()),
        circuit_gate(2, RxGate(√0.2), (1, 4)),
        circuit_gate((1,5), SwapGate(), (2,4)),
        circuit_gate(2, TGate()),
        circuit_gate(2, PhaseShiftGate(0.2π)),
        circuit_gate(2, SdagGate()),
        circuit_gate(3, RotationGate(0.1π, [1, 0, 0])),
        circuit_gate(4, TdagGate()),
        circuit_gate(1, RyGate(1.4π)),
        circuit_gate(3, RzGate(0.4π)),
        circuit_gate(1,2, SwapGate()),
        circuit_gate(4, SGate()),
        circuit_gate((3,5), SwapGate(), 4),
        circuit_gate((3,), Y, 2),
        circuit_gate((5,), Z, 1),
        circuit_gate((4,), X)
    ]

    cgc_refstring =
        "\n" *
        "    5 ——————————————x—————————————————————————————x——————————[Z ]———————\n" *
        "                    |                             |           |         \n" *
        "    4 ————————•—————•————————————————[T†]——[S ]———•————————————————[X ]—\n" *
        "              |     |                             |           |         \n" *
        "    3 —[H ]——————————————————————————[Rθ]——[Rz]———x————[Y ]—————————————\n" *
        "              |     |                                   |     |         \n" *
        "    2 ———————[Rx]———•————[T ]——[Pϕ]——[S†]———x———————————•———————————————\n" *
        "              |     |                       |                 |         \n" *
        "    1 ————————•—————x————————————————[Ry]———x—————————————————•—————————\n"

    io = IOBuffer()
    show(io, cgc)
    @test String(take!(io)) == cgc_refstring

end

@testset ExtendedTestSet "empty basic view" begin
    cgc = CircuitGate[]

    cgc_refstring = "[]\n"

    io = IOBuffer()
    show(io, cgc)
    @test String(take!(io)) == cgc_refstring
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "view circuit" begin

    N = 5
    cgc = Circuit{N}([
        Moment(
        [circuit_gate(3, HadamardGate()),
        circuit_gate(2, RxGate(√0.2), (1, 4))]
        ),
        Moment(circuit_gate((1,5), SwapGate(), (2,4))),
        Moment([circuit_gate(2, PhaseShiftGate(0.2π)),
        circuit_gate(3, RotationGate(0.1π, [1, 0, 0])),
        circuit_gate(1, RyGate(1.4π))]),
        Moment([circuit_gate(1,2, SwapGate()),
        circuit_gate((3,5), SwapGate(), 4)]),
    ])

    cgc_refstring =
        "\n" *
        "    5 ——————————————x———————————x———\n" *
        "                    |           |   \n" *
        "    4 ————————•—————•———————————•———\n" *
        "              |     |           |   \n" *
        "    3 —[H ]——————————————[Rθ]———x———\n" *
        "              |     |               \n" *
        "    2 ———————[Rx]———•————[Pϕ]———x———\n" *
        "              |     |           |   \n" *
        "    1 ————————•—————x————[Ry]———x———\n"

    io = IOBuffer()
    show(io, cgc)
    @test String(take!(io)) == cgc_refstring

end

@testset ExtendedTestSet "view empty circuit" begin
    N = 5
    cgc = Circuit{N}()
    cgc_refstring = "[]\n"

    io = IOBuffer()
    show(io, cgc)
    @test String(take!(io)) == cgc_refstring
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "view moment" begin

    N = 5
    m = Moment([circuit_gate(2, PhaseShiftGate(0.2π)),
        circuit_gate(3, RotationGate(0.1π, [1, 0, 0])),
        circuit_gate(1, RyGate(1.4π))])

    m_refstring =
        "\n" *
        "    3 ...|—[Rθ]—|...\n" *
        "         |      |   \n" *
        "    2 ...|—[Pϕ]—|...\n" *
        "         |      |   \n" *
        "    1 ...|—[Ry]—|...\n" *
        "\n"
    io = IOBuffer()
    show(io, [m])
    @test String(take!(io)) == m_refstring

    m = Moment([circuit_gate(5, PhaseShiftGate(0.2π)),
        circuit_gate(3, RotationGate(0.1π, [1, 0, 0])),
        circuit_gate(1, RyGate(1.4π)),
        circuit_gate((4), HadamardGate(), (2))])

    m_refstring =
        "\n" *
        "    5 ...|—[Pϕ]———————|...\n" *
        "         |            |   \n" *
        "    4 ...|———————[H ]—|...\n" *
        "         |        |   |   \n" *
        "    3 ...|—[Rθ]———————|...\n" *
        "         |        |   |   \n" *
        "    2 ...|————————•———|...\n" *
        "         |            |   \n" *
        "    1 ...|—[Ry]———————|...\n" *
        "\n"

    io = IOBuffer()
    show(io, [m])
    @test String(take!(io)) == m_refstring

    m = [Moment([circuit_gate(5, PhaseShiftGate(0.2π)),
        circuit_gate(3, RotationGate(0.1π, [1, 0, 0])),
        circuit_gate(1, RyGate(1.4π)),
        circuit_gate((4), HadamardGate(), (2))]),
        Moment([circuit_gate(5, PhaseShiftGate(0.2π)),
        circuit_gate(3, RotationGate(0.1π, [1, 0, 0])),
        circuit_gate(1, RyGate(1.4π)),
        circuit_gate((4), HadamardGate(), (2))]),
        ]

    m_refstring =
        "Vector{Moment}[\n" *
        "    5 ...|—[Pϕ]———————|...\n" *
        "         |            |   \n" *
        "    4 ...|———————[H ]—|...\n" *
        "         |        |   |   \n" *
        "    3 ...|—[Rθ]———————|...\n" *
        "         |        |   |   \n" *
        "    2 ...|————————•———|...\n" *
        "         |            |   \n" *
        "    1 ...|—[Ry]———————|...\n" *
        "\n" *
        "\n" *
        "    5 ...|—[Pϕ]———————|...\n" *
        "         |            |   \n" *
        "    4 ...|———————[H ]—|...\n" *
        "         |        |   |   \n" *
        "    3 ...|—[Rθ]———————|...\n" *
        "         |        |   |   \n" *
        "    2 ...|————————•———|...\n" *
        "         |            |   \n" *
        "    1 ...|—[Ry]———————|...\n" *
        "\n]"

    io = IOBuffer()
    show(io, [m])
    @test String(take!(io)) == m_refstring

    m = Moment([circuit_gate(5, PhaseShiftGate(0.2π)),
        circuit_gate(3, RotationGate(0.1π, [1, 0, 0])),
        circuit_gate((2), HadamardGate(), 4),
        circuit_gate(1, RyGate(1.4π))])
    m_refstring =
    "CircuitGate{1, PhaseShiftGate}((5,), PhaseShiftGate([0.6283185307179586]))\n" *
    "CircuitGate{1, RotationGate}((3,), RotationGate([0.3141592653589793, 0.0, 0.0]))\n" *
    "CircuitGate{2, ControlledGate{HadamardGate}}((2, 4), ControlledGate{HadamardGate}(HadamardGate(), 1))\n" *
    "CircuitGate{1, RyGate}((1,), RyGate([4.39822971502571]))\n"

    io = IOBuffer()
    show(io, m)
    @test String(take!(io)) == m_refstring
end

@testset ExtendedTestSet "view empty moment" begin
    N = 5
    cgc = CircuitGate[]
    m = Moment(cgc)
    m_refstring = ""

    io = IOBuffer()
    show(io, m)
    @test String(take!(io)) == m_refstring
end

@testset ExtendedTestSet "view empty moments" begin
    N = 5
    cgc = CircuitGate[]
    m = [Moment(cgc)]
    m_refstring = ""

    io = IOBuffer()
    show(io, m)
    @test String(take!(io)) == m_refstring
end


text:
module Zoo

using Peacock

"""
    wu_ep(x, y, angles, Rs, d1s, d2s, ep_bg, ep_cyl)

Returns the permittivity `ep(x,y)` for photonic crystals similar to Wu et al, 2015.
"""
function wu_ep(x, y, angles, Rs, d1s, d2s, ep_bg, ep_cyl)
    if length(Rs) == 1
        Rs = fill(Rs[1], length(angles))
    end
    if length(d1s) == 1
        d1s = fill(d1s[1], length(angles))
    end
    if length(d2s) == 1
        d2s = fill(d2s[1], length(angles))
    end
    @assert length(Rs) == length(d1s) == length(d2s) == length(angles)
    # Triangular lattice sites
    lx, ly = 1, 1*sqrt(3)
    sites = [
        (0, 0),
        (-lx, 0),
        (+lx, 0),
        (-lx/2, -ly/2),
        (+lx/2, -ly/2),
        (+lx/2, +ly/2),
        (-lx/2, +ly/2),
    ]
    # Hexagonal ring of cylinders at each site
    for (x0,y0) in sites
        for (R,d1,d2,th) in zip(Rs,d1s,d2s,angles)
            x_, y_ = [cosd(th) sind(th); -sind(th) cosd(th)] * [(x-x0); (y-y0)]
            x_ = (x_-R) / (d1/2)
            y_ = (y_-0) / (d2/2)
            if x_^2 + y_^2 <= 1
                return ep_cyl
            end
        end
    end
    return ep_bg
end

"""
    make_wu_generic(Rs, d1s, d2s, polarisation, cutoff;
                angles=[60n for n in 0:5], ep_bg=1, ep_cyl=11.7)

Returns a `NamedTuple` of `geometry`, `solver`, `polarisation`, and the `G`,
`M`, `K` high-symmetry points for photonic crystals similar to Wu et al, 2015.
"""
function make_wu_generic(Rs, d1s, d2s, polarisation, cutoff; angles=[60n for n in 0:5], ep_bg=1, ep_cyl=11.7)
    epf(x,y) = wu_ep(x, y, angles, Rs, d1s, d2s, ep_bg, ep_cyl)
    muf(x,y) = 1.0
    geometry = Geometry(epf, muf, -60, +60, 1/500, 1/500)
    solver = Solver(geometry, cutoff)
    G = BrillouinZoneCoordinate(  0,   0, "Γ")
    M = BrillouinZoneCoordinate(  0, 1/2, "M")
    K = BrillouinZoneCoordinate(1/3, 1/3, "K")
    # Return as a NamedTuple so variables can be easily unpacked using Parameters.jl
    return (; geometry=geometry, solver=solver,
            G=G, M=M, K=K, polarisation=polarisation)
end

"""
    make_wu(a0_div_R, cutoff)

Reproduces "Scheme for Achieving a Topological Photonic Crystal
by Using Dielectric Material", Wu et al, 2015.

Returns a `NamedTuple` of `geometry`, `solver`, `polarisation`, and the
`G`, `M`, `K` high-symmetry points.
"""
function make_wu(a0_div_R, cutoff)
    # Geometry
    R = 1 / a0_div_R # assuming a0 = 1
    d = 2R/3
    return make_wu_generic(R, d, d, TM, cutoff)
end

"""
    make_wu_topo(cutoff::Int)

Reproduces the trivial crystal (a0_div_R=3.125)from  "Scheme for Achieving a
Topological Photonic Crystal by Using Dielectric Material", Wu et al, 2015.

Returns a `NamedTuple` of `geometry`, `solver`, `polarisation`, and the
`G`, `M`, `K` high-symmetry points.
"""
function make_wu_triv(cutoff::Int)
    return make_wu(3.125, cutoff)
end

"""
    make_wu_topo(cutoff::Int)

Reproduces the topological crystal (`a0_div_R=2.9`) from "Scheme for Achieving a
Topological Photonic Crystal by Using Dielectric Material", Wu et al, 2015.

Returns a `NamedTuple` of `geometry`, `solver`, `polarisation`, and the
`G`, `M`, `K` high-symmetry points.
"""
function make_wu_topo(cutoff::Int)
    return make_wu(2.9, cutoff)
end

"""
    make_wu_primitive(cutoff)

Reproduces the primitive unit cell of "Scheme for Achieving a
Topological Photonic Crystal by Using Dielectric Material", Wu et al, 2015

Returns a `NamedTuple` of `geometry`, `solver`, `polarisation`, and the
`G`, `M`, `K` high-symmetry points.
"""
function make_wu_primitive(cutoff)
    # Geometry
    R = 1 / 3 # assuming a0 = 1
    d = 2R/3
    x = 1 / sqrt(3)
    a1 = [x*cosd(-30), x*sind(-30)]
    a2 = [x*cosd(+30), x*sind(+30)]
    epf(x,y) = wu_ep(x, y, 0:60:300, R, d, d, 1, 11.7)
    muf(x,y) = 1.0
    g = Geometry(epf, muf, a1, a2, 1/500, 1/500)
    solver = Solver(cg, TM)
    # BZ symmetry points
    G = BrillouinZoneCoordinate(  0,   0, "Γ")
    M = BrillouinZoneCoordinate(  0, 1/2, "M")
    K = BrillouinZoneCoordinate(-1/3, 1/3, "K")
    return g, solver, G, M, K
end

"""
    make_dePaz(d1::Real, d2::Real, cutoff::Int; R::Real=1/3)

Reproduces "Engineering fragile topology in photonic crystals:
Topological quantum chemistry of light", Blanco de Paz et al, 2019

Returns a `NamedTuple` of `geometry`, `solver`, `polarisation`, and the
`G`, `M`, `K` high-symmetry points.
"""
function make_dePaz(d1::Real, d2::Real, cutoff::Int; R::Real=1/3)
    return make_wu_generic(R, d1, d2, TM, cutoff)
end

"""
    make_dePaz_triv(cutoff::Int)

Reproduces the trivial (`d1=0.52`, `d2=0.31`) crystal from "Engineering fragile
topology in photonic crystals: Topological quantum chemistry of light",
Blanco de Paz et al, 2019

Returns a `NamedTuple` of `geometry`, `solver`, `polarisation`, and the
`G`, `M`, `K` high-symmetry points.
"""
function make_dePaz_triv(cutoff::Int)
    return make_dePaz(0.52, 0.31, cutoff)
end

"""
    make_dePaz_frag(cutoff::Int)

Reproduces the trivial (`d1=0.4`, `d2=0.13`) crystal from "Engineering fragile
topology in photonic crystals: Topological quantum chemistry of light",
Blanco de Paz et al, 2019

Returns a `NamedTuple` of `geometry`, `solver`, `polarisation`, and the
`G`, `M`, `K` high-symmetry points.
"""
function make_dePaz_frag(cutoff::Int)
    return make_dePaz(0.4, 0.13, cutoff)
end

"""
    make_dePaz_OAL(cutoff::Int)

Reproduces the obstructed atomic limit (`d1=0.4`, `d2=0.61`) crystal from
"Engineering fragile topology in photonic crystals: Topological quantum
chemistry of light", Blanco de Paz et al, 2019

Returns a `NamedTuple` of `geometry`, `solver`, `polarisation`, and the
`G`, `M`, `K` high-symmetry points.
"""
function make_dePaz_OAL(cutoff::Int)
    return make_dePaz(0.4, 0.61, cutoff)
end

"""
    make_zhu(r, l, cutoff; angles=[60n for n in 0:5])

Reproduces "Topological transitions in continuously deformed photonic crystals"
Zhu et al, 2018.

- `r` is the radius of the cylindrical air holes.
- `l` is the radius of the hexagonal ring that the air holes form

Note that 'l' in Zhu et al is the same as 'R' in wu_ep.

Returns a `NamedTuple` of `geometry`, `solver`, `polarisation`, and the
`G`, `M`, `K` high-symmetry points.
"""
function make_zhu(r, l, cutoff; angles=[60n for n in 0:5])
    return make_wu_generic(l, 2r, 2r, TE, cutoff, ep_cyl=1, ep_bg=10.2)
end

"""
    make_xu(r1, r2, cutoff; wedges=[(0,360)])

Reproduces "Accidental degeneracy in photonic bands and topological phase
transitions in two-dimensional core-shell dielectric photonic crystals",
Xu et al, 2016.

The rings can be split using `wedges = [(start,stop), (start,stop)...]`.

Returns a `NamedTuple` of `geometry`, `solver`, `polarisation`, and the
`G`, `M`, `K` high-symmetry points.
"""
function make_xu(r1, r2, cutoff; wedges=[(0,360)])
    function epf(x,y)
        lx, ly = 1, 1*sqrt(3)
        sites = [
            (    0,     0),
            (-lx/2, -ly/2),
            (+lx/2, -ly/2),
            (+lx/2, +ly/2),
            (-lx/2, +ly/2),
        ]
        for (x0,y0) in sites
            x_, y_ = x-x0, y-y0
            if r1^2 <= x_^2 + y_^2 <= r2^2
                θ = atan(y_,x_) * 180/pi
                for (θ1,θ2) in wedges
                    dθ = mod(θ - θ1, 360)
                    if dθ <= (θ2-θ1)
                        return 12.0
                    end
                end
            end
        end
        return 1.0
    end
    muf(x,y) = 1.0
    geometry = Geometry(epf, muf, -60, +60, 1/500, 1/500)
    solver = Solver(geometry, cutoff)
    G = BrillouinZoneCoordinate(  0,   0, "Γ")
    M = BrillouinZoneCoordinate(  0, 1/2, "M")
    K = BrillouinZoneCoordinate(1/3, 1/3, "K")
    return (; geometry=geometry, solver=solver,
            G=G, M=M, K=K, polarisation=TM)
end

export make_wu, make_wu_triv, make_wu_topo, make_wu_primitive
export make_dePaz, make_dePaz_triv, make_dePaz_frag, make_dePaz_OAL
export make_zhu
export make_xu

end # module


text:
# This was copied from Debugger.jl and then modified

using TerminalRegressionTests, Rebugger, Revise
using HeaderREPLs, REPL
using Test

function run_terminal_test(cmd, validation, commands)
    dirpath = joinpath(@__DIR__, "ui", "v$(VERSION.major).$(VERSION.minor)")
    isdir(dirpath) || mkpath(dirpath)
    filepath = joinpath(dirpath, validation)
    TerminalRegressionTests.automated_test(filepath, commands) do emuterm
    # TerminalRegressionTests.create_automated_test(filepath, commands) do emuterm
        main_repl = REPL.LineEditREPL(emuterm, true)
        main_repl.interface = REPL.setup_interface(main_repl)
        main_repl.specialdisplay = REPL.REPLDisplay(main_repl)
        main_repl.mistate = REPL.LineEdit.init_state(REPL.terminal(main_repl), main_repl.interface)
        iprompt, eprompt = Rebugger.rebugrepl_init(main_repl, true)
        repl = iprompt.repl
        s = repl.mistate
        s.current_mode = iprompt
        repl.t = emuterm
        REPL.LineEdit.edit_clear(s)
        REPL.LineEdit.edit_insert(s, cmd)
        Rebugger.interpret(s)
    end
end

Revise.track(Base)  # just to get Info printing out of the way

CTRL_C = "\x3"
EOT = "\x4"
UP_ARROW = "\e[A"

run_terminal_test("gcd(10, 20)",
                  "gcd.multiout",
                  ['\n'])
run_terminal_test("__gcdval__ = gcd(10, 20);",
                  "gcdsc.multiout",
                  ['\n'])
@test __gcdval__ == gcd(10, 20)


text:
@testset "Statistics" begin
  @testset "Ensemble" begin
    grid = CartesianGrid(3,3)
    reals = Dict(:value => [i*ones(nelements(grid)) for i in 1:3])
    ensemble = Ensemble(grid, reals)

    # mean
    mean2D = mean(ensemble)
    @test mean2D[:value] == 2.0*ones(nelements(mean2D))
    @test domain(mean2D) == ensemble.domain

    # variance
    var2D = var(ensemble)
    @test var2D[:value] == 1.0*ones(nelements(var2D))
    @test domain(var2D) == ensemble.domain

    # quantile (scalar)
    p = 0.5
    quant2D = quantile(ensemble, p)
    @test quant2D[:value] == 2.0*ones(nelements(quant2D))
    @test domain(quant2D) == ensemble.domain

    # quantile (vector)
    ps = [0.0, 0.5, 1.0]
    quants2D = quantile(ensemble, ps)
    @test quants2D[2][:value] == quant2D[:value]
  end

  @testset "Data" begin
    # load data with bias towards large values (gold mine)
    sdata = georef(CSV.File(joinpath(datadir,"clustered.csv")), (:x,:y))

    # spatial mean
    μn = mean(sdata[:Au])
    μs = mean(sdata, :Au)
    @test abs(μn - 0.5) > abs(μs - 0.5)
    @test mean(sdata)[:Au] ≈ μs

    # spatial variance
    σn = var(sdata[:Au])
    σs = var(sdata, :Au)
    @test σn ≤ σs
    @test var(sdata)[:Au] ≈ σs

    # spatial quantile
    qn = quantile(sdata[:Au], 0.5)
    qs = quantile(sdata, :Au, 0.5)
    @test qn ≥ qs
    @test quantile(sdata, 0.5)[:Au] ≈ qs
  end
end


text:
import RDatasets
import DataFrames: DataFrame
using LowRankModels

# pick a data set
df = RDatasets.dataset("psych", "msq")

# make a GLRM on the whole dataframe using type imputation
auto_glrm, labels = GLRM(df, 3)
# fit!(auto_glrm)  - this doesn't work yet, because the data contains some trivial columns (ordinals with <2 levels)

# now we'll try it without type imputation
# we'll just fit four of the columns, to try out all four data types
dd = DataFrame([df[s] for s in [:TOD, :Scale, :Vigorous, :Wakeful]])
dd[!,end] = (dd[:,end].==1)
datatypes = [:real, :cat, :ord, :bool]

# fit it!
glrm = GLRM(dd, 2, datatypes)

println("initializing")
init_svd!(glrm)

println("fitting")
X, Y, ch = fit!(glrm)

# print results
println(ch.objective)

println("imputing")
impute(glrm)

println("crossvalidating")
cross_validate(glrm, do_obs_check=false, init=init_svd!)


text:
#
# Circular correlation statistic of Jupp and Mardia
#
# See Jupp, P. E., & Mardia, K. V. (1980). A General Correlation
# Coefficient for Directional Data and Related Regression Problems.
# Biometrika, 67(1), 163–173. doi:10.2307/2335329
#
# The formula in terms of correlation coefficients given in the text
# is slightly wrong. The correct formula is given in:
#
# Mardia, K. V., & Jupp, P. E. (2009). Directional Statistics. John
# Wiley & Sons, p. 249.
immutable JuppMardiaR{Normalized} <: NormalizedPairwiseStatistic{Normalized} end
JuppMardiaR() = JuppMardiaR{false}()
Base.eltype{T<:Real}(::JuppMardiaR, X::AbstractArray{Complex{T}}) = T

function meandiff!(out, X)
    for j = 1:size(X, 2)
        # Compute mean resultant
        m = zero(eltype(X))
        @simd for i = 1:size(X, 1)
            @inbounds m += X[i, j]
        end
        m /= size(X, 1)

        # Compute difference from mean
        @simd for i = 1:size(X, 1)
            @inbounds z = X[i, j] - m
            @inbounds out[i, 2j-1] = real(z)
            @inbounds out[i, 2j] = imag(z)
        end
    end
    out
end

immutable JuppMardiaRWorkX{T<:Real}
    workX::Matrix{T}
    covwork::Matrix{T}

    JuppMardiaRWorkX(w1, w2) = new(w1, w2)
    JuppMardiaRWorkX(w1, w2, w3) = new(w1, w2, w3)
end
allocwork{T<:Real}(t::JuppMardiaR{true}, X::AbstractVecOrMat{Complex{T}}) =
    (Array(T, size(X, 1), size(X, 2)*2),  Array(T, size(X, 2)*2, size(X, 2)*2))
function computestat!{T<:Real}(t::JuppMardiaR{true}, out::AbstractMatrix{T}, work::Tuple{Matrix{T}, Matrix{T}},
                               X::AbstractVecOrMat{Complex{T}})
    chkinput(out, X)
    workX, covwork = work

    # Correlation between all components
    r = At_mul_A!(covwork, meandiff!(workX, X))
    cov2coh!(r)

    for j = 1:nchannels(X)
        c1 = 2j-1
        s1 = 2j
        for i = 1:j-1
            c2 = 2i-1
            s2 = 2i
            @inbounds out[i, j] = ((r[c2, c1]^2 + r[s2, c1]^2 + r[c2, s1]^2 + r[s2, s1]^2) +
                                   2*(r[c2, c1]*r[s2, s1] + r[s2, c1]*r[c2, s1])*r[c1, s1]*r[c2, s2] -
                                   2*(r[c2, c1]*r[s2, c1] + r[c2, s1]*r[s2, s1])*r[c2, s2] -
                                   2*(r[c2, c1]*r[c2, s1] + r[s2, c1]*r[s2, s1])*r[c1, s1]) /
                                  ((1-r[c1, s1]^2)*(1-r[c2, s2]^2))
        end
        @inbounds out[j, j] = 2
    end
    out
end

immutable JuppMardiaRWorkXY{T<:Real}
    workX::Matrix{T}
    workY::Matrix{T}
    covwork::Matrix{T}
    invsqrtsumX::Matrix{T}
    invsqrtsumY::Matrix{T}
end
function allocwork{T<:Real}(t::JuppMardiaR{true}, X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}})
    JuppMardiaRWorkXY{T}(Array(T, size(X, 1), size(X, 2)*2),
                         Array(T, size(Y, 1), size(Y, 2)*2),
                         Array(T, size(X, 2)*2, size(Y, 2)*2),
                         Array(T, 1, size(X, 2)*2),
                         Array(T, 1, size(Y, 2)*2))
end
function reimcor!(workX, invsqrtsumX, n)
    for j = 1:div(size(workX, 2), 2)
        v = zero(eltype(workX))
        @simd for i = 1:size(workX, 1)
            @inbounds v += workX[i, 2j-1] * workX[i, 2j]
        end
        @inbounds workX[j] = v * invsqrtsumX[2j-1] * invsqrtsumX[2j]
    end
    workX
end
function computestat!{T<:Real}(t::JuppMardiaR{true}, out::AbstractMatrix{T}, work::JuppMardiaRWorkXY{T},
                               X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}})
    chkinput(out, X, Y)

    # Compute mean phases
    workX = meandiff!(work.workX, X)
    workY = meandiff!(work.workY, Y)

    # Correlation between signals
    r = At_mul_B!(work.covwork, workX, workY)
    cov2coh!(r, workX, workY, work.invsqrtsumX, work.invsqrtsumY, r)

    # Correlation between real and imaginary components of each signal
    reimcor!(workX, work.invsqrtsumX, size(X, 1))
    reimcor!(workY, work.invsqrtsumY, size(Y, 1))

    for j = 1:nchannels(Y)
        c1 = 2j-1
        s1 = 2j
        for i = 1:nchannels(X)
            c2 = 2i-1
            s2 = 2i
            @inbounds out[i, j] = ((r[c2, c1]^2 + r[s2, c1]^2 + r[c2, s1]^2 + r[s2, s1]^2) +
                                   2*(r[c2, c1]*r[s2, s1] + r[s2, c1]*r[c2, s1])*workY[j]*workX[i] -
                                   2*(r[c2, c1]*r[s2, c1] + r[c2, s1]*r[s2, s1])*workX[i] -
                                   2*(r[c2, c1]*r[c2, s1] + r[s2, c1]*r[s2, s1])*workY[j]) /
                                  ((1-workY[j]^2)*(1-workX[i]^2))
        end
    end
    out
end


text:
# # Sampling from a GP
#
# ## Preliminary steps

# ### Loading necessary packages

using Plots
using AugmentedGaussianProcesses
using Distributions
using LinearAlgebra

# ### Loading the banana dataset from OpenML
kernel = SqExponentialKernel()
x = range(0, 10, length=50)
K = kernelmatrix(kernel, x)
f = rand(MvNormal(K + 1e-8I)) # Sample a random GP
y = rand.(Bernoulli.(AGP.logistic.(f)))
y_sign = sign.(y .- 0.5)

# ### We create a function to visualize the data

function plot_data(x, y; size=(300,500))
    Plots.scatter(x,
                y,
                alpha=0.2,
                markerstrokewidth=0.0,
                lab="",
                size=size
            )
end
plot_data(x, y; size = (500, 500))

# ### Run the variational gaussian process approximation
@info "Running full model"
mfull = VGP(x, y_sign,
            kernel,
            LogisticLikelihood(),
            AnalyticVI(),
            optimiser = false
            )
@time train!(mfull, 5)

# ### We can also create a sampling based model
@info "Sampling from model"
mmcmc = MCGP(x, y,
            kernel,
            LogisticLikelihood(),
            GibbsSampling(),
            optimiser = false
            )
m = mmcmc
@time samples = sample(mmcmc, 1000)

# ### We can now visualize the results of both models

# ### We first plot the latent function f (truth, the VI estimate, the samples)
p1 = plot(x, f, label="true f")
plot!(x, samples, label="", color=:black, alpha=0.02, lab="")
plot!(x, mean(mfull[1]), ribbon=sqrt.(var(mfull[1])), label="VI")
# ### And we can also plot the predictions vs the data
p2 = plot_data(x, y; size=(600,400))
μ_vi, σ_vi = proba_y(mfull, x)
plot!(x, μ_vi; ribbon=σ_vi, label="VI")
μ_mcmc, σ_mcmc = proba_y(mmcmc, x)
plot!(x, μ_mcmc; ribbon=σ_mcmc, label="MCMC")


text:
using Distributed
@everywhere using Mamba

## Data
kidney = Dict{Symbol, Any}(
  :t => permutedims(reshape(
    [8, 16, 23, NaN, 22, 28, 447, 318, 30, 12, 24, 245, 7, 9, 511, 30, 53, 196,
     15, 154, 7, 333, 141, NaN, 96, 38, NaN, NaN, 536, NaN, 17, NaN, 185, 177,
     292, 114, NaN, NaN, 15, NaN, 152, 562, 402, NaN, 13, 66, 39, NaN, 12, 40,
     NaN, 201, 132, 156, 34, 30, 2, 25, 130, 26, 27, 58, NaN, 43, 152, 30, 190,
     NaN, 119, 8, NaN, NaN, NaN, 78, 63, NaN],
    2, 38)),
  :tcensor => permutedims(reshape(
    [0, 0, 0, 13, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8, 0,
     0, 149, 70, 0, 25, 0, 4, 0, 0, 0, 0, 22, 159, 0, 108, 0, 0, 0, 24, 0, 0, 0,
     46, 0, 0, 113, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 5, 0, 0, 54,
     16, 6, 0, 0, 8],
    2, 38)),
  :age => permutedims(reshape(
    [28, 28, 48, 48, 32, 32, 31, 32, 10, 10, 16, 17, 51, 51, 55, 56, 69, 69, 51,
     52, 44, 44, 34, 34, 35, 35, 42, 42, 17, 17, 60, 60, 60, 60, 43, 44, 53, 53,
     44, 44, 46, 47, 30, 30, 62, 63, 42, 43, 43, 43, 57, 58, 10, 10, 52, 52, 53,
     53, 54, 54, 56, 56, 50, 51, 57, 57, 44, 45, 22, 22, 42, 42, 52, 52, 60, 60],
    2, 38)),
  :sex => [0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1,
           1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0],
  :disease => [1, 2, 1, 1, 1, 1, 2, 2, 3, 2, 3, 1, 3, 3, 1, 3, 1, 1, 2, 1, 4,
               1, 3, 3, 3, 3, 2, 3, 2, 2, 3, 3, 4, 2, 1, 1, 4, 4],
  :N => 38,
  :M => 2
)

kidney[:Dx] = Int[
  kidney[:disease][i] == j ? 1 : 0
  for i in 1:38, j in 2:4
]

## Model Specification
model = Model(

  t = Stochastic(2,
    (alpha, beta_age, age, beta_sex, sex, Dx, beta_Dx, b, r, tcensor, N, M) ->
      begin
        beta_dis = Dx * beta_Dx
        UnivariateDistribution[(
          mu = alpha + beta_age * age[i, j] + beta_sex * sex[i] + beta_dis[i] +
               b[i];
          lambda = exp(-mu / r);
          0 < lambda < Inf ?
            Truncated(Weibull(r, lambda), tcensor[i, j], Inf) :
            Uniform(0, Inf)) for i in 1:N, j in 1:M
        ]
      end,
    false
  ),

  b = Stochastic(1,
    s2 -> Normal(0, sqrt(s2)),
    false
  ),

  s2 = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  ),

  alpha = Stochastic(
    () -> Normal(0, 100)
  ),

  beta_age = Stochastic(
    () -> Normal(0, 100)
  ),

  beta_sex = Stochastic(
    () -> Normal(0, 100)
  ),

  beta_Dx = Stochastic(1,
    () -> Normal(0, 100)
  ),

  r = Stochastic(
    () -> Gamma(1, 1000)
  )

)

## Initial Values
inits = [
  Dict(:t => kidney[:t], :alpha => 0, :beta_age => 0, :beta_sex => 0,
       :beta_Dx => zeros(3), :s2 => 3, :r => 1.0, :b => zeros(kidney[:N])),
  Dict(:t => kidney[:t], :alpha => 1, :beta_age => -1, :beta_sex => 1,
       :beta_Dx => ones(3), :s2 => 1, :r => 1.5, :b => zeros(kidney[:N]))
]

## Sampling Scheme
scheme = [MISS(:t),
          Slice([:alpha, :beta_age, :beta_sex, :beta_Dx], 0.1),
          Slice(:b, 0.01),
          Slice(:s2, 0.1),
          Slice(:r, 0.001)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, kidney, inits, 20000, burnin=2500, thin=2, chains=2)
describe(sim)