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

@test integrate((b*tan(c+d*x))^-1//2*(A+C*cot(c+d*x)^2), x) == :(-2//3*C*b*d^-1*(b*tan(c+d*x))^-3//2+(1/2)*2^(1/2)*b^-1//2*d^-1*(A+-1C)*arctan(1+2^(1/2)*b^-1//2*(b*tan(c+d*x))^(1/2))+-1//2*2^(1/2)*b^-1//2*d^-1*(A+-1C)*arctan(1+-1*2^(1/2)*b^-1//2*(b*tan(c+d*x))^(1/2))+-1//4*2^(1/2)*b^-1//2*d^-1*(A+-1C)*log(b^(1/2)+b^(1/2)*tan(c+d*x)+-1*2^(1/2)*(b*tan(c+d*x))^(1/2))+1//4*2^(1/2)*b^-1//2*d^-1*(A+-1C)*log(b^(1/2)+2^(1/2)*(b*tan(c+d*x))^(1/2)+b^(1/2)*tan(c+d*x)))
@test integrate(a+b*cot(c+d*x)^2, x) == :(a*x+-1*b*x+-1*b*d^-1*cot(c+d*x))
@test integrate((a+b*cot(c+d*x)^2)^2, x) == :(x*(a+-1b)^2+-1//3*b^2*d^-1*cot(c+d*x)^3+-1*b*d^-1*(-1b+2a)*cot(c+d*x))
@test integrate((a+b*cot(c+d*x)^2)^3, x) == :(x*(a+-1b)^3+-1//5*b^3*d^-1*cot(c+d*x)^5+-1*b*d^-1*(b^2+3*a^2+-3*a*b)*cot(c+d*x)+-1//3*b^2*d^-1*cot(c+d*x)^3*(-1b+3a))
@test integrate((a+b*cot(c+d*x)^2)^-1, x) == :(x*(a+-1b)^-1+a^-1//2*b^(1/2)*d^-1*(a+-1b)^-1*arctan(a^-1//2*b^(1/2)*cot(c+d*x)))
@test integrate((a+b*cot(c+d*x)^2)^-2, x) == :(x*(a+-1b)^-2+(1/2)*b*a^-1*d^-1*(a+-1b)^-1*(a+b*cot(c+d*x)^2)^-1*cot(c+d*x)+(1/2)*a^-3//2*b^(1/2)*d^-1*(a+-1b)^-2*(-1b+3a)*arctan(a^-1//2*b^(1/2)*cot(c+d*x)))
@test integrate((a+b*cot(c+d*x)^2)^-3, x) == :(x*(a+-1b)^-3+1//4*b*a^-1*d^-1*(a+-1b)^-1*(a+b*cot(c+d*x)^2)^-2*cot(c+d*x)+1//8*a^-5//2*b^(1/2)*d^-1*(a+-1b)^-3*(3*b^2+15*a^2+-10*a*b)*arctan(a^-1//2*b^(1/2)*cot(c+d*x))+1//8*b*a^-2*d^-1*(a+-1b)^-2*(a+b*cot(c+d*x)^2)^-1*(-3b+7a)*cot(c+d*x))
@test integrate((1+cot(x)^2)^3//2, x) == :(-1//2*arcsinh(cot(x))+-1//2*(csc(x)^2)^(1/2)*cot(x))
@test integrate((1+cot(x)^2)^(1/2), x) == :(-1*arcsinh(cot(x)))
@test integrate((1+cot(x)^2)^-1//2, x) == :(-1*(csc(x)^2)^-1//2*cot(x))
@test integrate((-1+-1*cot(x)^2)^3//2, x) == :(-1//2*arctan((-1*csc(x)^2)^-1//2*cot(x))+(1/2)*(-1*csc(x)^2)^(1/2)*cot(x))
@test integrate((-1+-1*cot(x)^2)^(1/2), x) == :(arctan((-1*csc(x)^2)^-1//2*cot(x)))
@test integrate((-1+-1*cot(x)^2)^-1//2, x) == :(-1*(-1*csc(x)^2)^-1//2*cot(x))
@test integrate((a+a*cot(x)^2)^-1//2*cot(x)^3, x) == :(-1*(a*csc(x)^2)^-1//2+-1*a^-1*(a*csc(x)^2)^(1/2))
@test integrate((a+a*cot(x)^2)^-1//2*cot(x)^2, x) == :((a*csc(x)^2)^-1//2*cot(x)+-1*(a*csc(x)^2)^-1//2*arctanh(cos(x))*csc(x))
@test integrate((a+a*cot(x)^2)^-1//2*cot(x), x) == :((a*csc(x)^2)^-1//2)
@test integrate((a+a*cot(x)^2)^-1//2*tan(x), x) == :(-1*(a*csc(x)^2)^-1//2+a^-1//2*arctanh(a^-1//2*(a*csc(x)^2)^(1/2)))
@test integrate((a+a*cot(x)^2)^-1//2*tan(x)^2, x) == :((a*csc(x)^2)^-1//2*cot(x)+(a*csc(x)^2)^-1//2*csc(x)*sec(x))
@test integrate((a+b*cot(x)^2)^(1/2)*cot(x)^3, x) == :((a+b*cot(x)^2)^(1/2)+-1*(a+-1b)^(1/2)*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2))+-1//3*b^-1*(a+b*cot(x)^2)^3//2)
@test integrate((a+b*cot(x)^2)^(1/2)*cot(x), x) == :(-1*(a+b*cot(x)^2)^(1/2)+(a+-1b)^(1/2)*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2)))
@test integrate((a+b*cot(x)^2)^(1/2)*tan(x), x) == :(a^(1/2)*arctanh(a^-1//2*(a+b*cot(x)^2)^(1/2))+-1*(a+-1b)^(1/2)*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2)))
@test integrate((a+b*cot(x)^2)^(1/2)*cot(x)^2, x) == :((a+-1b)^(1/2)*arctan((a+-1b)^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x))+-1//2*(a+b*cot(x)^2)^(1/2)*cot(x)+-1//2*b^-1//2*(a+-2b)*arctanh(b^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x)))
@test integrate((a+b*cot(x)^2)^(1/2), x) == :(-1*b^(1/2)*arctanh(b^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x))+-1*(a+-1b)^(1/2)*arctan((a+-1b)^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x)))
@test integrate((a+b*cot(x)^2)^(1/2)*tan(x)^2, x) == :((a+-1b)^(1/2)*arctan((a+-1b)^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x))+(a+b*cot(x)^2)^(1/2)*tan(x))
@test integrate((a+b*cot(x)^2)^(1/2)*tan(x)^4, x) == :(-1*(a+-1b)^(1/2)*arctan((a+-1b)^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x))+1//3*(a+b*cot(x)^2)^(1/2)*tan(x)^3+-1//3*a^-1*(a+b*cot(x)^2)^(1/2)*(-1b+3a)*tan(x))
@test integrate((a+b*cot(x)^2)^3//2*cot(x)^3, x) == :(1//3*(a+b*cot(x)^2)^3//2+(a+b*cot(x)^2)^(1/2)*(a+-1b)+-1*(a+-1b)^3//2*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2))+-1//5*b^-1*(a+b*cot(x)^2)^5//2)
@test integrate((a+b*cot(x)^2)^3//2*cot(x)^2, x) == :((a+-1b)^3//2*arctan((a+-1b)^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x))+-1*(a+b*cot(x)^2)^(1/2)*(-1//2*b+5//8*a)*cot(x)+-1//4*b*(a+b*cot(x)^2)^(1/2)*cot(x)^3+-1//8*b^-1//2*(3*a^2+8*b^2+-12*a*b)*arctanh(b^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x)))
@test integrate((a+b*cot(x)^2)^3//2*cot(x), x) == :(-1//3*(a+b*cot(x)^2)^3//2+(a+-1b)^3//2*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2))+-1*(a+b*cot(x)^2)^(1/2)*(a+-1b))
@test integrate((a+b*cot(x)^2)^3//2*tan(x), x) == :(a^3//2*arctanh(a^-1//2*(a+b*cot(x)^2)^(1/2))+-1*b*(a+b*cot(x)^2)^(1/2)+-1*(a+-1b)^3//2*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2)))
@test integrate((a+b*cot(x)^2)^3//2*tan(x)^2, x) == :((a+-1b)^3//2*arctan((a+-1b)^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x))+-1*b^3//2*arctanh(b^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x))+a*(a+b*cot(x)^2)^(1/2)*tan(x))
@test integrate((a+b*cot(c+d*x)^2)^5//2, x) == :(-1*d^-1*(a+-1b)^5//2*arctan((a+-1b)^(1/2)*(a+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))+-1//4*b*d^-1*(a+b*cot(c+d*x)^2)^3//2*cot(c+d*x)+-1//8*b^(1/2)*d^-1*(8*b^2+15*a^2+-20*a*b)*arctanh(b^(1/2)*(a+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))+-1//8*b*d^-1*(a+b*cot(c+d*x)^2)^(1/2)*(-4b+7a)*cot(c+d*x))
@test integrate((a+b*cot(c+d*x)^2)^3//2, x) == :(-1*d^-1*(a+-1b)^3//2*arctan((a+-1b)^(1/2)*(a+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))+-1//2*b*d^-1*(a+b*cot(c+d*x)^2)^(1/2)*cot(c+d*x)+-1//2*b^(1/2)*d^-1*(-2b+3a)*arctanh(b^(1/2)*(a+b*cot(c+d*x)^2)^-1//2*cot(c+d*x)))
@test integrate((a+b*cot(c+d*x)^2)^(1/2), x) == :(-1*b^(1/2)*d^-1*arctanh(b^(1/2)*(a+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))+-1*d^-1*(a+-1b)^(1/2)*arctan((a+-1b)^(1/2)*(a+b*cot(c+d*x)^2)^-1//2*cot(c+d*x)))
@test integrate((a+b*cot(c+d*x)^2)^-1//2, x) == :(-1*d^-1*(a+-1b)^-1//2*arctan((a+-1b)^(1/2)*(a+b*cot(c+d*x)^2)^-1//2*cot(c+d*x)))
@test integrate((a+b*cot(c+d*x)^2)^-3//2, x) == :(-1*d^-1*(a+-1b)^-3//2*arctan((a+-1b)^(1/2)*(a+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))+b*a^-1*d^-1*(a+-1b)^-1*(a+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))
@test integrate((a+b*cot(c+d*x)^2)^-5//2, x) == :(-1*d^-1*(a+-1b)^-5//2*arctan((a+-1b)^(1/2)*(a+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))+1//3*b*a^-1*d^-1*(a+-1b)^-1*(a+b*cot(c+d*x)^2)^-3//2*cot(c+d*x)+1//3*b*a^-2*d^-1*(a+-1b)^-2*(a+b*cot(c+d*x)^2)^-1//2*(-2b+5a)*cot(c+d*x))
@test integrate((a+b*cot(c+d*x)^2)^-7//2, x) == :(-1*d^-1*(a+-1b)^-7//2*arctan((a+-1b)^(1/2)*(a+b*cot(c+d*x)^2)^-1//2*cot(c+d*x))+1//5*b*a^-1*d^-1*(a+-1b)^-1*(a+b*cot(c+d*x)^2)^-5//2*cot(c+d*x)+1//15*b*a^-3*d^-1*(a+-1b)^-3*(a+b*cot(c+d*x)^2)^-1//2*(8*b^2+33*a^2+-26*a*b)*cot(c+d*x)+1//15*b*a^-2*d^-1*(a+-1b)^-2*(a+b*cot(c+d*x)^2)^-3//2*(-4b+9a)*cot(c+d*x))
@test integrate((1+-1*cot(x)^2)^3//2, x) == :(5//2*arcsin(cot(x))+(1/2)*(1+-1*cot(x)^2)^(1/2)*cot(x)+-2*2^(1/2)*arctan(2^(1/2)*(1+-1*cot(x)^2)^-1//2*cot(x)))
@test integrate((1+-1*cot(x)^2)^(1/2), x) == :(-1*2^(1/2)*arctan(2^(1/2)*(1+-1*cot(x)^2)^-1//2*cot(x))+arcsin(cot(x)))
@test integrate((1+-1*cot(x)^2)^-1//2, x) == :(-1//2*2^(1/2)*arctan(2^(1/2)*(1+-1*cot(x)^2)^-1//2*cot(x)))
@test integrate((-1+cot(x)^2)^3//2, x) == :(5//2*arctanh((-1+cot(x)^2)^-1//2*cot(x))+-2*2^(1/2)*arctanh(2^(1/2)*(-1+cot(x)^2)^-1//2*cot(x))+-1//2*(-1+cot(x)^2)^(1/2)*cot(x))
@test integrate((-1+cot(x)^2)^(1/2), x) == :(-1*arctanh((-1+cot(x)^2)^-1//2*cot(x))+2^(1/2)*arctanh(2^(1/2)*(-1+cot(x)^2)^-1//2*cot(x)))
@test integrate((-1+cot(x)^2)^-1//2, x) == :(-1//2*2^(1/2)*arctanh(2^(1/2)*(-1+cot(x)^2)^-1//2*cot(x)))
@test integrate((a+b*cot(x)^2)^-1//2*cot(x)^3, x) == :(-1*b^-1*(a+b*cot(x)^2)^(1/2)+-1*(a+-1b)^-1//2*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2)))
@test integrate((a+b*cot(x)^2)^-1//2*cot(x)^2, x) == :((a+-1b)^-1//2*arctan((a+-1b)^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x))+-1*b^-1//2*arctanh(b^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x)))
@test integrate((a+b*cot(x)^2)^-1//2*cot(x), x) == :((a+-1b)^-1//2*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2)))
@test integrate((a+b*cot(x)^2)^-1//2*tan(x), x) == :(a^-1//2*arctanh(a^-1//2*(a+b*cot(x)^2)^(1/2))+-1*(a+-1b)^-1//2*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2)))
@test integrate((a+b*cot(x)^2)^-1//2*tan(x)^2, x) == :((a+-1b)^-1//2*arctan((a+-1b)^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x))+a^-1*(a+b*cot(x)^2)^(1/2)*tan(x))
@test integrate((a+b*cot(x)^2)^-3//2*cot(x)^3, x) == :(-1*(a+-1b)^-3//2*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2))+a*b^-1*(a+-1b)^-1*(a+b*cot(x)^2)^-1//2)
@test integrate((a+b*cot(x)^2)^-3//2*cot(x)^2, x) == :((a+-1b)^-3//2*arctan((a+-1b)^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x))+-1*(a+-1b)^-1*(a+b*cot(x)^2)^-1//2*cot(x))
@test integrate((a+b*cot(x)^2)^-3//2*cot(x), x) == :((a+-1b)^-3//2*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2))+-1*(a+-1b)^-1*(a+b*cot(x)^2)^-1//2)
@test integrate((a+b*cot(x)^2)^-3//2*tan(x), x) == :(a^-3//2*arctanh(a^-1//2*(a+b*cot(x)^2)^(1/2))+-1*(a+-1b)^-3//2*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2))+b*a^-1*(a+-1b)^-1*(a+b*cot(x)^2)^-1//2)
@test integrate((a+b*cot(x)^2)^-3//2*tan(x)^2, x) == :((a+-1b)^-3//2*arctan((a+-1b)^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x))+b*a^-1*(a+-1b)^-1*(a+b*cot(x)^2)^-1//2*tan(x)+a^-2*(a+-1b)^-1*(a+b*cot(x)^2)^(1/2)*(a+-2b)*tan(x))
@test integrate((a+b*cot(x)^2)^-5//2*cot(x)^3, x) == :((a+-1b)^-2*(a+b*cot(x)^2)^-1//2+-1*(a+-1b)^-5//2*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2))+a*b^-1*(a+b*cot(x)^2)^-3//2*(-3b+3a)^-1)
@test integrate((a+b*cot(x)^2)^-5//2*cot(x)^2, x) == :((a+-1b)^-5//2*arctan((a+-1b)^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x))+-1*(a+b*cot(x)^2)^-3//2*(-3b+3a)^-1*cot(x)+-1//3*a^-1*(a+-1b)^-2*(a+b*cot(x)^2)^-1//2*(b+2a)*cot(x))
@test integrate((a+b*cot(x)^2)^-5//2*cot(x), x) == :((a+-1b)^-5//2*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2))+-1*(a+-1b)^-2*(a+b*cot(x)^2)^-1//2+-1*(a+b*cot(x)^2)^-3//2*(-3b+3a)^-1)
@test integrate((a+b*cot(x)^2)^-5//2*tan(x), x) == :(a^-5//2*arctanh(a^-1//2*(a+b*cot(x)^2)^(1/2))+-1*(a+-1b)^-5//2*arctanh((a+-1b)^-1//2*(a+b*cot(x)^2)^(1/2))+1//3*b*a^-1*(a+-1b)^-1*(a+b*cot(x)^2)^-3//2+b*a^-2*(a+-1b)^-2*(a+b*cot(x)^2)^-1//2*(-1b+2a))
@test integrate((a+b*cot(x)^2)^-5//2*tan(x)^2, x) == :((a+-1b)^-5//2*arctan((a+-1b)^(1/2)*(a+b*cot(x)^2)^-1//2*cot(x))+1//3*b*a^-1*(a+-1b)^-1*(a+b*cot(x)^2)^-3//2*tan(x)+1//3*b*a^-2*(a+-1b)^-2*(a+b*cot(x)^2)^-1//2*(-4b+7a)*tan(x)+1//3*a^-3*(a+-1b)^-2*(a+b*cot(x)^2)^(1/2)*(a+-4b)*(-2b+3a)*tan(x))
@test integrate((1+cot(x)^3)^-1, x) == :((1/2)*x+(1/2)*log(sin(x))+-1//6*log(1+cot(x))+1//3*log(1+cot(x)^2+-1*cot(x)))
@test integrate((a+b*cot(x)^4)^(1/2)*cot(x), x) == :(-1//2*(a+b*cot(x)^4)^(1/2)+(1/2)*b^(1/2)*arctanh(b^(1/2)*(a+b*cot(x)^4)^-1//2*cot(x)^2)+(1/2)*(a+b)^(1/2)*arctanh((a+b)^-1//2*(a+b*cot(x)^4)^-1//2*(a+-1*b*cot(x)^2)))
@test integrate((a+b*cot(x)^4)^3//2*cot(x), x) == :(-1//6*(a+b*cot(x)^4)^3//2+(1/2)*(a+b)^3//2*arctanh((a+b)^-1//2*(a+b*cot(x)^4)^-1//2*(a+-1*b*cot(x)^2))+-1*(a+b*cot(x)^4)^(1/2)*((1/2)*a+(1/2)*b+-1//4*b*cot(x)^2)+1//4*b^(1/2)*(2b+3a)*arctanh(b^(1/2)*(a+b*cot(x)^4)^-1//2*cot(x)^2))
@test integrate((a+b*cot(x)^4)^-1//2*cot(x), x) == :((1/2)*(a+b)^-1//2*arctanh((a+b)^-1//2*(a+b*cot(x)^4)^-1//2*(a+-1*b*cot(x)^2)))
@test integrate((a+b*cot(x)^4)^-3//2*cot(x), x) == :((1/2)*(a+b)^-3//2*arctanh((a+b)^-1//2*(a+b*cot(x)^4)^-1//2*(a+-1*b*cot(x)^2))+-1//2*a^-1*(a+b)^-1*(a+b*cot(x)^4)^-1//2*(a+b*cot(x)^2))
@test integrate((a+b*cot(x)^4)^-5//2*cot(x), x) == :((1/2)*(a+b)^-5//2*arctanh((a+b)^-1//2*(a+b*cot(x)^4)^-1//2*(a+-1*b*cot(x)^2))+-1//6*a^-1*(a+b)^-1*(a+b*cot(x)^4)^-3//2*(a+b*cot(x)^2)+-1//6*a^-2*(a+b)^-2*(a+b*cot(x)^4)^-1//2*(3*a^2+b*cot(x)^2*(2b+5a)))

max_stars_repo_path: test/4 Trig functions/4.4 Cotangent/4.4.7 (d trig)^m (a+b (c cot)^n)^p.jl


text:
import Pkg

function _generate_packagecompiler_install_content(config::Config)
    packagecompiler_installation_command = config.packagecompiler_installation_command
    lines = String[
        "empty!(Base.DEPOT_PATH)",
        "pushfirst!(Base.DEPOT_PATH, \"/opt/simplecontainergenerator_containers/packagecompiler_depot\")",
        "import Pkg",
        "$(packagecompiler_installation_command)",
    ]
    content = join(lines, "\n") * "\n"
    return content
end

ext: jl


text:
## B-splines, recursive definition
B0(u::Vector{Float64}, i::Integer) = x::Real -> (u[i] <= x < u[i+1]) ? float(1) : float(0) ## non-continuous
B(p::Integer,u::Vector{Float64},i::Integer) = (p==0) ? B0(u,i) : x -> float(((x-u[i])/(u[i+p]-u[i]))*B(p-1,u,i)(x) + ((u[i+p+1]-x)/(u[i+p+1]-u[i+1]))*B(p-1,u,i+1)(x))

totalNumNodes(numNodes::Integer, p::Integer) = numNodes + 2p ## adding p nodes on each side
numOfBasisElements(numNodes::Integer,p::Integer) = numNodes + p - 1 ## number of non-zero spline functions

## returns an index-function of basis-functions, which takes as input i in 1:k
## creates evenly-spaced nodes
function splineBasisFuns(p::Integer, range::(Real,Real), numNodes::Integer)
    intervalLen = (range[2]-range[1])/(numNodes-1)
    u = linspace(range[1] - p*intervalLen, range[2]+ p*intervalLen, totalNumNodes(numNodes,p))
    i -> B(p,u,i)
end

## B-splines with pre-specific nodes
## WARNING: all functions in this basis will be 0 at the leftmost and rightmost nodes, and the same is true of their first p-1 derivatives.
function splineBasisFuns(p::Integer, nodes::Vector{Float64})
    i::Integer -> B(p,nodes,i)
end

cubicSplineBasisFuns(range::(Real,Real), numNodes::Integer) = splineBasisFuns(3, range, numNodes)
quadraticSplineBasisFuns(range::(Real,Real), numNodes::Integer) = splineBasisFuns(2, range, numNodes)
linearSplineBasisFuns(range::(Real,Real), numNodes::Integer) = splineBasisFuns(1, range, numNodes)

## BFE "Basis Function Expansion"

## evenly-spaced nodes
function splineBFE(p::Integer, x::Vector{Float64}, range::(Real,Real), numNodes::Integer)
    funs = splineBasisFuns(p, range, numNodes)
    num = numOfBasisElements(numNodes,p)
    Float64[funs(i)(xj) for xj in x, i in 1:num]
end

## pre-specified nodes
function splineBFE(p::Integer, x::Vector{Float64}, nodes::Vector{Float64})
    funs = splineBasisFuns(p, nodes)
    n = length(nodes) - p - 1
    Float64[funs(i)(xj) for xj in x, i in 1:n]
end

cubicSplineBFE(x::Vector{Float64}, range::(Real,Real), numNodes::Integer) =
    splineBFE(3, x, range, numNodes)
quadraticSplineBFE(x::Vector{Float64}, range::(Real,Real), numNodes::Integer) =
    splineBFE(2, x, range, numNodes)
linearSplineBFE(x::Vector{Float64}, range::(Real,Real), numNodes::Integer) =
    splineBFE(1, x, range, numNodes)

cubicSplineBFE(x::Vector{Float64}, nodes::Vector{Float64}) = splineBFE(3, x, nodes)
quadraticSplineBFE(x::Vector{Float64}, nodes::Vector{Float64}) = splineBFE(2, x, nodes)
linearSplineBFE(x::Vector{Float64}, nodes::Vector{Float64}) = splineBFE(1, x, nodes)


text:
module UnitfulCosmo
# based heavily on Mason Protter's NaturallyUnitful.jl

using Reexport
@reexport using Unitful
@reexport using UnitfulAstro
import PhysicalConstants.CODATA2018: SpeedOfLightInVacuum, ReducedPlanckConstant,
    BoltzmannConstant, NewtonianConstantOfGravitation, VacuumElectricPermittivity

const c  = float(SpeedOfLightInVacuum)
const ħ  = float(ReducedPlanckConstant)
const G = float(NewtonianConstantOfGravitation)
const kB = float(BoltzmannConstant)
const ε0 = float(VacuumElectricPermittivity)

setDimPow(D::Dict, x::Unitful.Dimension{T}) where {T} = D[T] = x.power

function dimDict(x::Unitful.Dimensions{T}) where {T}
    D = Dict(:Length => 0,
             :Mass => 0,
             :Temperature => 0,
             :Time => 0,
             :Current => 0)
    for t in T
        setDimPow(D,t)
    end
    return D
end
dimDict(x::Unitful.Quantity{T}) where {T} = dimDict(dimension(x))
dimDict(x) = Dict(:Length => 0, :Mass => 0, :Temperature => 0, :Time => 0, :Current => 0)

"""
    mpc(q; base=u"Mpc")

Convert `q` to Mpc units based on the units specified by `base`. If `base` is
unspecified, `mpc` will default to `Mpc`. You could choose some other length unit
if you wanted.
"""
mpc(q; base=u"Mpc") = _mpc(base, q)

function _mpc(base::Unitful.FreeUnits{T, Unitful.𝐋, U}, q) where {T, U}
    D = dimDict(q)
    (α,β,γ,δ,ϕ) = (D[:Length], D[:Mass], D[:Temperature], D[:Time], D[:Current])
    uconvert(base^(α+δ-β-γ-ϕ),
            q * (c/ħ)^β * (kB/(c*ħ))^γ * c^δ * (c^3*ε0*ħ)^(-ϕ/2))
end

function _mpc(base::Unitful.FreeUnits{T, U, V}, q) where {T, U, V}
    throw(ArgumentError(
        """mpc(q; base)` where `base` has dimensions `$(U)` has not yet been implemented. Please use a base with dimensions of `length` and then use `unmpc` to convert to your desired units."""
        )
    )
end

"""
    unmpc(targetUnit, q)

Convert a quantity `q` to units specified by `targetUnits` automatically inserting whatever
factors  of `ħ`, `c`, `ϵ₀` and `kb` are required to make the conversion work.
"""
function unmpc(targetUnit, q)
    natTarget = mpc(1targetUnit)
    natQ = mpc(q)
    ratio = natQ/natTarget
    if typeof(ratio |> dimension) == Unitful.Dimensions{()}
        (ratio)targetUnit
    else
        throw(Unitful.DimensionError(q,targetUnit))
    end
end

"""
    planck(q)

Strip units by inserting whatever `ħ`, `c`, `G` and `kb` are required to make the conversion work.
"""
function planck(q)
    D = dimDict(q)
    (α,β,γ,δ,ϕ) = (D[:Length], D[:Mass], D[:Temperature], D[:Time], D[:Current])
    return uconvert(Unitful.NoUnits,
        q *
        (G * ħ / c^3)^(-α/2) *  # length
        (c * ħ / G)^(-β/2) *    # mass
        (c^5 * ħ / (G * kB^2))^(-γ/2) *  # temperature
        (G * ħ / c^5)^(ϕ/2 - δ/2) *  # time and current
        (4π * ε0 * ħ * c)^(-ϕ/2))  # charge
end

"""
    unplanck(targetUnit, q)

Convert a quantity `q` to units specified by `targetUnits` automatically inserting whatever
factors  of `ħ`, `c`, `G` and `kb` are required to make the conversion work.
"""
function unplanck(targetUnit, q)
    if unit(q) != unit(1)
        throw(Unitful.DimensionError(q,Unitful.NoUnits))
    end
    natTarget = planck(1targetUnit)
    natQ = planck(q)
    ratio = natQ/natTarget
    (ratio)targetUnit
end

export planck, unplanck
export mpc, unmpc

end


text:
"Nearest interpolation type"
abstract type Nearest <: InterpMode end

function interp(xi::Vector{T}, yi::VecOrMat{T}, ::Type{Nearest};
                copy=true, sorted=false) where T
    _xi, _yi = _prepare_input(xi, yi, copy, sorted)

    function f(x)
        idx = _searchsortednearest(_xi, x)
        return ndims(_yi) == 1 ? _yi[idx] : _yi[:, idx]
    end

    return f
end

function _searchsortednearest(a,x)
   idx = searchsortedfirst(a,x)
   idx==1 && return idx
   idx>length(a) && return length(a)
   a[idx]==x && return idx
   if abs(a[idx]-x) < abs(a[idx-1]-x)
      return idx
   else
      return idx-1
   end
end


text:
println("Loading libraries...")
using NeRCA
using DrWatson
using ProgressMeter

struct RecoFile{T}
    filepath::AbstractString
    _fobj

    function RecoFile{T}(filepath::AbstractString) where T
        if isfile(filepath)
            @warn "Reconstruction file '$filepath' present, overwriting..."
        end
        fobj = open(filepath, "w")
        write(fobj, join(["event_id", "du", "Q", "n_triggered_doms", "n_hits", "n_triggered_hits", fieldnames(T)...], ",") * "\n")
        new(filepath, fobj)
    end
end

close(f::RecoFile) = close(f._fobj)

function Base.write(f::RecoFile, event_id, du, Q, n_triggered_doms, n_hits, n_triggered_hits, s::SingleDUParams)
    write(f._fobj, "$event_id,$du,$Q,$n_triggered_doms,$n_hits,$n_triggered_hits,")
    write(f._fobj, join([getfield(s, field) for field in fieldnames(typeof(s))], ","))
    write(f._fobj, "\n")
end

if length(ARGS) < 3
    println("Usage: ./single_du_fit.jl OUTPATH DETX ROOTFILE")
    exit(1)
end

function main()
    println("Starting reconstruction.")

    sparams = SingleDURecoParams(max_iter=200)

    outpath = ARGS[1]
    detx = ARGS[2]
    rootfile = ARGS[3]
    outfile = joinpath(outpath, basename(rootfile) * ".ROyFit.csv")

    mkpath(outpath)
    recofile = RecoFile{SingleDUParams}(outfile)

    calib = Calibration(detx)
    f = NeRCA.OnlineFile(rootfile)

    event_shits = NeRCA.read_snapshot_hits(f)
    event_thits = NeRCA.read_triggered_hits(f)
    n_events = length(event_shits)
    println("$n_events events found")

    @showprogress for (event_id, shits, thits) in zip(1:n_events, event_shits, event_thits))
        hits = calibrate(calib, NeRCA.combine(shits, thits))

        triggered_hits = triggered(hits)

        dus = sort(unique(map(h->h.du, hits)))
        triggered_dus = sort(unique(map(h->h.du, triggered_hits)))
        n_dus = length(dus)
        n_triggered_dus = length(triggered_dus)

        for du in dus
            du_hits = filter(h->h.du == du, hits)
            triggered_du_hits = triggered(du_hits)
            n_triggered_hits = length(triggered_du_hits)
            n_triggered_doms = length(unique(h->h.dom_id, triggered_du_hits))
            if n_triggered_hits == 0
                continue
            end
            fit = NeRCA.single_du_fit(du_hits, sparams)
            write(recofile, event_id, du, fit.Q, n_triggered_doms, length(du_hits), n_triggered_hits, fit.sdp)
        end
    end
end

main()


text:
using Metalhead, Test
using Flux

@testset "ViT" begin
    for mode in [:small, :base, :large] # :tiny, #,:huge, :giant, :gigantic]
        m = ViT(mode)
        @test size(m(x_256)) == (1000, 1)
        @test gradtest(m, x_256)
        GC.safepoint()
        GC.gc()
    end
end


text:
using SpecialFunctions

# from datetime import date
# User-provided date:
#
stosym = "wmt"
exyear = 2019
exmonth = 5
exday = 16
stopr = 92.68
strike = 94
callBid = 0.97
callAsk = 1.01
rfir = 0.0225
#
#  Initialize key variables below
#  (This is not actually required. This is your Prof's habit).
#
daystd = 1
cipd = 0.0
d1 = 0.0
durvol = 0.0
cumd1 = 0.0
cumd2 = 0.0
newcp = 0.0
tempcp = 0.0
timedecay = 0.0
tdprice = 0.0
callpr = 0.0
#
# Calculating strike value (PEG), which is somewhere between BB and BA.
# PEG coefficient is explained in lecture.
peg_coef = 0.60
spread = callAsk - callBid
callpr = callBid + (peg_coef*spread)
# Calculating days to expiry:
tnow = Dates.today()
expiry = Dates.Date(exyear, exmonth, exday)
days = Int(expiry - tnow)
#  Our method for calculating the cumulative standard normal distribution:
function csnd(dval)
	return (1.0 + erf(dval/sqrt(2.0)))/2.0
end
#
#  Below we calculate implied daily volatility using a recursive method to
#  converge to an answer.  First we make an initializing calculation, then
#  we converge inside the while loop. The convergence method shown is Prof E's
#  "Divide and Conquer", which is relatively efficient, usually using only 8
#  (+/- 2) steps. The student should recognize that the cumd1 is the delta.
#  The implied daily volatility that we are seeking is "cipd" below.
#
target = callpr
precision = 1e-4
count = 0
low = 0.0
high = 1.0
cipd = (high+low)/2
d1 = log(stopr/strike)+((rfir/365)+(cipd^2)/2)*days
durvol = cipd*sqrt(days)
cumd1 = csnd(d1/durvol)
cumd2 = csnd((d1/durvol) - durvol)
discount = exp(-rfir*days/365)
tempcp = (stopr*cumd1)-(strike*discount*cumd2)
while tempcp<=(target-precision) || tempcp>=(target+precision)
	if tempcp >= (target+precision)
		high = cipd
	else
		low = cipd
    end
	cipd = (high+low)/2
	d1 = log(stopr/strike)+((rfir/365)+(cipd^2)/2)*days
	durvol = cipd*sqrt(days)
	cumd1 = csnd(d1/durvol)
	cumd2 = csnd((d1/durvol) - durvol)
	discount = exp(-rfir*days/365)
	tempcp = (stopr*cumd1)-(strike*discount*cumd2)
	count +=1
end
#	Below we calculate one day time decay using our new value for volatility
days = days - daystd
d1 = log(stopr/strike)+((rfir/365)+(cipd^2)/2)*days
durvol = cipd*sqrt(days)
cumd1 = csnd(d1/durvol)
cumd2 = csnd((d1/durvol) - durvol)
discount = exp(-rfir*days/365)
newcp = (stopr*cumd1)-(strike*discount*cumd2)
timedecay = callpr - newcp

println( "")
println( "Days to expiry: ", days)
println( "Stock price: ", stopr)
println( "Strike price: ", strike)
println( "Call ASK: ", callAsk)
println( "Call BID: ", callBid)
println( "Call price (PEG): ", callpr)
println( "The Delta:", cumd1)
println( "One day time decay:", timedecay)
println( "The call's implied volatility: ", cipd)


text:
function timestep()
  dt = dt_max
  for i in active_particle
    @assert sph[i].P > 0

    s = index_list[i]
    vmax = 0
    for t = (s - Nngb) : (s + Nngb)
      t = order_periodic(t)
      k = index_order[t]
      vsig = cs(i) + cs(k) - 3 * eij(i, k) * vij(i, k)
      if vsig > vmax
        vmax = vsig
      end
    end
    dt_i = cfl * sph[i].hsml / vmax
    @assert dt_i > 0
    dt = min(dt, dt_i)
  end
  return dt
end

function update_time(t, dt)
  global t += dt
end


text:
function __init__()
	register(
		DataDep(
			"shapes2D",
			"""
			Dataset: 2D shapes
			Authors: Loic Matthey and Irina Higgins and Demis Hassabis and Alexander Lerchner
			Website: https://github.com/deepmind/dsprites-dataset
			""",
			[
				"https://github.com/deepmind/dsprites-dataset/raw/master/dsprites_ndarray_co1sh3sc6or40x32y32_64x64.npz"
			],
			"857dc001287e8efa226e0025fac6b29e5745c621b380b0d3363e4e33f4e28fa2"
		))
end

"""
	symmetric_angle(shape, angle)

Returns a normalized rotation angle that is dependent on the input shape and the number
of its symmetry axes.
"""
function symmetric_angle(shape, angle)
    if shape == 1 # square
        nsym = 4
    elseif shape == 2 # oval
        nsym = 2
    elseif shape == 3 # heart
        nsym = 1
    else
        error("Unknown shape code $shape.")
    end
    symfac = 360/nsym
    return (angle%symfac)/symfac
end
function _preprocess_shapes2D()
	@info "Preprocessing 2D shapes data..."
	f = joinpath(datadep"shapes2D", "dsprites_ndarray_co1sh3sc6or40x32y32_64x64.npz")
	np = pyimport("numpy")
	data = np.load(f, allow_pickle=true, encoding="latin1")

	# image data
	x = get(data, "imgs")
	n = size(x,1)
	imax = 400000
	X1 = similar(x, 64, 64, 1, imax)
	X2 = similar(x, 64, 64, 1, n-imax)
	for i in 1:imax
		X1[:,:,1,i] .= x[i,:,:]
	end
	for i in imax+1:n
		X2[:,:,1,(i-imax)] .= x[i,:,:]
	end
	outf1 = datadir("shapes2D/shapes_part1.bson")
	outf2 = datadir("shapes2D/shapes_part2.bson")
	save(outf1, Dict(:imgs=>X1))
	save(outf2, Dict(:imgs=>X2))

	# labels
	factor_names = ["color", "shape", "scale", "orientation", "posX", "posY"]
	y = get(data, "latents_values")
	y[:,4] .= y[:,4]/(2*pi)*360 #rads to degrees
	df = DataFrame(y)
	map(x->rename!(df, x[1]=>Symbol(x[2])), zip(names(df), factor_names))
	# now add rotation factor that is different for each shape based on its symmetry
	df[!,:normalized_orientation] = symmetric_angle.(df[!,:shape], df[!,:orientation])
	outf = datadir("shapes2D/labels.csv")
	CSV.write(outf, df)
	@info "Done."
end
function load_shapes2D()
	# load X
	infs = map(x->datadir("shapes2D/shapes_part$x.bson"), [1,2])
	if !all(isfile.(infs))
		_preprocess_shapes2D()
	end
	Xs = map(load, infs)
	X = Float32.(cat([X[:imgs] for X in Xs]..., dims=4))

	# load labels
	f = datadir("shapes2D/labels.csv")
	if !isfile(f)
		_preprocess_shapes2D()
	end
	y = CSV.read(f, DataFrame)

	return X,y
end

max_stars_repo_name: vitskvara/HierarchicalAD

max_issues_repo_name: vitskvara/HierarchicalAD

max_forks_repo_name: vitskvara/HierarchicalAD


text:
@testset "unwind function" begin
    iterators = ([1, 2], ["a","b"], ["x", "y", "z"])
    @test MLJTuning.unwind(iterators...) ==
        [1  "a"  "x";
         2  "a"  "x";
         1  "b"  "x";
         2  "b"  "x";
         1  "a"  "y";
         2  "a"  "y";
         1  "b"  "y";
         2  "b"  "y";
         1  "a"  "z";
         2  "a"  "z";
         1  "b"  "z";
         2  "b"  "z"]
end

true


text:
real_rep(a::AbstractArray{Complex{T}, N}) where {T, N} = reinterpret(reshape, T, a)
#PtrArray(Ptr{T}(pointer(a)), (StaticInt(2), size(a)...))

@inline function _matmul!(_C::AbstractMatrix{Complex{T}}, _A::AbstractMatrix{Complex{U}}, _B::AbstractMatrix{Complex{V}},
                         α=One(), β=Zero(), nthread::Nothing=nothing, MKN=nothing, contig_axis=nothing) where {T,U,V}
    C, A, B =  real_rep.((_C, _A, _B))

    η = ifelse(ArrayInterface.is_lazy_conjugate(_A), StaticInt(-1), StaticInt(1))
    θ = ifelse(ArrayInterface.is_lazy_conjugate(_B), StaticInt(-1), StaticInt(1))
    (+ᶻ, -ᶻ) = ifelse(ArrayInterface.is_lazy_conjugate(_C), (-, +), (+, -))
    ηθ = η*θ

    @avxt for n ∈ indices((C, B), 3), m ∈ indices((C, A), 2)
        Cmn_re = zero(T)
        Cmn_im = zero(T)
        for k ∈ indices((A, B), (3, 2))
            Cmn_re +=     A[1, m, k] * B[1, k, n] - ηθ * A[2, m, k] * B[2, k, n]
            Cmn_im += θ * A[1, m, k] * B[2, k, n] +  η * A[2, m, k] * B[1, k, n]
        end
        C[1,m,n] = (real(α) * Cmn_re -ᶻ imag(α) * Cmn_im) + (real(β) * C[1,m,n] -ᶻ imag(β) * C[2,m,n])
        C[2,m,n] = (imag(α) * Cmn_re +ᶻ real(α) * Cmn_im) + (imag(β) * C[1,m,n] +ᶻ real(β) * C[2,m,n])
    end
    _C
end

@inline function _matmul!(_C::AbstractMatrix{Complex{T}}, A::AbstractMatrix{U}, _B::AbstractMatrix{Complex{V}},
                         α=One(), β=Zero(), nthread::Nothing=nothing, MKN=nothing, contig_axis=nothing) where {T,U,V}
    C, B = real_rep.((_C, _B))

    θ = ifelse(ArrayInterface.is_lazy_conjugate(_B), StaticInt(-1), StaticInt(1))
    (+ᶻ, -ᶻ) = ifelse(ArrayInterface.is_lazy_conjugate(_C), (-, +), (+, -))

    @avxt for n ∈ indices((C, B), 3), m ∈ indices((C, A), (2, 1))
        Cmn_re = zero(T)
        Cmn_im = zero(T)
        for k ∈ indices((A, B), (2, 2))
            Cmn_re +=     A[m, k] * B[1, k, n]
            Cmn_im += θ * A[m, k] * B[2, k, n]
        end
        C[1,m,n] = (real(α) * Cmn_re -ᶻ imag(α) * Cmn_im) + (real(β) * C[1,m,n] -ᶻ imag(β) * C[2,m,n])
        C[2,m,n] = (imag(α) * Cmn_re +ᶻ real(α) * Cmn_im) + (imag(β) * C[1,m,n] +ᶻ real(β) * C[2,m,n])
    end
    _C
end

@inline function _matmul!(_C::AbstractMatrix{Complex{T}}, _A::AbstractMatrix{Complex{U}}, B::AbstractMatrix{V},
                         α=One(), β=Zero(), nthread::Nothing=nothing, MKN=nothing, contig_axis=nothing) where {T,U,V}
    C, A = real_rep.((_C, _A))

    η = ifelse(ArrayInterface.is_lazy_conjugate(_A), StaticInt(-1), StaticInt(1))
    (+ᶻ, -ᶻ) = ifelse(ArrayInterface.is_lazy_conjugate(_C), (-, +), (+, -))

    @avxt for n ∈ indices((C, B), (3, 2)), m ∈ indices((C, A), 2)
        Cmn_re = zero(T)
        Cmn_im = zero(T)
        for k ∈ indices((A, B), (3, 1))
            Cmn_re +=     A[1, m, k] * B[k, n]
            Cmn_im += η * A[2, m, k] * B[k, n]
        end
        C[1,m,n] = (real(α) * Cmn_re -ᶻ imag(α) * Cmn_im) + (real(β) * C[1,m,n] -ᶻ imag(β) * C[2,m,n])
        C[2,m,n] = (imag(α) * Cmn_re +ᶻ real(α) * Cmn_im) + (imag(β) * C[1,m,n] +ᶻ real(β) * C[2,m,n])
    end
    _C
end

@inline function _matmul_serial!(_C::AbstractMatrix{Complex{T}}, _A::AbstractMatrix{Complex{U}}, _B::AbstractMatrix{Complex{V}},
                         α=One(), β=Zero(), MKN=nothing, contig_axis=nothing) where {T,U,V}
    C, A, B = real_rep.((_C, _A, _B))

    η = ifelse(ArrayInterface.is_lazy_conjugate(_A), StaticInt(-1), StaticInt(1))
    θ = ifelse(ArrayInterface.is_lazy_conjugate(_B), StaticInt(-1), StaticInt(1))
    (+ᶻ, -ᶻ) = ifelse(ArrayInterface.is_lazy_conjugate(_C), (-, +), (+, -))
    ηθ = η*θ
    @avxt for n ∈ indices((C, B), 3), m ∈ indices((C, A), 2)
        Cmn_re = zero(T)
        Cmn_im = zero(T)
        for k ∈ indices((A, B), (3, 2))
            Cmn_re +=     A[1, m, k] * B[1, k, n] - ηθ * A[2, m, k] * B[2, k, n]
            Cmn_im += θ * A[1, m, k] * B[2, k, n] +  η * A[2, m, k] * B[1, k, n]
        end
        C[1,m,n] = (real(α) * Cmn_re -ᶻ imag(α) * Cmn_im) + (real(β) * C[1,m,n] -ᶻ imag(β) * C[2,m,n])
        C[2,m,n] = (imag(α) * Cmn_re +ᶻ real(α) * Cmn_im) + (imag(β) * C[1,m,n] +ᶻ real(β) * C[2,m,n])
    end
    _C
end

@inline function _matmul_serial!(_C::AbstractMatrix{Complex{T}}, A::AbstractMatrix{U}, _B::AbstractMatrix{Complex{V}},
                         α=One(), β=Zero(), MKN=nothing, contig_axis=nothing) where {T,U,V}
    C, B = real_rep.((_C, _B))

    θ = ifelse(ArrayInterface.is_lazy_conjugate(_B), StaticInt(-1), StaticInt(1))
    (+ᶻ, -ᶻ) = ifelse(ArrayInterface.is_lazy_conjugate(_C), (-, +), (+, -))

    @avx for n ∈ indices((C, B), 3), m ∈ indices((C, A), (2, 1))
        Cmn_re = zero(T)
        Cmn_im = zero(T)
        for k ∈ indices((A, B), (2, 2))
            Cmn_re +=     A[m, k] * B[1, k, n]
            Cmn_im += θ * A[m, k] * B[2, k, n]
        end
        C[1,m,n] = (real(α) * Cmn_re -ᶻ imag(α) * Cmn_im) + (real(β) * C[1,m,n] -ᶻ imag(β) * C[2,m,n])
        C[2,m,n] = (imag(α) * Cmn_re +ᶻ real(α) * Cmn_im) + (imag(β) * C[1,m,n] +ᶻ real(β) * C[2,m,n])
    end
    _C
end

@inline function _matmul_serial!(_C::AbstractMatrix{Complex{T}}, _A::AbstractMatrix{Complex{U}}, B::AbstractMatrix{V},
                         α=One(), β=Zero(), MKN=nothing, contig_axis=nothing) where {T,U,V}
    C, A = real_rep.((_C, _A))

    η = ifelse(ArrayInterface.is_lazy_conjugate(_A), StaticInt(-1), StaticInt(1))
    (+ᶻ, -ᶻ) = ifelse(ArrayInterface.is_lazy_conjugate(_C), (-, +), (+, -))

    @avx for n ∈ indices((C, B), (3, 2)), m ∈ indices((C, A), 2)
        Cmn_re = zero(T)
        Cmn_im = zero(T)
        for k ∈ indices((A, B), (3, 1))
            Cmn_re +=     A[1, m, k] * B[k, n]
            Cmn_im += η * A[2, m, k] * B[k, n]
        end
        C[1,m,n] = (real(α) * Cmn_re -ᶻ imag(α) * Cmn_im) + (real(β) * C[1,m,n] -ᶻ imag(β) * C[2,m,n])
        C[2,m,n] = (imag(α) * Cmn_re +ᶻ real(α) * Cmn_im) + (imag(β) * C[1,m,n] +ᶻ real(β) * C[2,m,n])
    end
    _C
end


text: sumaverage(arr) = map(mean,arr) |> sum |> floor

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
"""
struct ArrayModel{T <: MillFunction} <: MillModel
m::T
end

use a Chain, Dense, or any other function on an ArrayNode
"""
struct ArrayModel{T} <: MillModel
    m::T
end

Flux.@functor ArrayModel

(m::ArrayModel)(x::ArrayNode) = mapdata(x -> m.m(x), x)

modelprint(io::IO, m::ArrayModel; pad=[], s="", tr=false) = paddedprint(io, "ArrayModel(", m.m, ")$(tr_repr(s, tr))")

function HiddenLayerModel(m::ArrayModel, x::ArrayNode, k::Int)
	os = Flux.activations(m.m, x.data)
	layers = map(x -> Dense(size(x,1), k), os)
	ArrayModel(layers), ArrayNode(os[end])
end

function mapactivations(hm::ArrayModel, x::ArrayNode, m::ArrayModel)
	os = Flux.activations(m.m, x.data)
	hx = mapfoldl((mx) -> mx[1](mx[2]),+,zip(hm.m, os))
	(ArrayNode(hx), ArrayNode(os[end]))
end

Flux.activations(::typeof(identity), x::Array{Float32,2}) = (x,)


text:
#!/usr/bin/env julia

#=
Copyright (c) 2016, Matthew Morris

Permission to use, copy, modify, and/or distribute this software for any
purpose with or without fee is hereby granted, provided that the above
copyright notice and this permission notice appear in all copies.

THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
=#

fizzbuzz(num::Int, mod3::Bool, mod5::Bool) = num
fizzbuzz(num::Int, mod3::Int, mod5::Bool) = "fizz"
fizzbuzz(num::Int, mod3::Bool, mod5::Int) = "buzz"
fizzbuzz(num::Int, mod3::Int, mod5::Int) = "fizzbuzz"

for n = 1:100
    println(fizzbuzz(n, n % 3 == 0 ? 1 : false, n % 5 == 0 ? 1 : false))
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct FlowLogProperties <: SwaggerModel
    storageId::Any # spec type: Union{ Nothing, String } # spec name: storageId
    enabled::Any # spec type: Union{ Nothing, Bool } # spec name: enabled
    retentionPolicy::Any # spec type: Union{ Nothing, RetentionPolicyParameters } # spec name: retentionPolicy
    format::Any # spec type: Union{ Nothing, FlowLogFormatParameters } # spec name: format

    function FlowLogProperties(;storageId=nothing, enabled=nothing, retentionPolicy=nothing, format=nothing)
        o = new()
        validate_property(FlowLogProperties, Symbol("storageId"), storageId)
        setfield!(o, Symbol("storageId"), storageId)
        validate_property(FlowLogProperties, Symbol("enabled"), enabled)
        setfield!(o, Symbol("enabled"), enabled)
        validate_property(FlowLogProperties, Symbol("retentionPolicy"), retentionPolicy)
        setfield!(o, Symbol("retentionPolicy"), retentionPolicy)
        validate_property(FlowLogProperties, Symbol("format"), format)
        setfield!(o, Symbol("format"), format)
        o
    end
end # type FlowLogProperties

const _property_map_FlowLogProperties = Dict{Symbol,Symbol}(Symbol("storageId")=>Symbol("storageId"), Symbol("enabled")=>Symbol("enabled"), Symbol("retentionPolicy")=>Symbol("retentionPolicy"), Symbol("format")=>Symbol("format"))
const _property_types_FlowLogProperties = Dict{Symbol,String}(Symbol("storageId")=>"String", Symbol("enabled")=>"Bool", Symbol("retentionPolicy")=>"RetentionPolicyParameters", Symbol("format")=>"FlowLogFormatParameters")
Base.propertynames(::Type{ FlowLogProperties }) = collect(keys(_property_map_FlowLogProperties))
Swagger.property_type(::Type{ FlowLogProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_FlowLogProperties[name]))}
Swagger.field_name(::Type{ FlowLogProperties }, property_name::Symbol) =  _property_map_FlowLogProperties[property_name]

function check_required(o::FlowLogProperties)
    (getproperty(o, Symbol("storageId")) === nothing) && (return false)
    (getproperty(o, Symbol("enabled")) === nothing) && (return false)
    true
end

function validate_property(::Type{ FlowLogProperties }, name::Symbol, val)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct VirtualMachineExtensionImagesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets a virtual machine extension image.
Param: location::String (required)
Param: publisherName::String (required)
Param: type::String (required)
Param: version::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualMachineExtensionImage
"""
function _swaggerinternal_virtualMachineExtensionImagesGet(_api::VirtualMachineExtensionImagesApi, location::String, publisherName::String, type::String, version::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualMachineExtensionImage, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/publishers/{publisherName}/artifacttypes/vmextension/types/{type}/versions/{version}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "publisherName", publisherName)  # type String
    Swagger.set_param(_ctx.path, "type", type)  # type String
    Swagger.set_param(_ctx.path, "version", version)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualMachineExtensionImagesGet(_api::VirtualMachineExtensionImagesApi, location::String, publisherName::String, type::String, version::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineExtensionImagesGet(_api, location, publisherName, type, version, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualMachineExtensionImagesGet(_api::VirtualMachineExtensionImagesApi, response_stream::Channel, location::String, publisherName::String, type::String, version::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineExtensionImagesGet(_api, location, publisherName, type, version, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets a list of virtual machine extension image types.
Param: location::String (required)
Param: publisherName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Vector{VirtualMachineExtensionImage}
"""
function _swaggerinternal_virtualMachineExtensionImagesListTypes(_api::VirtualMachineExtensionImagesApi, location::String, publisherName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Vector{VirtualMachineExtensionImage}, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/publishers/{publisherName}/artifacttypes/vmextension/types", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "publisherName", publisherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualMachineExtensionImagesListTypes(_api::VirtualMachineExtensionImagesApi, location::String, publisherName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineExtensionImagesListTypes(_api, location, publisherName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualMachineExtensionImagesListTypes(_api::VirtualMachineExtensionImagesApi, response_stream::Channel, location::String, publisherName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineExtensionImagesListTypes(_api, location, publisherName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets a list of virtual machine extension image versions.
Param: location::String (required)
Param: publisherName::String (required)
Param: type::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: filter::String
Param: top::Int32
Param: orderby::String
Return: Vector{VirtualMachineExtensionImage}
"""
function _swaggerinternal_virtualMachineExtensionImagesListVersions(_api::VirtualMachineExtensionImagesApi, location::String, publisherName::String, type::String, api_version::String, subscriptionId::String; filter=nothing, top=nothing, orderby=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Vector{VirtualMachineExtensionImage}, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/publishers/{publisherName}/artifacttypes/vmextension/types/{type}/versions", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "publisherName", publisherName)  # type String
    Swagger.set_param(_ctx.path, "type", type)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$filter", filter)  # type String
    Swagger.set_param(_ctx.query, "\$top", top)  # type Int32
    Swagger.set_param(_ctx.query, "\$orderby", orderby)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualMachineExtensionImagesListVersions(_api::VirtualMachineExtensionImagesApi, location::String, publisherName::String, type::String, api_version::String, subscriptionId::String; filter=nothing, top=nothing, orderby=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineExtensionImagesListVersions(_api, location, publisherName, type, api_version, subscriptionId; filter=filter, top=top, orderby=orderby, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualMachineExtensionImagesListVersions(_api::VirtualMachineExtensionImagesApi, response_stream::Channel, location::String, publisherName::String, type::String, api_version::String, subscriptionId::String; filter=nothing, top=nothing, orderby=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_virtualMachineExtensionImagesListVersions(_api, location, publisherName, type, api_version, subscriptionId; filter=filter, top=top, orderby=orderby, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export virtualMachineExtensionImagesGet, virtualMachineExtensionImagesListTypes, virtualMachineExtensionImagesListVersions


text:
using DataValues
using Test

@testset "DataValueArray: Broadcast" begin

A1 = rand(Int, 10)
M1 = rand(Bool, 10)
n = rand(2:5)
dims = [ rand(2:5) for i in 1:n]
A2 = rand(10, dims...)
M2 = rand(Bool, 10, dims...)
C2 = Array{Float64}(undef, 10, dims...)
i = rand(2:5)
A3 = rand(10, [dims; i]...)
M3 = rand(Bool, 10, [dims; i]...)
C3 = Array{Float64}(undef, 10, [dims; i]...)

m1 = broadcast((x,y)->x, M1, M2)
m2 = broadcast((x,y)->x, M2, M3)
R2 = reshape(Bool[ (x->(m1[x] | M2[x]))(i) for i in 1:length(M2) ], size(M2))
R3 = reshape(Bool[ (x->(m2[x] | M3[x]))(i) for i in 1:length(M3) ], size(M3))
r2 = broadcast((x,y)->x, R2, R3)
Q3 = reshape(Bool[ (x->(r2[x] | R3[x]))(i) for i in 1:length(R3) ], size(M3))

U1 = DataValueArray(A1)
U2 = DataValueArray(A2)
U3 = DataValueArray(A3)
V1 = DataValueArray(A1, M1)
V2 = DataValueArray(A2, M2)
V3 = DataValueArray(A3, M3)
Z2 = DataValueArray{Float64}(10, dims...)
Z3 = DataValueArray{Float64}(10, [dims; i]...)

f() = 5
f(x) = 5 * x
f(x, y) = x * y
f(x, y, z) = x * y * z

for (dests, arrays, nullablearrays, mask) in
    ( ((C2, Z2), (A1, A2), (U1, U2), ()),
      ((C3, Z3), (A2, A3), (U2, U3), ()),
      ((C3, Z3), (A1, A2, A3), (U1, U2, U3), ()),

      ((C2, Z2), (A1, A2), (V1, V2), (R2,)),
      ((C3, Z3), (A2, A3), (V2, V3), (R3,)),
      ((C3, Z3), (A1, A2, A3), (V1, V2, V3), (Q3,)),
)

    # Base.broadcast!(f, B::DataValueArray, As::DataValueArray...)
    broadcast!(f, dests[1], arrays...)
    broadcast!(f, dests[2], nullablearrays...)
    @test isequal(dests[2], DataValueArray(dests[1], mask...))

    # Base.broadcast(f, As::DataValueArray...)
    D = broadcast(f, arrays...)
    X = broadcast(f, nullablearrays...)
    @test isequal(X, DataValueArray(D, mask...))
end

# Base.broadcast!(f, X::DataValueArray)
for (array, nullablearray, mask) in
    ( (A1, U1, ()), (A2, U2, ()), (A3, U3, ()),
      (A1, V1, (M1,)), (A2, V2, (M2,)), (A3, V3, (M3,)),
)
    broadcast!(f, array)
    broadcast!(f, nullablearray)
    @test isequal(nullablearray, DataValueArray(array, mask...))
end

# test broadcasted arithmetic operators
A = rand(10)
X1 = DataValueArray(A)
n = rand(2:5)
dims = rand(2:5, n)
B = rand(Float64, 10, dims...)
X2 = DataValueArray(B)
M = rand(Bool, 10, dims...)
Y = DataValueArray(B, M)

A = rand(Bool, 100)
B = rand(Bool, 100)
M1 = rand(Bool, 100)
M2 = rand(Bool, 100)
X = DataValueArray(A, M1)
Y = DataValueArray(B, M2)
# TODO This NullableArray test doesn't use 3VL
# @test isequal(broadcast(&, X, Y), DataValueArray(A .& B, M1 .| M2))
# @test isequal(broadcast(|, X, Y), DataValueArray(A .| B, M1 .| M2))

end # module


text:
@block AnshulSinghvi "UMAP" begin

    @cell "Basic UMAP visualization" [2d, scatter, umap, ml] begin
        using MLDatasets, UMAP
        mnist_x = MNIST.convert2features(MNIST.traintensor(Float64))
        mnist_y = MNIST.trainlabels(1:size(mnist_x, 2));
        res_jl = UMAP.umap(mnist_x; n_neighbors=10, min_dist=0.001, n_epochs=200)
        scatter(res_jl[1, :], res_jl[2, :]; color = mnist_y, colormap = :magma) # as in repo
    end

    @cell "3D UMAP visualization" [3d, scatter, umap, ml] begin
        using MLDatasets, UMAP
        mnist_x = MNIST.convert2features(MNIST.traintensor(Float64))
        mnist_y = MNIST.trainlabels(1:size(mnist_x, 2));
        res_jl = UMAP.umap(mnist_x, 3; n_neighbors=10, min_dist=0.001, n_epochs=200)
        scatter(eachrow(res)...; color = mnist_y, colormap = :magma) # 3D, Makie style!
    end

    # @cell "3D UMAP visualization with mouse highlighting" [3d, scatter, umap, ml] begin
    #     using MLDatasets, UMAP
    #     mnist_x = MNIST.convert2features(MNIST.traintensor(Float64))
    #     mnist_y = MNIST.trainlabels(1:size(mnist_x, 2));
    #     res_jl = UMAP.umap(mnist_x, 3; n_neighbors=10, min_dist=0.001, n_epochs=200)
    #     sc = scatter(eachrow(res)...; color = mnist_y, colormap = :magma) # 3D, Makie style!
    #     ju
    # end
end


text:
using JuMP, Cbc

"""
    kantorovich_distance(ν, μ)

Calculate the kantorovich distance between two probability distributions.

Thanks to this [hero] (https://stla.github.io/stlapblog/posts/KantorovichWithJulia.html)

```jldoctest
julia> mu = [1/7, 2/7, 4/7];
julia> nu = [1/4, 1/4, 1/2];
julia> kantorovich_distance(mu, nu)
0.10714285714285714
```
"""
function kantorovich_distance(ν, μ)
    n = length(ν)
    m = Model(solver=CbcSolver())
    @variable(m, p[1:n, 1:n] >= 0)
    @objective(m, Min, sum(p[i, j] for i in 1:n for j in 1:n if i != j))
    for k in 1:n
        @constraint(m, sum(p[k, :]) == μ[k])
        @constraint(m, sum(p[:, k]) == ν[k])
    end
    solve(m)
    return getobjectivevalue(m)
end


text:
using FEAST
using Test

@testset "FEAST.jl" begin
    # Write your tests here.
end


text:
export ScoreSpherical

"""
    ScoreSpherical <: ScoringRule

 Calculate Spherical Scores for a set of outcomes and predictions

## References
* https://en.wikipedia.org/wiki/Scoring_rule
*

## Parameters
* `normalise::false`: if true, normalise scores so that for binary choice
     + S( [0.5,0.5], 1) = 0.0
     + S( [0.0,1.0], 1) = 1.0

```
"""
struct ScoreSpherical <: ScoringRule
    normalise::Bool
end
ScoreSpherical(; normalise::Bool=false) = ScoreSpherical(normalise)

function score_direction(srule::ScoreSpherical)
    return 1
end

function scoring_function(srule::ScoreSpherical, predicted_probabilities::Array{Float64,1}, outcome::Int)
    if abs( sum(predicted_probabilities) - 1.0 ) > 1.0e-2
        error("probabilities must sum to 1")
    end
    result = predicted_probabilities[outcome] / sqrt( sum( predicted_probabilities.^2 ) )
    if srule.normalise
        x = 2/(2 - sqrt(2))
        y = 1 - x
        return result*x + y
    else
        return result
    end
end


text:
"""
Generate a binary tree of tiles of depth ```n```.
The tree is generated by recursvely splitting one tile into two subtiles, altrernating the orientation of the split.
The top-level tile will have ```size(tile) = (1,1)```

Additional arguments:
    - ```orientation::Orientation```: orientation of the top-level split. Defaults to ```H```.

"""
function generate_tree(n::Int; orientation = H)
    if n == 1
        return TileNode(1, 1)
    else
        tile1 = generate_tree(n-1; orientation = orientation == H ? V : H)

        if orientation == H
            h = height(tile1)
            b1 = base(tile1)
            b2 = rand(Float64)
            b = b2 / (b1 + b2)
            rescale!((b1/(b1 + b2),1),tile1)
        else
            b = base(tile1)
            h1 = height(tile1)
            h2 = rand(Float64)
            h = h2 / (h1 + h2)
            rescale!((1,h1/(h1 + h2)),tile1)
        end

        tile2 = TileNode(b,h)

        if rand(Float64) < 0.5
            return TileLeaf((tile1,tile2),orientation)
        else
            return TileLeaf((tile2,tile1),orientation)
        end
    end
end

"""
Draw ```b*h``` grid of tiles generated by ```generate_tree(n)```, and save it to ```file```.

See [```Luxor``` documentation](http://juliagraphics.github.io/Luxor.jl/stable/) for a list of supported file formats.

Optional arguments:
    - ```size```: size in pixel of a single square of the grid. Defaults to ```500```;
    - ```palette```: the palette used to color tiles. Defaults to ```palette_StampingItOut```.

"""
function paint_tree((b,h), n, file; size=500, palette=palette_StampingItOut)
    Drawing(b*size,h*size,file)
    setline(1)
    for i in 0:b-1, j in 0:h-1
        start = (size*i,size*j)
        orient = rand(Float64) < 0.5 ? H : V
        tile = rescale!((size, size),generate_tree(n; orientation = orient))
        draw(tile; origin=start, palette=palette)
    end
    finish()
end


text:
function assign_scenario(s::BasePriceScenario, model::JuMP.Model, config::OptionsV2016, params::ParametersV2016, vars::VariablesV2016)
    # We add the first solve since our run is infeasible without it.
    JuMP.solve(model);
    setvalue(params.ψ₂, 0.0);
    JuMP.solve(model);

    photel = getvalue(vars.CPRICE);

    setvalue(params.ψ₂, config.ψ₂);
    for i in 1:config.N
        if i <= config.tnopol
            setupperbound(vars.CPRICE[i], max(photel[i],params.cpricebase[i]));
        end
    end
end

function assign_scenario(s::OptimalPriceScenario, model::JuMP.Model, config::OptionsV2016, params::ParametersV2016, vars::VariablesV2016)
    setupperbound(vars.μ[1], config.μ₀);
end

function assign_scenario(s::Scenario, model::JuMP.Model, config::OptionsV2016, params::ParametersV2016, vars::VariablesV2016)
    error("$(s) is not a valid scenario for v2016R beta");
end


text:
"""
Based partly on code in LegacyStrings that used to be part of Julia
Licensed under MIT License, see LICENSE.md

(Written by Scott P. Jones in series of PRs contributed to the Julia project in 2015)
"""
module StrErrors

## Error messages for Unicode / UTF support

const SHORT =
  "invalid UTF-8 sequence starting at index <<1>> (0x<<2>>) missing one or more continuation bytes"
const CONT =
  "invalid UTF-8 sequence starting at index <<1>> (0x<<2>>) is not a continuation byte"
const LONG =
  "invalid UTF-8 sequence, overlong encoding starting at index <<1>> (0x<<2>>)"
const NOT_LEAD =
  "not a leading Unicode surrogate code unit at index <<1>> (0x<<2>>)"
const NOT_TRAIL =
  "not a trailing Unicode surrogate code unit at index <<1>> (0x<<2>>)"
const NOT_SURROGATE =
  "not a valid Unicode surrogate code unit at index <<1>> (0x<<2>>)"
const MISSING_SURROGATE =
  "missing trailing Unicode surrogate code unit after index <<1>> (0x<<2>>)"
const INVALID =
  "invalid Unicode character starting at index <<1>> (0x<<2>> > 0x10ffff)"
const SURROGATE =
  "surrogate encoding not allowed in UTF-8 or UTF-32, at index <<1>> (0x<<2>>)"
const ODD_BYTES_16 =
  "UTF16String can't have odd number of bytes <<1>>"
const ODD_BYTES_32 =
  "UTF32String must have multiple of 4 bytes <<1>>"
const INVALID_ASCII =
  "invalid ASCII character at index <<1>> (0x<<2>> > 0x7f)"
const INVALID_LATIN1 =
  "invalid Latin1 character at index <<1>> (0x<<2>> > 0xff)"
const INVALID_CHAR =
  "invalid Unicode character (0x<<2>> > 0x10ffff)"
const INVALID_8 =
  "invalid UTF-8 data"
const INVALID_16 =
  "invalid UTF-16 data"
const INVALID_UCS2 =
  "invalid UCS-2 character (surrogate present)"
const INVALID_INDEX =
  "invalid character index <<1>>"
const DECOMPOSE_COMPOSE =
  "only one of decompose or compose may be true"
const COMPAT_STRIPMARK =
  "compat or stripmark true requires compose or decompose true"
const NL_CONVERSION =
  "only one newline conversion may be specified"
const NORMALIZE =
  " is not one of :NFC, :NFD, :NFKC, :NFKD"

end # module StrErrors

@static if isdefined(Base, :UnicodeError)
const StringError = UnicodeError
else

struct StringError <: Exception
    errmsg::String           # Error message
    errpos::Int32            # Position of invalid character
    errchr::UInt32           # Invalid character
    StringError(msg, pos, chr) = new(msg, pos%Int32, chr%UInt32)
end

_repmsg(msg, pos, chr) =
    replace(replace(msg, "<<1>>" => string(pos)), "<<2>>" =>  outhex(chr))
Base.show(io::IO, exc::StringError) =
    print(io, "StringError: ", _repmsg(exc.errmsg, exc.errpos, exc.errchr))
end
(::Type{StringError})(msg, pos=0%Int32) = StringError(msg, pos%Int32, 0%UInt32)

@noinline boundserr(s, pos)      = throw(BoundsError(s, pos))
@noinline strerror(err)          = throw(StringError(err))
@noinline strerror(err, pos, ch) = throw(StringError(err, pos, ch))
@noinline strerror(err, v)       = strerror(string(":", v, err))
@noinline nulerr()               = strerror("cannot convert NULL to string")
@noinline neginderr(s, n)        = strerror("Index ($n) must be non negative")
@noinline codepoint_error(T, v)  = strerror(string("Invalid CodePoint: ", T, " 0x", outhex(v)))
@noinline argerror(startpos, endpos) =
    strerror(string("End position ", endpos, " is less than start position (", startpos, ")"))

@noinline ascii_err()    = throw(ArgumentError("Not a valid ASCII string"))
@noinline ncharerr(n)    = throw(ArgumentError(string("nchar (", n, ") must be not be negative")))
@noinline repeaterr(cnt) = throw(ArgumentError("repeat count $cnt must be >= 0"))

@static isdefined(Base, :string_index_err) && (const index_error = Base.string_index_err)


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct CurrentModeControl <: InnerControl
        kpc::Float64
        kic::Float64
        kffv::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
    end

Parameters of an inner loop PI current control using based on Purba, Dhople, Jafarpour, Bullo and Johnson.
"Reduced-order Structure-preserving Model for Parallel-connected Three-phase Grid-tied Inverters."
2017 IEEE 18th Workshop on Control and Modeling for Power Electronics (COMPEL): 1-7.

# Arguments
- `kpc::Float64`: Current controller proportional gain, validation range: `(0, nothing)`
- `kic::Float64`: Current controller integral gain, validation range: `(0, nothing)`
- `kffv::Float64`: Gain to enable feed-forward gain of voltage., validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the CurrentModeControl model are:
	γd_ic: d-axis integrator state of the PI current controller,
	γq_ic: q-axis integrator state of the PI current controller
- `n_states::Int`: CurrentControl has 2 states
"""
mutable struct CurrentModeControl <: InnerControl
    "Current controller proportional gain"
    kpc::Float64
    "Current controller integral gain"
    kic::Float64
    "Gain to enable feed-forward gain of voltage."
    kffv::Float64
    ext::Dict{String, Any}
    "The states of the CurrentModeControl model are:
	γd_ic: d-axis integrator state of the PI current controller,
	γq_ic: q-axis integrator state of the PI current controller"
    states::Vector{Symbol}
    "CurrentControl has 2 states"
    n_states::Int
end

function CurrentModeControl(kpc, kic, kffv, ext=Dict{String, Any}(), )
    CurrentModeControl(kpc, kic, kffv, ext, [:γd_ic, :γq_ic], 2, )
end

function CurrentModeControl(; kpc, kic, kffv, ext=Dict{String, Any}(), states=[:γd_ic, :γq_ic], n_states=2, )
    CurrentModeControl(kpc, kic, kffv, ext, states, n_states, )
end

# Constructor for demo purposes; non-functional.
function CurrentModeControl(::Nothing)
    CurrentModeControl(;
        kpc=0,
        kic=0,
        kffv=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`CurrentModeControl`](@ref) `kpc`."""
get_kpc(value::CurrentModeControl) = value.kpc
"""Get [`CurrentModeControl`](@ref) `kic`."""
get_kic(value::CurrentModeControl) = value.kic
"""Get [`CurrentModeControl`](@ref) `kffv`."""
get_kffv(value::CurrentModeControl) = value.kffv
"""Get [`CurrentModeControl`](@ref) `ext`."""
get_ext(value::CurrentModeControl) = value.ext
"""Get [`CurrentModeControl`](@ref) `states`."""
get_states(value::CurrentModeControl) = value.states
"""Get [`CurrentModeControl`](@ref) `n_states`."""
get_n_states(value::CurrentModeControl) = value.n_states

"""Set [`CurrentModeControl`](@ref) `kpc`."""
set_kpc!(value::CurrentModeControl, val) = value.kpc = val
"""Set [`CurrentModeControl`](@ref) `kic`."""
set_kic!(value::CurrentModeControl, val) = value.kic = val
"""Set [`CurrentModeControl`](@ref) `kffv`."""
set_kffv!(value::CurrentModeControl, val) = value.kffv = val
"""Set [`CurrentModeControl`](@ref) `ext`."""
set_ext!(value::CurrentModeControl, val) = value.ext = val


text:
function change_point_test_prep(N=10000, λ=N/10)
    L = @SMatrix [1. 0.;
                  0. 1.]
    Σ = @SMatrix [10^(-5) 0.;
         0. 10^(-5)]

    θ₀ = [10.0, -8.0, 15.0, 0.0, 3.0]

    param = :complexConjug
    # Target law
    P˟ = BSI.FitzhughDiffusion(param, θ₀...)

    # Auxiliary law
    t₀ = 1.0
    T = 2.0
    x0 = ℝ{2}(-0.5, 2.25)
    xT = ℝ{2}(1.0, 0.0)
    P̃ = BSI.FitzhughDiffusionAux(param, θ₀..., t₀, L*x0, T, L*xT)

    τ(t₀,T) = (x) ->  t₀ + (x-t₀) * (2-(x-t₀)/(T-t₀))
    dt = (T-t₀)/N
    tt = τ(t₀,T).(t₀:dt:T)

    P₁ = BSI.GuidPropBridge(eltype(x0), tt, P˟, P̃, L, L*x0, Σ;
                        change_pt=BSI.NoChangePt(), solver=BSI.Vern7())

    P₂ = BSI.GuidPropBridge(eltype(x0), tt, P˟, P̃, L, L*x0, Σ;
                        change_pt=BSI.SimpleChangePt(λ), solver=BSI.Vern7())
    P₁, P₂
end

@testset "change point between ODE solvers" begin
    N = 10000
    P₁, P₂ = change_point_test_prep(N)
    @testset "comparing H[$i]" for i in 1:div(N,20):N
        @test P₁.H[i] ≈ P₂.H[i]
    end
    @testset "comparing Hν[$i]" for i in 1:div(N,20):N
        @test P₁.Hν[i] ≈ P₂.Hν[i]
    end
    @testset "comparing c[$i]" for i in 1:div(N,20):N
        @test P₁.c[i] ≈ P₂.c[i]
    end
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

struct RawSMPS{T <: AbstractFloat}
    name::String
    tim::RawTime
    cor::RawCor{T}
    sto::RawStoch

    function RawSMPS(tim::RawTime, cor::RawCor{T}, sto::RawStoch) where T <: AbstractFloat
        tim.name == cor.name && tim.name == sto.name || error("Inconsistent names of SMPS files.")
        return new{T}(tim.name, tim, cor, sto)
    end
end

function Base.read(io::IO, ::Type{RawSMPS})
    # Get filepath
    m = match(r"<file (.*).smps>", io.name)
    m === nothing && error("SMPS specification is malformed. Correct filename usage: /path/to/smps/files.smps")
    path = only(m.captures)
    # Parse tim file
    timfile = path * ".tim"
    tim = parse_tim(timfile)
    # Parse cor file
    corfile = path * ".cor"
    cor = parse_cor(corfile)
    # Parse sto file
    stofile = path * ".sto"
    sto = parse_sto(tim, cor, stofile)
    return RawSMPS(tim, cor, sto)
end

function Base.read(io::IO, ::Type{RawSMPS{T}}) where T <: AbstractFloat
    # Get filepath
    m = match(r"<file (.*).smps>", io.name)
    m === nothing && error("SMPS specification is malformed. Correct filename usage: /path/to/smps/files.smps")
    path = only(m.captures)
    # Parse tim file
    timfile = path * ".tim"
    tim = parse_tim(T, timfile)
    # Parse cor file
    corfile = path * ".cor"
    cor = parse_cor(T, corfile)
    # Parse sto file
    stofile = path * ".sto"
    sto = parse_sto(T, tim, cor, io)
    return RawSMPS(tim, cor, sto)
end


text:
using BinDeps,Compat

@BinDeps.setup

libmetis = library_dependency("libmetis", aliases=["libmetis5"])

@windows_only begin
    using WinRPM
    provides(WinRPM.RPM, "metis", libmetis, os = :Windows)
end

@osx_only begin
    if Pkg.installed("Homebrew") === nothing
        error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")  end
    using Homebrew
    provides(Homebrew.HB, "metis", libmetis, os = :Darwin)
end

provides(AptGet, "libmetis5", libmetis)

provides(Yum, "metis-5.1.0", libmetis)

provides(Sources, URI("http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/metis-5.1.0.tar.gz"), libmetis)

metisname = "metis-5.1.0"

depsdir = BinDeps.depsdir(libmetis)
prefix = joinpath(depsdir,"usr")
srcdir = joinpath(depsdir,"src",metisname)

provides(SimpleBuild,
         (@build_steps begin
             GetSources(libmetis)
             (@build_steps begin
	         ChangeDirectory(srcdir)
	         (@build_steps begin
	             `make config shared=1 prefix=$prefix`
	             `make`
	             `make install`
                 end)
	     end)
         end), [libmetis], os=:Unix)

@BinDeps.install @compat Dict(:libmetis => :libmetis)


text:
cd("C:\\Users\\maxim\\.julia\\dev\\sem\\test\\comparisons\\factor_simulation")

include("factor_functions.jl")

config = DataFrame(CSV.File("config_factor.csv"))

data_vec = read_files("data", get_data_paths(config))
par_vec = read_files("parest", get_data_paths(config))
start_vec = read_files("start", get_data_paths(config))

##############################################
models = gen_models(config, data_vec, start_vec)

par_order = [model[2] for model in models]

models = [model[1] for model in models]

fits = get_fits(models)
##############################################

check_solution(fits, par_vec, par_order)

##############################################

benchmarks = benchmark_models(models)

##############################################

start_val = models[1].imply.start_val
grad = similar(start_val)
grad2 = similar(start_val)

@btime models[1](start_val, grad)

FiniteDiff.finite_difference_gradient!(grad2, models[2], start_val)

isapprox(grad, grad2)

diff_ana =
    SemAnalyticDiff(
        :LD_LBFGS,
        nothing,
        (grad_fiml,))

using Cthulhu

start_val = convert(Vector{Float64}, start_vec[1].est[par_order[1]])

@descend models[1](start_val)

grad = similar(start_val)

models[1](start_val, grad)

using FiniteDiff

grad2 = FiniteDiff.finite_difference_gradient(models[1], start_val)

grad ≈ grad2

@descend models[1](start_val, grad)

### NLopt

sol = sem.sem_fit_nlopt(models[1])

maximum(abs.(sol[2] .- par_vec[1].est[par_order[1]]))

using BenchmarkTools

@benchmark sem.sem_fit_nlopt(models[1])

sol2 = sem_fit_nlopt(models[2])

@time sol2 = sem_fit_nlopt(models[2])

maximum(abs.(sol2[2] .- par_vec[2].est[par_order[2]]))


text:
# Methods to implement DataStreams.jl sources and sinks (simple way to enable I/O).

# Table

# Sink implementation
function Table(???)
end


text:
using DifferentialEquations, LinearAlgebra
using Plots
gr()

const N0 =  2.5000013629442974e25
const M0 = 1.5264213305629408e16
# number of air molecules inside the fiber is 1.5264213305629408e16
#=
ozone from
http://www.columbia.edu/itc/chemistry/chem-c2407/hw/ozone_kinetics.pdf
equations 5 & 6
=#

function gaussianPulse(;τ= τ, FWHM=FWHM, t0=t0)
    res0 = sqrt.((exp.(-((τ.-t0).^2.0) ./ (FWHM^2.0))))
end

function createTimeGrid(time_window; npts = 2^14)
    dt = time_window/npts
    τ = LinRange(0, time_window, npts)
    τ = collect(τ)
end

t=createTimeGrid(200)
FWHM=10
gauss(τ) = gaussianPulse(τ= τ, FWHM=FWHM, t0=10) .+ gaussianPulse(τ= τ, FWHM=FWHM, t0=60) .+ gaussianPulse(τ= τ, FWHM=FWHM, t0=110) .+
            gaussianPulse(τ= τ, FWHM=FWHM, t0=160)
# plot(t, gauss.(t))

function fun(du,u,p,t)
    # k1,k2,k3,k4 = p
    k1 = 1.03E-14*gauss(t)
    k3 = 8.58E-2*gauss(t)

    # k1 = 1.03E-14
    k2 = 5.60E-34
    # k3 = 8.58E-2
    # k4 = 8.34E-15
    k4=8.34E-19#sice the reactions are slow compared to the time scale of the ultrafast pulses
    # A = N0*0.50*1e-6*1e-3
    M = M0
    # print(t)
    #=
    u[1] = O, u[2] = O3, u[3] = O2, du[4] = dM = 0
    =#
    du[1] = 2*k1*u[3]-k2*u[1]*u[3]*M+k3*u[2]-k4*u[1]*u[2]
    du[2] = k2*u[1]*u[3]*M - k3*u[2]-k4*u[1]*u[2]
    du[3] = -k1*u[3] - k2*u[3]*M*u[1] + k3*u[2] + 2*k4*u[1]*u[2]
end

# k1 = 1.03E-14
# k2 = 5.60E-34
# k3 = 8.58E-2
# k4 = 8.34E-15
# # A = N0*0.50*1e-6*1e-3
# M = N0*0.80*1e-6*1e-3
p = []

tspan = (0.0, 200)
u0 = [1e7,1e7, N0*0.20*1e-6*1e-3]

prob1 = ODEProblem(fun, u0, tspan, p)
sol = solve(prob1, Tsit5(),reltol = 1e-10, abstol=1e-16)

plot(sol,vars=(0,1))
plot!(sol,vars=(0,2))
# plot!(sol,vars=(0,3))

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev


text:
module MicrobiomeAnalysis

# using SummarizedExperiments: assay, rowdata, coldata
# using MultiAssayExperiments: MultiAssayExperiment
using Distances: pairwise, BrayCurtis, HellingerDist, Jaccard
using MultivariateStats: fit, MDS
using FdeSolver: FDEsolver
using Plots: plot
using Statistics: mean, std
using Microbiome: CommunityProfile, abundances, features, featurenames, samples, metadata, set!, Taxon, MicrobiomeSample, GeneFunction
using DataFrames: rename!, select!, sort, dropmissing
using Random, CSV

import Microbiome: shannon, ginisimpson, shannon!, ginisimpson!, braycurtis, jaccard, hellinger, pcoa, CommunityProfile
import Base: size, log10
import DataFrames: nrow, ncol, transform, transform!, dropmissing
import SummarizedExperiments: SummarizedExperiment

using Reexport
@reexport using SummarizedExperiments, MultiAssayExperiments
@reexport using DataStructures: OrderedDict
@reexport using DataFrames: DataFrame

include("alpha.jl")
export(shannon)
export(ginisimpson)
export(shannon!)
export(ginisimpson!)

include("beta.jl")
export(braycurtis)
export(jaccard)
export(hellinger)
export(pcoa)

include("wrappers.jl")
export(LVmodel)
export(abundance_plot)
export(size)
export(nrow)
export(ncol)

include("importers.jl")
export(import_se_from_csv)
export(import_mae_from_csv)

include("converters.jl")
export(SummarizedExperiment)
export(CommunityProfile)

include("artifacts.jl")
export(HintikkaXOData)
export(OKeefeDSData)

include("transform.jl")
export(transform)
export(transform!)
export(log10)
export(pa)
export(relabund)
export(ztransform)
export(clr)

include("sorters.jl")
export(select_top_taxa)
export(dropmissing)

end


text:
using Test, Random

function padnum(m,x)
    digs = floor(Int, log10(m))
    digs_x = floor(Int, log10(x))
    v = fill(' ', digs-digs_x)
    for d in digits(x)[end:-1:1] push!(v, '0' + d) end
    String(v)
end

Random.seed!(1)
@testset verbose = true "FFTA" begin
    @testset verbose = true "1D" begin
        @testset verbose = false "Complex" begin
            include("onedim/complex_forward.jl")
            include("onedim/complex_backward.jl")
            x = rand(ComplexF64, 100)
            y = fft(x)
            x2 = bfft(y)/length(x)
            @test x ≈ x2 atol=1e-12
        end
        @testset verbose = false "Real" begin
            include("onedim/real_forward.jl")
            include("onedim/real_backward.jl")
            x = rand(Float64, 100)
            y = fft(x)
            x2 = bfft(y)/length(x)
            @test x ≈ x2 atol=1e-12
        end
    end
    @testset verbose = false "2D" begin
        @testset verbose = true "Complex" begin
            include("twodim/complex_forward.jl")
            include("twodim/complex_backward.jl")
            x = rand(ComplexF64, 100, 100)
            y = fft(x)
            x2 = bfft(y)/length(x)
            @test x ≈ x2
        end
        @testset verbose = true "Real" begin
            include("twodim/real_forward.jl")
            include("twodim/real_backward.jl")
            x = rand(Float64, 100, 100)
            y = fft(x)
            x2 = bfft(y)/length(x)
            @test x ≈ x2
        end
    end
end


text:
export AbstractAdditiveRungeKutta
export LowStorageVariant, NaiveVariant
export AdditiveRungeKutta
export ARK2GiraldoKellyConstantinescu
export ARK548L2SA2KennedyCarpenter, ARK437L2SA1KennedyCarpenter

# Naive formulation that uses equation 3.8 from Giraldo, Kelly, and
# Constantinescu (2013) directly.  Seems to cut the number of solver iterations
# by half but requires Nstages - 1 additional storage.
struct NaiveVariant end
additional_storage(::NaiveVariant, Q, Nstages) =
    (Lstages = ntuple(i -> similar(Q), Nstages),)

# Formulation that does things exactly as in Giraldo, Kelly, and Constantinescu
# (2013).  Uses only one additional vector of storage regardless of the number
# of stages.
struct LowStorageVariant end
additional_storage(::LowStorageVariant, Q, Nstages) = (Qtt = similar(Q),)

abstract type AbstractAdditiveRungeKutta <: AbstractODESolver end

"""
    AdditiveRungeKutta(f, l, backward_euler_solver, RKAe, RKAi, RKB, RKC, Q;
                       split_explicit_implicit, variant, dt, t0 = 0)

This is a time stepping object for implicit-explicit time stepping of a
decomposed differential equation. When `split_explicit_implicit == false`
the equation is assumed to be decomposed as

```math
  \\dot{Q} = [l(Q, t)] + [f(Q, t) - l(Q, t)]
```

where `Q` is the state, `f` is the full tendency and
`l` is the chosen implicit operator. When `split_explicit_implicit == true`
the assumed decomposition is

```math
  \\dot{Q} = [l(Q, t)] + [f(Q, t)]
```

where `f` is now only the nonlinear tendency. For both decompositions the implicit
operator `l` is integrated implicitly whereas the remaining part is integrated
explicitly. Other arguments are the required time step size `dt` and the
optional initial time `t0`. The resulting backward Euler type systems are solved
using the provided `backward_euler_solver`. This time stepping object is
intended to be passed to the `solve!` command.

The constructor builds an additive Runge--Kutta scheme based on the provided
`RKAe`, `RKAi`, `RKB` and `RKC` coefficient arrays.  Additionally `variant`
specifies which of the analytically equivalent but numerically different
formulations of the scheme is used.

The available concrete implementations are:

  - [`ARK2GiraldoKellyConstantinescu`](@ref)
  - [`ARK548L2SA2KennedyCarpenter`](@ref)
  - [`ARK437L2SA1KennedyCarpenter`](@ref)
"""
mutable struct AdditiveRungeKutta{T, RT, AT, BE, V, VS, Nstages, Nstages_sq} <:
               AbstractAdditiveRungeKutta
    "time step"
    dt::RT
    "time"
    t::RT
    "rhs function"
    rhs!
    "rhs linear operator"
    rhs_implicit!
    "backwark Euler solver"
    besolver!::BE
    "Storage for solution during the AdditiveRungeKutta update"
    Qstages::NTuple{Nstages, AT}
    "Storage for RHS during the AdditiveRungeKutta update"
    Rstages::NTuple{Nstages, AT}
    "Storage for the linear solver rhs vector"
    Qhat::AT
    "RK coefficient matrix A for the explicit scheme"
    RKA_explicit::SArray{NTuple{2, Nstages}, RT, 2, Nstages_sq}
    "RK coefficient matrix A for the implicit scheme"
    RKA_implicit::SArray{NTuple{2, Nstages}, RT, 2, Nstages_sq}
    "RK coefficient vector B (rhs add in scaling)"
    RKB::SArray{Tuple{Nstages}, RT, 1, Nstages}
    "RK coefficient vector C (time scaling)"
    RKC::SArray{Tuple{Nstages}, RT, 1, Nstages}
    split_explicit_implicit::Bool
    "Variant of the ARK scheme"
    variant::V
    "Storage dependent on the variant of the ARK scheme"
    variant_storage::VS

    function AdditiveRungeKutta(
        rhs!,
        rhs_implicit!,
        backward_euler_solver,
        RKA_explicit,
        RKA_implicit,
        RKB,
        RKC,
        split_explicit_implicit,
        variant,
        Q::AT;
        dt = nothing,
        t0 = 0,
    ) where {AT <: AbstractArray}

        @assert dt != nothing

        T = eltype(Q)
        RT = real(T)

        Nstages = length(RKB)

        Qstages = (Q, ntuple(i -> similar(Q), Nstages - 1)...)
        Rstages = ntuple(i -> similar(Q), Nstages)
        Qhat = similar(Q)

        V = typeof(variant)
        variant_storage = additional_storage(variant, Q, Nstages)
        VS = typeof(variant_storage)

        # The code throughout assumes SDIRK implicit tableau so we assert that
        # here.
        for is in 2:Nstages
            @assert RKA_implicit[is, is] ≈ RKA_implicit[2, 2]
        end

        α = dt * RKA_implicit[2, 2]
        besolver! = setup_backward_Euler_solver(
            backward_euler_solver,
            Q,
            α,
            rhs_implicit!,
        )
        @assert besolver! isa AbstractBackwardEulerSolver
        @assert besolver! isa LinBESolver || variant isa NaiveVariant
        BE = typeof(besolver!)

        new{T, RT, AT, BE, V, VS, Nstages, Nstages^2}(
            RT(dt),
            RT(t0),
            rhs!,
            rhs_implicit!,
            besolver!,
            Qstages,
            Rstages,
            Qhat,
            RKA_explicit,
            RKA_implicit,
            RKB,
            RKC,
            split_explicit_implicit,
            variant,
            variant_storage,
        )
    end
end

# this will only work for iterative solves
# direct solvers use prefactorization
function updatedt!(ark::AdditiveRungeKutta, dt)
    @assert Δt_is_adjustable(ark.besolver!)
    ark.dt = dt
    α = dt * ark.RKA_implicit[2, 2]
    update_backward_Euler_solver!(ark.besolver!, ark.Qstages[1], α)
end

function dostep!(
    Q,
    ark::AdditiveRungeKutta,
    p,
    time,
    slow_δ = nothing,
    slow_rv_dQ = nothing,
    slow_scaling = nothing,
)
    dostep!(Q, ark, ark.variant, p, time, slow_δ, slow_rv_dQ, slow_scaling)
end

function dostep!(
    Q,
    ark::AdditiveRungeKutta,
    variant::NaiveVariant,
    p,
    time::Real,
    slow_δ = nothing,
    slow_rv_dQ = nothing,
    slow_scaling = nothing,
)
    dt = ark.dt

    besolver! = ark.besolver!
    RKA_explicit, RKA_implicit = ark.RKA_explicit, ark.RKA_implicit
    RKB, RKC = ark.RKB, ark.RKC
    rhs!, rhs_implicit! = ark.rhs!, ark.rhs_implicit!
    Qstages, Rstages = ark.Qstages, ark.Rstages
    Qhat = ark.Qhat
    split_explicit_implicit = ark.split_explicit_implicit
    Lstages = ark.variant_storage.Lstages

    rv_Q = realview(Q)
    rv_Qstages = realview.(Qstages)
    rv_Lstages = realview.(Lstages)
    rv_Rstages = realview.(Rstages)
    rv_Qhat = realview(Qhat)

    Nstages = length(RKB)

    groupsize = 256

    # calculate the rhs at first stage to initialize the stage loop
    rhs!(Rstages[1], Qstages[1], p, time + RKC[1] * dt, increment = false)

    rhs_implicit!(
        Lstages[1],
        Qstages[1],
        p,
        time + RKC[1] * dt,
        increment = false,
    )

    # note that it is important that this loop does not modify Q!
    for istage in 2:Nstages
        stagetime = time + RKC[istage] * dt

        # this kernel also initializes Qstages[istage] with an initial guess
        # for the linear solver
        event = Event(device(Q))
        event = stage_update!(device(Q), groupsize)(
            variant,
            rv_Q,
            rv_Qstages,
            rv_Lstages,
            rv_Rstages,
            rv_Qhat,
            RKA_explicit,
            RKA_implicit,
            dt,
            Val(istage),
            Val(split_explicit_implicit),
            slow_δ,
            slow_rv_dQ;
            ndrange = length(rv_Q),
            dependencies = (event,),
        )
        wait(device(Q), event)

        # solves
        # Qs = Qhat + dt * RKA_implicit[istage, istage] * rhs_implicit!(Qs)
        α = dt * RKA_implicit[istage, istage]
        besolver!(Qstages[istage], Qhat, α, p, stagetime)

        rhs!(Rstages[istage], Qstages[istage], p, stagetime, increment = false)
        rhs_implicit!(
            Lstages[istage],
            Qstages[istage],
            p,
            stagetime,
            increment = false,
        )
    end

    # compose the final solution
    event = Event(device(Q))
    event = solution_update!(device(Q), groupsize)(
        variant,
        rv_Q,
        rv_Lstages,
        rv_Rstages,
        RKB,
        dt,
        Val(Nstages),
        Val(split_explicit_implicit),
        slow_δ,
        slow_rv_dQ,
        slow_scaling;
        ndrange = length(rv_Q),
        dependencies = (event,),
    )
    wait(device(Q), event)
end

function dostep!(
    Q,
    ark::AdditiveRungeKutta,
    variant::LowStorageVariant,
    p,
    time::Real,
    slow_δ = nothing,
    slow_rv_dQ = nothing,
    slow_scaling = nothing,
)
    dt = ark.dt

    besolver! = ark.besolver!
    RKA_explicit, RKA_implicit = ark.RKA_explicit, ark.RKA_implicit
    RKB, RKC = ark.RKB, ark.RKC
    rhs!, rhs_implicit! = ark.rhs!, ark.rhs_implicit!
    Qstages, Rstages = ark.Qstages, ark.Rstages
    Qhat = ark.Qhat
    split_explicit_implicit = ark.split_explicit_implicit
    Qtt = ark.variant_storage.Qtt

    rv_Q = realview(Q)
    rv_Qstages = realview.(Qstages)
    rv_Rstages = realview.(Rstages)
    rv_Qhat = realview(Qhat)
    rv_Qtt = realview(Qtt)

    Nstages = length(RKB)

    groupsize = 256

    # calculate the rhs at first stage to initialize the stage loop
    rhs!(Rstages[1], Qstages[1], p, time + RKC[1] * dt, increment = false)

    # note that it is important that this loop does not modify Q!
    for istage in 2:Nstages
        stagetime = time + RKC[istage] * dt

        # this kernel also initializes Qtt for the linear solver
        event = Event(device(Q))
        event = stage_update!(device(Q), groupsize)(
            variant,
            rv_Q,
            rv_Qstages,
            rv_Rstages,
            rv_Qhat,
            rv_Qtt,
            RKA_explicit,
            RKA_implicit,
            dt,
            Val(istage),
            Val(split_explicit_implicit),
            slow_δ,
            slow_rv_dQ;
            ndrange = length(rv_Q),
            dependencies = (event,),
        )
        wait(device(Q), event)

        # solves
        # Q_tt = Qhat + dt * RKA_implicit[istage, istage] * rhs_implicit!(Q_tt)
        α = dt * RKA_implicit[istage, istage]
        besolver!(Qtt, Qhat, α, p, stagetime)

        # update Qstages
        Qstages[istage] .+= Qtt

        rhs!(Rstages[istage], Qstages[istage], p, stagetime, increment = false)
    end

    if split_explicit_implicit
        for istage in 1:Nstages
            stagetime = time + RKC[istage] * dt
            rhs_implicit!(
                Rstages[istage],
                Qstages[istage],
                p,
                stagetime,
                increment = true,
            )
        end
    end

    # compose the final solution
    event = Event(device(Q))
    event = solution_update!(device(Q), groupsize)(
        variant,
        rv_Q,
        rv_Rstages,
        RKB,
        dt,
        Val(Nstages),
        slow_δ,
        slow_rv_dQ,
        slow_scaling;
        ndrange = length(rv_Q),
        dependencies = (event,),
    )
    wait(device(Q), event)
end

@kernel function stage_update!(
    ::NaiveVariant,
    Q,
    Qstages,
    Lstages,
    Rstages,
    Qhat,
    RKA_explicit,
    RKA_implicit,
    dt,
    ::Val{is},
    ::Val{split_explicit_implicit},
    slow_δ,
    slow_dQ,
) where {is, split_explicit_implicit}
    i = @index(Global, Linear)
    @inbounds begin
        Qhat_i = Q[i]
        Qstages_is_i = Q[i]

        if slow_δ !== nothing
            Rstages[is - 1][i] += slow_δ * slow_dQ[i]
        end

        @unroll for js in 1:(is - 1)
            R_explicit = dt * RKA_explicit[is, js] * Rstages[js][i]
            L_explicit = dt * RKA_explicit[is, js] * Lstages[js][i]
            L_implicit = dt * RKA_implicit[is, js] * Lstages[js][i]
            Qhat_i += (R_explicit + L_implicit)
            Qstages_is_i += R_explicit
            if split_explicit_implicit
                Qstages_is_i += L_explicit
            else
                Qhat_i -= L_explicit
            end
        end
        Qstages[is][i] = Qstages_is_i
        Qhat[i] = Qhat_i
    end
end

@kernel function stage_update!(
    ::LowStorageVariant,
    Q,
    Qstages,
    Rstages,
    Qhat,
    Qtt,
    RKA_explicit,
    RKA_implicit,
    dt,
    ::Val{is},
    ::Val{split_explicit_implicit},
    slow_δ,
    slow_dQ,
) where {is, split_explicit_implicit}
    i = @index(Global, Linear)
    @inbounds begin
        Qhat_i = Q[i]
        Qstages_is_i = -zero(eltype(Q))

        if slow_δ !== nothing
            Rstages[is - 1][i] += slow_δ * slow_dQ[i]
        end

        @unroll for js in 1:(is - 1)
            if split_explicit_implicit
                rkcoeff = RKA_implicit[is, js] / RKA_implicit[is, is]
            else
                rkcoeff =
                    (RKA_implicit[is, js] - RKA_explicit[is, js]) /
                    RKA_implicit[is, is]
            end
            commonterm = rkcoeff * Qstages[js][i]
            Qhat_i += commonterm + dt * RKA_explicit[is, js] * Rstages[js][i]
            Qstages_is_i -= commonterm
        end
        Qstages[is][i] = Qstages_is_i
        Qhat[i] = Qhat_i
        Qtt[i] = Qhat_i
    end
end

@kernel function solution_update!(
    ::NaiveVariant,
    Q,
    Lstages,
    Rstages,
    RKB,
    dt,
    ::Val{Nstages},
    ::Val{split_explicit_implicit},
    slow_δ,
    slow_dQ,
    slow_scaling,
) where {Nstages, split_explicit_implicit}
    i = @index(Global, Linear)
    @inbounds begin
        if slow_δ !== nothing
            Rstages[Nstages][i] += slow_δ * slow_dQ[i]
        end
        if slow_scaling !== nothing
            slow_dQ[i] *= slow_scaling
        end

        @unroll for is in 1:Nstages
            Q[i] += RKB[is] * dt * Rstages[is][i]
            if split_explicit_implicit
                Q[i] += RKB[is] * dt * Lstages[is][i]
            end
        end
    end
end

@kernel function solution_update!(
    ::LowStorageVariant,
    Q,
    Rstages,
    RKB,
    dt,
    ::Val{Nstages},
    slow_δ,
    slow_dQ,
    slow_scaling,
) where {Nstages}
    i = @index(Global, Linear)
    @inbounds begin
        if slow_δ !== nothing
            Rstages[Nstages][i] += slow_δ * slow_dQ[i]
        end
        if slow_scaling !== nothing
            slow_dQ[i] *= slow_scaling
        end

        @unroll for is in 1:Nstages
            Q[i] += RKB[is] * dt * Rstages[is][i]
        end
    end
end

"""
    ARK2GiraldoKellyConstantinescu(f, l, backward_euler_solver, Q; dt, t0,
                                   split_explicit_implicit, variant, paperversion)

This function returns an [`AdditiveRungeKutta`](@ref) time stepping object,
see the documentation of [`AdditiveRungeKutta`](@ref) for arguments definitions.
This time stepping object is intended to be passed to the `solve!` command.

`paperversion=true` uses the coefficients from the paper, `paperversion=false`
uses coefficients that make the scheme (much) more stable but less accurate

This uses the second-order-accurate 3-stage additive Runge--Kutta scheme of
Giraldo, Kelly and Constantinescu (2013).

### References
    @article{giraldo2013implicit,
      title={Implicit-explicit formulations of a three-dimensional
             nonhydrostatic unified model of the atmosphere ({NUMA})},
      author={Giraldo, Francis X and Kelly, James F and Constantinescu, Emil M},
      journal={SIAM Journal on Scientific Computing},
      volume={35},
      number={5},
      pages={B1162--B1194},
      year={2013},
      publisher={SIAM}
    }
"""
function ARK2GiraldoKellyConstantinescu(
    F,
    L,
    backward_euler_solver,
    Q::AT;
    dt = nothing,
    t0 = 0,
    split_explicit_implicit = false,
    variant = LowStorageVariant(),
    paperversion = false,
) where {AT <: AbstractArray}

    @assert dt != nothing

    T = eltype(Q)
    RT = real(T)

    a32 = RT(paperversion ? (3 + 2 * sqrt(2)) / 6 : 1 // 2)
    RKA_explicit = [
        RT(0) RT(0) RT(0)
        RT(2 - sqrt(2)) RT(0) RT(0)
        RT(1 - a32) RT(a32) RT(0)
    ]

    RKA_implicit = [
        RT(0) RT(0) RT(0)
        RT(1 - 1 / sqrt(2)) RT(1 - 1 / sqrt(2)) RT(0)
        RT(1 / (2 * sqrt(2))) RT(1 / (2 * sqrt(2))) RT(1 - 1 / sqrt(2))
    ]

    RKB = [RT(1 / (2 * sqrt(2))), RT(1 / (2 * sqrt(2))), RT(1 - 1 / sqrt(2))]
    RKC = [RT(0), RT(2 - sqrt(2)), RT(1)]

    Nstages = length(RKB)

    AdditiveRungeKutta(
        F,
        L,
        backward_euler_solver,
        RKA_explicit,
        RKA_implicit,
        RKB,
        RKC,
        split_explicit_implicit,
        variant,
        Q;
        dt = dt,
        t0 = t0,
    )
end

"""
    ARK548L2SA2KennedyCarpenter(f, l, backward_euler_solver, Q; dt, t0,
                                split_explicit_implicit, variant)

This function returns an [`AdditiveRungeKutta`](@ref) time stepping object,
see the documentation of [`AdditiveRungeKutta`](@ref) for arguments definitions.
This time stepping object is intended to be passed to the `solve!` command.

This uses the fifth-order-accurate 8-stage additive Runge--Kutta scheme of
Kennedy and Carpenter (2013).

### References

    @article{kennedy2019higher,
      title={Higher-order additive Runge--Kutta schemes for ordinary
             differential equations},
      author={Kennedy, Christopher A and Carpenter, Mark H},
      journal={Applied Numerical Mathematics},
      volume={136},
      pages={183--205},
      year={2019},
      publisher={Elsevier}
    }
"""
function ARK548L2SA2KennedyCarpenter(
    F,
    L,
    backward_euler_solver,
    Q::AT;
    dt = nothing,
    t0 = 0,
    split_explicit_implicit = false,
    variant = LowStorageVariant(),
) where {AT <: AbstractArray}

    @assert dt != nothing

    T = eltype(Q)
    RT = real(T)

    Nstages = 8
    gamma = RT(2 // 9)

    # declared as Arrays for mutability, later these will be converted to static
    # arrays
    RKA_explicit = zeros(RT, Nstages, Nstages)
    RKA_implicit = zeros(RT, Nstages, Nstages)
    RKB = zeros(RT, Nstages)
    RKC = zeros(RT, Nstages)

    # the main diagonal
    for is in 2:Nstages
        RKA_implicit[is, is] = gamma
    end

    RKA_implicit[3, 2] = RT(2366667076620 // 8822750406821)
    RKA_implicit[4, 2] = RT(-257962897183 // 4451812247028)
    RKA_implicit[4, 3] = RT(128530224461 // 14379561246022)
    RKA_implicit[5, 2] = RT(-486229321650 // 11227943450093)
    RKA_implicit[5, 3] = RT(-225633144460 // 6633558740617)
    RKA_implicit[5, 4] = RT(1741320951451 // 6824444397158)
    RKA_implicit[6, 2] = RT(621307788657 // 4714163060173)
    RKA_implicit[6, 3] = RT(-125196015625 // 3866852212004)
    RKA_implicit[6, 4] = RT(940440206406 // 7593089888465)
    RKA_implicit[6, 5] = RT(961109811699 // 6734810228204)
    RKA_implicit[7, 2] = RT(2036305566805 // 6583108094622)
    RKA_implicit[7, 3] = RT(-3039402635899 // 4450598839912)
    RKA_implicit[7, 4] = RT(-1829510709469 // 31102090912115)
    RKA_implicit[7, 5] = RT(-286320471013 // 6931253422520)
    RKA_implicit[7, 6] = RT(8651533662697 // 9642993110008)

    RKA_explicit[3, 1] = RT(1 // 9)
    RKA_explicit[3, 2] = RT(1183333538310 // 1827251437969)
    RKA_explicit[4, 1] = RT(895379019517 // 9750411845327)
    RKA_explicit[4, 2] = RT(477606656805 // 13473228687314)
    RKA_explicit[4, 3] = RT(-112564739183 // 9373365219272)
    RKA_explicit[5, 1] = RT(-4458043123994 // 13015289567637)
    RKA_explicit[5, 2] = RT(-2500665203865 // 9342069639922)
    RKA_explicit[5, 3] = RT(983347055801 // 8893519644487)
    RKA_explicit[5, 4] = RT(2185051477207 // 2551468980502)
    RKA_explicit[6, 1] = RT(-167316361917 // 17121522574472)
    RKA_explicit[6, 2] = RT(1605541814917 // 7619724128744)
    RKA_explicit[6, 3] = RT(991021770328 // 13052792161721)
    RKA_explicit[6, 4] = RT(2342280609577 // 11279663441611)
    RKA_explicit[6, 5] = RT(3012424348531 // 12792462456678)
    RKA_explicit[7, 1] = RT(6680998715867 // 14310383562358)
    RKA_explicit[7, 2] = RT(5029118570809 // 3897454228471)
    RKA_explicit[7, 3] = RT(2415062538259 // 6382199904604)
    RKA_explicit[7, 4] = RT(-3924368632305 // 6964820224454)
    RKA_explicit[7, 5] = RT(-4331110370267 // 15021686902756)
    RKA_explicit[7, 6] = RT(-3944303808049 // 11994238218192)
    RKA_explicit[8, 1] = RT(2193717860234 // 3570523412979)
    RKA_explicit[8, 2] = RKA_explicit[8, 1]
    RKA_explicit[8, 3] = RT(5952760925747 // 18750164281544)
    RKA_explicit[8, 4] = RT(-4412967128996 // 6196664114337)
    RKA_explicit[8, 5] = RT(4151782504231 // 36106512998704)
    RKA_explicit[8, 6] = RT(572599549169 // 6265429158920)
    RKA_explicit[8, 7] = RT(-457874356192 // 11306498036315)

    RKB[2] = 0
    RKB[3] = RT(3517720773327 // 20256071687669)
    RKB[4] = RT(4569610470461 // 17934693873752)
    RKB[5] = RT(2819471173109 // 11655438449929)
    RKB[6] = RT(3296210113763 // 10722700128969)
    RKB[7] = RT(-1142099968913 // 5710983926999)

    RKC[2] = RT(4 // 9)
    RKC[3] = RT(6456083330201 // 8509243623797)
    RKC[4] = RT(1632083962415 // 14158861528103)
    RKC[5] = RT(6365430648612 // 17842476412687)
    RKC[6] = RT(18 // 25)
    RKC[7] = RT(191 // 200)

    for is in 2:Nstages
        RKA_implicit[is, 1] = RKA_implicit[is, 2]
    end

    for is in 1:(Nstages - 1)
        RKA_implicit[Nstages, is] = RKB[is]
    end

    RKB[1] = RKB[2]
    RKB[8] = gamma

    RKA_explicit[2, 1] = RKC[2]
    RKA_explicit[Nstages, 1] = RKA_explicit[Nstages, 2]

    RKC[1] = 0
    RKC[Nstages] = 1

    # conversion to static arrays
    RKA_explicit = SMatrix{Nstages, Nstages}(RKA_explicit)
    RKA_implicit = SMatrix{Nstages, Nstages}(RKA_implicit)
    RKB = SVector{Nstages}(RKB)
    RKC = SVector{Nstages}(RKC)

    ark = AdditiveRungeKutta(
        F,
        L,
        backward_euler_solver,
        RKA_explicit,
        RKA_implicit,
        RKB,
        RKC,
        split_explicit_implicit,
        variant,
        Q;
        dt = dt,
        t0 = t0,
    )
end

"""
    ARK437L2SA1KennedyCarpenter(f, l, backward_euler_solver, Q; dt, t0,
                                split_explicit_implicit, variant)

This function returns an [`AdditiveRungeKutta`](@ref) time stepping object,
see the documentation of [`AdditiveRungeKutta`](@ref) for arguments definitions.
This time stepping object is intended to be passed to the `solve!` command.

This uses the fourth-order-accurate 7-stage additive Runge--Kutta scheme of
Kennedy and Carpenter (2013).

### References
    @article{kennedy2019higher,
      title={Higher-order additive Runge--Kutta schemes for ordinary
             differential equations},
      author={Kennedy, Christopher A and Carpenter, Mark H},
      journal={Applied Numerical Mathematics},
      volume={136},
      pages={183--205},
      year={2019},
      publisher={Elsevier}
    }
"""
function ARK437L2SA1KennedyCarpenter(
    F,
    L,
    backward_euler_solver,
    Q::AT;
    dt = nothing,
    t0 = 0,
    split_explicit_implicit = false,
    variant = LowStorageVariant(),
) where {AT <: AbstractArray}

    @assert dt != nothing

    T = eltype(Q)
    RT = real(T)

    Nstages = 7
    gamma = RT(1235 // 10000)

    # declared as Arrays for mutability, later these will be converted to static
    # arrays
    RKA_explicit = zeros(RT, Nstages, Nstages)
    RKA_implicit = zeros(RT, Nstages, Nstages)
    RKB = zeros(RT, Nstages)
    RKC = zeros(RT, Nstages)

    # the main diagonal
    for is in 2:Nstages
        RKA_implicit[is, is] = gamma
    end

    RKA_implicit[3, 2] = RT(624185399699 // 4186980696204)
    RKA_implicit[4, 2] = RT(1258591069120 // 10082082980243)
    RKA_implicit[4, 3] = RT(-322722984531 // 8455138723562)
    RKA_implicit[5, 2] = RT(-436103496990 // 5971407786587)
    RKA_implicit[5, 3] = RT(-2689175662187 // 11046760208243)
    RKA_implicit[5, 4] = RT(4431412449334 // 12995360898505)
    RKA_implicit[6, 2] = RT(-2207373168298 // 14430576638973)
    RKA_implicit[6, 3] = RT(242511121179 // 3358618340039)
    RKA_implicit[6, 4] = RT(3145666661981 // 7780404714551)
    RKA_implicit[6, 5] = RT(5882073923981 // 14490790706663)
    RKA_implicit[7, 2] = 0
    RKA_implicit[7, 3] = RT(9164257142617 // 17756377923965)
    RKA_implicit[7, 4] = RT(-10812980402763 // 74029279521829)
    RKA_implicit[7, 5] = RT(1335994250573 // 5691609445217)
    RKA_implicit[7, 6] = RT(2273837961795 // 8368240463276)

    RKA_explicit[3, 1] = RT(247 // 4000)
    RKA_explicit[3, 2] = RT(2694949928731 // 7487940209513)
    RKA_explicit[4, 1] = RT(464650059369 // 8764239774964)
    RKA_explicit[4, 2] = RT(878889893998 // 2444806327765)
    RKA_explicit[4, 3] = RT(-952945855348 // 12294611323341)
    RKA_explicit[5, 1] = RT(476636172619 // 8159180917465)
    RKA_explicit[5, 2] = RT(-1271469283451 // 7793814740893)
    RKA_explicit[5, 3] = RT(-859560642026 // 4356155882851)
    RKA_explicit[5, 4] = RT(1723805262919 // 4571918432560)
    RKA_explicit[6, 1] = RT(6338158500785 // 11769362343261)
    RKA_explicit[6, 2] = RT(-4970555480458 // 10924838743837)
    RKA_explicit[6, 3] = RT(3326578051521 // 2647936831840)
    RKA_explicit[6, 4] = RT(-880713585975 // 1841400956686)
    RKA_explicit[6, 5] = RT(-1428733748635 // 8843423958496)
    RKA_explicit[7, 2] = RT(760814592956 // 3276306540349)
    RKA_explicit[7, 3] = RT(-47223648122716 // 6934462133451)
    RKA_explicit[7, 4] = RT(71187472546993 // 9669769126921)
    RKA_explicit[7, 5] = RT(-13330509492149 // 9695768672337)
    RKA_explicit[7, 6] = RT(11565764226357 // 8513123442827)

    RKB[2] = 0
    RKB[3] = RT(9164257142617 // 17756377923965)
    RKB[4] = RT(-10812980402763 // 74029279521829)
    RKB[5] = RT(1335994250573 // 5691609445217)
    RKB[6] = RT(2273837961795 // 8368240463276)
    RKB[7] = RT(247 // 2000)

    RKC[2] = RT(247 // 1000)
    RKC[3] = RT(4276536705230 // 10142255878289)
    RKC[4] = RT(67 // 200)
    RKC[5] = RT(3 // 40)
    RKC[6] = RT(7 // 10)

    for is in 2:Nstages
        RKA_implicit[is, 1] = RKA_implicit[is, 2]
    end

    for is in 1:(Nstages - 1)
        RKA_implicit[Nstages, is] = RKB[is]
    end

    RKB[1] = RKB[2]

    RKA_explicit[2, 1] = RKC[2]
    RKA_explicit[Nstages, 1] = RKA_explicit[Nstages, 2]

    RKC[1] = 0
    RKC[Nstages] = 1

    # conversion to static arrays
    RKA_explicit = SMatrix{Nstages, Nstages}(RKA_explicit)
    RKA_implicit = SMatrix{Nstages, Nstages}(RKA_implicit)
    RKB = SVector{Nstages}(RKB)
    RKC = SVector{Nstages}(RKC)

    ark = AdditiveRungeKutta(
        F,
        L,
        backward_euler_solver,
        RKA_explicit,
        RKA_implicit,
        RKB,
        RKC,
        split_explicit_implicit,
        variant,
        Q;
        dt = dt,
        t0 = t0,
    )
end


text:
infeasible = @stochastic_model begin
    @stage 1 begin
        @decision(model, x₁ >= 0)
        @decision(model, x₂ >= 0)
        @objective(model, Min, 3*x₁ + 2*x₂)
    end
    @stage 2 begin
        @uncertain ξ₁ ξ₂
        @recourse(model, 0.8*ξ₁ <= y₁ <= ξ₁)
        @recourse(model, 0.8*ξ₂ <= y₂ <= ξ₂)
        @objective(model, Min, -15*y₁ - 12*y₂)
        @constraint(model, 3*y₁ + 2*y₂ <= x₁)
        @constraint(model, 2*y₁ + 5*y₂ <= x₂)
    end
end

ξ₁ = @scenario ξ₁ = 6. ξ₂ = 8. probability = 0.5
ξ₂ = @scenario ξ₁ = 4. ξ₂ = 4. probability = 0.5

infeasible_res = SPResult([27.2,41.6], Dict(1 => [4.8, 6.4], 2 => [4., 4.]), 36.4, 9.2, 27.2, Inf, 9.2, Inf)
push!(problems, (infeasible, [ξ₁,ξ₂], infeasible_res, "Infeasible"))