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

@test integrate((d+e*x^2)^4*(a+b*arcsinh(c*x)), x) == :(x*d^4*(a+b*arcsinh(c*x))+1//9*e^4*x^9*(a+b*arcsinh(c*x))+-1//81*b*c^-9*e^4*(1+c^2*x^2)^9//2+-1//315*b*c^-9*(1+c^2*x^2)^(1/2)*(35*e^4+315*c^8*d^4+-420*e*c^6*d^3+-180*d*c^2*e^3+378*c^4*d^2*e^2)+4//3*e*d^3*x^3*(a+b*arcsinh(c*x))+4//7*d*e^3*x^7*(a+b*arcsinh(c*x))+6//5*d^2*e^2*x^5*(a+b*arcsinh(c*x))+-4//441*b*c^-9*e^3*(1+c^2*x^2)^7//2*(-7*e+9*d*c^2)+-4//945*b*e*c^-9*(1+c^2*x^2)^3//2*(-35*e^3+105*c^6*d^3+-189*e*c^4*d^2+135*d*c^2*e^2)+-2//525*b*c^-9*e^2*(1+c^2*x^2)^5//2*(35*e^2+63*c^4*d^2+-90*d*e*c^2))
@test integrate((d+e*x^2)^3*(a+b*arcsinh(c*x)), x) == :(x*d^3*(a+b*arcsinh(c*x))+1//7*e^3*x^7*(a+b*arcsinh(c*x))+e*d^2*x^3*(a+b*arcsinh(c*x))+-1//35*b*c^-7*(1+c^2*x^2)^(1/2)*(-5*e^3+35*c^6*d^3+-35*e*c^4*d^2+21*d*c^2*e^2)+-1//49*b*c^-7*e^3*(1+c^2*x^2)^7//2+3//5*d*e^2*x^5*(a+b*arcsinh(c*x))+-3//175*b*c^-7*e^2*(1+c^2*x^2)^5//2*(-5*e+7*d*c^2)+-1//105*b*e*c^-7*(1+c^2*x^2)^3//2*(15*e^2+35*c^4*d^2+-42*d*e*c^2))
@test integrate((d+e*x^2)^2*(a+b*arcsinh(c*x)), x) == :(x*d^2*(a+b*arcsinh(c*x))+1//5*e^2*x^5*(a+b*arcsinh(c*x))+-1//15*b*c^-5*(1+c^2*x^2)^(1/2)*(3*e^2+15*c^4*d^2+-10*d*e*c^2)+-1//25*b*c^-5*e^2*(1+c^2*x^2)^5//2+2//3*d*e*x^3*(a+b*arcsinh(c*x))+-2//45*b*e*c^-5*(1+c^2*x^2)^3//2*(-3*e+5*d*c^2))
@test integrate((a+b*arcsinh(c*x))*(d+e*x^2), x) == :(d*x*(a+b*arcsinh(c*x))+1//3*e*x^3*(a+b*arcsinh(c*x))+-1//3*b*c^-3*(1+c^2*x^2)^(1/2)*(-1*e+3*d*c^2)+-1//9*b*e*c^-3*(1+c^2*x^2)^3//2)
@test integrate(a+b*arcsinh(c*x), x) == :(a*x+b*x*arcsinh(c*x)+-1*b*c^-1*(1+c^2*x^2)^(1/2))
@test integrate((d+e*x^2)^-1*(a+b*arcsinh(c*x)), x) == :((1/2)*b*e^-1//2*(-1d)^-1//2*Polylogarithms.polylog(2,e^(1/2)*((e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+(1/2)*b*e^-1//2*(-1d)^-1//2*Polylogarithms.polylog(2,e^(1/2)*(-1*(e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+(1/2)*e^-1//2*(-1d)^-1//2*(a+b*arcsinh(c*x))*log(1+-1*e^(1/2)*((e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+(1/2)*e^-1//2*(-1d)^-1//2*(a+b*arcsinh(c*x))*log(1+-1*e^(1/2)*(-1*(e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+-1//2*b*e^-1//2*(-1d)^-1//2*Polylogarithms.polylog(2,-1*e^(1/2)*((e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+-1//2*b*e^-1//2*(-1d)^-1//2*Polylogarithms.polylog(2,-1*e^(1/2)*(-1*(e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+-1//2*e^-1//2*(-1d)^-1//2*(a+b*arcsinh(c*x))*log(1+e^(1/2)*((e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+-1//2*e^-1//2*(-1d)^-1//2*(a+b*arcsinh(c*x))*log(1+e^(1/2)*(-1*(e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x))))
@test integrate((d+e*x^2)^-2*(a+b*arcsinh(c*x)), x) == :(-1//4*b*e^-1//2*(-1d)^-3//2*Polylogarithms.polylog(2,e^(1/2)*((e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+-1//4*b*e^-1//2*(-1d)^-3//2*Polylogarithms.polylog(2,e^(1/2)*(-1*(e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+-1//4*d^-1*e^-1//2*((-1d)^(1/2)+-1*x*e^(1/2))^-1*(a+b*arcsinh(c*x))+-1//4*e^-1//2*(-1d)^-3//2*(a+b*arcsinh(c*x))*log(1+-1*e^(1/2)*((e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+-1//4*e^-1//2*(-1d)^-3//2*(a+b*arcsinh(c*x))*log(1+-1*e^(1/2)*(-1*(e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+1//4*b*e^-1//2*(-1d)^-3//2*Polylogarithms.polylog(2,-1*e^(1/2)*((e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+1//4*b*e^-1//2*(-1d)^-3//2*Polylogarithms.polylog(2,-1*e^(1/2)*(-1*(e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+1//4*d^-1*e^-1//2*((-1d)^(1/2)+x*e^(1/2))^-1*(a+b*arcsinh(c*x))+1//4*e^-1//2*(-1d)^-3//2*(a+b*arcsinh(c*x))*log(1+e^(1/2)*((e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+1//4*e^-1//2*(-1d)^-3//2*(a+b*arcsinh(c*x))*log(1+e^(1/2)*(-1*(e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+-1//4*b*c*d^-1*e^-1//2*(-1*e+d*c^2)^-1//2*arctan((1+c^2*x^2)^-1//2*(-1*e+d*c^2)^-1//2*(e^(1/2)+x*c^2*(-1d)^(1/2)))+-1//4*b*c*d^-1*e^-1//2*(-1*e+d*c^2)^-1//2*arctan((1+c^2*x^2)^-1//2*(-1*e+d*c^2)^-1//2*(e^(1/2)+-1*x*c^2*(-1d)^(1/2))))
@test integrate((a+b*arcsinh(c*x))^2*(d+e*x^2)^3, x) == :(x*d^3*(a+b*arcsinh(c*x))^2+2*x*b^2*d^3+1//7*e^3*x^7*(a+b*arcsinh(c*x))^2+2//343*b^2*e^3*x^7+e*d^2*x^3*(a+b*arcsinh(c*x))^2+-32//245*x*b^2*c^-6*e^3+-12//1225*b^2*c^-2*e^3*x^5+2//9*e*b^2*d^2*x^3+3//5*d*e^2*x^5*(a+b*arcsinh(c*x))^2+6//125*d*b^2*e^2*x^5+16//735*b^2*c^-4*e^3*x^3+-2*b*c^-1*d^3*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+-8//75*d*b^2*c^-2*e^2*x^3+-4//3*e*x*b^2*c^-2*d^2+16//25*d*x*b^2*c^-4*e^2+32//245*b*c^-7*e^3*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+-16//25*b*d*c^-5*e^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+-16//245*b*c^-5*e^3*x^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+-2//49*b*c^-1*e^3*x^6*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+4//3*b*e*c^-3*d^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+12//245*b*c^-3*e^3*x^4*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+-6//25*b*d*c^-1*e^2*x^4*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+-2//3*b*e*c^-1*d^2*x^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+8//25*b*d*c^-3*e^2*x^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x)))
@test integrate((a+b*arcsinh(c*x))^2*(d+e*x^2)^2, x) == :(x*d^2*(a+b*arcsinh(c*x))^2+2*x*b^2*d^2+1//5*e^2*x^5*(a+b*arcsinh(c*x))^2+2//125*b^2*e^2*x^5+-8//225*b^2*c^-2*e^2*x^3+2//3*d*e*x^3*(a+b*arcsinh(c*x))^2+4//27*d*e*b^2*x^3+16//75*x*b^2*c^-4*e^2+-2*b*c^-1*d^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+-16//75*b*c^-5*e^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+-8//9*d*e*x*b^2*c^-2+-2//25*b*c^-1*e^2*x^4*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+8//9*b*d*e*c^-3*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+8//75*b*c^-3*e^2*x^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+-4//9*b*d*e*c^-1*x^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x)))
@test integrate((a+b*arcsinh(c*x))^2*(d+e*x^2), x) == :(d*x*(a+b*arcsinh(c*x))^2+2*d*x*b^2+1//3*e*x^3*(a+b*arcsinh(c*x))^2+2//27*e*b^2*x^3+-4//9*e*x*b^2*c^-2+-2*b*d*c^-1*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+4//9*b*e*c^-3*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))+-2//9*b*e*c^-1*x^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x)))
@test integrate((a+b*arcsinh(c*x))^2, x) == :(x*(a+b*arcsinh(c*x))^2+2*x*b^2+-2*b*c^-1*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x)))
@test integrate((a+b*arcsinh(c*x))^2*(d+e*x^2)^-1, x) == :(b^2*e^-1//2*(-1d)^-1//2*Polylogarithms.polylog(3,-1*e^(1/2)*((e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+b^2*e^-1//2*(-1d)^-1//2*Polylogarithms.polylog(3,-1*e^(1/2)*(-1*(e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+(1/2)*e^-1//2*(-1d)^-1//2*(a+b*arcsinh(c*x))^2*log(1+-1*e^(1/2)*((e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+(1/2)*e^-1//2*(-1d)^-1//2*(a+b*arcsinh(c*x))^2*log(1+-1*e^(1/2)*(-1*(e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+-1*b^2*e^-1//2*(-1d)^-1//2*Polylogarithms.polylog(3,e^(1/2)*((e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+-1*b^2*e^-1//2*(-1d)^-1//2*Polylogarithms.polylog(3,e^(1/2)*(-1*(e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+-1//2*e^-1//2*(-1d)^-1//2*(a+b*arcsinh(c*x))^2*log(1+e^(1/2)*((e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+-1//2*e^-1//2*(-1d)^-1//2*(a+b*arcsinh(c*x))^2*log(1+e^(1/2)*(-1*(e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+b*e^-1//2*(-1d)^-1//2*(a+b*arcsinh(c*x))*Polylogarithms.polylog(2,e^(1/2)*((e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+b*e^-1//2*(-1d)^-1//2*(a+b*arcsinh(c*x))*Polylogarithms.polylog(2,e^(1/2)*(-1*(e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+-1*b*e^-1//2*(-1d)^-1//2*(a+b*arcsinh(c*x))*Polylogarithms.polylog(2,-1*e^(1/2)*((e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x)))+-1*b*e^-1//2*(-1d)^-1//2*(a+b*arcsinh(c*x))*Polylogarithms.polylog(2,-1*e^(1/2)*(-1*(e+-1*d*c^2)^(1/2)+c*(-1d)^(1/2))^-1*exp(arcsinh(c*x))))
@test integrate((a+b*arcsinh(c*x))^-1*(d+e*x^2)^3, x) == :(b^-1*c^-1*d^3*cosh(a*b^-1)*coshintegral(b^-1*(a+b*arcsinh(c*x)))+-1*b^-1*c^-1*d^3*sinh(a*b^-1)*sinhintegral(b^-1*(a+b*arcsinh(c*x)))+-9//64*b^-1*c^-7*e^3*sinh(3*a*b^-1)*sinhintegral(3*arcsinh(c*x)+3*a*b^-1)+-5//64*b^-1*c^-7*e^3*cosh(a*b^-1)*coshintegral(a*b^-1+arcsinh(c*x))+-5//64*b^-1*c^-7*e^3*cosh(5*a*b^-1)*coshintegral(5*arcsinh(c*x)+5*a*b^-1)+-1//64*b^-1*c^-7*e^3*sinh(7*a*b^-1)*sinhintegral(7*arcsinh(c*x)+7*a*b^-1)+1//64*b^-1*c^-7*e^3*cosh(7*a*b^-1)*coshintegral(7*arcsinh(c*x)+7*a*b^-1)+5//64*b^-1*c^-7*e^3*sinh(a*b^-1)*sinhintegral(a*b^-1+arcsinh(c*x))+5//64*b^-1*c^-7*e^3*sinh(5*a*b^-1)*sinhintegral(5*arcsinh(c*x)+5*a*b^-1)+9//64*b^-1*c^-7*e^3*cosh(3*a*b^-1)*coshintegral(3*arcsinh(c*x)+3*a*b^-1)+-9//16*d*b^-1*c^-5*e^2*cosh(3*a*b^-1)*coshintegral(3*arcsinh(c*x)+3*a*b^-1)+-3//4*e*b^-1*c^-3*d^2*cosh(a*b^-1)*coshintegral(a*b^-1+arcsinh(c*x))+-3//4*e*b^-1*c^-3*d^2*sinh(3*a*b^-1)*sinhintegral(3*arcsinh(c*x)+3*a*b^-1)+-3//8*d*b^-1*c^-5*e^2*sinh(a*b^-1)*sinhintegral(a*b^-1+arcsinh(c*x))+-3//16*d*b^-1*c^-5*e^2*sinh(5*a*b^-1)*sinhintegral(5*arcsinh(c*x)+5*a*b^-1)+3//4*e*b^-1*c^-3*d^2*cosh(3*a*b^-1)*coshintegral(3*arcsinh(c*x)+3*a*b^-1)+3//4*e*b^-1*c^-3*d^2*sinh(a*b^-1)*sinhintegral(a*b^-1+arcsinh(c*x))+3//8*d*b^-1*c^-5*e^2*cosh(a*b^-1)*coshintegral(a*b^-1+arcsinh(c*x))+3//16*d*b^-1*c^-5*e^2*cosh(5*a*b^-1)*coshintegral(5*arcsinh(c*x)+5*a*b^-1)+9//16*d*b^-1*c^-5*e^2*sinh(3*a*b^-1)*sinhintegral(3*arcsinh(c*x)+3*a*b^-1))
@test integrate((a+b*arcsinh(c*x))^-1*(d+e*x^2)^2, x) == :(b^-1*c^-1*d^2*cosh(a*b^-1)*coshintegral(b^-1*(a+b*arcsinh(c*x)))+-1*b^-1*c^-1*d^2*sinh(a*b^-1)*sinhintegral(b^-1*(a+b*arcsinh(c*x)))+-3//16*b^-1*c^-5*e^2*cosh(3*a*b^-1)*coshintegral(3*arcsinh(c*x)+3*a*b^-1)+-1//8*b^-1*c^-5*e^2*sinh(a*b^-1)*sinhintegral(a*b^-1+arcsinh(c*x))+-1//16*b^-1*c^-5*e^2*sinh(5*a*b^-1)*sinhintegral(5*arcsinh(c*x)+5*a*b^-1)+1//8*b^-1*c^-5*e^2*cosh(a*b^-1)*coshintegral(a*b^-1+arcsinh(c*x))+1//16*b^-1*c^-5*e^2*cosh(5*a*b^-1)*coshintegral(5*arcsinh(c*x)+5*a*b^-1)+3//16*b^-1*c^-5*e^2*sinh(3*a*b^-1)*sinhintegral(3*arcsinh(c*x)+3*a*b^-1)+(1/2)*d*e*b^-1*c^-3*cosh(3*a*b^-1)*coshintegral(3*arcsinh(c*x)+3*a*b^-1)+(1/2)*d*e*b^-1*c^-3*sinh(a*b^-1)*sinhintegral(a*b^-1+arcsinh(c*x))+-1//2*d*e*b^-1*c^-3*cosh(a*b^-1)*coshintegral(a*b^-1+arcsinh(c*x))+-1//2*d*e*b^-1*c^-3*sinh(3*a*b^-1)*sinhintegral(3*arcsinh(c*x)+3*a*b^-1))
@test integrate((a+b*arcsinh(c*x))^-1*(d+e*x^2), x) == :(d*b^-1*c^-1*cosh(a*b^-1)*coshintegral(b^-1*(a+b*arcsinh(c*x)))+-1*d*b^-1*c^-1*sinh(a*b^-1)*sinhintegral(b^-1*(a+b*arcsinh(c*x)))+-1//4*e*b^-1*c^-3*cosh(a*b^-1)*coshintegral(a*b^-1+arcsinh(c*x))+-1//4*e*b^-1*c^-3*sinh(3*a*b^-1)*sinhintegral(3*arcsinh(c*x)+3*a*b^-1)+1//4*e*b^-1*c^-3*cosh(3*a*b^-1)*coshintegral(3*arcsinh(c*x)+3*a*b^-1)+1//4*e*b^-1*c^-3*sinh(a*b^-1)*sinhintegral(a*b^-1+arcsinh(c*x)))
@test integrate((a+b*arcsinh(c*x))^-1, x) == :(b^-1*c^-1*cosh(a*b^-1)*coshintegral(b^-1*(a+b*arcsinh(c*x)))+-1*b^-1*c^-1*sinh(a*b^-1)*sinhintegral(b^-1*(a+b*arcsinh(c*x))))
@test integrate((a+b*arcsinh(c*x))^-2*(d+e*x^2)^2, x) == :(b^-2*c^-1*d^2*cosh(a*b^-1)*sinhintegral(a*b^-1+arcsinh(c*x))+-1*b^-1*c^-1*d^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))^-1+-1*b^-2*c^-1*d^2*coshintegral(a*b^-1+arcsinh(c*x))*sinh(a*b^-1)+-9//16*b^-2*c^-5*e^2*cosh(3*a*b^-1)*sinhintegral(3*arcsinh(c*x)+3*a*b^-1)+-5//16*b^-2*c^-5*e^2*coshintegral(5*arcsinh(c*x)+5*a*b^-1)*sinh(5*a*b^-1)+-1//8*b^-2*c^-5*e^2*coshintegral(a*b^-1+arcsinh(c*x))*sinh(a*b^-1)+1//8*b^-2*c^-5*e^2*cosh(a*b^-1)*sinhintegral(a*b^-1+arcsinh(c*x))+5//16*b^-2*c^-5*e^2*cosh(5*a*b^-1)*sinhintegral(5*arcsinh(c*x)+5*a*b^-1)+9//16*b^-2*c^-5*e^2*coshintegral(3*arcsinh(c*x)+3*a*b^-1)*sinh(3*a*b^-1)+(1/2)*d*e*b^-2*c^-3*coshintegral(a*b^-1+arcsinh(c*x))*sinh(a*b^-1)+-1*b^-1*c^-1*e^2*x^4*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))^-1+-3//2*d*e*b^-2*c^-3*coshintegral(3*arcsinh(c*x)+3*a*b^-1)*sinh(3*a*b^-1)+-1//2*d*e*b^-2*c^-3*cosh(a*b^-1)*sinhintegral(a*b^-1+arcsinh(c*x))+3//2*d*e*b^-2*c^-3*cosh(3*a*b^-1)*sinhintegral(3*arcsinh(c*x)+3*a*b^-1)+-2*d*e*b^-1*c^-1*x^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))^-1)
@test integrate((a+b*arcsinh(c*x))^-2*(d+e*x^2), x) == :(d*b^-2*c^-1*cosh(a*b^-1)*sinhintegral(a*b^-1+arcsinh(c*x))+-1*d*b^-1*c^-1*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))^-1+-1*d*b^-2*c^-1*coshintegral(a*b^-1+arcsinh(c*x))*sinh(a*b^-1)+-3//4*e*b^-2*c^-3*coshintegral(3*arcsinh(c*x)+3*a*b^-1)*sinh(3*a*b^-1)+-1//4*e*b^-2*c^-3*cosh(a*b^-1)*sinhintegral(a*b^-1+arcsinh(c*x))+1//4*e*b^-2*c^-3*coshintegral(a*b^-1+arcsinh(c*x))*sinh(a*b^-1)+3//4*e*b^-2*c^-3*cosh(3*a*b^-1)*sinhintegral(3*arcsinh(c*x)+3*a*b^-1)+-1*e*b^-1*c^-1*x^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))^-1)
@test integrate((a+b*arcsinh(c*x))^-2, x) == :(b^-2*c^-1*cosh(a*b^-1)*sinhintegral(a*b^-1+arcsinh(c*x))+-1*b^-1*c^-1*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))^-1+-1*b^-2*c^-1*coshintegral(a*b^-1+arcsinh(c*x))*sinh(a*b^-1))
@test integrate((a+b*arcsinh(c*x))^(1/2)*(d+e*x^2)^2, x) == :(x*d^2*(a+b*arcsinh(c*x))^(1/2)+1//5*e^2*x^5*(a+b*arcsinh(c*x))^(1/2)+2//3*d*e*x^3*(a+b*arcsinh(c*x))^(1/2)+-1//4*pi^(1/2)*b^(1/2)*c^-1*d^2*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+-1//32*pi^(1/2)*b^(1/2)*c^-5*e^2*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+1//4*pi^(1/2)*b^(1/2)*c^-1*d^2*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+1//32*pi^(1/2)*b^(1/2)*c^-5*e^2*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+-1//8*d*e*pi^(1/2)*b^(1/2)*c^-3*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+-1//192*3^(1/2)*pi^(1/2)*b^(1/2)*c^-5*e^2*erf(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(3*a*b^-1)+-1//1600*5^(1/2)*pi^(1/2)*b^(1/2)*c^-5*e^2*erfi(5^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-5*a*b^-1)+1//8*d*e*pi^(1/2)*b^(1/2)*c^-3*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+1//192*3^(1/2)*pi^(1/2)*b^(1/2)*c^-5*e^2*erfi(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-3*a*b^-1)+1//1600*5^(1/2)*pi^(1/2)*b^(1/2)*c^-5*e^2*erf(5^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(5*a*b^-1)+-1//72*d*e*3^(1/2)*pi^(1/2)*b^(1/2)*c^-3*erfi(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-3*a*b^-1)+1//72*d*e*3^(1/2)*pi^(1/2)*b^(1/2)*c^-3*erf(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(3*a*b^-1))
@test integrate((a+b*arcsinh(c*x))^(1/2)*(d+e*x^2), x) == :(d*x*(a+b*arcsinh(c*x))^(1/2)+1//3*e*x^3*(a+b*arcsinh(c*x))^(1/2)+-1//4*d*pi^(1/2)*b^(1/2)*c^-1*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+-1//16*e*pi^(1/2)*b^(1/2)*c^-3*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+1//4*d*pi^(1/2)*b^(1/2)*c^-1*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+1//16*e*pi^(1/2)*b^(1/2)*c^-3*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+-1//144*e*3^(1/2)*pi^(1/2)*b^(1/2)*c^-3*erfi(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-3*a*b^-1)+1//144*e*3^(1/2)*pi^(1/2)*b^(1/2)*c^-3*erf(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(3*a*b^-1))
@test integrate((a+b*arcsinh(c*x))^(1/2), x) == :(x*(a+b*arcsinh(c*x))^(1/2)+-1//4*pi^(1/2)*b^(1/2)*c^-1*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+1//4*pi^(1/2)*b^(1/2)*c^-1*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1))
@test integrate((a+b*arcsinh(c*x))^3//2*(d+e*x^2), x) == :(d*x*(a+b*arcsinh(c*x))^3//2+1//3*e*x^3*(a+b*arcsinh(c*x))^3//2+-3//2*b*d*c^-1*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))^(1/2)+1//3*b*e*c^-3*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))^(1/2)+-3//32*e*pi^(1/2)*b^3//2*c^-3*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+-3//32*e*pi^(1/2)*b^3//2*c^-3*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+-1//6*b*e*c^-1*x^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))^(1/2)+3//8*d*pi^(1/2)*b^3//2*c^-1*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+3//8*d*pi^(1/2)*b^3//2*c^-1*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+1//288*e*3^(1/2)*pi^(1/2)*b^3//2*c^-3*erf(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(3*a*b^-1)+1//288*e*3^(1/2)*pi^(1/2)*b^3//2*c^-3*erfi(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-3*a*b^-1))
@test integrate((a+b*arcsinh(c*x))^3//2, x) == :(x*(a+b*arcsinh(c*x))^3//2+-3//2*b*c^-1*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))^(1/2)+3//8*pi^(1/2)*b^3//2*c^-1*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+3//8*pi^(1/2)*b^3//2*c^-1*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1))
@test integrate((a+b*arcsinh(c*x))^-1//2*(d+e*x^2)^2, x) == :((1/2)*pi^(1/2)*b^-1//2*c^-1*d^2*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+(1/2)*pi^(1/2)*b^-1//2*c^-1*d^2*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+1//16*pi^(1/2)*b^-1//2*c^-5*e^2*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+1//16*pi^(1/2)*b^-1//2*c^-5*e^2*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+-1//4*d*e*pi^(1/2)*b^-1//2*c^-3*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+-1//4*d*e*pi^(1/2)*b^-1//2*c^-3*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+-1//32*3^(1/2)*pi^(1/2)*b^-1//2*c^-5*e^2*erf(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(3*a*b^-1)+-1//32*3^(1/2)*pi^(1/2)*b^-1//2*c^-5*e^2*erfi(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-3*a*b^-1)+1//160*5^(1/2)*pi^(1/2)*b^-1//2*c^-5*e^2*erf(5^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(5*a*b^-1)+1//160*5^(1/2)*pi^(1/2)*b^-1//2*c^-5*e^2*erfi(5^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-5*a*b^-1)+1//12*d*e*3^(1/2)*pi^(1/2)*b^-1//2*c^-3*erf(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(3*a*b^-1)+1//12*d*e*3^(1/2)*pi^(1/2)*b^-1//2*c^-3*erfi(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-3*a*b^-1))
@test integrate((a+b*arcsinh(c*x))^-1//2*(d+e*x^2), x) == :((1/2)*d*pi^(1/2)*b^-1//2*c^-1*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+(1/2)*d*pi^(1/2)*b^-1//2*c^-1*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+-1//8*e*pi^(1/2)*b^-1//2*c^-3*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+-1//8*e*pi^(1/2)*b^-1//2*c^-3*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+1//24*e*3^(1/2)*pi^(1/2)*b^-1//2*c^-3*erf(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(3*a*b^-1)+1//24*e*3^(1/2)*pi^(1/2)*b^-1//2*c^-3*erfi(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-3*a*b^-1))
@test integrate((a+b*arcsinh(c*x))^-1//2, x) == :((1/2)*pi^(1/2)*b^-1//2*c^-1*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+(1/2)*pi^(1/2)*b^-1//2*c^-1*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1))
@test integrate((a+b*arcsinh(c*x))^-3//2*(d+e*x^2), x) == :(-2*d*b^-1*c^-1*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))^-1//2+d*pi^(1/2)*b^-3//2*c^-1*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+-1*d*pi^(1/2)*b^-3//2*c^-1*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+-2*e*b^-1*c^-1*x^2*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))^-1//2+-1//4*e*pi^(1/2)*b^-3//2*c^-3*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+1//4*e*pi^(1/2)*b^-3//2*c^-3*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1)+-1//4*e*3^(1/2)*pi^(1/2)*b^-3//2*c^-3*erf(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(3*a*b^-1)+1//4*e*3^(1/2)*pi^(1/2)*b^-3//2*c^-3*erfi(3^(1/2)*b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-3*a*b^-1))
@test integrate((a+b*arcsinh(c*x))^-3//2, x) == :(-2*b^-1*c^-1*(1+c^2*x^2)^(1/2)*(a+b*arcsinh(c*x))^-1//2+pi^(1/2)*b^-3//2*c^-1*erfi(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(-1*a*b^-1)+-1*pi^(1/2)*b^-3//2*c^-1*erf(b^-1//2*(a+b*arcsinh(c*x))^(1/2))*exp(a*b^-1))
@test integrate((d+e*x^2)^-3//2*(a+b*arcsinh(c*x)), x) == :(x*d^-1*(d+e*x^2)^-1//2*(a+b*arcsinh(c*x))+-1*b*d^-1*e^-1//2*arctanh(c^-1*e^(1/2)*(1+c^2*x^2)^(1/2)*(d+e*x^2)^-1//2))
@test integrate((d+e*x^2)^-5//2*(a+b*arcsinh(c*x)), x) == :(-2//3*b*d^-2*e^-1//2*arctanh(c^-1*e^(1/2)*(1+c^2*x^2)^(1/2)*(d+e*x^2)^-1//2)+1//3*x*d^-1*(d+e*x^2)^-3//2*(a+b*arcsinh(c*x))+2//3*x*d^-2*(d+e*x^2)^-1//2*(a+b*arcsinh(c*x))+-1//3*b*c*d^-1*(1+c^2*x^2)^(1/2)*(d+e*x^2)^-1//2*(-1*e+d*c^2)^-1)
@test integrate((d+e*x^2)^-7//2*(a+b*arcsinh(c*x)), x) == :(-8//15*b*d^-3*e^-1//2*arctanh(c^-1*e^(1/2)*(1+c^2*x^2)^(1/2)*(d+e*x^2)^-1//2)+1//5*x*d^-1*(d+e*x^2)^-5//2*(a+b*arcsinh(c*x))+4//15*x*d^-2*(d+e*x^2)^-3//2*(a+b*arcsinh(c*x))+8//15*x*d^-3*(d+e*x^2)^-1//2*(a+b*arcsinh(c*x))+-1//15*b*c*d^-1*(1+c^2*x^2)^(1/2)*(d+e*x^2)^-3//2*(-1*e+d*c^2)^-1+-2//15*b*c*d^-2*(1+c^2*x^2)^(1/2)*(d+e*x^2)^-1//2*(-1*e+d*c^2)^-2*(-2*e+3*d*c^2))

max_stars_repo_path: test/7 Inverse hyperbolic functions/7.1 Inverse hyperbolic sine/7.1.4b (f x)^m (d+e x^2)^p (a+b arcsinh(c x))^n.jl


text:
#--------------------------------------------------------------------
# Note: some of this code was copied/adapted from JuliaIO/ProtoBuf.jl which has the following license:
#
#   The ProtoBuf.jl package is licensed under the MIT "Expat" License:
#
#   Copyright (c) 2014: Tanmay Mohapatra.
#
#   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.
#--------------------------------------------------------------------

module Serialization

export Proto, wiretype, jtype, is_primitive, writeproto, defaultval

const MSB = 0x80
const MASK7 = 0x7f
const MASK8 = 0xff
const MASK3 = 0x07

const WIRETYP_VARINT   = 0
const WIRETYP_64BIT    = 1
const WIRETYP_LENDELIM = 2
const WIRETYP_GRPSTART = 3   # deprecated
const WIRETYP_GRPEND   = 4   # deprecated
const WIRETYP_32BIT    = 5

"""
The abstract type from which all generated protobuf structs extend.
"""
abstract type Proto end

proto_types(::Type{Int32})                            = [:int32, :sint32, :enum, :sfixed32]
proto_types(::Type{Int64})                            = [:int64, :sint64, :sfixed64]
proto_types(::Type{UInt32})                           = [:uint32, :fixed32]
proto_types(::Type{UInt64})                           = [:uint64, :fixed64]
proto_types(::Type{Bool})                             = [:bool]
proto_types(::Type{Float64})                          = [:double]
proto_types(::Type{Float32})                          = [:float]
proto_types(::Type{T}) where {T<:AbstractString}      = [:string]
proto_types(::Type{Vector{UInt8}})                    = [:bytes]
proto_types(::Type{Dict{K,V}}) where {K,V}            = [:map]
proto_types(::Type)                                   = [:obj]
proto_types(::Type{Vector{T}}) where {T}              = proto_types(T)

proto_type(::Type{T}) where {T}                       = proto_types(T)[1]

const WIRETYPES = Dict{Symbol, Int64}(
    :int32          => WIRETYP_VARINT,
    :int64          => WIRETYP_VARINT,
    :uint32         => WIRETYP_VARINT,
    :uint64         => WIRETYP_VARINT,
    :sint32         => WIRETYP_VARINT,
    :sint64         => WIRETYP_VARINT,
    :bool           => WIRETYP_VARINT,
    :enum           => WIRETYP_VARINT,
    :fixed64        => WIRETYP_64BIT,
    :sfixed64       => WIRETYP_64BIT,
    :double         => WIRETYP_64BIT,
    :string         => WIRETYP_LENDELIM,
    :bytes          => WIRETYP_LENDELIM,
    :map            => WIRETYP_LENDELIM,
    :fixed32        => WIRETYP_32BIT,
    :sfixed32       => WIRETYP_32BIT,
    :float          => WIRETYP_32BIT
)

wiretype(s::AbstractString) = wiretype(Symbol(s))
wiretype(s::Symbol) = haskey(WIRETYPES, s) ? WIRETYPES[s] : WIRETYP_LENDELIM

const JTYPES = Dict{Symbol, Symbol}(
    :int32    => :Int32,
    :int64    => :Int64,
    :uint32   => :UInt32,
    :uint64   => :UInt64,
    :sint32   => :Int32,
    :sint64   => :Int64,
    :bool     => :Bool,
    :enum     => :Int32,
    :fixed64  => :UInt64,
    :sfixed64 => :Int64,
    :double   => :Float64,
    :string   => :AbstractString,
    :bytes    => Symbol("Vector{UInt8}"),
    # :map      => :Dict,
    :fixed32  => :UInt32,
    :sfixed32 => :Int32,
    :float    => :Float32
)

jtype(::Missing) = :Missing
jtype(s::AbstractString) = jtype(Symbol(s))
jtype(s::Symbol) = haskey(JTYPES, s) ? JTYPES[s] : s

is_primitive(_type) = haskey(JTYPES, Symbol(_type))

defaultval(::Type{T}) where {T<:Number}             = [zero(T)]
defaultval(::Type{T}) where {T<:AbstractString}     = [convert(T,"")]
defaultval(::Type{Bool})                            = [false]
defaultval(::Type{Vector{T}}) where {T}             = Any[T[]]
defaultval(::Type{Dict{K,V}}) where {K,V}           = [Dict{K,V}()]
defaultval(::Type)                                  = []

function _write_uleb(io::IO, x::T) where T <: Integer
    nw = 0
    cont = true
    while cont
        byte = x & MASK7
        if (x >>>= 7) != 0
            byte |= MSB
        else
            cont = false
        end
        nw += write(io, UInt8(byte))
    end
    nw
end

# max number of 7bit blocks for reading n bytes
# d,r = divrem(sizeof(T)*8, 7)
# (r > 0) && (d += 1)
const _max_n = [2, 3, 4, 5, 6, 7, 8, 10]

function _read_uleb_base(io::IO, ::Type{T}) where T <: Integer
    res = zero(T)
    n = 0
    byte = UInt8(MSB)
    while (byte & MSB) != 0
        byte = read(io, UInt8)
        res |= (convert(T, byte & MASK7) << (7*n))
        n += 1
    end
    n, res
end

function _read_uleb(io::IO, ::Type{Int32})
    n, res = _read_uleb_base(io, Int32)

    # negative int32 are encoded in 10 bytes (ref: https://developers.google.com/protocol-buffers/docs/encoding)
    # > if you use int32 or int64 as the type for a negative number, the resulting varint is always ten bytes long
    #
    # but Julia can be tolerant like the C protobuf implementation (unlike python)

    if n > _max_n[sizeof(res < 0 ? Int64 : Int32)]
        @debug("overflow reading Int32. returning 0")
        return Int32(0)
    end

    res
end

function _read_uleb(io::IO, ::Type{T}) where T <: Integer
    n, res = _read_uleb_base(io, T)
    # in case of overflow, consider it as missing field and return default value
    if n > _max_n[sizeof(T)]
        @debug("overflow reading integer type. returning 0", T)
        return zero(T)
    end
    res
end

function _write_zigzag(io::IO, x::T) where T <: Integer
    nbits = 8*sizeof(x)
    zx = (x << 1) ⊻ (x >> (nbits-1))
    _write_uleb(io, zx)
end

function _read_zigzag(io::IO, ::Type{T}) where T <: Integer
    zx = _read_uleb(io, UInt64)
    # result is positive if zx is even
    convert(T, iseven(zx) ? (zx >>> 1) : -signed((zx+1) >>> 1))
end

##
# read and write field keys
_write_key(io::IO, fld::Int, wiretyp::Int) = _write_uleb(io, (fld << 3) | wiretyp)
function _read_key(io::IO)
    key = _read_uleb(io, UInt64)
    wiretyp = key & MASK3
    fld = key >>> 3
    (fld, wiretyp)
end

##
# read and write field values
write_varint(io::IO, x::T) where {T <: Integer} = _write_uleb(io, x)
write_varint(io::IO, x::Int32) = _write_uleb(io, x < 0 ? Int64(x) : x)
write_bool(io::IO, x::Bool) = _write_uleb(io, x ? 1 : 0)
write_svarint(io::IO, x::T) where {T <: Integer} = _write_zigzag(io, x)

read_varint(io::IO, ::Type{T}) where {T <: Integer} = _read_uleb(io, T)
read_bool(io::IO) = Bool(_read_uleb(io, UInt64))
read_bool(io::IO, ::Type{Bool}) = read_bool(io)
read_svarint(io::IO, ::Type{T}) where {T <: Integer} = _read_zigzag(io, T)

write_fixed(io::IO, x::UInt32) = _write_fixed(io, x)
write_fixed(io::IO, x::UInt64) = _write_fixed(io, x)
write_fixed(io::IO, x::Int32) = _write_fixed(io, reinterpret(UInt32, x))
write_fixed(io::IO, x::Int64) = _write_fixed(io, reinterpret(UInt64, x))
write_fixed(io::IO, x::Float32) = _write_fixed(io, reinterpret(UInt32, x))
write_fixed(io::IO, x::Float64) = _write_fixed(io, reinterpret(UInt64, x))
function _write_fixed(io::IO, ux::T) where T <: Unsigned
    N = sizeof(ux)
    for n in 1:N
        write(io, UInt8(ux & MASK8))
        ux >>>= 8
    end
    N
end

read_fixed(io::IO, typ::Type{UInt32}) = _read_fixed(io, convert(UInt32,0), 4)
read_fixed(io::IO, typ::Type{UInt64}) = _read_fixed(io, convert(UInt64,0), 8)
read_fixed(io::IO, typ::Type{Int32}) = reinterpret(Int32, _read_fixed(io, convert(UInt32,0), 4))
read_fixed(io::IO, typ::Type{Int64}) = reinterpret(Int64, _read_fixed(io, convert(UInt64,0), 8))
read_fixed(io::IO, typ::Type{Float32}) = reinterpret(Float32, _read_fixed(io, convert(UInt32,0), 4))
read_fixed(io::IO, typ::Type{Float64}) = reinterpret(Float64, _read_fixed(io, convert(UInt64,0), 8))
function _read_fixed(io::IO, ret::T, N::Int) where T <: Unsigned
    for n in 0:(N-1)
        byte = convert(T, read(io, UInt8))
        ret |= (byte << (8*n))
    end
    ret
end

function write_bytes(io::IO, data::Vector{UInt8})
    n = _write_uleb(io, sizeof(data))
    n += write(io, data)
    n
end

function read_bytes(io::IO)
    n = _read_uleb(io, UInt64)
    data = Vector{UInt8}(undef, n)
    read!(io, data)
    data
end
read_bytes(io::IO, ::Type{Vector{UInt8}}) = read_bytes(io)

write_string(io::IO, x::AbstractString) = write_string(io, String(x))
write_string(io::IO, x::String) = write_bytes(io, @static isdefined(Base, :codeunits) ? unsafe_wrap(Vector{UInt8}, x) : Vector{UInt8}(x))

read_string(io::IO) = String(read_bytes(io))
read_string(io::IO, ::Type{T}) where {T <: AbstractString} = convert(T, read_string(io))

function write_map(io::IO, fldnum::Int, dict::Dict)
    dmeta = mapentry_meta(typeof(dict))
    iob = IOBuffer()

    n = 0
    for key in keys(dict)
        @debug("write_map", key)
        val = dict[key]
        writeproto(iob, key, dmeta.ordered[1])
        @debug("write_map", val)
        writeproto(iob, val, dmeta.ordered[2])
        n += _write_key(io, fldnum, WIRETYP_LENDELIM)
        n += write_bytes(io, take!(iob))
    end
    n
end

function read_map(io, dict::Dict{K,V}) where {K,V}
    iob = IOBuffer(read_bytes(io))

    dmeta = mapentry_meta(Dict{K,V})
    key_wtyp, key_wfn, key_rfn, key_jtyp = WIRETYPES[dmeta.numdict[1].ptyp]
    val_wtyp, val_wfn, val_rfn, val_jtyp = WIRETYPES[dmeta.numdict[2].ptyp]
    key_val = Vector{Union{K,V}}(undef, 2)

    while !eof(iob)
        fldnum, wiretyp = _read_key(iob)
        @debug("reading map", fldnum)

        fldnum = Int(fldnum)
        attrib = dmeta.numdict[fldnum]

        if fldnum == 1
            key_val[1] = read_field(iob, nothing, attrib, key_wtyp, K)
        elseif fldnum == 2
            key_val[2] = read_field(iob, nothing, attrib, val_wtyp, V)
        else
            skip_field(iob, wiretyp)
        end
    end
    @debug("read map", key=key_val[1], val=key_val[2])
    dict[key_val[1]] = key_val[2]
    dict
end

# const WIRETYPES = Dict{Symbol,Tuple}(
#     :int32          => (WIRETYP_VARINT,     write_varint,  read_varint,   Int32),
#     :int64          => (WIRETYP_VARINT,     write_varint,  read_varint,   Int64),
#     :uint32         => (WIRETYP_VARINT,     write_varint,  read_varint,   UInt32),
#     :uint64         => (WIRETYP_VARINT,     write_varint,  read_varint,   UInt64),
#     :sint32         => (WIRETYP_VARINT,     write_svarint, read_svarint,  Int32),
#     :sint64         => (WIRETYP_VARINT,     write_svarint, read_svarint,  Int64),
#     :bool           => (WIRETYP_VARINT,     write_bool,    read_bool,     Bool),
#     :enum           => (WIRETYP_VARINT,     write_varint,  read_varint,   Int32),

#     :fixed64        => (WIRETYP_64BIT,      write_fixed,   read_fixed,    UInt64),
#     :sfixed64       => (WIRETYP_64BIT,      write_fixed,   read_fixed,    Int64),
#     :double         => (WIRETYP_64BIT,      write_fixed,   read_fixed,    Float64),

#     :string         => (WIRETYP_LENDELIM,   write_string,  read_string,   AbstractString),
#     :bytes          => (WIRETYP_LENDELIM,   write_bytes,   read_bytes,    Vector{UInt8}),
#     # :obj            => (WIRETYP_LENDELIM,   writeproto,    readproto,     Any),
#     :map            => (WIRETYP_LENDELIM,   write_map,     read_map,      Dict),

#     :fixed32        => (WIRETYP_32BIT,      write_fixed,   read_fixed,    UInt32),
#     :sfixed32       => (WIRETYP_32BIT,      write_fixed,   read_fixed,    Int32),
#     :float          => (WIRETYP_32BIT,      write_fixed,   read_fixed,    Float32)
# )

writefield(io::IO, idx::Int64, val, ::Val{:int32})    = writefield(io, idx, val, WIRETYP_VARINT,     write_varint)
writefield(io::IO, idx::Int64, val, ::Val{:int64})    = writefield(io, idx, val, WIRETYP_VARINT,     write_varint)
writefield(io::IO, idx::Int64, val, ::Val{:uint32})   = writefield(io, idx, val, WIRETYP_VARINT,     write_varint)
writefield(io::IO, idx::Int64, val, ::Val{:uint64})   = writefield(io, idx, val, WIRETYP_VARINT,     write_varint)
writefield(io::IO, idx::Int64, val, ::Val{:sint32})   = writefield(io, idx, val, WIRETYP_VARINT,     write_svarint)
writefield(io::IO, idx::Int64, val, ::Val{:sint64})   = writefield(io, idx, val, WIRETYP_VARINT,     write_svarint)
writefield(io::IO, idx::Int64, val, ::Val{:bool})     = writefield(io, idx, val, WIRETYP_VARINT,     write_bool)
writefield(io::IO, idx::Int64, val, ::Val{:enum})     = writefield(io, idx, val, WIRETYP_VARINT,     write_varint)

writefield(io::IO, idx::Int64, val, ::Val{:fixed64})  = writefield(io, idx, val, WIRETYP_64BIT,      write_fixed)
writefield(io::IO, idx::Int64, val, ::Val{:sfixed64}) = writefield(io, idx, val, WIRETYP_64BIT,      write_fixed)
writefield(io::IO, idx::Int64, val, ::Val{:double})   = writefield(io, idx, val, WIRETYP_64BIT,      write_fixed)

writefield(io::IO, idx::Int64, val, ::Val{:string})   = writefield(io, idx, val, WIRETYP_LENDELIM,   write_string)
writefield(io::IO, idx::Int64, val, ::Val{:bytes})    = writefield(io, idx, val, WIRETYP_LENDELIM,   write_bytes)
# writefield(io::IO, idx::Int64, val, ::Val{:obj})      = writefield(io, idx, val, WIRETYP_LENDELIM,   writeproto)
writefield(io::IO, idx::Int64, val, ::Val{:map})      = writefield(io, idx, val, WIRETYP_LENDELIM,   write_map)

writefield(io::IO, idx::Int64, val, ::Val{:fixed32})  = writefield(io, idx, val, WIRETYP_32BIT,      write_fixed)
writefield(io::IO, idx::Int64, val, ::Val{:sfixed32}) = writefield(io, idx, val, WIRETYP_32BIT,      write_fixed)
writefield(io::IO, idx::Int64, val, ::Val{:float})    = writefield(io, idx, val, WIRETYP_32BIT,      write_fixed)

function writefield(io::IO, idx::Int64, val, wiretype::Int64, writefn::Function)
	n = 0
	n += _write_key(io, idx, wiretype)
	n += writefn(io, val)
	n
end

"""
This serializes a Proto object to an IOStream.
"""
function writeproto end

# this is a helper function to write to a temp buffer
function writeproto(x)
	io = IOBuffer()
	writeproto(io, x)
	take!(io)
end

readfield(io::IO, ::Val{:int32})    = read_varint(io, Int32)
readfield(io::IO, ::Val{:int64})    = read_varint(io, Int64)
readfield(io::IO, ::Val{:uint32})   = read_varint(io, UInt32)
readfield(io::IO, ::Val{:uint64})   = read_varint(io, UInt64)
readfield(io::IO, ::Val{:sint32})   = read_svarint(io, Int32)
readfield(io::IO, ::Val{:sint64})   = read_svarint(io, Int64)
readfield(io::IO, ::Val{:bool})     = read_bool(io, Bool)
readfield(io::IO, ::Val{:enum})     = read_varint(io, Int32)

readfield(io::IO, ::Val{:fixed64})  = read_fixed(io, UInt64)
readfield(io::IO, ::Val{:sfixed64}) = read_fixed(io, Int64)
readfield(io::IO, ::Val{:double})   = read_fixed(io, Float64)

readfield(io::IO, ::Val{:string})   = read_string(io, AbstractString)
readfield(io::IO, ::Val{:bytes})    = read_bytes(io, Vector{UInt8})
# readfield(io::IO, ::Val{:obj})      = readproto(io, Any)
readfield(io::IO, ::Val{:map})      = read_map(io, Dict)

readfield(io::IO, ::Val{:fixed32})  = read_fixed(io, UInt32)
readfield(io::IO, ::Val{:sfixed32}) = read_fixed(io, Int32)
readfield(io::IO, ::Val{:float})    = read_fixed(io, Float32)

end # module

ext: jl


text:
"""
QQQQ
"""
function hypoexp_test()
    mean_desired = 5.7
    scv_desired = 0.23
    d = hypo_exp_init(mean_desired, scv_desired)
    (mean_desired ≈ mean(d)) || (return false)
    (scv_desired ≈ scv(d)) || (return false)
    return true
end

function hyperexp_test()
    # mean_desired = 2.7
    # scv_desired = 1.5
    # d = hyper_exp_init(mean_desired, scv_desired)
    # (mean_desired ≈ mean(d)) || (return false)
    # (scv_desired ≈ scv(d)) || (return false)
    return true
end


text:
using LoopVectorization
using Plots
using UnPack

using StartUpDG
using StartUpDG.ExplicitTimestepUtils

N = 3
K1D = 16
T = .5
CFL = 1.

# init ref element and mesh
rd = init_reference_tri(N)
VX,VY,EToV = uniform_tri_mesh(K1D)
md = init_DG_mesh(VX,VY,EToV,rd)

# Make domain periodic
LX,LY = (x->maximum(x)-minimum(x)).((VX,VY)) # find lengths of domain
build_periodic_boundary_maps!(md,rd,LX,LY)

#####
##### Define initial coefficients and time-stepping
#####

# "Define the initial conditions by interpolation"
@unpack x,y = md
p = @. exp(-100*(x^2+y^2))
u = @. 0*x
v = @. 0*x

# "Define function to evaluate the RHS: Q = (p,u,v)"
function rhs(Q, rd::RefElemData, md::MeshData)

    @unpack Dr,Ds,LIFT,Vf = rd
    @unpack rxJ,sxJ,ryJ,syJ,J = md
    @unpack nxJ,nyJ,sJ = md
    @unpack mapP,mapB = md

    Qf = (x->Vf*x).(Q)
    QP = (x->x[mapP]).(Qf)
    dp,du,dv = QP.-Qf

    tau = .5
    dun   = @. (du*nxJ + dv*nyJ)/sJ
    pflux = @. .5*(du*nxJ + dv*nyJ) - tau*dp*sJ
    uflux = @. .5*dp*nxJ - tau*dun*nxJ
    vflux = @. .5*dp*nyJ - tau*dun*nyJ

    pr,ur,vr = (x->Dr*x).(Q)
    ps,us,vs = (x->Ds*x).(Q)
    dpdx = @. rxJ*pr + sxJ*ps
    dpdy = @. ryJ*pr + syJ*ps
    dudx = @. rxJ*ur + sxJ*us
    dvdy = @. ryJ*vr + syJ*vs
    rhsp = dudx + dvdy + LIFT*pflux
    rhsu = dpdx        + LIFT*uflux
    rhsv = dpdy        + LIFT*vflux

    return map(x->-x./J,(rhsp,rhsu,rhsv))
end

function rhs_avx(Q, rd::RefElemData, md::MeshData)

    @unpack Dr,Ds,LIFT,Vf = rd
    @unpack rxJ,sxJ,ryJ,syJ,J = md
    @unpack nxJ,nyJ,sJ = md
    @unpack mapP,mapB = md

    Qf = (x->Vf*x).(Q)
    QP = (x->x[mapP]).(Qf)
    dp,du,dv = QP.-Qf

    tau = .5
    dun   = @avx @. (du*nxJ + dv*nyJ)/sJ
    pflux = @avx @. .5*(du*nxJ + dv*nyJ) - tau*dp*sJ
    uflux = @avx @. .5*dp*nxJ - tau*dun*nxJ
    vflux = @avx @. .5*dp*nyJ - tau*dun*nyJ

    pr,ur,vr = (x->Dr*x).(Q)
    ps,us,vs = (x->Ds*x).(Q)
    dpdx = @avx @. rxJ*pr + sxJ*ps
    dpdy = @avx @. ryJ*pr + syJ*ps
    dudx = @avx @. rxJ*ur + sxJ*us
    dvdy = @avx @. ryJ*vr + syJ*vs
    rhsp = dudx + dvdy + LIFT*pflux
    rhsu = dpdx        + LIFT*uflux
    rhsv = dpdy        + LIFT*vflux

    return (x->@avx -x./J).((rhsp,rhsu,rhsv))
end

# "Time integration coefficients"
rk4a,rk4b,rk4c = ck45()
CN = (N+1)*(N+2)/2  # estimated trace constant
dt = CFL * 2 / (CN*K1D)
Nsteps = convert(Int,ceil(T/dt))
dt = T/Nsteps

# "Perform time-stepping"
Q = (p,u,v) # use mutable arrays for Q, resQ for scoping
resQ = zero.(Q)
for i = 1:Nsteps
    for INTRK = 1:5
        rhsQ = rhs_avx(Q,rd,md)
        bcopy!.(resQ, @. rk4a[INTRK]*resQ + dt*rhsQ)
        bcopy!.(Q,@. Q + rk4b[INTRK]*resQ)
    end

    if i%10==0 || i==Nsteps
        println("Time step $i out of $Nsteps")
    end
end

# # dopri initialization
# rka,rkE,rkc = dp56()
# PIparams = init_PI_controller(5)
# Qtmp = similar.(Q)
# rhsQrk = ntuple(x->zero.(Q),length(rkE))
# prevErrEst = nothing
# rhsQ = rhs(Q,rd,md)
# bcopy!.(rhsQrk[1],rhsQ) # initialize DOPRI rhs (FSAL property)
#
# t = 0.0
# i = 0
# dthist = [dt]
#
# while t < T
#     for INTRK = 2:7
#         k = zero.(Qtmp)
#         for s = 1:INTRK-1
#             bcopy!.(k, @. k + rka[INTRK,s]*rhsQrk[s])
#         end
#         bcopy!.(Qtmp, @. Q + dt*k)
#         rhsQ = rhs_avx(Qtmp,rd,md)
#         bcopy!.(rhsQrk[INTRK],rhsQ)
#     end
#
#     global t,dt,i,prevErrEst # modify compute_adaptive_dt!(...prevErrEst)
#     accept_step, dtnew, prevErrEst = compute_adaptive_dt(
#                                         Q,rhsQrk,dt,rkE,PIparams,prevErrEst)
#     if accept_step
#         t += dt
#         bcopy!.(Q, Qtmp)
#         bcopy!.(rhsQrk[1], rhsQrk[7]) # use FSAL property
#         push!(dthist,dt) # store dt history
#     end
#     dt = min(T-t,dtnew)
#     i = i + 1  # number of total steps attempted
#
#     if i%10==0
#         println("i = $i, t = $t, dt = $dtnew, errEst = $prevErrEst")
#     end
# end

# plotting nodes
gr(aspect_ratio=1,legend=false)
using NodesAndModes.Tri
rp,sp = equi_nodes_2D(10)
Vp = vandermonde_2D(N,rp,sp)/rd.VDM
scatter((x->Vp*x).((x,y,Q[1])),zcolor=Vp*Q[1], msw=0,ms=2,cam=(0,90))

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
push!(LOAD_PATH, "../src")
using Base.Test
using PumiInterface
using ODLCommonTools
using SummationByParts
using PdePumiInterface

function not_isapprox(args...; kwargs...)
  return !isapprox(args...; kwargs...)
end

include("defs.jl")
include("common_functions.jl")
include("test_funcs.jl")
include("test_adapt.jl")
include("test_parallel_types.jl")

@testset "----- Testing 4 process PDEPumiInterface3DG -----" begin
  degree = 1
  Tsbp = Float64
  sbp = getTetSBPOmega(degree=degree)
#  sbp = TetSBP{Tsbp}(degree=degree, internal=true)
  ref_verts = sbp.vtx
  interp_op = SummationByParts.buildinterpolation(sbp, ref_verts.')
  face_verts = SummationByParts.SymCubatures.getfacevertexindices(sbp.cub)
  topo = ElementTopology{3}(face_verts)
  sbpface = TetFace{Tsbp}(degree, sbp.cub, ref_verts)

  dmg_name = ".null"
  smb_name = "pcube10.smb"

  opts = PdePumiInterface.get_defaults()
  opts["dmg_name"] = dmg_name
  opts["smb_name"] = smb_name
  opts["use_linear_metrics"] = true
  opts["coloring_distance"] = 2
  opts["order"] = degree
  opts["numBC"] = 1
  opts["BC1"] = [0,1,2,3,4,5]

#  interp_op = eye(4)
  mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo)
#  mesh = PumiMeshDG3{Float64, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)

  # check coloring
  @test  mesh.maxColors  < 18
  @test  mesh.numColors  < mesh.maxColors+1
  @test  mesh.numColors  > 6

  for i=1:mesh.numEl
    el_i = mesh.elements[i]
    cnt = apf.countBridgeAdjacent(mesh.m_ptr, el_i, mesh.dim-1, mesh.dim)
    @test  countnz(mesh.pertNeighborEls[i, :])  > cnt
  end

  # check mapping interpolation
  for i=1:mesh.npeers

    dxidx_peer = mesh.dxidx_sharedface[i]
    jac_peer = mesh.jac_sharedface[i]
    for j=1:mesh.peer_face_counts[i]
      iface_i = mesh.shared_interfaces[i][j]
      el_j = iface_i.elementL
      dxidx_el = mesh.dxidx[:, :, 1, el_j]
      jac_el = mesh.jac[:, el_j]
      jac_face = jac_peer[:, j]

      for k=1:mesh.numNodesPerFace
        dxidx_face = dxidx_peer[:, :, k, j]
        for p=1:3
          for n=1:3
            @test isapprox( dxidx_face[n, p], dxidx_el[n, p]) atol=1e-13
          end
        end
        @test isapprox( jac_face[k], jac_el[k]) atol=1e-13
      end
    end   # end loop over peer faces
  end  # end loop over peers

  test_initSendToOwner(mesh)

  mesh_c = PumiMeshDG3(Complex128, sbp, opts, sbpface, topo)
  mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo)
  compare_meshes(mesh, mesh_c)

  # note: this trashes the mesh arrays
  test_parallel_metrics(mesh, sbp, opts)

  # this allocates a ton of memory - don't run
#  test_metrics_rev(mesh, mesh_c, sbp)

  test_adapt_3d(parallel=true)

  # test curvilinear metrics
  opts["use_linear_metrics"] = false
  mesh_c = PumiMeshDG3(Complex128, sbp, opts, sbpface, topo)
  mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo)

  test_metrics_rev_1d(mesh_c, sbp, opts)
  #=
  #increment = mesh.dim*mesh.dim*mesh.numNodesPerElement
  increment=mesh.dim*mesh.coord_numNodesPerElement
  idx = 1
  for i=1:increment:length(mesh.vert_coords)
    println("testing element ", idx)
    test_metrics_rev_1d(mesh_c, sbp, opts, i)
    idx += 1
  end
  =#

  ctor = (T, sbp, opts, sbpface) -> PumiMeshDG3(T, sbp, opts, sbpface, topo)
  test_coords_file(opts, sbp, sbpface, ctor)
end

MPI.Barrier(MPI.COMM_WORLD)
if MPI.Initialized()
  MPI.Finalize()
end


text:
@doc raw"""
    TabularMDP(T, R, discount)

Specify a discrete MDP in tabular form.

- The transition matrix is 3 dimensional of size ``|\mathcal{S}|\times |\mathcal{A}| \times |\mathcal{S}|``, then `T[sj, a, si]` corresponds to the probability of ending in `sj` while taking action `a` in `si`.
- The reward matrix is 2 dimensional of size ``|\mathcal{S}|\times |\mathcal{A}|``, where `R[s, a]` is the reward obtained when taking action `a` in state `s`.
"""
mutable struct TabularMDP <: MDP{Int64, Int64}
    T::Array{Float64, 3} # SPxAxS
    R::Matrix{Float64} # SxA
    discount::Float64
    InitialTable::Array{Float64}
end

@doc raw"""
    TabularPOMDP(T, R, O, discount)

Specify a discrete POMDP in tabular form.

- The transition matrix is 3 dimensional of size ``|\mathcal{S}|\times |\mathcal{A}| \times |\mathcal{S}|``, then `T[sj, a, si]` corresponds to the probability of ending in `sj` while taking action `a` in `si`.
- The observation matrix is also 3 dimensional of size ``|\mathcal{O}| \times |\mathcal{A}| \times |\mathcal{S}|``, `O[o, a, sp]` represents the probability of observing `o` in in state `sp` and action `a`.
- The reward matrix is 2 dimensional of size ``|\mathcal{S}|\times |\mathcal{A}|``, where `R[s, a]` is the reward obtained when taking action `a` in state `s`.
"""
mutable struct TabularPOMDP <: POMDP{Int64, Int64, Int64}
    T::Array{Float64, 3} # SPxAxS
    R::Matrix{Float64} # SxA
    O::Array{Float64, 3} # OxAxSP
    discount::Float64
    InitialTable::Array{Float64}
end

const TabularProblem = Union{TabularMDP, TabularPOMDP}

# Distribution Type and methods
# XXX: this should be replaced with Categorical when https://github.com/JuliaStats/Distributions.jl/issues/743 is fixed
struct DiscreteDistribution{P<:AbstractVector{Float64}}
    p::P
end

support(d::DiscreteDistribution) = 1:length(d.p)

pdf(d::DiscreteDistribution, sp::Int64) = d.p[sp] # T(s', a, s)

rand(rng::AbstractRNG, d::DiscreteDistribution) = sample(rng, Weights(d.p))

# MDP and POMDP common methods
states(p::TabularProblem) = 1:size(p.T, 1)
actions(p::TabularProblem) = 1:size(p.T, 2)

stateindex(::TabularProblem, s::Int64) = s
actionindex(::TabularProblem, a::Int64) = a

discount(p::TabularProblem) = p.discount

transition(p::TabularProblem, s::Int64, a::Int64) = DiscreteDistribution(view(p.T, :, a, s))

reward(prob::TabularProblem, s::Int64, a::Int64) = prob.R[s, a]

initialstate(p::TabularProblem) = DiscreteDistribution(p.InitialTable)

# POMDP only methods
n_observations(p::TabularProblem) = size(p.O, 1)

observations(p::TabularPOMDP) = 1:n_observations(p)

observation(p::TabularPOMDP, a::Int64, sp::Int64) = DiscreteDistribution(view(p.O, :, a, sp))

obsindex(p::TabularPOMDP, o::Int64) = o

# deprecated in POMDPs v0.9
POMDPs.initialstate_distribution(p::TabularProblem) = DiscreteDistribution(p.InitialTable)


text:
#=
function vϕ(model::EoSModel,P,T,z;phase=:unknown)
    v = volume(model,P,T,z,phase = phase)
    vx = volume(model,P,T,z)

    n = sum(z)
    p = pressure(model,v,T,z)
    @show p
    lnZ = log(p*v/(n*R̄*T))
    fun(z) = a_res(model,v,T,z)*n
    ϕ = ForwardDiff.gradient(fun,z)
    return v,ϕ .- lnZ
end

function flash_tpd(xi, lnϕ,x0,lnϕ0)
    f(xiᵢ,lnϕᵢ,x0ᵢ,lnϕ0ᵢ) = xiᵢ*(log(xiᵢ) + lnϕᵢ-log(x0ᵢ)- lnϕ0ᵢ)
    return mapreduce(f,+,xi, lnϕ,x0,lnϕ0)
end

function x0_PT_flash(model::EoSModel,T,P,x0)
    #check each T with T_scale, if treshold is over, replace Pi with inf

 _0 = zero(P+T+first(x0))
 pure = split_model(model)
 crit = crit_pure.(pure)

 T_c = [tup[1] for tup in crit]
 V_c = [tup[3] for tup in crit]
 nan = _0/_0
 sat_nan = (nan,nan,nan)
 replaceP = ifelse.(T_c .< T,true,false)

 sat = [if !replaceP[i] sat_pure(pure[i],T) else sat_nan end for i in 1:length(pure)]

 P_sat = [tup[1] for tup in sat]
 Ki = zeros(eltype(x0),length(x0))
 for i in 1:length(x0)
     if !replaceP[i]
         Ki[i] = P_sat[i][1]/P
     else
         Ki[i] = pressure(pure[i],V_c[i],T)/P
     end
 end

return Ki
end

function PT_flash(model::EoSModel,P,T,z,K0=nothing)

    _0 = zero(P+T+first(z))
    _1 = one(_0)
    if K0 === nothing
        K0 = x0_PT_flash(model,T,P,z)
    end

    β = rr_vle_vapor_fraction(K0,z)
    v0 = volume(model,P,T,z)
    if β <= 0  #bubble point assumption, initial phase supposed
        β = _0
        β0 = _0
        xil = copy(z)
        xiv = rr_flash_vapor(K0,z,_0)

    elseif β >= 0 #dew point assumption
        β = _1
        β0 = _1
        xil = rr_flash_liquid(K0,z,_1)
        xiv = copy(z)
    else #two phase assumption
        xil = rr_flash_liquid(K0,z,β)
        xiv = rr_flash_vapor(K0,z,β)
    end

    logk = similar(xil)
    k = similar(xil)
    v0,logϕ0= vϕ(model,P,T,z)
    vv,logϕv = vϕ(model,P,T,xiv,phase=:v)
    vl,logϕl = vϕ(model,P,T,xil,phase=:l)
    stable_phase = false
    @show K0
    #Rachford Rice sucessive substitution, 5 iters to improve the result a bit
    for i = 1:5
        logk = logϕl - logϕv
        k .= exp.(logk)
        β = rr_vle_vapor_fraction(k,z)
        @show k,β
        if !(_0 <= β <= _1)
            stable_phase = true
            break
        end
        xil = rr_flash_liquid(k,z,β)
        xiv = rr_flash_vapor(k,z,β)
        vv,logϕv = vϕ(model,P,T,xiv,phase=:v)

        vl,logϕl = vϕ(model,P,T,xil,phase=:l)
        if !isfinite(vl+vv)
            stable_phase = true
            break
        end
    end
    #TPD analysis
    if stable_phase == false
        tpdl = flash_tpd(xil, logϕl,z,logϕ0)
        tpdv = flash_tpd(xiv, logϕv,z,logϕ0)
        ΔG = (1 - β) * tpdl + β * tpdv
        if all(>(0),(ΔG,tpdv,tpdl))
            stable_phase = true # more detailed stability phase analisis
        end
    end
    yy = copy(xiv)
    prepend!(yy,log10.([vl,vv]))
    @show yy
    #if TPD analysis worked, pass to optimization.
    #if not, try to check if a phase could be made stable.
    return xil, xiv
end

=#


text:
############################################################
# judiWeights ##############################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: June 2019

# Updated by Ziyi Yin (ziyi.yin@gatech.edu), Nov 2020

export judiWeights, judiWeightsException, subsample

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

# structure for seismic data as an abstract vector
mutable struct judiWeights{vDT<:Number} <: joAbstractLinearOperator{vDT,vDT}
    name::String
    m::Integer
    n::Integer
    nsrc::Integer
    weights
end

mutable struct judiWeightsException <: Exception
    msg :: String
end

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

## outer constructors

"""
    judiWeights
        name::String
        m::Integer
        n::Integer
        nsrc::Integer
        weights
Abstract vector for weighting an extended source, which is injected at every grid point, as weighted by this vector.
Constructors
============
Construct vector cell array of weights. The `weights` keyword\\
can also be a single (non-cell) array, in which case the weights are the same for all source positions:
    judiWeights(weights; nsrc=1)
"""
function judiWeights(weights::Array{T, N}; nsrc=1, vDT::DataType=Float32) where {T, N}
    (T != vDT) && (weights = convert(Array{vDT},weights))
    # length of vector
    n = 1
    m = prod(size(weights))*nsrc
    weightsCell = [deepcopy(weights) for j=1:nsrc]
    return judiWeights{vDT}("Extended source weights",m,n,nsrc,weightsCell)
end

# constructor if weights are passed as a cell array
function judiWeights(weights::Array{Array{T, N},1}; vDT::DataType=Float32) where {T, N}
    nsrc = length(weights)
    (T != vDT) &&  (weights = convert.(Array{vDT, N}, weights))
    # length of vector
    n = 1
    m = prod(size(weights[1]))*nsrc
    return judiWeights{vDT}("Extended source weights",m,n,nsrc,weights)
end

############################################################
## overloaded Base functions

# conj(jo)
conj(a::judiWeights{vDT}) where vDT =
    judiWeights{vDT}("conj("*a.name*")",a.m,a.n,a.nsrc,a.weights)

# transpose(jo)
transpose(a::judiWeights{vDT}) where vDT =
    judiWeights{vDT}(""*a.name*".'",a.n,a.m,a.nsrc,a.weights)

# adjoint(jo)
adjoint(a::judiWeights{vDT}) where vDT =
        judiWeights{vDT}(""*a.name*".'",a.n,a.m,a.nsrc,a.weights)

##########################################################
# *(joLinearFunction, judiWeights)
function *(A::joLinearFunction{ADDT,ARDT},v::judiWeights{avDT}) where {ADDT, ARDT, avDT}
    A.n == size(v,1) || throw(judiWeightsException("Shape mismatch: A:$(size(A)), v: $(size(v))"))
    jo_check_type_match(ADDT,avDT,join(["DDT for *(joLinearFunction,judiWeights):",A.name,typeof(A),avDT]," / "))
    # Evaluate as mat-mat over the weights
    n = size(v.weights[1])
    V = A.fop(vcat([vec(v.weights[i]) for i=1:v.nsrc]...))
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(joLinearFunction,judiWeights):",A.name,typeof(A),eltype(V)]," / "))
    m = length(V)
    return judiWeights{avDT}("Extended source weights", m, 1, v.nsrc, [reshape(V[:, i], n) for i=1:length(v.weights)])
end

# *(joLinearOperator, judiWeights)
function *(A::joAbstractLinearOperator{ADDT,ARDT}, v::judiWeights{avDT}) where {ADDT, ARDT, avDT}
    A.name == "joDirac" && return avDT(1) .* v
    A.name == "(N*joDirac)" && return avDT(A.fop.a) .* v
    A.name == "adjoint(joDirac)" && return avDT(1) .* v
    A.name == "adjoint((N*joDirac))" && return avDT(A.fop.a) .* v
    return mulJ(A, v)
end

function mulJ(A::joAbstractLinearOperator{ADDT,ARDT}, v::judiWeights{avDT}) where {ADDT, ARDT, avDT}
    A.n == size(v,1) || throw(judiWeightsException("Shape mismatch: A:$(size(A)), v: $(size(v))"))
    jo_check_type_match(ADDT,avDT,join(["DDT for *(joLinearFunction,judiWeights):",A.name,typeof(A),avDT]," / "))
    # Evaluate as mat-mat over the weights
    n = size(v.weights[1])
    try
        # Mul may be defined for judiWeights
        V = A.fop(v)
        jo_check_type_match(ARDT,eltype(V),join(["RDT from *(joLinearFunction,judiWeights):",A.name,typeof(A),eltype(V)]," / "))
        return V
    catch e
        V = A.fop(vcat([vec(v.weights[i]) for i=1:v.nsrc]...))
        jo_check_type_match(ARDT,eltype(V),join(["RDT from *(joLinearFunction,judiWeights):",A.name,typeof(A),eltype(V)]," / "))
        return V
    end
end

# *(joCoreBlock, judiWeights)
function *(A::joCoreBlock{ADDT,ARDT}, v::judiWeights{avDT}) where {ADDT, ARDT, avDT}
    A.n == size(v,1) || throw(judiWeightsException("Shape mismatch: A:$(size(A)), v: $(size(v))"))
    jo_check_type_match(ADDT,avDT,join(["DDT for *(joLinearFunction,judiWeights):",A.name,typeof(A),avDT]," / "))
    # Evaluate as mat-mat over the weights
    V = collect(A.fop[i]*v for i=1:length(A.fop))
    [jo_check_type_match(ARDT,eltype(V[i]),join(["RDT from *(joLinearFunction,judiWeights):",A.fop[i].name,typeof(A.fop[i]),eltype(V)]," / ")) for i=1:length(V)]
    return vcat(V...)
end

# vcat
function vcat(ai::Vararg{judiWeights{avDT}, N}) where {avDT, N}
    N == 1 && (return ai[1])
    N > 2 && (return vcat(ai[1], vcat(ai[2:end]...)))
    a, b = ai
    m = a.m + b.m
    n = 1
    nsrc = a.nsrc + b.nsrc
    weights = vcat(a.weights, b.weights)
    return judiWeights{avDT}(a.name,m,n,nsrc,weights)
end

# dot product
function dot(a::judiWeights{avDT}, b::judiWeights{bvDT}) where {avDT, bvDT}
# Dot product for data containers
    size(a) == size(b) || throw(judiWeightsException("dimension mismatch"))
    dotprod = 0f0
    for j=1:a.nsrc
        dotprod += dot(vec(a.weights[j]), vec(b.weights[j]))
    end
    return dotprod
end

# norm
function norm(a::judiWeights{avDT}, p::Real=2) where avDT
    if p == Inf
        return max([maximum(abs.(a.weights[i])) for i=1:a.nsrc]...)
    end
    x = 0.f0
    for j=1:a.nsrc
        x += sum(abs.(vec(a.weights[j])).^p)
    end
    return x^(1.f0/p)
end

#maximum
maximum(a::judiWeights{avDT}) where avDT =   max([maximum(a.weights[i]) for i=1:a.nsrc]...)

#minimum
minimum(a::judiWeights{avDT}) where avDT =   min([minimum(a.weights[i]) for i=1:a.nsrc]...)

# abs
function abs(a::judiWeights{avDT}) where avDT
    b = deepcopy(a)
    for j=1:a.nsrc
        b.weights[j] = abs.(a.weights[j])
    end
    return b
end

# Subsample weights container
"""
    subsample(x,source_numbers)
Subsample seismic weights vectors or matrix-free linear operators and extract the entries that correspond\\
to the shot positions defined by `source_numbers`. Works for inputs of type `judiWeights`, `judiModeling`, \\
`judiProjection`, `judiJacobian`, `Geometry`, `judiRHS`, `judiPDE`, `judiPDEfull`.
Examples
========
(1) Extract 2 shots from `judiWeights` vector:
    dsub = subsample(dobs,[1,2])
(2) Extract geometry for shot location 100:
    geometry_sub = subsample(dobs.geometry,100)
(3) Extract Jacobian for shots 10 and 20:
    Jsub = subsample(J,[10,20])
"""
function subsample(a::judiWeights{avDT},srcnum) where avDT
    srcnum > a.nsrc ? sub = 0 : sub = judiWeights(a.weights[srcnum];vDT=avDT)
    return sub
end

length(x::judiWeights) = x.m

getindex(x::judiWeights, a::Integer) = subsample(x,a)

setindex!(x::judiWeights, y, i) = x.weights[i][:] = y

firstindex(x::judiWeights) = 1

lastindex(x::judiWeights) = x.nsrc

axes(x::judiWeights) = Base.OneTo(x.nsrc)

ndims(x::judiWeights) = length(size(x))

copy(x::judiWeights) = judiWeights{Float32}("Extended source weights", x.m, x.n, x.nsrc, x.weights)

similar(x::judiWeights) = judiWeights{Float32}("Extended source weights", x.m, x.n, x.nsrc, 0f0 * x.weights)

similar(x::judiWeights, element_type::DataType, dims::Union{AbstractUnitRange, Integer}...) = similar(x)

isfinite(x::judiWeights) = all(all(isfinite.(x.weights[i])) for i=1:length(x.weights))

iterate(S::judiWeights, state::Integer=1) = state > length(S.weights) ? nothing : (S.weights[state], state+1)

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

BroadcastStyle(::Type{judiWeights}) = Base.Broadcast.DefaultArrayStyle{1}()

ndims(::Type{judiWeights{Float32}}) = 1

### +/- ####
broadcasted(::typeof(+), x::judiWeights, y::judiWeights) = x + y
broadcasted(::typeof(-), x::judiWeights, y::judiWeights) = x - y

broadcasted(::typeof(+), x::judiWeights, y::Number) = x + y
broadcasted(::typeof(-), x::judiWeights, y::Number) = x - y

broadcasted(::typeof(+), x::Number, y::judiWeights) = x + y
broadcasted(::typeof(-), x::Number, y::judiWeights) = x - y

### * ####
function broadcasted(::typeof(*), x::judiWeights, y::judiWeights)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    z = deepcopy(x)
    for j=1:length(x.weights)
        z.weights[j] = x.weights[j] .* y.weights[j]
    end
    return z
end

function broadcasted(::typeof(*), x::judiWeights, y::Number)
    z = deepcopy(x)
    for j=1:length(x.weights)
        z.weights[j] .*= y
    end
    return z
end

broadcasted(::typeof(*), x::Number, y::judiWeights) = broadcasted(*, y, x)

### / ####
function broadcasted(::typeof(/), x::judiWeights, y::judiWeights)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    z = deepcopy(x)
    for j=1:length(x.weights)
        z.weights[j] = x.weights[j] ./ y.weights[j]
    end
    return z
end

broadcasted(::typeof(/), x::judiWeights, y::Number) = broadcasted(*, x, 1/y)

# Materialize for broadcasting
function materialize!(x::judiWeights, y::judiWeights)
    for j=1:length(x.weights)
        x.weights[j] .= y.weights[j]
    end
    return x
end

function broadcast!(identity, x::judiWeights, y::judiWeights)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    copy!(x,y)
end

function broadcasted(identity, x::judiWeights)
    return x
end

function copy!(x::judiWeights, y::judiWeights)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    for j=1:x.nsrc
        x.weights[j] = y.weights[j]
    end
end

function isapprox(x::judiWeights, y::judiWeights; rtol::Real=sqrt(eps()), atol::Real=0)
    isapprox(x.weights, y.weights; rtol=rtol, atol=atol)
end

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

function A_mul_B!(x::judiWeights, F::Union{joAbstractLinearOperator, joLinearFunction}, y::judiWeights)
    F.m == size(y, 1) ? z = adjoint(F)*y : z = F*y
    for j=1:length(x.weights)
        x.weights[j] .= z.weights[j]
    end
end

function A_mul_B!(x::judiWeights, F::Union{joAbstractLinearOperator, joLinearFunction}, y::Array)
    F.m == size(y, 1) ? z = adjoint(F)*y : z = F*y
    for j=1:length(x.weights)
        x.weights[j] .= z.weights[j]
    end
end

function A_mul_B!(x::Array, F::Union{joAbstractLinearOperator, joLinearFunction}, y::judiWeights)
    F.m == size(y, 1) ? x[:] .= adjoint(F)*y : x[:] .= F*y
end

mul!(x::judiWeights, F::Union{joAbstractLinearOperator, joLinearFunction}, y::judiWeights) = A_mul_B!(x, F, y)
mul!(x::judiWeights, F::Union{joAbstractLinearOperator, joLinearFunction}, y::Array) = A_mul_B!(x, F, y)
mul!(x::Array, F::Union{joAbstractLinearOperator, joLinearFunction}, y::judiWeights) = A_mul_B!(x, F, y)


text:
using LinearAlgebra

# linear combination of constant term, sin(k*π*(x-a)/(b-a)) (for k = 1, 2, ...) and cos(k*π*(x-a)/(b-a)) (for k = 1, 2, ...) terms on the domain (a,b)

struct SinCosSeries{T,S<:Real}
    coeff0::Base.RefValue{T}
    sincoeffs::Vector{T}
    coscoeffs::Vector{T}
    domain::Tuple{S,S}
end
function SinCosSeries(coeff0::T,
                        sincoeffs::Vector{T},
                        coscoeffs::Vector{T},
                        domain::Tuple{Real,Real} = (0,1)) where T
    d = promote(domain...)
    SinCosSeries{T,typeof(d[1])}(Ref(coeff0), sincoeffs, coscoeffs, d)
end

Base.eltype(F::SinCosSeries) = eltype(typeof(F))
Base.eltype(::Type{<:SinCosSeries{T}}) where {T} = eltype(T)

numcosmodes(F::SinCosSeries) = length(F.coscoeffs)
numsinmodes(F::SinCosSeries) = length(F.sincoeffs)

domain(F::SinCosSeries) = F.domain

function Base.similar(F::SinCosSeries, ::Type{T}) where T
    G0 = similar(F[0], T)
    return SinCosSeries(G0, fill(G0, 0), fill(G0, 0))
end

function Base.getindex(F::SinCosSeries, k::Integer)
    if k == 0
        return F.coeff0[]
    elseif isodd(k)
        l = (k+1)>>1
        if l > numsinmodes(F)
            return zero(F.coeff0[])
        else
            return F.sincoeffs[l]
        end
    else
        l = k>>1
        if l > numcosmodes(F)
            return zero(F.coeff0[])
        else
            return F.coscoeffs[l]
        end
    end
end
function Base.setindex!(F::SinCosSeries, v, k::Integer)
    if k == 0
        F.coeff0[] = v
    elseif isodd(k)
        l = (k+1)>>1
        if l > numsinmodes(F)
            setnumsinmodes!(F, l)
        end
        return F.sincoeffs[l] = v
    else
        l = k>>1
        if l > numcosmodes(F)
            setnumcosmodes!(F, l)
        end
        return F.coscoeffs[l] = v
    end
end
function Base.getindex(F::SinCosSeries, l::Integer, s::Symbol)
    if s == :s || s ==:S
        if l > numsinmodes(F)
            return zero(F.coeff0[])
        else
            return F.sincoeffs[l]
        end
    elseif s == :c || s == :C
        if l > numcosmodes(F)
            return zero(F.coeff0[])
        elseif l == 0
            return F.coeff0[]
        else
            return F.coscoeffs[l]
        end
    else
        return BoundsError(F, (k,s))
    end
end
function Base.setindex!(F::SinCosSeries, v, l::Integer, s::Symbol)
    if s == :s || s == :S
        if l > numsinmodes(F)
            setnumsinmodes!(F, l)
        end
        return F.sincoeffs[l] = v
    elseif s == :c || s == :C
        if l > numcosmodes(F)
            setnumcosmodes!(F, l)
        end
        if l == 0
            return F.coeff0[] = v
        else
            return F.coscoeffs[l] = v
        end
    else
        return BoundsError(F, (k,s))
    end
end

eachindex(F::SinCosSeries) =
    Iterators.flatten((0:1:0, 2*(1:numsinmodes(F)) .- 1, 2*(1:numcosmodes(F))))

function (F::SinCosSeries)(x)
    a, b = domain(F)
    a <= x <= b || throw(DomainError())
    θ = π*(x-a)/(b-a)
    f = F[0]
    for k = 1:numsinmodes(F)
        f += F[k,:s]*sin(k*θ)
    end
    for k = 1:numcosmodes(F)
        f += F[k,:c]*cos(k*θ)
    end
    return f
end

function trunc!(F::SinCosSeries;
                    Kmax::Integer = max(numsinmodes(F), numcosmodes(F)), tol::Real = 0)
    Ksintol = findlast(x->norm(x)>=tol, F.sincoeffs)
    Ksinmax = min(Kmax, Ksintol === nothing ? 0 : Ksintol)
    if Ksinmax < numsinmodes(F)
        resize!(F.sincoeffs, Ksinmax)
    end
    Kcostol = findlast(x->norm(x)>=tol, F.coscoeffs)
    Kcosmax = min(Kmax, Kcostol === nothing ? 0 : Kcostol)
    if Kcosmax < numcosmodes(F)
        resize!(F.coscoeffs, Kcosmax)
    end
    return F
end

function setnumsinmodes!(F::SinCosSeries, K::Int)
    while numsinmodes(F) < K
        push!(F.sincoeffs, zero(F[0]))
    end
    if numsinmodes(F) > K
        resize!(F.sincoeffs, K)
    end
    return F
end
function setnumcosmodes!(F::SinCosSeries, K::Int)
    while numcosmodes(F) < K
        push!(F.coscoeffs, zero(F[0]))
    end
    if numcosmodes(F) > K
        resize!(F.coscoeffs, K)
    end
    return F
end

function Base.copyto!(Fdst::SinCosSeries, Fsrc::SinCosSeries)
    domain(Fdst) == domain(Fsrc) || throw(DomainError())
    setnumsinmodes!(Fdst, numsinmodes(Fsrc))
    setnumcosmodes!(Fdst, numcosmodes(Fsrc))
    @inbounds for i in eachindex(Fsrc)
        Fdst[i] = Fsrc[i]
    end
    return Fdst
end

function Base.fill!(F::SinCosSeries{<:Number}, α)
    @inbounds for i in eachindex(F)
        F[i] = α
    end
    return F
end
function LinearAlgebra.rmul!(F::SinCosSeries{<:Number}, α)
    @inbounds for i in eachindex(F)
        F[i] *= α
    end
    return F
end
function LinearAlgebra.lmul!(α, F::SinCosSeries{<:Number})
    @inbounds for i in eachindex(F)
        F[i] *= α
    end
    return F
end
function Base.fill!(F::SinCosSeries{<:AbstractArray}, α)
    @inbounds for i in eachindex(F)
        fill!(F[i], α)
    end
    return F
end
function LinearAlgebra.rmul!(F::SinCosSeries{<:AbstractArray}, α)
    @inbounds for i in eachindex(F)
        rmul!(F[i], α)
    end
    return F
end
function LinearAlgebra.lmul!(α, F::SinCosSeries{<:AbstractArray})
    @inbounds for i in eachindex(F)
        lmul!(α, F[i])
    end
    return F
end

function LinearAlgebra.mul!(Fdst::SinCosSeries, α, Fsrc::SinCosSeries)
    domain(Fdst) == domain(Fsrc) || throw(DomainError())
    setnumsinmodes!(Fdst, numsinmodes(Fsrc))
    setnumcosmodes!(Fdst, numcosmodes(Fsrc))
    @inbounds for i in eachindex(Fsrc)
        Fdst[i] = α * Fsrc[i]
    end
    return Fdst
end
function LinearAlgebra.mul!(Fdst::SinCosSeries, Fsrc::SinCosSeries, α)
    domain(Fdst) == domain(Fsrc) || throw(DomainError())
    setnumsinmodes!(Fdst, numsinmodes(Fsrc))
    setnumcosmodes!(Fdst, numcosmodes(Fsrc))
    @inbounds for i in eachindex(Fsrc)
        Fdst[i] = Fsrc[i] * α
    end
    return Fdst
end

function LinearAlgebra.axpy!(α, Fx::SinCosSeries, Fy::SinCosSeries)
    domain(Fy) == domain(Fx) || throw(DomainError())
    setnumsinmodes!(Fy, max(numsinmodes(Fx), numsinmodes(Fy)))
    setnumcosmodes!(Fy, max(numcosmodes(Fx), numcosmodes(Fy)))
    for i in eachindex(Fx)
        Fy[i] += α*Fx[i]
    end
    return Fy
end
function LinearAlgebra.axpby!(α, Fx::SinCosSeries, β, Fy::SinCosSeries)
    domain(Fy) == domain(Fx) || throw(DomainError())
    for i in numsinmodes(Fx)+1:numsinmodes(Fy)
        Fy[i, :s] *= β
    end
    for i in numcosmodes(Fx)+1:numcosmodes(Fy)
        Fy[i, :c] *= β
    end
    setnumsinmodes!(Fy, max(numsinmodes(Fx), numsinmodes(Fy)))
    setnumcosmodes!(Fy, max(numcosmodes(Fx), numcosmodes(Fy)))
    for i in eachindex(Fx)
        Fy[i] = α*Fx[i] + β*Fy[i]
    end
    return Fy
end

Base.zero(F::SinCosSeries{T}) where T = SinCosSeries(zero(F[0]), T[], T[], domain(F))
function Base.one(F::SinCosSeries)
    G = zero(F)
    G[0] = one(G[0])
    return G
end

function LinearAlgebra.dot(F1::SinCosSeries, F2::SinCosSeries)
    domain(F1) == domain(F2) || throw(DomainError())
    a, b = domain(F1)
    L = (b-a)/1
    Lpi = L/π
    s = L*dot(F1[0], F2[0])
    for k = 1:2:numsinmodes(F1)
        s += (2*Lpi/k)*dot(F1[k,:s], F2[0])
    end
    for l = 1:2:numsinmodes(F2)
        s += (2*Lpi/l)*dot(F1[0], F2[l,:s])
    end
    for k = 1:min(numcosmodes(F1), numcosmodes(F2))
        s += L/2*dot(F1[k,:c], F2[k,:c])
    end
    for k = 1:min(numsinmodes(F1), numsinmodes(F2))
        s += L/2*dot(F1[k,:s], F2[k,:s])
    end
    for k = 1:numcosmodes(F1)
        for l = 1:numsinmodes(F2)
            if isodd(k+l)
                s += Lpi*dot(F1[k,:c], F2[l,:s])*(2*l)/(l^2 - k^2)
            end
        end
    end
    for k = 1:numsinmodes(F1)
        for l = 1:numcosmodes(F2)
            if isodd(k+l)
                s += Lpi*dot(F1[k,:s], F2[l,:c])*(2*k)/(k^2 - l^2)
            end
        end
    end
    return s
end
function LinearAlgebra.norm(F::SinCosSeries)
    a, b = domain(F)
    L = (b-a)/1
    Lpi = L/π
    s = L*norm(F[0])^2
    for k = 1:2:numsinmodes(F)
        s += (4*Lpi/k)*real(dot(F[k,:s], F[0]))
    end
    for k = 1:numcosmodes(F)
        s += L/2*norm(F[k,:c])^2
    end
    for k = 1:numsinmodes(F)
        s += L/2*norm(F[k,:s])^2
    end
    for k = 1:numcosmodes(F)
        for l = 1:numsinmodes(F)
            if isodd(k+l)
                s += Lpi*real(dot(F[k,:c], F[l,:s]))*(4*l)/(l^2 - k^2)
            end
        end
    end
    return sqrt(s)
end

function LinearAlgebra.isapprox(x::SinCosSeries, y::SinCosSeries;
                                atol::Real=0,
                                rtol::Real=Base.rtoldefault(eltype(x), eltype(y), atol))
    return norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))
end

function Base.:+(F1::SinCosSeries, F2::SinCosSeries)
    domain(F1) == domain(F2) || throw(DomainError())
    Ks = min(numsinmodes(F1), numsinmodes(F2))
    sincoeffs = [F1[k,:s]+F2[k,:s] for k = 1:Ks]
    for k = Ks+1:numsinmodes(F1)
        push!(sincoeffs, (+1)*F1[k, :s])
    end
    for k = Ks+1:numsinmodes(F2)
        push!(sincoeffs, (+1)*F2[k, :s])
    end
    Kc = min(numcosmodes(F1), numcosmodes(F2))
    coscoeffs = [F1[k,:c]+F2[k,:c] for k = 1:Kc]
    for k = Kc+1:numcosmodes(F1)
        push!(coscoeffs, (+1)*F1[k, :c])
    end
    for k = Kc+1:numcosmodes(F2)
        push!(coscoeffs, (+1)*F2[k, :c])
    end
    return SinCosSeries(F1[0]+F2[0], sincoeffs, coscoeffs, domain(F1))
end

function Base.:-(F1::SinCosSeries, F2::SinCosSeries)
    domain(F1) == domain(F2) || throw(DomainError())
    Ks = min(numsinmodes(F1), numsinmodes(F2))
    sincoeffs = [F1[k,:s]-F2[k,:s] for k = 1:Ks]
    for k = Ks+1:numsinmodes(F1)
        push!(sincoeffs, (+1)*F1[k, :s])
    end
    for k = Ks+1:numsinmodes(F2)
        push!(sincoeffs, (-1)*F2[k, :s])
    end
    Kc = min(numcosmodes(F1), numcosmodes(F2))
    coscoeffs = [F1[k,:c]-F2[k,:c] for k = 1:Kc]
    for k = Kc+1:numcosmodes(F1)
        push!(coscoeffs, (+1)*F1[k, :c])
    end
    for k = Kc+1:numcosmodes(F2)
        push!(coscoeffs, (-1)*F1[k, :c])
    end
    return SinCosSeries(F1[0]-F2[0], sincoeffs, coscoeffs, domain(F1))
end

Base.:-(F::SinCosSeries) = SinCosSeries(-F[0], -F.sincoeffs, -F.coscoeffs)

Base.:*(F::SinCosSeries, a::Number) = SinCosSeries(F[0]*a, F.sincoeffs*a, F.coscoeffs*a)
Base.:*(a::Number, F::SinCosSeries) = SinCosSeries(a*F[0], a*F.sincoeffs, a*F.coscoeffs)
Base.:/(F::SinCosSeries, a::Number) = SinCosSeries(F[0]/a, F.sincoeffs/a, F.coscoeffs/a)
Base.:\(a::Number, F::SinCosSeries) = SinCosSeries(a\F[0], a\F.sincoeffs, a\F.coscoeffs)

Base.conj(F::SinCosSeries) =
    SinCosSeries(conj(F[0]), conj.(F.sincoeffs), conj.(F.coscoeffs), domain(F))
Base.adjoint(F::SinCosSeries) =
    SinCosSeries(adjoint(F[0]), adjoint.(F.sincoeffs), adjoint.(F.coscoeffs), domain(F))
Base.transpose(F::SinCosSeries) =
    SinCosSeries(transpose(F[0]),
                    transpose.(F.sincoeffs),
                    transpose.(F.coscoeffs),
                    domain(F))

function Base.:*(F1::SinCosSeries, F2::SinCosSeries; Kmax = -1)
    domain(F1) == domain(F2) || throw(DomainError())
    F0 = zero(F1[0]) * zero(F2[0])
    Ks = max(numsinmodes(F1)+numcosmodes(F2), numcosmodes(F1)+numsinmodes(F2))
    Kc = max(numsinmodes(F1)+numsinmodes(F2), numcosmodes(F1)+numcosmodes(F2))
    if Kmax >= 0
        Ks = min(Ks, Kmax)
        Kc = min(Kc, Kmax)
    end
    sincoeffs = [copy(F0) for k = 1:Ks]
    coscoeffs = [copy(F0) for k = 1:Kc]
    F = SinCosSeries(F0, sincoeffs, coscoeffs, domain(F1))
    return mul!(F, F1, F2, true, true)
end
function LinearAlgebra.mul!(F::SinCosSeries, F1::SinCosSeries, F2::SinCosSeries,
                                α = true, β = false)
    domain(F) == domain(F1) == domain(F2) || throw(DomainError())
    if β != true
        rmul!(F, β)
    end
    for k in eachindex(F1)
        f1k = F1[k]
        for l in eachindex(F2)
            f2l = F2[l]
            for (p,s) in SimpsonIterator(k,l)
                if (iseven(p) && p <= 2*numcosmodes(F)) ||
                    (isodd(p) && p <= 2*numsinmodes(F)-1)
                    fp = F[p]
                    if fp isa AbstractArray
                        if f1k isa AbstractArray && f2l isa AbstractArray
                                mul!(fp, f1k, f2l, s * α, true)
                        elseif f2l isa Number
                            axpy!((s * α) * f2l, f1k, fp)
                        elseif f1k isa Number
                            axpy!((s * α) * f1k, f2l, fp)
                        else
                            fp .+= (f1k * f2l) .* (s * α)
                        end
                    else
                        F[p] += (f1k * f2l) * (s * α)
                    end
                end
            end
        end
    end
    return F
end

function differentiate(F::SinCosSeries)
    a, b = domain(F)
    πL = π/(b-a)
    sincoeffs = [-πL*k*F[k,:c] for k in 1:numcosmodes(F)]
    coscoeffs = [+πL*k*F[k,:s] for k in 1:numsinmodes(F)]
    return SinCosSeries(zero(F[0]), sincoeffs, coscoeffs, domain(F))
end
function integrate(F::SinCosSeries, interval::Tuple{Real,Real} = domain(F))
    @assert interval == domain(F)
    a, b = interval
    Lπ = (b-a)/π
    s = ((b-a)/1)*F[0]
    for k = 1:2:numsinmodes(F)
        s += (2*Lπ/k)*F[k,:s]
    end
    return s
end

# Note that k == 0 => constant, iseven(k) => cos((k/2)*pi*x/L, isodd(l) => sin((l+1)/2*pi*x/L)
struct SimpsonIterator
    k::Int
    l::Int
end
function Base.iterate(it::SimpsonIterator, state = 0)
    k, l = it.k, it.l
    if k == 0
        if state == 0
            return l=>1//1, 1
        else
            return nothing
        end
    elseif l == 0
        if state == 0
            return k=>1//1, 1
        else
            return nothing
        end
    elseif isodd(k) && isodd(l) # sin * sin
        if state == 0
            return (k+l+2)=>-1//2, 1
        elseif state == 1
            return abs(k-l)=>1//2, 2
        else
            return nothing
        end
    elseif iseven(k) && iseven(l) # cos * cos
        if state == 0
            return k+l=>1//2, 1
        elseif state == 1
            return abs(k-l)=>1//2, 2
        else
            return nothing
        end
    elseif iseven(k) && isodd(l) # cos * sin
        if state == 0
            return k+l=>1//2, 1
        elseif state == 1 && k != l+1
            return (k < l ? (l-k)=>1//2 : (k-l-2)=>-1//2), 2
        else
            return nothing
        end
    else # sin * cos
        if state == 0
            return k+l=>1//2, 1
        elseif state == 1 && k+1 != l
            return (k+1 > l ? (k-l)=>1//2 : (l-k-2)=>-1//2), 2
        else
            return nothing
        end
    end
end


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ bb6cae46-4108-11eb-2309-3d11f13a19fa
begin
	using Pkg
	Pkg.activate(mktempdir())
end

# ╔═╡ d8b334ae-4108-11eb-2546-99e75efc18d6
begin
	Pkg.add(["Plots", "PlutoUI","StatsBase", "DataStructures", "Unitful", "PhysicalConstants", "UnitfulEquivalences", "DataFrames", "CSV", "Dates", "TimeZones", "StatsPlots", "ImageFiltering", "OffsetArrays"])

	using Plots
	plotly()
	using PlutoUI
	using StatsBase
	using DataStructures
	using Unitful
	using PhysicalConstants.CODATA2018: c_0, h, ħ
	using UnitfulEquivalences
	using DataFrames
	using CSV
	using TimeZones
	using Dates
	using StatsPlots
	using ImageFiltering
	using OffsetArrays
end

# ╔═╡ e6b122ce-5548-11eb-3980-53dfd8e5897c
begin
	frmt = dateformat"yyyy-mm-dd HH:MM:SS"
	column_types = Dict(:ResultTime=>DateTime, :AnalysisInsertTime=>DateTime)
	df = CSV.read("/home/kris/Downloads/opendata_covid19_test_results.csv", DataFrame, types=column_types, dateformat=frmt)
	dropmissing!(df)
end

# ╔═╡ eaab22ca-5b2c-11eb-1eb6-cf73054ee98e
begin
	inimarv_raw = CSV.read("/home/kris/Downloads/RV0222U_20012021163956877.csv", DataFrame)
	inimarv = Dict(eachrow(inimarv_raw[!, [:Maakond, :Value]]));
end

# ╔═╡ 6fcdea12-35f3-4648-8ba0-cfd6e89547ef
counties = collect(setdiff(keys(inimarv), ["..Tallinn", "..Tartu linn", "Maakond teadmata"]))

# ╔═╡ 2a795f4a-6bca-11eb-2b3f-b385b41d4759
@bind selectedCounty MultiSelect(counties)

# ╔═╡ ed6d5c68-5b1f-11eb-134b-57fcc502469c
selected_df = filter(x -> x[:County] ∈ selectedCounty, df);

# ╔═╡ 2c31d1ca-5b20-11eb-1f1b-db947a36b2c6
selected_df_p = filter(x -> x[:ResultValue] == "P", selected_df);

# ╔═╡ fc91c29a-6bca-11eb-1cd7-132e515bcc53
all = combine(groupby(selected_df_p, [:StatisticsDate, :County]), nrow => :count);

# ╔═╡ 7b2f2674-5b22-11eb-225c-cbd01d9b1058
begin
	smooth_range = -7:0
	n = size(smooth_range)[1]
	kernel = OffsetArray(fill(1/n, n), smooth_range)
end;

# ╔═╡ baf09934-5b20-11eb-3283-d9b7bdfb5a90
allSelected = transform(combine(groupby(selected_df_p, [:StatisticsDate]), nrow => :count), (arr -> imfilter(arr[!, :count], kernel)));

# ╔═╡ dc1fdc56-6bcd-11eb-32c7-f3515450573c
dates = allSelected[!, :StatisticsDate];

# ╔═╡ f9b7501b-abb7-4998-bf5a-787cba8355d6
@df allSelected plot(:StatisticsDate, [:count, :x1], ticks=:native, legend=:topleft)

# ╔═╡ 74295a2c-6bc9-11eb-1474-9befdff2e2bf
next = transform( groupby(all, [:County]), (arr -> imfilter(arr[!, :count], kernel) / inimarv[arr[1,:County]] * 10000))

# ╔═╡ 080ecf02-5549-11eb-06cd-b10c7ed4e1e4
groups = groupby(df, [:ResultValue]);

# ╔═╡ 0f8d9246-6bd2-11eb-3c89-7dcd26ee5f0e
tartu_p = filter(x -> x[:County] ∈ selectedCounty && x[:ResultValue] == "P" , df);

# ╔═╡ 742c4eaa-7fe5-11eb-37e5-b96d6d7d86cd
@df next plot(:StatisticsDate, [:x1], group=:County, title="New infections per 10000", legend=:topleft)

# ╔═╡ d49f6f5e-5cb5-11eb-14e1-45d31ab53e89
function week_multi(week_p)
	week_multiplier = zeros(size(week_p)[1])
	for i = 2:size(week_p)[1]
		prev = week_p[i-1, :count_in_week]
		current = week_p[i, :count_in_week]
		week_multiplier[i] = current/prev
	end
	week_multiplier
end;

# ╔═╡ d212b59e-810d-11eb-3b97-5dd770b3c949
function multi_smooth(next) # Not used, but it's a nice function and I spent time on it, so it should stay for now
	weeks = reshape(next[1:size(next, 1)÷7*7], 7, :)
	sums = sum(weeks, dims=1)
	res = map(zip(sums[:], sums[2:end])) do (prev, cur)
		cur/prev
	end
	prepend!(res, 1)
	append!(res, 1)
	repeat(res, inner=7)[1:size(next,1)]
end;

# ╔═╡ ef39fc56-554b-11eb-3082-77349ed217d4
counts = combine(groups, nrow => :count)

# ╔═╡ c85e156e-554a-11eb-218b-df5eabe2d8ab
@df counts bar(:ResultValue, :count)

# ╔═╡ 5123d73c-6c9b-11eb-0231-675c2b98d0d3
yearweek(date) = (year(date) - 2020)*52 + dayofyear(date) ÷ 7

# ╔═╡ 15f92d78-5b26-11eb-0050-27f2434d1f75
week_p = combine(groupby(transform(tartu_p, :StatisticsDate => ByRow(x -> yearweek(x)) => :week), [:week, :County]), nrow => :count_in_week);

# ╔═╡ dcba4866-5b26-11eb-3e27-8fc9dd4bae77
@df week_p plot(:week, [:count_in_week], group=:County, title="New infections in week")

# ╔═╡ 1889b0c6-6c97-11eb-1577-3ba3e8940484
transform!(groupby(week_p, [:County]), week_multi);

# ╔═╡ d7a78b2e-5cb7-11eb-1967-1b3dc0a388bd
begin
	@df week_p plot(:week, [:x1], group=:County, title="New infections multiplier in county", )
	hline!([1])
end

# ╔═╡ 446f8889-a549-4fd2-b3e4-c09af9961405
no_county = transform(combine(groupby(week_p, [:week]), :count_in_week => sum => :count_in_week), week_multi);

# ╔═╡ 65004f86-3c67-43a6-bb01-16df98e58148
begin
	begin
		@df no_county plot(:week, [:x1], title="Combined new infections multiplier")
		hline!([1])
	end
end

# ╔═╡ 6eccfc72-cf1f-4259-942d-5898c83f4df7
no_county

# ╔═╡ 15499246-5cb6-11eb-2683-d34604bf5914
size(no_county)

# ╔═╡ d8ed8d5c-6c9b-11eb-2215-2fb6b80d30e2
countmap(week.(df[!, :StatisticsDate]))

# ╔═╡ d5b5cf78-5a87-11eb-1d85-176632857548
function mylinspace(d1::DateTime, d2::DateTime, n::Int)
    Δ = d2 - d1
    T = typeof(Δ)
    δ = T(round(Int, Dates.value(Δ)/(n - 1)))
    d2 = d1 + δ*(n - 1)
    return d1:δ:d2
end

# ╔═╡ 7f5be60a-5b23-11eb-3d24-d17513534e16
ticks = mylinspace(DateTime(2020,03,01), DateTime(2021,03,06), 30)

# ╔═╡ 620a967e-5b2f-11eb-3dc5-6108d315124b
(..)(x::AbstractDict,i...) = getindex.(Ref(x), i...)

# ╔═╡ 9e3c3714-5b2d-11eb-2a50-413ab7ed65d1
df[!, :population] = inimarv.. df[!, :County];

# ╔═╡ dd2bc6ca-5a87-11eb-1cfd-9bbf310b48fa
md"""
# Functions in physics
1. In the figure below you see a graph of some function. It is a square function.
"""

# ╔═╡ 9b1d7a9a-5a8b-11eb-3721-510c5ee70087
md"""
### 2. Physical example

This function could be the function of velocity with regards to time with a constant _jerk_.

$$v(t) = \frac{jt^2}{2} + a_0t + v_0$$
$$1 \frac{m}{s} = \frac{1 \frac{m}{s^3} \cdot (1 s)^2}{2} + 1 \frac{m}{s^2} \cdot 1s + 1 \frac{m}{s}$$

* _$v$_ velocity - it's the first deriative of position, showing the change of change of postion.
* __$v_0$__ starting velocity
* __$a_0$__ starting acceleration. acceleration is the first deriative of velocity and the second deriative of position.
* __$j$__ jerk - the third deriative of position. constantly changing velocity.
* __$t$__ time
"""

# ╔═╡ 15a2336e-5a91-11eb-3f7b-a3e55c5a3be4
md"""
### 3. Integral

$$v(t) = \frac{jt^2}{2} + a_0t + v_0$$
$$ds = v\ dt =  \left( \frac{jt^2}{2} + a_0t + v_0 \right) dt$$
$$\int_{s_0}^s ds = \int_0^t \left( \frac{jt^2}{2} + a_0t + v_0 \right) dt$$
$$\int_{s_0}^s ds = \int_0^t \left( \frac{jt^2}{2} + a_0t + v_0 \right) dt$$
$$s-s_0 = \frac{jt^3}{6} + \frac{a_0t^2}{2} + v_0t $$
$$s = \frac{jt^3}{6} + \frac{a_0t^2}{2} + v_0t + s_0$$

I chose to derive this integral with the limits of $(s_0, s)$ to get the positional displacement in familiar terms. (where $s_0$ is the starting positional displacement and $s$ is the final displacement)

The integral of the constant jerk velocity function with regards to time is the positional displacement physical quantity.
"""

# ╔═╡ 935814f4-5a91-11eb-2dc7-a9db013f9d95
md"""

### 4. Property of Nature

The square function describes the concept of compound growth. How even a linearly changing variable, when summed up over time, changes polynominally. This can happen when measuring speed and the jerk increases linearly.

"""

# ╔═╡ 02b9ba5c-5a8b-11eb-39b1-0d23e4c12f0b
square(a,x,c) = a*x^2 + x

# ╔═╡ 4d0b47b8-5a88-11eb-094a-1d816b3d3a9b
x = 0:1:50

# ╔═╡ d78890c6-5a88-11eb-3601-13e50420b603
@bind j_ Slider(0:10, show_value=true)

# ╔═╡ e7510d12-5a88-11eb-0536-8de677458800
@bind a_0 Slider(-50:50, show_value=true)

# ╔═╡ ee2ef6ca-5a89-11eb-0b85-0d59759a8bac
@bind s_0 Slider(-500:500, show_value=true)

# ╔═╡ de09158c-5a89-11eb-23d3-adc6acc42782
@bind v_0 Slider(-200:200, show_value=true)

# ╔═╡ 5b8a2d9a-5a88-11eb-1193-998aca8aa0a0
jitter(j, a_0, v_0, s_0, t) = j*t^3 / 6 + a_0*t^2/2 + v_0*t + s_0

# ╔═╡ b0a9669c-5a88-11eb-34ee-cd73f0892f7f
jitter(j,t) = jitter(j, 0, 0, 0, t)

# ╔═╡ b2ddbdf0-5a88-11eb-3dba-03b51da9048e
s = jitter.(j_,a_0, v_0, s_0, x)

# ╔═╡ 84c21a6c-4109-11eb-166e-151c460d28c8
md"""
## E2.1.1

In region AB car covered the distance of 25 m, in region BC the distance of 15 m and in CD the distance of 20 m. The acceleration of the car in the region AB was 2.5 m/s2, in the region BC 0 m/s2 and in the region CD -5 m/s2. During the whole trip the car was affected by the drag force of 2 kN. Motor's traction force in region AB was 7 kN
; in region BC 8) Answer
2 kN

. In region CD driver hit the brakes. Braking affected the car with the force of 9) Answer
10 kN

. The work done by motor in the region AB was 10) Answer
14 kJ

, in the region BC 11) Answer
2 kJ

. In the end of the region AB motor reached the power of 12) Answer
7 kW

 and in the region BC the power of 13) Answer
2 kW

.
"""

# ╔═╡ 50b8ef74-5591-11eb-26bf-9feb4bba0ca3
v_A = 10u"m/s"

# ╔═╡ 7b147c4e-559f-11eb-37ce-075bf3998087
v_B = 15u"m/s"

# ╔═╡ f690e870-5591-11eb-0a58-6faf52de3fec
t_AB = 1u"s"

# ╔═╡ 57c9516e-5591-11eb-226e-f163d5df253d
v_C = 15u"m/s"

# ╔═╡ 070ea700-5592-11eb-1334-abb3cb13262a
t_BC = 2u"s"

# ╔═╡ 5def440c-5591-11eb-05b9-0dfd6d2952c8
v_D = 5u"m/s"

# ╔═╡ 18c7119e-5592-11eb-1afb-e71d7e7b9ca8
t_CD = 2u"s"

# ╔═╡ b29c4136-55a1-11eb-0f19-85789ee1a471

# ╔═╡ 7f1987a6-55a1-11eb-2e64-631c13345a8f

# ╔═╡ ed11f462-55a0-11eb-364b-6b9fd7282d48

# ╔═╡ 3d797c44-55a0-11eb-2490-d98ffcf5f006
drag = 6u"kN"

# ╔═╡ 61e0edf6-55a0-11eb-0d43-69315ade5c2e
mass = 4000u"kg"

# ╔═╡ 71d73800-55a0-11eb-28aa-895ff4f0a36b
force(a) = uconvert(u"kN", mass*a + drag)

# ╔═╡ a76941d0-559f-11eb-2412-bbdefbbfeb04
distance(v1, v2, t) = (v1 + v2) / 2 * t

# ╔═╡ 6edc69c2-5591-11eb-0b05-b35d24abe8aa
s_AB = distance(v_A, v_B, t_AB)

# ╔═╡ c135bbca-559f-11eb-21d0-81f51845afcb
s_BC = distance(v_B, v_C, t_BC)

# ╔═╡ ec9ee372-559f-11eb-3e88-3f4a9cd4c0ec
s_CD = distance(v_C, v_D, t_CD)

# ╔═╡ fe32a588-559f-11eb-25c8-25c3fdb3fb3f
acc(v1, v2, t) = (v2-v1) / t

# ╔═╡ 3ad405cc-55a0-11eb-169a-43cfdaa2d3af
a_AB = acc(v_A, v_B, t_AB)

# ╔═╡ 3be46f9c-55a0-11eb-1e67-f9c641ab7eaf
f_AB = force(a_AB)

# ╔═╡ 42b15790-55a0-11eb-0a84-a79553b87993
a_BC = acc(v_B, v_C, t_BC)

# ╔═╡ cc0c493c-55a0-11eb-37ee-dff68db5868b
f_BC = force(a_BC)

# ╔═╡ 1f54ed0c-55a0-11eb-14bd-01ed4f5743e7
a_CD = acc(v_C, v_D, t_CD)

# ╔═╡ ecfdecf4-55a0-11eb-0427-d9eacdea896c
f_CD = force(a_CD)

# ╔═╡ 045ab740-55a1-11eb-0751-a59eb196863e
work(f, s) = uconvert(u"kJ", f*s)

# ╔═╡ c98e09ac-55a0-11eb-29d3-b12c978f1e57
w_AB = work(f_AB, s_AB)

# ╔═╡ cc6a6710-55a0-11eb-1da6-45a6bb3bb173
w_BC = work(f_BC, s_BC)

# ╔═╡ 962f7d2e-55a1-11eb-3299-215a5c9504c5
power(w, t) = uconvert(u"kW", w/t)

# ╔═╡ 3198fbb0-55a1-11eb-016e-712a7d8f07f1
p_AB = power(w_AB, t_AB)

# ╔═╡ 74fa7748-410c-11eb-0688-a37608618c96
md"""
## E2

"""

# ╔═╡ fdae38fe-410c-11eb-18dd-9bae51fb4e00
t2 = 14u"s"

# ╔═╡ fdbe0894-410c-11eb-0ace-3d269b2942aa
grav = 10u"m/s^2"

# ╔═╡ 220ec448-410d-11eb-17c3-d1682a1588c9
v_2 = grav *t2

# ╔═╡ 4b4e5710-55a3-11eb-2f5d-3df5543e480c
x2 = 0u"s":0.5u"s":30u"s"

# ╔═╡ 5a8d3296-55a3-11eb-3d1c-65245f677262
vs2 = map(x -> v_2 - x*grav, x2)

# ╔═╡ 75ef05aa-55a3-11eb-1c89-fd0fb600608f
plot(ustrip(x2),ustrip(vs2))

# ╔═╡ 22b3b39e-55a4-11eb-3d8e-a18fe4e18872
hs2 = cumsum(vs2 .* 0.5u"s")

# ╔═╡ 6e1b79c0-55a4-11eb-2d3f-296762a8604f
plot(ustrip(x2),ustrip(hs2))

# ╔═╡ 7c3d0eb0-55a4-11eb-3ed0-81c97785c7f0
maximum(ustrip(hs2))

# ╔═╡ 7f1fb2fe-55a4-11eb-1b91-3bd39c4872fe

# ╔═╡ 3bb15c6c-410d-11eb-1b1b-4b13708ceb70
md"""
## E3

"""

# ╔═╡ 0f68fa64-55a5-11eb-24ed-fd2b3fc9ed68
v_e = v_0 + ((a_0 + a_1) / 2)*t3

# ╔═╡ 528ef22e-55a5-11eb-2b03-3d92a4b9c3c6
j = 0.5u"m/s^3"

# ╔═╡ 3de00658-55a5-11eb-102e-3103b223001d
s_e = j*t3^3 / 6 + a_0*t3^2/2 + v_0*t3

# ╔═╡ c094cb3a-4117-11eb-168a-df762ad5739e
md"""
## E4

"""

# ╔═╡ a837ef5c-55a5-11eb-06f5-cb27048bc9a4
t4 = 3u"s"

# ╔═╡ c5cecd9e-55a5-11eb-343a-a3c29893be1f
d4 = 0.4u"s^-1"

# ╔═╡ cff2f4c4-55a5-11eb-32c2-1db64f8ea677
one = d4^-1

# ╔═╡ 74480530-55a6-11eb-0a7b-8b30ed458234
two = d4 * t4

# ╔═╡ b0f10e8e-55a6-11eb-34a1-755dceb63e08
thr = 1.44 * one

# ╔═╡ ef0bdbcc-55a6-11eb-114f-7d43ab25736a

# ╔═╡ 7fec90b4-411b-11eb-04db-a913940a9f22
md"""
## E5

The time counting is started when a mass-spring oscillator passes the equilibrium state. The oscillator stops at the distance of 5 cm from equilibrium state. The system makes one oscillation during time interval of 1.57 seconds. The  amplitude of the oscillator is 1)

5.00 cm

, oscillation period is 2)

1.57 s

 and angular frequency is 3)

4.00 rad/s

. The maximal value of acceleration is 4)

80.00 cm/s2

. Assume that at the start of clock, the deviation is increasing in the positive direction. In this case, the pendulum will reach its maximum negative acceleration 5)

0.39

 seconds after the clock started. At this moment, pendulum will pass the phase angle of 6)

1.57

 radians. The pendulum will reach its maximum positive acceleration 7)

1.18

 seconds after the clock started. At this moment, pendulum will pass the phase angle of 8)

4.71

 radians. The maximal value of velocity is 9)

20.00 cm/s

. The oscillator will reach the maximum negative velocity 10)

0.78

 seconds after the clock started. At this moment, the oscillator will pass the phase angle of 11)

3.14

 radians.

[NB: give answers with precision of 2 characters after dot]
[NB: case insensitive!]

"""

# ╔═╡ 46652b7e-4122-11eb-1576-29e34caf1b03
d = 4u"cm"

# ╔═╡ 78e0e354-4122-11eb-1a11-0b1d36a38b70
T = 3.14u"s"

# ╔═╡ 46798970-4122-11eb-31a8-ab967162b1ae
ω = 2*pi / T

# ╔═╡ 451edba2-4122-11eb-0ea9-b3e5a7021bd1
md"""
# E5

The time counting is started when a mass-spring oscillator passes the equilibrium state. The oscillator stops at the distance of 5 cm from equilibrium state. The system makes one oscillation during time interval of 1.57 seconds. The  amplitude of the oscillator is 1)

$d

, oscillation period is 2)

$T

 and angular frequency is 3)

$ω rad/s

. The maximal value of acceleration is 4)

80.00 cm/s2

. Assume that at the start of clock, the deviation is increasing in the positive direction. In this case, the pendulum will reach its maximum negative acceleration 5)

0.39

 seconds after the clock started. At this moment, pendulum will pass the phase angle of 6)

1.57

 radians. The pendulum will reach its maximum positive acceleration 7)

1.18

 seconds after the clock started. At this moment, pendulum will pass the phase angle of 8)

4.71

 radians. The maximal value of velocity is 9)

20.00 cm/s

. The oscillator will reach the maximum negative velocity 10)

0.78

 seconds after the clock started. At this moment, the oscillator will pass the phase angle of 11)

3.14

 radians.

"""

# ╔═╡ 20c55224-55a7-11eb-2490-010a8c63cd4b
v5 = ω*d/2

# ╔═╡ aa066c66-4126-11eb-1909-713aeda74edc
speed6 = 314u"m/s"

# ╔═╡ c4375a50-4126-11eb-2667-7b693682e9a8
freq6 = 125u"Hz"

# ╔═╡ aa323b5c-4126-11eb-18fa-a5aa100d0678
period6 = uconvert(u"ms", 1/freq6)

# ╔═╡ aa487444-4126-11eb-306c-1d98f391cfc6
wavelength6 = uconvert(u"cm", speed6*period6)

# ╔═╡ aa708326-4126-11eb-2129-df2f5666aa93
wave_n = uconvert(u"m^-1", 2*pi/wavelength6)

# ╔═╡ c98691d4-4125-11eb-3e71-d1ffa4fc0f3a
md"""
## E6

Sound plane wave propagates in the cold air with the speed of 314 m/s and frequency 500 Hz. Sound wave oscillation period is 1)
 $period6
 ms, wavelength 2)

$wavelength6

 cm and wave number is 3)

$wave_n

 m-1.  At the start of time counting, the deviation of the sound wave source was minimal (-A). After 3.5 ms from the start and at the distance of 15.7 cm from the wave source the deviation is 4)

maximal

(maximal or A, minimal or –A or zero)."

"""

# ╔═╡ Cell order:
# ╠═bb6cae46-4108-11eb-2309-3d11f13a19fa
# ╠═d8b334ae-4108-11eb-2546-99e75efc18d6
# ╠═e6b122ce-5548-11eb-3980-53dfd8e5897c
# ╠═eaab22ca-5b2c-11eb-1eb6-cf73054ee98e
# ╠═9e3c3714-5b2d-11eb-2a50-413ab7ed65d1
# ╠═6fcdea12-35f3-4648-8ba0-cfd6e89547ef
# ╠═2a795f4a-6bca-11eb-2b3f-b385b41d4759
# ╠═ed6d5c68-5b1f-11eb-134b-57fcc502469c
# ╠═2c31d1ca-5b20-11eb-1f1b-db947a36b2c6
# ╠═baf09934-5b20-11eb-3283-d9b7bdfb5a90
# ╠═dc1fdc56-6bcd-11eb-32c7-f3515450573c
# ╠═f9b7501b-abb7-4998-bf5a-787cba8355d6
# ╠═fc91c29a-6bca-11eb-1cd7-132e515bcc53
# ╠═74295a2c-6bc9-11eb-1474-9befdff2e2bf
# ╠═7b2f2674-5b22-11eb-225c-cbd01d9b1058
# ╠═080ecf02-5549-11eb-06cd-b10c7ed4e1e4
# ╠═0f8d9246-6bd2-11eb-3c89-7dcd26ee5f0e
# ╠═742c4eaa-7fe5-11eb-37e5-b96d6d7d86cd
# ╠═15f92d78-5b26-11eb-0050-27f2434d1f75
# ╠═dcba4866-5b26-11eb-3e27-8fc9dd4bae77
# ╠═d49f6f5e-5cb5-11eb-14e1-45d31ab53e89
# ╠═d212b59e-810d-11eb-3b97-5dd770b3c949
# ╠═1889b0c6-6c97-11eb-1577-3ba3e8940484
# ╠═d7a78b2e-5cb7-11eb-1967-1b3dc0a388bd
# ╠═65004f86-3c67-43a6-bb01-16df98e58148
# ╠═446f8889-a549-4fd2-b3e4-c09af9961405
# ╠═6eccfc72-cf1f-4259-942d-5898c83f4df7
# ╠═15499246-5cb6-11eb-2683-d34604bf5914
# ╠═ef39fc56-554b-11eb-3082-77349ed217d4
# ╠═c85e156e-554a-11eb-218b-df5eabe2d8ab
# ╠═7f5be60a-5b23-11eb-3d24-d17513534e16
# ╠═5123d73c-6c9b-11eb-0231-675c2b98d0d3
# ╠═d8ed8d5c-6c9b-11eb-2215-2fb6b80d30e2
# ╠═d5b5cf78-5a87-11eb-1d85-176632857548
# ╠═620a967e-5b2f-11eb-3dc5-6108d315124b
# ╟─dd2bc6ca-5a87-11eb-1cfd-9bbf310b48fa
# ╟─9b1d7a9a-5a8b-11eb-3721-510c5ee70087
# ╠═15a2336e-5a91-11eb-3f7b-a3e55c5a3be4
# ╠═935814f4-5a91-11eb-2dc7-a9db013f9d95
# ╠═02b9ba5c-5a8b-11eb-39b1-0d23e4c12f0b
# ╠═4d0b47b8-5a88-11eb-094a-1d816b3d3a9b
# ╠═d78890c6-5a88-11eb-3601-13e50420b603
# ╠═e7510d12-5a88-11eb-0536-8de677458800
# ╠═ee2ef6ca-5a89-11eb-0b85-0d59759a8bac
# ╠═de09158c-5a89-11eb-23d3-adc6acc42782
# ╟─b2ddbdf0-5a88-11eb-3dba-03b51da9048e
# ╠═5b8a2d9a-5a88-11eb-1193-998aca8aa0a0
# ╠═b0a9669c-5a88-11eb-34ee-cd73f0892f7f
# ╠═84c21a6c-4109-11eb-166e-151c460d28c8
# ╠═50b8ef74-5591-11eb-26bf-9feb4bba0ca3
# ╠═7b147c4e-559f-11eb-37ce-075bf3998087
# ╠═f690e870-5591-11eb-0a58-6faf52de3fec
# ╠═57c9516e-5591-11eb-226e-f163d5df253d
# ╠═070ea700-5592-11eb-1334-abb3cb13262a
# ╠═5def440c-5591-11eb-05b9-0dfd6d2952c8
# ╠═18c7119e-5592-11eb-1afb-e71d7e7b9ca8
# ╠═6edc69c2-5591-11eb-0b05-b35d24abe8aa
# ╠═3ad405cc-55a0-11eb-169a-43cfdaa2d3af
# ╠═3be46f9c-55a0-11eb-1e67-f9c641ab7eaf
# ╠═c98e09ac-55a0-11eb-29d3-b12c978f1e57
# ╠═3198fbb0-55a1-11eb-016e-712a7d8f07f1
# ╠═b29c4136-55a1-11eb-0f19-85789ee1a471
# ╠═c135bbca-559f-11eb-21d0-81f51845afcb
# ╠═42b15790-55a0-11eb-0a84-a79553b87993
# ╠═cc0c493c-55a0-11eb-37ee-dff68db5868b
# ╠═cc6a6710-55a0-11eb-1da6-45a6bb3bb173
# ╠═7f1987a6-55a1-11eb-2e64-631c13345a8f
# ╠═ec9ee372-559f-11eb-3e88-3f4a9cd4c0ec
# ╠═1f54ed0c-55a0-11eb-14bd-01ed4f5743e7
# ╠═ecfdecf4-55a0-11eb-0427-d9eacdea896c
# ╠═ed11f462-55a0-11eb-364b-6b9fd7282d48
# ╠═3d797c44-55a0-11eb-2490-d98ffcf5f006
# ╠═61e0edf6-55a0-11eb-0d43-69315ade5c2e
# ╠═71d73800-55a0-11eb-28aa-895ff4f0a36b
# ╠═a76941d0-559f-11eb-2412-bbdefbbfeb04
# ╠═fe32a588-559f-11eb-25c8-25c3fdb3fb3f
# ╠═045ab740-55a1-11eb-0751-a59eb196863e
# ╠═962f7d2e-55a1-11eb-3299-215a5c9504c5
# ╠═74fa7748-410c-11eb-0688-a37608618c96
# ╠═fdae38fe-410c-11eb-18dd-9bae51fb4e00
# ╠═fdbe0894-410c-11eb-0ace-3d269b2942aa
# ╠═220ec448-410d-11eb-17c3-d1682a1588c9
# ╠═4b4e5710-55a3-11eb-2f5d-3df5543e480c
# ╠═5a8d3296-55a3-11eb-3d1c-65245f677262
# ╠═75ef05aa-55a3-11eb-1c89-fd0fb600608f
# ╠═22b3b39e-55a4-11eb-3d8e-a18fe4e18872
# ╠═6e1b79c0-55a4-11eb-2d3f-296762a8604f
# ╠═7c3d0eb0-55a4-11eb-3ed0-81c97785c7f0
# ╠═7f1fb2fe-55a4-11eb-1b91-3bd39c4872fe
# ╠═3bb15c6c-410d-11eb-1b1b-4b13708ceb70
# ╠═0f68fa64-55a5-11eb-24ed-fd2b3fc9ed68
# ╠═528ef22e-55a5-11eb-2b03-3d92a4b9c3c6
# ╠═3de00658-55a5-11eb-102e-3103b223001d
# ╠═c094cb3a-4117-11eb-168a-df762ad5739e
# ╠═a837ef5c-55a5-11eb-06f5-cb27048bc9a4
# ╠═c5cecd9e-55a5-11eb-343a-a3c29893be1f
# ╠═cff2f4c4-55a5-11eb-32c2-1db64f8ea677
# ╠═74480530-55a6-11eb-0a7b-8b30ed458234
# ╠═b0f10e8e-55a6-11eb-34a1-755dceb63e08
# ╠═ef0bdbcc-55a6-11eb-114f-7d43ab25736a
# ╠═7fec90b4-411b-11eb-04db-a913940a9f22
# ╠═451edba2-4122-11eb-0ea9-b3e5a7021bd1
# ╠═46652b7e-4122-11eb-1576-29e34caf1b03
# ╠═78e0e354-4122-11eb-1a11-0b1d36a38b70
# ╠═46798970-4122-11eb-31a8-ab967162b1ae
# ╠═20c55224-55a7-11eb-2490-010a8c63cd4b
# ╠═c98691d4-4125-11eb-3e71-d1ffa4fc0f3a
# ╠═aa066c66-4126-11eb-1909-713aeda74edc
# ╠═c4375a50-4126-11eb-2667-7b693682e9a8
# ╠═aa323b5c-4126-11eb-18fa-a5aa100d0678
# ╠═aa487444-4126-11eb-306c-1d98f391cfc6
# ╠═aa708326-4126-11eb-2129-df2f5666aa93

max_stars_repo_path: exam1 (copy).jl

max_issues_repo_path: exam1 (copy).jl

max_forks_repo_path: exam1 (copy).jl


text:
module TryCooperative
import PatternDispatch.Nodes.julia_signature_of
using PatternDispatch, PatternDispatch.PartialOrder, PatternDispatch.Dispatch
using PatternDispatch.DecisionTree, PatternDispatch.Patterns

@pattern f(1)      = 42
@pattern f(x::Int) = x
@pattern f(x)      = 0

mt = PatternDispatch.method_tables[f]

nodes = subDAGof(mt.top)
methods = methodsof(mt)
actual_methods = filter(m->(m != nomethod), methods)
hullTs = Set{Tuple}(Tuple[m.hullT for m in actual_methods])

ltT(S,T) = !(S==T) && (S<:T)
function dominated_wrt(node::MethodNode, cut::Intension)
    any([intentof(child) & cut == intentof(node) & cut for child in node.gt])
end

for hullT in hullTs
    top = copyDAG(mt.top)

    # filter out too specific methods
    keep = Set{Method}(filter(m->!ltT(m.hullT, hullT), methods))
    @assert !isempty(keep)
    raw_filter!(top, keep)

    # filter out non-questions
    hull = intension(hullT)
    top = simplify!(top, hull)

    @show hullT
    @show tuple([node.value.sig for node in subDAGof(top)]...)
    @show top.value.sig
    println()
end

end # module


text:
# a newspaper headline:
str = "The Gold and Blue Loses a Bit of Its Luster"
println(str)
nchars = length(str)
println("The headline counts $nchars characters") # 43
str2 = replace(str, "Blue", "Red")
# strings are immutable
println(str) # The Gold and Blue Loses a Bit of Its Luster
println(str2)
println("Here are the characters at position 25 to 30:")
subs = str[25:30]
print("-$(lowercase(subs))-") # "-a bit -"
println("Here are all the characters:")
for c in str
    println(c)
end
arr = split(str,' ')
show(arr) # ["The","Gold","and","Blue","Loses","a","Bit","of","Its","Luster"]
nwords = length(arr)
println("The headline counts $nwords words") # 10
println("Here are all the words:")
for word in arr
    println(word)
end
arr[4] = "Red"
show(arr) # arrays are mutable
println("Convert back to a sentence:")
nstr = join(arr, ' ')
println(nstr) # The Gold and Red Loses a Bit of Its Luster
# working with arrays:
println("arrays: calculates the sum, mean and standard deviation ")
arr = [1:100]
typeof(arr) #>
println(sum(arr)) #> 5050
println(mean(arr)) #> 50.5
println(std(arr)) #> 29.011491975882016

max_stars_repo_path: Module 1/Chapter02/strings_arrays.jl

max_issues_repo_path: Module 1/Chapter02/strings_arrays.jl

max_forks_repo_path: Module 1/Chapter02/strings_arrays.jl


text:
# if isinstalled("OhMyREPL")
#     @eval using OhMyREPL
#     enable_autocomplete_brackets(false)
# else
#     warn("'OhMyREPL' is not installed")
# end

atreplinit() do repl
    try
        @eval using Revise
        @async Revise.wait_steal_repl_backend()
        @info "Using Revise"
    catch
        @warn "Loading Revise failed!"
    end
end

@info "Startup done."

max_stars_repo_name: phipsgabler/dotfiles

max_issues_repo_name: phipsgabler/dotfiles

max_forks_repo_name: phipsgabler/dotfiles


text:
using Jtac

G = MetaTac

chain = @chain G begin
          @stack stack_input = true Conv(512, relu, stride = 3) Conv(1024, relu)
	  Dense(1024, relu)
          Dropout(0.3)
	  Dense(512, relu)
          Dropout(0.3)
          Dense(256, relu)
          Dropout(0.3)
          Dense(256, relu)
        end

#chain = @chain G Conv(2048) Dense(1024) Dense(512) Dense(256) Dense(128)
model = NeuralModel(G, chain) |> to_gpu |> Async
player = MCTSPlayer(model, power = 250)
#player = MCTSPlayer(power = 50)

max_stars_repo_name: roSievers/jtac


text:
using ExamplePackage
using Test

 @test foo() == 0.11388071406436832
 @test foo(1, 1.5) == 0.2731856314283442
 @test_broken foo(1, 0) # tells us this is broken
 @test_broken @inferred foo(3) # tells us we're type-unstable