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

@test integrate(x^m*sinhintegral(b*x), x) == :(x^(1+m)*(1+m)^-1*sinhintegral(b*x)+-1//2*b^-1*x^m*(b*x)^(-1m)*(1+m)^-1*SpecialFunctions.gamma(1+m,b*x)+-1//2*b^-1*x^m*(-1*b*x)^(-1m)*(1+m)^-1*SpecialFunctions.gamma(1+m,-1*b*x))
@test integrate(x^3*sinhintegral(b*x), x) == :(1//4*x^4*sinhintegral(b*x)+3//2*b^-4*sinh(b*x)+-3//2*x*b^-3*cosh(b*x)+-1//4*b^-1*x^3*cosh(b*x)+3//4*b^-2*x^2*sinh(b*x))
@test integrate(x^2*sinhintegral(b*x), x) == :(-2//3*b^-3*cosh(b*x)+1//3*x^3*sinhintegral(b*x)+-1//3*b^-1*x^2*cosh(b*x)+2//3*x*b^-2*sinh(b*x))
@test integrate(x*sinhintegral(b*x), x) == :((1/2)*b^-2*sinh(b*x)+(1/2)*x^2*sinhintegral(b*x)+-1//2*x*b^-1*cosh(b*x))
@test integrate(sinhintegral(b*x), x) == :(x*sinhintegral(b*x)+-1*b^-1*cosh(b*x))
@test integrate(x^-2*sinhintegral(b*x), x) == :(b*coshintegral(b*x)+-1*x^-1*sinh(b*x)+-1*x^-1*sinhintegral(b*x))
@test integrate(x^-3*sinhintegral(b*x), x) == :(-1//2*x^-2*sinhintegral(b*x)+-1//4*x^-2*sinh(b*x)+1//4*b^2*sinhintegral(b*x)+-1//4*b*x^-1*cosh(b*x))
@test integrate(x^3*sinhintegral(b*x)^2, x) == :((1/2)*b^-2*x^2+2*b^-4*sinh(b*x)^2+-3//2*b^-4*coshintegral(2*b*x)+1//4*x^4*sinhintegral(b*x)^2+3//2*b^-4*log(x)+3*b^-4*sinh(b*x)*sinhintegral(b*x)+1//4*b^-2*x^2*sinh(b*x)^2+-1*x*b^-3*cosh(b*x)*sinh(b*x)+-3*x*b^-3*cosh(b*x)*sinhintegral(b*x)+-1//2*b^-1*x^3*cosh(b*x)*sinhintegral(b*x)+3//2*b^-2*x^2*sinh(b*x)*sinhintegral(b*x))
@test integrate(x^2*sinhintegral(b*x)^2, x) == :(1//3*x^3*sinhintegral(b*x)^2+2//3*b^-3*sinhintegral(2*b*x)+5//6*x*b^-2+-5//6*b^-3*cosh(b*x)*sinh(b*x)+-4//3*b^-3*cosh(b*x)*sinhintegral(b*x)+1//3*x*b^-2*sinh(b*x)^2+-2//3*b^-1*x^2*cosh(b*x)*sinhintegral(b*x)+4//3*x*b^-2*sinh(b*x)*sinhintegral(b*x))
@test integrate(x*sinhintegral(b*x)^2, x) == :((1/2)*b^-2*sinh(b*x)^2+(1/2)*b^-2*log(x)+(1/2)*x^2*sinhintegral(b*x)^2+-1//2*b^-2*coshintegral(2*b*x)+b^-2*sinh(b*x)*sinhintegral(b*x)+-1*x*b^-1*cosh(b*x)*sinhintegral(b*x))
@test integrate(sinhintegral(b*x)^2, x) == :(x*sinhintegral(b*x)^2+b^-1*sinhintegral(2*b*x)+-2*b^-1*cosh(b*x)*sinhintegral(b*x))
@test integrate(x^3*sinhintegral(a+b*x), x) == :(1//4*x^4*sinhintegral(a+b*x)+3//2*b^-4*sinh(a+b*x)+(1/2)*a*b^-4*cosh(a+b*x)+-3//2*x*b^-3*cosh(a+b*x)+-1//4*a^4*b^-4*sinhintegral(a+b*x)+-1//4*b^-1*x^3*cosh(a+b*x)+1//4*a^2*b^-4*sinh(a+b*x)+1//4*a^3*b^-4*cosh(a+b*x)+3//4*b^-2*x^2*sinh(a+b*x)+-1//2*a*x*b^-3*sinh(a+b*x)+-1//4*x*a^2*b^-3*cosh(a+b*x)+1//4*a*b^-2*x^2*cosh(a+b*x))
@test integrate(x^2*sinhintegral(a+b*x), x) == :(-2//3*b^-3*cosh(a+b*x)+1//3*x^3*sinhintegral(a+b*x)+-1//3*a*b^-3*sinh(a+b*x)+-1//3*a^2*b^-3*cosh(a+b*x)+-1//3*b^-1*x^2*cosh(a+b*x)+1//3*a^3*b^-3*sinhintegral(a+b*x)+2//3*x*b^-2*sinh(a+b*x)+1//3*a*x*b^-2*cosh(a+b*x))
@test integrate(x*sinhintegral(a+b*x), x) == :((1/2)*b^-2*sinh(a+b*x)+(1/2)*x^2*sinhintegral(a+b*x)+(1/2)*a*b^-2*cosh(a+b*x)+-1//2*x*b^-1*cosh(a+b*x)+-1//2*a^2*b^-2*sinhintegral(a+b*x))
@test integrate(sinhintegral(a+b*x), x) == :(-1*b^-1*cosh(a+b*x)+b^-1*(a+b*x)*sinhintegral(a+b*x))
@test integrate(x^-2*sinhintegral(a+b*x), x) == :(-1*x^-1*sinhintegral(a+b*x)+-1*b*a^-1*sinhintegral(a+b*x)+b*a^-1*cosh(a)*sinhintegral(b*x)+b*a^-1*coshintegral(b*x)*sinh(a))
@test integrate(x^-3*sinhintegral(a+b*x), x) == :(-1//2*x^-2*sinhintegral(a+b*x)+(1/2)*a^-2*b^2*sinhintegral(a+b*x)+(1/2)*a^-1*b^2*cosh(a)*coshintegral(b*x)+(1/2)*a^-1*b^2*sinh(a)*sinhintegral(b*x)+-1//2*b*a^-1*x^-1*sinh(a+b*x)+-1//2*a^-2*b^2*cosh(a)*sinhintegral(b*x)+-1//2*a^-2*b^2*coshintegral(b*x)*sinh(a))
@test integrate(x^2*sinhintegral(a+b*x)^2, x) == :(-1//12*b^-3*sinh(2a+2*b*x)+2//3*x*b^-2+2//3*b^-3*sinhintegral(2a+2*b*x)+a*b^-3*coshintegral(2a+2*b*x)+a^2*b^-3*sinhintegral(2a+2*b*x)+-1*a*b^-3*log(a+b*x)+-4//3*b^-3*cosh(a+b*x)*sinhintegral(a+b*x)+-2//3*b^-3*cosh(a+b*x)*sinh(a+b*x)+-1//3*a*b^-3*cosh(2a+2*b*x)+1//6*x*b^-2*cosh(2a+2*b*x)+-2//3*a*b^-3*sinh(a+b*x)*sinhintegral(a+b*x)+-2//3*a^2*b^-3*cosh(a+b*x)*sinhintegral(a+b*x)+-2//3*b^-1*x^2*cosh(a+b*x)*sinhintegral(a+b*x)+1//3*a^2*b^-3*sinhintegral(a+b*x)^2*(a+b*x)+1//3*b^-1*x^2*sinhintegral(a+b*x)^2*(a+b*x)+4//3*x*b^-2*sinh(a+b*x)*sinhintegral(a+b*x)+-1//3*a*x*b^-2*sinhintegral(a+b*x)^2*(a+b*x)+2//3*a*x*b^-2*cosh(a+b*x)*sinhintegral(a+b*x))
@test integrate(x*sinhintegral(a+b*x)^2, x) == :((1/2)*b^-2*log(a+b*x)+-1//2*b^-2*coshintegral(2a+2*b*x)+1//4*b^-2*cosh(2a+2*b*x)+b^-2*sinh(a+b*x)*sinhintegral(a+b*x)+-1*a*b^-2*sinhintegral(2a+2*b*x)+a*b^-2*cosh(a+b*x)*sinhintegral(a+b*x)+(1/2)*x*b^-1*sinhintegral(a+b*x)^2*(a+b*x)+-1*x*b^-1*cosh(a+b*x)*sinhintegral(a+b*x)+-1//2*a*b^-2*sinhintegral(a+b*x)^2*(a+b*x))
@test integrate(sinhintegral(a+b*x)^2, x) == :(b^-1*sinhintegral(2a+2*b*x)+b^-1*sinhintegral(a+b*x)^2*(a+b*x)+-2*b^-1*cosh(a+b*x)*sinhintegral(a+b*x))
@test integrate(x^2*sinhintegral(d*(a+b*log(c*x^n))), x) == :(1//3*x^3*sinhintegral(d*(a+b*log(c*x^n)))+-1//6*x^3*(c*x^n)^(-3*n^-1)*exp(-3*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(3+b*d*n)*(a+b*log(c*x^n)))+1//6*x^3*(c*x^n)^(-3*n^-1)*exp(-3*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(3+-1*b*d*n)*(a+b*log(c*x^n))))
@test integrate(x*sinhintegral(d*(a+b*log(c*x^n))), x) == :((1/2)*x^2*sinhintegral(d*(a+b*log(c*x^n)))+-1//4*x^2*(c*x^n)^(-2*n^-1)*exp(-2*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(2+b*d*n)*(a+b*log(c*x^n)))+1//4*x^2*(c*x^n)^(-2*n^-1)*exp(-2*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(2+-1*b*d*n)*(a+b*log(c*x^n))))
@test integrate(sinhintegral(d*(a+b*log(c*x^n))), x) == :(x*sinhintegral(d*(a+b*log(c*x^n)))+(1/2)*x*(c*x^n)^(-1*n^-1)*exp(-1*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(1+-1*b*d*n)*(a+b*log(c*x^n)))+-1//2*x*(c*x^n)^(-1*n^-1)*exp(-1*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(1+b*d*n)*(a+b*log(c*x^n))))
@test integrate(x^-1*sinhintegral(d*(a+b*log(c*x^n))), x) == :(b^-1*n^-1*(a+b*log(c*x^n))*sinhintegral(d*(a+b*log(c*x^n)))+-1*b^-1*d^-1*n^-1*cosh(d*(a+b*log(c*x^n))))
@test integrate(x^-2*sinhintegral(d*(a+b*log(c*x^n))), x) == :(-1*x^-1*sinhintegral(d*(a+b*log(c*x^n)))+(1/2)*x^-1*(c*x^n)^(n^-1)*exp(a*b^-1*n^-1)*expintegralei(-1*b^-1*n^-1*(1+-1*b*d*n)*(a+b*log(c*x^n)))+-1//2*x^-1*(c*x^n)^(n^-1)*exp(a*b^-1*n^-1)*expintegralei(-1*b^-1*n^-1*(1+b*d*n)*(a+b*log(c*x^n))))
@test integrate(x^-3*sinhintegral(d*(a+b*log(c*x^n))), x) == :(-1//2*x^-2*sinhintegral(d*(a+b*log(c*x^n)))+-1//4*x^-2*(c*x^n)^(2*n^-1)*exp(2*a*b^-1*n^-1)*expintegralei(-1*b^-1*n^-1*(2+b*d*n)*(a+b*log(c*x^n)))+1//4*x^-2*(c*x^n)^(2*n^-1)*exp(2*a*b^-1*n^-1)*expintegralei(-1*b^-1*n^-1*(2+-1*b*d*n)*(a+b*log(c*x^n))))
@test integrate((e*x)^m*sinhintegral(d*(a+b*log(c*x^n))), x) == :(e^-1*(e*x)^(1+m)*(1+m)^-1*sinhintegral(d*(a+b*log(c*x^n)))+(1/2)*x*(c*x^n)^(-1*n^-1*(1+m))*(e*x)^m*(1+m)^-1*exp(-1*a*b^-1*n^-1*(1+m))*expintegralei(b^-1*n^-1*(a+b*log(c*x^n))*(1+m+-1*b*d*n))+-1//2*x*(c*x^n)^(-1*n^-1*(1+m))*(e*x)^m*(1+m)^-1*exp(-1*a*b^-1*n^-1*(1+m))*expintegralei(b^-1*n^-1*(a+b*log(c*x^n))*(1+m+b*d*n)))
@test integrate(x^-3*sinh(b*x)*sinhintegral(b*x), x) == :(b^2*coshintegral(2*b*x)+-1//4*x^-2*sinh(b*x)^2+1//4*b^2*sinhintegral(b*x)^2+-1//2*x^-2*sinh(b*x)*sinhintegral(b*x)+-1//4*b*x^-1*sinh(2*b*x)+-1//2*b*x^-1*cosh(b*x)*sinh(b*x)+-1//2*b*x^-1*cosh(b*x)*sinhintegral(b*x))
@test integrate(x^-1*sinh(b*x)*sinhintegral(b*x), x) == :((1/2)*sinhintegral(b*x)^2)
@test integrate(sinh(b*x)*sinhintegral(b*x), x) == :(-1//2*b^-1*sinhintegral(2*b*x)+b^-1*cosh(b*x)*sinhintegral(b*x))
@test integrate(x*sinh(b*x)*sinhintegral(b*x), x) == :((1/2)*b^-2*coshintegral(2*b*x)+-1//2*b^-2*sinh(b*x)^2+-1//2*b^-2*log(x)+-1*b^-2*sinh(b*x)*sinhintegral(b*x)+x*b^-1*cosh(b*x)*sinhintegral(b*x))
@test integrate(x^2*sinh(b*x)*sinhintegral(b*x), x) == :(-1*b^-3*sinhintegral(2*b*x)+-5//4*x*b^-2+2*b^-3*cosh(b*x)*sinhintegral(b*x)+-1//2*x*b^-2*sinh(b*x)^2+5//4*b^-3*cosh(b*x)*sinh(b*x)+b^-1*x^2*cosh(b*x)*sinhintegral(b*x)+-2*x*b^-2*sinh(b*x)*sinhintegral(b*x))
@test integrate(x^3*sinh(b*x)*sinhintegral(b*x), x) == :(-1*b^-2*x^2+-4*b^-4*sinh(b*x)^2+-3*b^-4*log(x)+3*b^-4*coshintegral(2*b*x)+-6*b^-4*sinh(b*x)*sinhintegral(b*x)+-1//2*b^-2*x^2*sinh(b*x)^2+b^-1*x^3*cosh(b*x)*sinhintegral(b*x)+-3*b^-2*x^2*sinh(b*x)*sinhintegral(b*x)+2*x*b^-3*cosh(b*x)*sinh(b*x)+6*x*b^-3*cosh(b*x)*sinhintegral(b*x))
@test integrate(x^-2*cosh(b*x)*sinhintegral(b*x), x) == :(b*coshintegral(2*b*x)+(1/2)*b*sinhintegral(b*x)^2+-1//2*x^-1*sinh(2*b*x)+-1*x^-1*cosh(b*x)*sinhintegral(b*x))
@test integrate(cosh(b*x)*sinhintegral(b*x), x) == :((1/2)*b^-1*log(x)+-1//2*b^-1*coshintegral(2*b*x)+b^-1*sinh(b*x)*sinhintegral(b*x))
@test integrate(x*cosh(b*x)*sinhintegral(b*x), x) == :((1/2)*x*b^-1+(1/2)*b^-2*sinhintegral(2*b*x)+-1*b^-2*cosh(b*x)*sinhintegral(b*x)+-1//2*b^-2*cosh(b*x)*sinh(b*x)+x*b^-1*sinh(b*x)*sinhintegral(b*x))
@test integrate(x^2*cosh(b*x)*sinhintegral(b*x), x) == :(b^-3*log(x)+-1*b^-3*coshintegral(2*b*x)+1//4*b^-1*x^2+5//4*b^-3*sinh(b*x)^2+2*b^-3*sinh(b*x)*sinhintegral(b*x)+b^-1*x^2*sinh(b*x)*sinhintegral(b*x)+-2*x*b^-2*cosh(b*x)*sinhintegral(b*x)+-1//2*x*b^-2*cosh(b*x)*sinh(b*x))
@test integrate(x^3*cosh(b*x)*sinhintegral(b*x), x) == :(3*b^-4*sinhintegral(2*b*x)+4*x*b^-3+1//6*b^-1*x^3+-6*b^-4*cosh(b*x)*sinhintegral(b*x)+-4*b^-4*cosh(b*x)*sinh(b*x)+2*x*b^-3*sinh(b*x)^2+b^-1*x^3*sinh(b*x)*sinhintegral(b*x)+-3*b^-2*x^2*cosh(b*x)*sinhintegral(b*x)+6*x*b^-3*sinh(b*x)*sinhintegral(b*x)+-1//2*b^-2*x^2*cosh(b*x)*sinh(b*x))
@test integrate(sinh(5x)*sinhintegral(2x), x) == :(-1//10*sinhintegral(7x)+1//10*sinhintegral(3x)+1//5*cosh(5x)*sinhintegral(2x))
@test integrate(cosh(5x)*sinhintegral(2x), x) == :(-1//10*coshintegral(7x)+1//10*coshintegral(3x)+1//5*sinh(5x)*sinhintegral(2x))
@test integrate(x^2*sinh(a+b*x)*sinhintegral(a+b*x), x) == :(-1*x*b^-2+-1*b^-3*sinhintegral(2a+2*b*x)+1//8*b^-3*sinh(2a+2*b*x)+a*b^-3*log(a+b*x)+b^-3*cosh(a+b*x)*sinh(a+b*x)+-1*a*b^-3*coshintegral(2a+2*b*x)+2*b^-3*cosh(a+b*x)*sinhintegral(a+b*x)+-1//2*a^2*b^-3*sinhintegral(2a+2*b*x)+-1//4*x*b^-2*cosh(2a+2*b*x)+1//4*a*b^-3*cosh(2a+2*b*x)+b^-1*x^2*cosh(a+b*x)*sinhintegral(a+b*x)+-2*x*b^-2*sinh(a+b*x)*sinhintegral(a+b*x))
@test integrate(x*sinh(a+b*x)*sinhintegral(a+b*x), x) == :((1/2)*b^-2*coshintegral(2a+2*b*x)+-1//2*b^-2*log(a+b*x)+-1//4*b^-2*cosh(2a+2*b*x)+(1/2)*a*b^-2*sinhintegral(2a+2*b*x)+-1*b^-2*sinh(a+b*x)*sinhintegral(a+b*x)+x*b^-1*cosh(a+b*x)*sinhintegral(a+b*x))
@test integrate(sinh(a+b*x)*sinhintegral(a+b*x), x) == :(-1//2*b^-1*sinhintegral(2a+2*b*x)+b^-1*cosh(a+b*x)*sinhintegral(a+b*x))
@test integrate(x^2*cosh(a+b*x)*sinhintegral(a+b*x), x) == :(b^-3*log(a+b*x)+(1/2)*b^-3*cosh(2a+2*b*x)+-1*b^-3*coshintegral(2a+2*b*x)+1//4*b^-1*x^2+1//4*b^-3*sinh(a+b*x)^2+(1/2)*a^2*b^-3*log(a+b*x)+-1*a*b^-3*sinhintegral(2a+2*b*x)+2*b^-3*sinh(a+b*x)*sinhintegral(a+b*x)+-1//2*a*x*b^-2+-1//2*a^2*b^-3*coshintegral(2a+2*b*x)+b^-1*x^2*sinh(a+b*x)*sinhintegral(a+b*x)+(1/2)*a*b^-3*cosh(a+b*x)*sinh(a+b*x)+-2*x*b^-2*cosh(a+b*x)*sinhintegral(a+b*x)+-1//2*x*b^-2*cosh(a+b*x)*sinh(a+b*x))
@test integrate(x*cosh(a+b*x)*sinhintegral(a+b*x), x) == :((1/2)*x*b^-1+(1/2)*b^-2*sinhintegral(2a+2*b*x)+(1/2)*a*b^-2*coshintegral(2a+2*b*x)+-1*b^-2*cosh(a+b*x)*sinhintegral(a+b*x)+-1//2*a*b^-2*log(a+b*x)+-1//2*b^-2*cosh(a+b*x)*sinh(a+b*x)+x*b^-1*sinh(a+b*x)*sinhintegral(a+b*x))
@test integrate(cosh(a+b*x)*sinhintegral(a+b*x), x) == :((1/2)*b^-1*log(a+b*x)+-1//2*b^-1*coshintegral(2a+2*b*x)+b^-1*sinh(a+b*x)*sinhintegral(a+b*x))
@test integrate(x*sinh(a+b*x)*sinhintegral(c+d*x), x) == :((1/2)*b^-1*(b+-1d)^-1*cosh(a+-1c+x*(b+-1d))+(1/2)*b^-2*cosh(a+-1*b*c*d^-1)*coshintegral(x*(b+d)+c*d^-1*(b+d))+(1/2)*b^-2*sinh(a+-1*b*c*d^-1)*sinhintegral(x*(b+d)+c*d^-1*(b+d))+-1*b^-2*sinh(a+b*x)*sinhintegral(c+d*x)+-1//2*b^-1*(b+d)^-1*cosh(a+c+x*(b+d))+-1//2*b^-2*cosh(a+-1*b*c*d^-1)*coshintegral(x*(b+-1d)+c*d^-1*(b+-1d))+-1//2*b^-2*sinh(a+-1*b*c*d^-1)*sinhintegral(x*(b+-1d)+c*d^-1*(b+-1d))+x*b^-1*cosh(a+b*x)*sinhintegral(c+d*x)+(1/2)*c*b^-1*d^-1*cosh(a+-1*b*c*d^-1)*sinhintegral(x*(b+d)+c*d^-1*(b+d))+(1/2)*c*b^-1*d^-1*coshintegral(x*(b+d)+c*d^-1*(b+d))*sinh(a+-1*b*c*d^-1)+-1//2*c*b^-1*d^-1*cosh(a+-1*b*c*d^-1)*sinhintegral(x*(b+-1d)+c*d^-1*(b+-1d))+-1//2*c*b^-1*d^-1*coshintegral(x*(b+-1d)+c*d^-1*(b+-1d))*sinh(a+-1*b*c*d^-1))
@test integrate(sinh(a+b*x)*sinhintegral(c+d*x), x) == :(b^-1*cosh(a+b*x)*sinhintegral(c+d*x)+(1/2)*b^-1*cosh(a+-1*b*c*d^-1)*sinhintegral(x*(b+-1d)+c*d^-1*(b+-1d))+(1/2)*b^-1*coshintegral(x*(b+-1d)+c*d^-1*(b+-1d))*sinh(a+-1*b*c*d^-1)+-1//2*b^-1*cosh(a+-1*b*c*d^-1)*sinhintegral(x*(b+d)+c*d^-1*(b+d))+-1//2*b^-1*coshintegral(x*(b+d)+c*d^-1*(b+d))*sinh(a+-1*b*c*d^-1))
@test integrate(x*cosh(a+b*x)*sinhintegral(c+d*x), x) == :((1/2)*b^-1*(b+-1d)^-1*sinh(a+-1c+x*(b+-1d))+(1/2)*b^-2*cosh(a+-1*b*c*d^-1)*sinhintegral(x*(b+d)+c*d^-1*(b+d))+(1/2)*b^-2*coshintegral(x*(b+d)+c*d^-1*(b+d))*sinh(a+-1*b*c*d^-1)+-1*b^-2*cosh(a+b*x)*sinhintegral(c+d*x)+-1//2*b^-1*(b+d)^-1*sinh(a+c+x*(b+d))+-1//2*b^-2*cosh(a+-1*b*c*d^-1)*sinhintegral(x*(b+-1d)+c*d^-1*(b+-1d))+-1//2*b^-2*coshintegral(x*(b+-1d)+c*d^-1*(b+-1d))*sinh(a+-1*b*c*d^-1)+x*b^-1*sinh(a+b*x)*sinhintegral(c+d*x)+(1/2)*c*b^-1*d^-1*cosh(a+-1*b*c*d^-1)*coshintegral(x*(b+d)+c*d^-1*(b+d))+(1/2)*c*b^-1*d^-1*sinh(a+-1*b*c*d^-1)*sinhintegral(x*(b+d)+c*d^-1*(b+d))+-1//2*c*b^-1*d^-1*cosh(a+-1*b*c*d^-1)*coshintegral(x*(b+-1d)+c*d^-1*(b+-1d))+-1//2*c*b^-1*d^-1*sinh(a+-1*b*c*d^-1)*sinhintegral(x*(b+-1d)+c*d^-1*(b+-1d)))
@test integrate(cosh(a+b*x)*sinhintegral(c+d*x), x) == :(b^-1*sinh(a+b*x)*sinhintegral(c+d*x)+(1/2)*b^-1*cosh(a+-1*b*c*d^-1)*coshintegral(x*(b+-1d)+c*d^-1*(b+-1d))+(1/2)*b^-1*sinh(a+-1*b*c*d^-1)*sinhintegral(x*(b+-1d)+c*d^-1*(b+-1d))+-1//2*b^-1*cosh(a+-1*b*c*d^-1)*coshintegral(x*(b+d)+c*d^-1*(b+d))+-1//2*b^-1*sinh(a+-1*b*c*d^-1)*sinhintegral(x*(b+d)+c*d^-1*(b+d)))
@test integrate(x^m*coshintegral(b*x), x) == :(x^(1+m)*(1+m)^-1*coshintegral(b*x)+(1/2)*b^-1*x^m*(b*x)^(-1m)*(1+m)^-1*SpecialFunctions.gamma(1+m,b*x)+-1//2*b^-1*x^m*(-1*b*x)^(-1m)*(1+m)^-1*SpecialFunctions.gamma(1+m,-1*b*x))
@test integrate(x^3*coshintegral(b*x), x) == :(1//4*x^4*coshintegral(b*x)+3//2*b^-4*cosh(b*x)+-3//2*x*b^-3*sinh(b*x)+-1//4*b^-1*x^3*sinh(b*x)+3//4*b^-2*x^2*cosh(b*x))
@test integrate(x^2*coshintegral(b*x), x) == :(-2//3*b^-3*sinh(b*x)+1//3*x^3*coshintegral(b*x)+-1//3*b^-1*x^2*sinh(b*x)+2//3*x*b^-2*cosh(b*x))
@test integrate(x*coshintegral(b*x), x) == :((1/2)*b^-2*cosh(b*x)+(1/2)*x^2*coshintegral(b*x)+-1//2*x*b^-1*sinh(b*x))
@test integrate(coshintegral(b*x), x) == :(x*coshintegral(b*x)+-1*b^-1*sinh(b*x))
@test integrate(x^-2*coshintegral(b*x), x) == :(b*sinhintegral(b*x)+-1*x^-1*cosh(b*x)+-1*x^-1*coshintegral(b*x))
@test integrate(x^-3*coshintegral(b*x), x) == :(-1//2*x^-2*coshintegral(b*x)+-1//4*x^-2*cosh(b*x)+1//4*b^2*coshintegral(b*x)+-1//4*b*x^-1*sinh(b*x))
@test integrate(x^3*coshintegral(b*x)^2, x) == :(-3//2*b^-4*coshintegral(2*b*x)+-3//2*b^-4*log(x)+-1//4*b^-2*x^2+1//4*x^4*coshintegral(b*x)^2+3//8*b^-4*cosh(b*x)^2+13//8*b^-4*sinh(b*x)^2+3*b^-4*cosh(b*x)*coshintegral(b*x)+1//4*b^-2*x^2*sinh(b*x)^2+-1*x*b^-3*cosh(b*x)*sinh(b*x)+-3*x*b^-3*coshintegral(b*x)*sinh(b*x)+-1//2*b^-1*x^3*coshintegral(b*x)*sinh(b*x)+3//2*b^-2*x^2*cosh(b*x)*coshintegral(b*x))
@test integrate(x^2*coshintegral(b*x)^2, x) == :(-1//2*x*b^-2+1//3*x^3*coshintegral(b*x)^2+2//3*b^-3*sinhintegral(2*b*x)+-5//6*b^-3*cosh(b*x)*sinh(b*x)+-4//3*b^-3*coshintegral(b*x)*sinh(b*x)+1//3*x*b^-2*sinh(b*x)^2+-2//3*b^-1*x^2*coshintegral(b*x)*sinh(b*x)+4//3*x*b^-2*cosh(b*x)*coshintegral(b*x))
@test integrate(x*coshintegral(b*x)^2, x) == :((1/2)*b^-2*sinh(b*x)^2+(1/2)*x^2*coshintegral(b*x)^2+-1//2*b^-2*coshintegral(2*b*x)+-1//2*b^-2*log(x)+b^-2*cosh(b*x)*coshintegral(b*x)+-1*x*b^-1*coshintegral(b*x)*sinh(b*x))
@test integrate(coshintegral(b*x)^2, x) == :(x*coshintegral(b*x)^2+b^-1*sinhintegral(2*b*x)+-2*b^-1*coshintegral(b*x)*sinh(b*x))
@test integrate(x^3*coshintegral(a+b*x), x) == :(1//4*x^4*coshintegral(a+b*x)+3//2*b^-4*cosh(a+b*x)+(1/2)*a*b^-4*sinh(a+b*x)+-3//2*x*b^-3*sinh(a+b*x)+-1//4*a^4*b^-4*coshintegral(a+b*x)+-1//4*b^-1*x^3*sinh(a+b*x)+1//4*a^2*b^-4*cosh(a+b*x)+1//4*a^3*b^-4*sinh(a+b*x)+3//4*b^-2*x^2*cosh(a+b*x)+-1//2*a*x*b^-3*cosh(a+b*x)+-1//4*x*a^2*b^-3*sinh(a+b*x)+1//4*a*b^-2*x^2*sinh(a+b*x))
@test integrate(x^2*coshintegral(a+b*x), x) == :(-2//3*b^-3*sinh(a+b*x)+1//3*x^3*coshintegral(a+b*x)+-1//3*a*b^-3*cosh(a+b*x)+-1//3*a^2*b^-3*sinh(a+b*x)+-1//3*b^-1*x^2*sinh(a+b*x)+1//3*a^3*b^-3*coshintegral(a+b*x)+2//3*x*b^-2*cosh(a+b*x)+1//3*a*x*b^-2*sinh(a+b*x))
@test integrate(x*coshintegral(a+b*x), x) == :((1/2)*b^-2*cosh(a+b*x)+(1/2)*x^2*coshintegral(a+b*x)+(1/2)*a*b^-2*sinh(a+b*x)+-1//2*x*b^-1*sinh(a+b*x)+-1//2*a^2*b^-2*coshintegral(a+b*x))
@test integrate(coshintegral(a+b*x), x) == :(-1*b^-1*sinh(a+b*x)+b^-1*(a+b*x)*coshintegral(a+b*x))
@test integrate(x^-2*coshintegral(a+b*x), x) == :(-1*x^-1*coshintegral(a+b*x)+-1*b*a^-1*coshintegral(a+b*x)+b*a^-1*cosh(a)*coshintegral(b*x)+b*a^-1*sinh(a)*sinhintegral(b*x))
@test integrate(x^-3*coshintegral(a+b*x), x) == :(-1//2*x^-2*coshintegral(a+b*x)+(1/2)*a^-2*b^2*coshintegral(a+b*x)+(1/2)*a^-1*b^2*cosh(a)*sinhintegral(b*x)+(1/2)*a^-1*b^2*coshintegral(b*x)*sinh(a)+-1//2*b*a^-1*x^-1*cosh(a+b*x)+-1//2*a^-2*b^2*cosh(a)*coshintegral(b*x)+-1//2*a^-2*b^2*sinh(a)*sinhintegral(b*x))
@test integrate(x^2*coshintegral(a+b*x)^2, x) == :(-2//3*x*b^-2+-1//12*b^-3*sinh(2a+2*b*x)+2//3*b^-3*sinhintegral(2a+2*b*x)+a*b^-3*coshintegral(2a+2*b*x)+a*b^-3*log(a+b*x)+a^2*b^-3*sinhintegral(2a+2*b*x)+-4//3*b^-3*coshintegral(a+b*x)*sinh(a+b*x)+-2//3*b^-3*cosh(a+b*x)*sinh(a+b*x)+-1//3*a*b^-3*cosh(2a+2*b*x)+1//6*x*b^-2*cosh(2a+2*b*x)+-2//3*a*b^-3*cosh(a+b*x)*coshintegral(a+b*x)+-2//3*a^2*b^-3*coshintegral(a+b*x)*sinh(a+b*x)+-2//3*b^-1*x^2*coshintegral(a+b*x)*sinh(a+b*x)+1//3*a^2*b^-3*coshintegral(a+b*x)^2*(a+b*x)+1//3*b^-1*x^2*coshintegral(a+b*x)^2*(a+b*x)+4//3*x*b^-2*cosh(a+b*x)*coshintegral(a+b*x)+-1//3*a*x*b^-2*coshintegral(a+b*x)^2*(a+b*x)+2//3*a*x*b^-2*coshintegral(a+b*x)*sinh(a+b*x))
@test integrate(x*coshintegral(a+b*x)^2, x) == :(-1//2*b^-2*coshintegral(2a+2*b*x)+-1//2*b^-2*log(a+b*x)+1//4*b^-2*cosh(2a+2*b*x)+b^-2*cosh(a+b*x)*coshintegral(a+b*x)+-1*a*b^-2*sinhintegral(2a+2*b*x)+a*b^-2*coshintegral(a+b*x)*sinh(a+b*x)+(1/2)*x*b^-1*coshintegral(a+b*x)^2*(a+b*x)+-1*x*b^-1*coshintegral(a+b*x)*sinh(a+b*x)+-1//2*a*b^-2*coshintegral(a+b*x)^2*(a+b*x))
@test integrate(coshintegral(a+b*x)^2, x) == :(b^-1*sinhintegral(2a+2*b*x)+b^-1*coshintegral(a+b*x)^2*(a+b*x)+-2*b^-1*coshintegral(a+b*x)*sinh(a+b*x))
@test integrate(x^2*coshintegral(d*(a+b*log(c*x^n))), x) == :(1//3*x^3*coshintegral(d*(a+b*log(c*x^n)))+-1//6*x^3*(c*x^n)^(-3*n^-1)*exp(-3*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(3+b*d*n)*(a+b*log(c*x^n)))+-1//6*x^3*(c*x^n)^(-3*n^-1)*exp(-3*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(3+-1*b*d*n)*(a+b*log(c*x^n))))
@test integrate(x*coshintegral(d*(a+b*log(c*x^n))), x) == :((1/2)*x^2*coshintegral(d*(a+b*log(c*x^n)))+-1//4*x^2*(c*x^n)^(-2*n^-1)*exp(-2*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(2+b*d*n)*(a+b*log(c*x^n)))+-1//4*x^2*(c*x^n)^(-2*n^-1)*exp(-2*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(2+-1*b*d*n)*(a+b*log(c*x^n))))
@test integrate(coshintegral(d*(a+b*log(c*x^n))), x) == :(x*coshintegral(d*(a+b*log(c*x^n)))+-1//2*x*(c*x^n)^(-1*n^-1)*exp(-1*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(1+b*d*n)*(a+b*log(c*x^n)))+-1//2*x*(c*x^n)^(-1*n^-1)*exp(-1*a*b^-1*n^-1)*expintegralei(b^-1*n^-1*(1+-1*b*d*n)*(a+b*log(c*x^n))))
@test integrate(x^-1*coshintegral(d*(a+b*log(c*x^n))), x) == :(b^-1*n^-1*(a+b*log(c*x^n))*coshintegral(d*(a+b*log(c*x^n)))+-1*b^-1*d^-1*n^-1*sinh(d*(a+b*log(c*x^n))))
@test integrate(x^-2*coshintegral(d*(a+b*log(c*x^n))), x) == :(-1*x^-1*coshintegral(d*(a+b*log(c*x^n)))+(1/2)*x^-1*(c*x^n)^(n^-1)*exp(a*b^-1*n^-1)*expintegralei(-1*b^-1*n^-1*(1+b*d*n)*(a+b*log(c*x^n)))+(1/2)*x^-1*(c*x^n)^(n^-1)*exp(a*b^-1*n^-1)*expintegralei(-1*b^-1*n^-1*(1+-1*b*d*n)*(a+b*log(c*x^n))))
@test integrate(x^-3*coshintegral(d*(a+b*log(c*x^n))), x) == :(-1//2*x^-2*coshintegral(d*(a+b*log(c*x^n)))+1//4*x^-2*(c*x^n)^(2*n^-1)*exp(2*a*b^-1*n^-1)*expintegralei(-1*b^-1*n^-1*(2+b*d*n)*(a+b*log(c*x^n)))+1//4*x^-2*(c*x^n)^(2*n^-1)*exp(2*a*b^-1*n^-1)*expintegralei(-1*b^-1*n^-1*(2+-1*b*d*n)*(a+b*log(c*x^n))))
@test integrate((e*x)^m*coshintegral(d*(a+b*log(c*x^n))), x) == :(e^-1*(e*x)^(1+m)*(1+m)^-1*coshintegral(d*(a+b*log(c*x^n)))+-1//2*x*(c*x^n)^(-1*n^-1*(1+m))*(e*x)^m*(1+m)^-1*exp(-1*a*b^-1*n^-1*(1+m))*expintegralei(b^-1*n^-1*(a+b*log(c*x^n))*(1+m+b*d*n))+-1//2*x*(c*x^n)^(-1*n^-1*(1+m))*(e*x)^m*(1+m)^-1*exp(-1*a*b^-1*n^-1*(1+m))*expintegralei(b^-1*n^-1*(a+b*log(c*x^n))*(1+m+-1*b*d*n)))
@test integrate(x^-3*cosh(b*x)*coshintegral(b*x), x) == :(b^2*coshintegral(2*b*x)+-1//4*x^-2*cosh(b*x)^2+1//4*b^2*coshintegral(b*x)^2+-1//2*x^-2*cosh(b*x)*coshintegral(b*x)+-1//4*b*x^-1*sinh(2*b*x)+-1//2*b*x^-1*cosh(b*x)*sinh(b*x)+-1//2*b*x^-1*coshintegral(b*x)*sinh(b*x))
@test integrate(x^-1*cosh(b*x)*coshintegral(b*x), x) == :((1/2)*coshintegral(b*x)^2)
@test integrate(cosh(b*x)*coshintegral(b*x), x) == :(-1//2*b^-1*sinhintegral(2*b*x)+b^-1*coshintegral(b*x)*sinh(b*x))
@test integrate(x*cosh(b*x)*coshintegral(b*x), x) == :((1/2)*b^-2*coshintegral(2*b*x)+(1/2)*b^-2*log(x)+-1//2*b^-2*sinh(b*x)^2+-1*b^-2*cosh(b*x)*coshintegral(b*x)+x*b^-1*coshintegral(b*x)*sinh(b*x))
@test integrate(x^2*cosh(b*x)*coshintegral(b*x), x) == :(-1*b^-3*sinhintegral(2*b*x)+3//4*x*b^-2+2*b^-3*coshintegral(b*x)*sinh(b*x)+-1//2*x*b^-2*sinh(b*x)^2+5//4*b^-3*cosh(b*x)*sinh(b*x)+b^-1*x^2*coshintegral(b*x)*sinh(b*x)+-2*x*b^-2*cosh(b*x)*coshintegral(b*x))
@test integrate(x^3*cosh(b*x)*coshintegral(b*x), x) == :((1/2)*b^-2*x^2+3*b^-4*coshintegral(2*b*x)+3*b^-4*log(x)+-13//4*b^-4*sinh(b*x)^2+-3//4*b^-4*cosh(b*x)^2+-6*b^-4*cosh(b*x)*coshintegral(b*x)+-1//2*b^-2*x^2*sinh(b*x)^2+b^-1*x^3*coshintegral(b*x)*sinh(b*x)+-3*b^-2*x^2*cosh(b*x)*coshintegral(b*x)+2*x*b^-3*cosh(b*x)*sinh(b*x)+6*x*b^-3*coshintegral(b*x)*sinh(b*x))
@test integrate(x^-2*coshintegral(b*x)*sinh(b*x), x) == :(b*coshintegral(2*b*x)+(1/2)*b*coshintegral(b*x)^2+-1//2*x^-1*sinh(2*b*x)+-1*x^-1*coshintegral(b*x)*sinh(b*x))
@test integrate(coshintegral(b*x)*sinh(b*x), x) == :(-1//2*b^-1*coshintegral(2*b*x)+-1//2*b^-1*log(x)+b^-1*cosh(b*x)*coshintegral(b*x))
@test integrate(x*coshintegral(b*x)*sinh(b*x), x) == :((1/2)*b^-2*sinhintegral(2*b*x)+-1//2*x*b^-1+-1*b^-2*coshintegral(b*x)*sinh(b*x)+-1//2*b^-2*cosh(b*x)*sinh(b*x)+x*b^-1*cosh(b*x)*coshintegral(b*x))
@test integrate(x^2*coshintegral(b*x)*sinh(b*x), x) == :(b^-3*sinh(b*x)^2+-1*b^-3*coshintegral(2*b*x)+-1*b^-3*log(x)+-1//4*b^-1*x^2+1//4*b^-3*cosh(b*x)^2+2*b^-3*cosh(b*x)*coshintegral(b*x)+b^-1*x^2*cosh(b*x)*coshintegral(b*x)+-2*x*b^-2*coshintegral(b*x)*sinh(b*x)+-1//2*x*b^-2*cosh(b*x)*sinh(b*x))
@test integrate(x^3*coshintegral(b*x)*sinh(b*x), x) == :(3*b^-4*sinhintegral(2*b*x)+-5//2*x*b^-3+-1//6*b^-1*x^3+(1/2)*x*b^-3*cosh(b*x)^2+-6*b^-4*coshintegral(b*x)*sinh(b*x)+-4*b^-4*cosh(b*x)*sinh(b*x)+3//2*x*b^-3*sinh(b*x)^2+b^-1*x^3*cosh(b*x)*coshintegral(b*x)+-3*b^-2*x^2*coshintegral(b*x)*sinh(b*x)+6*x*b^-3*cosh(b*x)*coshintegral(b*x)+-1//2*b^-2*x^2*cosh(b*x)*sinh(b*x))
@test integrate(coshintegral(2x)*sinh(5x), x) == :(-1//10*coshintegral(3x)+-1//10*coshintegral(7x)+1//5*cosh(5x)*coshintegral(2x))
@test integrate(cosh(5x)*coshintegral(2x), x) == :(-1//10*sinhintegral(3x)+-1//10*sinhintegral(7x)+1//5*coshintegral(2x)*sinh(5x))
@test integrate(x^2*coshintegral(a+b*x)*sinh(a+b*x), x) == :((1/2)*b^-3*cosh(2a+2*b*x)+-1*b^-3*coshintegral(2a+2*b*x)+-1*b^-3*log(a+b*x)+-1//4*b^-1*x^2+1//4*b^-3*cosh(a+b*x)^2+(1/2)*a*x*b^-2+-1*a*b^-3*sinhintegral(2a+2*b*x)+2*b^-3*cosh(a+b*x)*coshintegral(a+b*x)+-1//2*a^2*b^-3*coshintegral(2a+2*b*x)+-1//2*a^2*b^-3*log(a+b*x)+b^-1*x^2*cosh(a+b*x)*coshintegral(a+b*x)+(1/2)*a*b^-3*cosh(a+b*x)*sinh(a+b*x)+-2*x*b^-2*coshintegral(a+b*x)*sinh(a+b*x)+-1//2*x*b^-2*cosh(a+b*x)*sinh(a+b*x))
@test integrate(x*coshintegral(a+b*x)*sinh(a+b*x), x) == :((1/2)*b^-2*sinhintegral(2a+2*b*x)+-1//2*x*b^-1+(1/2)*a*b^-2*coshintegral(2a+2*b*x)+(1/2)*a*b^-2*log(a+b*x)+-1*b^-2*coshintegral(a+b*x)*sinh(a+b*x)+-1//2*b^-2*cosh(a+b*x)*sinh(a+b*x)+x*b^-1*cosh(a+b*x)*coshintegral(a+b*x))
@test integrate(coshintegral(a+b*x)*sinh(a+b*x), x) == :(-1//2*b^-1*coshintegral(2a+2*b*x)+-1//2*b^-1*log(a+b*x)+b^-1*cosh(a+b*x)*coshintegral(a+b*x))
@test integrate(x^2*cosh(a+b*x)*coshintegral(a+b*x), x) == :(x*b^-2+-1*b^-3*sinhintegral(2a+2*b*x)+1//8*b^-3*sinh(2a+2*b*x)+b^-3*cosh(a+b*x)*sinh(a+b*x)+-1*a*b^-3*coshintegral(2a+2*b*x)+-1*a*b^-3*log(a+b*x)+2*b^-3*coshintegral(a+b*x)*sinh(a+b*x)+-1//2*a^2*b^-3*sinhintegral(2a+2*b*x)+-1//4*x*b^-2*cosh(2a+2*b*x)+1//4*a*b^-3*cosh(2a+2*b*x)+b^-1*x^2*coshintegral(a+b*x)*sinh(a+b*x)+-2*x*b^-2*cosh(a+b*x)*coshintegral(a+b*x))
@test integrate(x*cosh(a+b*x)*coshintegral(a+b*x), x) == :((1/2)*b^-2*coshintegral(2a+2*b*x)+(1/2)*b^-2*log(a+b*x)+-1//4*b^-2*cosh(2a+2*b*x)+(1/2)*a*b^-2*sinhintegral(2a+2*b*x)+-1*b^-2*cosh(a+b*x)*coshintegral(a+b*x)+x*b^-1*coshintegral(a+b*x)*sinh(a+b*x))
@test integrate(cosh(a+b*x)*coshintegral(a+b*x), x) == :(-1//2*b^-1*sinhintegral(2a+2*b*x)+b^-1*coshintegral(a+b*x)*sinh(a+b*x))
@test integrate(x*coshintegral(c+d*x)*sinh(a+b*x), x) == :((1/2)*b^-2*cosh(a+-1*b*c*d^-1)*sinhintegral(x*(b+d)+c*d^-1*(b+d))+(1/2)*b^-2*cosh(a+-1*b*c*d^-1)*sinhintegral(x*(b+-1d)+c*d^-1*(b+-1d))+(1/2)*b^-2*coshintegral(x*(b+d)+c*d^-1*(b+d))*sinh(a+-1*b*c*d^-1)+(1/2)*b^-2*coshintegral(x*(b+-1d)+c*d^-1*(b+-1d))*sinh(a+-1*b*c*d^-1)+-1*b^-2*coshintegral(c+d*x)*sinh(a+b*x)+-1//2*b^-1*(b+d)^-1*sinh(a+c+x*(b+d))+-1//2*b^-1*(b+-1d)^-1*sinh(a+-1c+x*(b+-1d))+x*b^-1*cosh(a+b*x)*coshintegral(c+d*x)+(1/2)*c*b^-1*d^-1*cosh(a+-1*b*c*d^-1)*coshintegral(x*(b+d)+c*d^-1*(b+d))+(1/2)*c*b^-1*d^-1*cosh(a+-1*b*c*d^-1)*coshintegral(x*(b+-1d)+c*d^-1*(b+-1d))+(1/2)*c*b^-1*d^-1*sinh(a+-1*b*c*d^-1)*sinhintegral(x*(b+d)+c*d^-1*(b+d))+(1/2)*c*b^-1*d^-1*sinh(a+-1*b*c*d^-1)*sinhintegral(x*(b+-1d)+c*d^-1*(b+-1d)))
@test integrate(coshintegral(c+d*x)*sinh(a+b*x), x) == :(b^-1*cosh(a+b*x)*coshintegral(c+d*x)+-1//2*b^-1*cosh(a+-1*b*c*d^-1)*coshintegral(x*(b+d)+c*d^-1*(b+d))+-1//2*b^-1*cosh(a+-1*b*c*d^-1)*coshintegral(x*(b+-1d)+c*d^-1*(b+-1d))+-1//2*b^-1*sinh(a+-1*b*c*d^-1)*sinhintegral(x*(b+d)+c*d^-1*(b+d))+-1//2*b^-1*sinh(a+-1*b*c*d^-1)*sinhintegral(x*(b+-1d)+c*d^-1*(b+-1d)))
@test integrate(x*cosh(a+b*x)*coshintegral(c+d*x), x) == :((1/2)*b^-2*cosh(a+-1*b*c*d^-1)*coshintegral(x*(b+d)+c*d^-1*(b+d))+(1/2)*b^-2*cosh(a+-1*b*c*d^-1)*coshintegral(x*(b+-1d)+c*d^-1*(b+-1d))+(1/2)*b^-2*sinh(a+-1*b*c*d^-1)*sinhintegral(x*(b+d)+c*d^-1*(b+d))+(1/2)*b^-2*sinh(a+-1*b*c*d^-1)*sinhintegral(x*(b+-1d)+c*d^-1*(b+-1d))+-1*b^-2*cosh(a+b*x)*coshintegral(c+d*x)+-1//2*b^-1*(b+d)^-1*cosh(a+c+x*(b+d))+-1//2*b^-1*(b+-1d)^-1*cosh(a+-1c+x*(b+-1d))+x*b^-1*coshintegral(c+d*x)*sinh(a+b*x)+(1/2)*c*b^-1*d^-1*cosh(a+-1*b*c*d^-1)*sinhintegral(x*(b+d)+c*d^-1*(b+d))+(1/2)*c*b^-1*d^-1*cosh(a+-1*b*c*d^-1)*sinhintegral(x*(b+-1d)+c*d^-1*(b+-1d))+(1/2)*c*b^-1*d^-1*coshintegral(x*(b+d)+c*d^-1*(b+d))*sinh(a+-1*b*c*d^-1)+(1/2)*c*b^-1*d^-1*coshintegral(x*(b+-1d)+c*d^-1*(b+-1d))*sinh(a+-1*b*c*d^-1))
@test integrate(cosh(a+b*x)*coshintegral(c+d*x), x) == :(b^-1*coshintegral(c+d*x)*sinh(a+b*x)+-1//2*b^-1*cosh(a+-1*b*c*d^-1)*sinhintegral(x*(b+d)+c*d^-1*(b+d))+-1//2*b^-1*cosh(a+-1*b*c*d^-1)*sinhintegral(x*(b+-1d)+c*d^-1*(b+-1d))+-1//2*b^-1*coshintegral(x*(b+d)+c*d^-1*(b+d))*sinh(a+-1*b*c*d^-1)+-1//2*b^-1*coshintegral(x*(b+-1d)+c*d^-1*(b+-1d))*sinh(a+-1*b*c*d^-1))

max_stars_repo_path: test/8 Special functions/8.5 Hyperbolic integral functions.jl


text:
module LambdaFn
export @λ, @lf
using MacroTools: prewalk, @capture
using Printf

@enum BodyType Implicit Numbered Named
bodytype(sym) =
    sym == "_" ? Implicit :
    occursin(r"^_[0-9]+",sym) ? Numbered :
    occursin(r"^_[[:alnum:]_]+",sym) ? Named :
    missing

btype_str(bt,arg) = bt == Implicit ? "_" : string(arg)

macro λ(body)
    args = Vector{Symbol}()
    found_symbols = Dict{Symbol,Symbol}()
    btype::Union{Missing,BodyType} = missing

    body = prewalk(body) do expr
        if @capture(expr,@λ(body__))
            macroexpand(__module__,expr,recursive=true)
        elseif expr isa Symbol
            sym = string(expr)
            newbtype = bodytype(sym)
            if !ismissing(newbtype)
                if ismissing(btype)
                    btype = newbtype
                elseif btype != newbtype
                    error("Anonymous lambda function can have all implicit ",
                          "(_), all numbered (_1) or all named (_a) arguments.",
                          "You cannot mix these types. But we found two types",
                          " of symbols: `$(sym)` and ",
                          "`$(btype_str(btype,args[end]))`.")
                end
                if btype == Implicit
                    expr = gensym(string("_",length(args)+1))
                    push!(args,expr)
                elseif expr ∉ keys(found_symbols)
                    gens = gensym(expr)
                    found_symbols[expr] = gens

                    expr = gens
                    push!(args,gens)
                else
                    expr = found_symbols[expr]
                end

                expr
            else
                expr
            end
        else
            expr
        end
    end

    if ismissing(btype)
        quote
            function lambda()
                $(esc(body))
            end
        end
    else
        if btype == Numbered
            symbols = Dict((parse(Int,match(r"_([0-9]+)",string(x))[1]) => x
                            for x in args)...)
            if minimum(keys(symbols)) < 1
                error("Numbered arguments in lambda function must be ≥ 1.")
            end

            args = map(1:maximum(keys(symbols))) do num
                get(symbols,num,gensym(string("_",num)))
            end
        end

        quote
            function lambda($(map(esc,args)...))
                $(esc(body))
            end
        end
    end
end
@eval const $(Symbol("@lf")) = $(Symbol("@λ"))

end # module

ext: jl


text:
using Pkg

for p in ("ArgParse", "Knet", "AutoGrad", "Gym")
    if !haskey(Pkg.installed(),p)
        Pkg.add(p)
        if p == "Gym"
            ENV["GYM_ENVS"] = "atari:algorithmic:box2d:classic_control"
            Pkg.build("Gym")
        end
    end
end

"""
julia reinforce_continous.jl

This example implements the REINFORCE algorithm from
`Simple statistical gradient-following algorithms for
connectionist reinforcement learning.`,  Williams, Ronald J.
Machine learning, 8(3-4):229–256, 1992. This example also
demonstrates the usage of the `@zerograd` function for
stopping the gradient flow.
"""
module REINFORCE_CONTINUOUS

using Gym, ArgParse, Knet, AutoGrad, Statistics

function lrelu(x, leak=0.2)
    f1 = Float32(0.5 * (1 + leak))
    f2 = Float32(0.5 * (1 - leak))
    return f1 .* x .+ f2 .* abs.(x)
end

function predict_μ(w, ob)
    hidden = lrelu(w["w1"] * ob .+ w["b1"])
    linear = w["w2"] * hidden .+ w["b2"]
    return linear
end

function logpdf(μ, x; σ=Float32(1.0))
    fac = Float32(-log(sqrt(2pi)))
    r = (x.-μ) / σ
    return -r.* r.* Float32(0.5) .- Float32(log(σ)) .+ fac
end

function sample_action(μ; σ=1.0)
    μ = convert(Array{Float32}, μ)
    a = μ .+ randn(size(μ)) .* σ
end

AutoGrad.@zerograd sample_action(mu)

function play(w, ob)
    μ = predict_μ(w, ob)
    action = sample_action(μ)
    return action, μ
end

function play_episode(w, env, o)
    ob = reset!(env)
    rewards = Float32[]
    μs = Any[]
    actions = Array{Float32,1}()
    total = 0

    for t=1:env.spec.max_episode_steps
        ob_inp = convert(o["atype"], reshape(ob, size(ob, 1), 1))
        action, μ = play(w, ob_inp)
        push!(μs, μ)
        ob, reward, done, _ = step!(env, action.-1)

        total += reward[1]

        push!(rewards, reward[1])
        push!(actions, action[1])

        o["render"] && render(env)

        done && break
    end
    return μs, actions, rewards, total
end

function loss(w, env, o; totalr=nothing)
    μs, actions, rewards, total = play_episode(w, env, o)
    totalr[1] = total

    actions = convert(o["atype"], reshape(actions, size(w["w2"],1), length(actions)))
    discounted = discount(rewards; γ=o["gamma"])
    discounted = discounted .- mean(discounted)#mean R as a baseline
    discounted = discounted ./ (std(discounted) + 1e-6)

    discounted = convert(o["atype"], reshape(discounted, 1, size(actions, 2)))
    -sum((logpdf(hcat(μs...), actions) .* discounted)) / size(μs, 2)
end

lossgradient = grad(loss)

function init_params(input, hidden, output, atype)
    w = Dict()

    w["w1"] = xavier(hidden, input)
    w["b1"] = zeros(hidden, 1)
    w["w2"] = xavier(output, hidden)
    w["b2"] = zeros(output, 1)

    for k in keys(w)
        w[k] = convert(atype, w[k])
    end

    return w
end

function discount(rewards; γ=0.9)
    discounted = zeros(Float32, length(rewards), 1)
    discounted[end] = rewards[end]

    for i=(length(rewards)-1):-1:1
        discounted[i] = rewards[i] + γ * discounted[i+1]
    end
    return discounted
end

function train!(w, opts, env, o)
    totalr = [0.0]
    g = lossgradient(w, env, o; totalr=totalr)
    update!(w, g, opts)
    return totalr[1]
end

function main(args=ARGS)
    s = ArgParseSettings()
    s.description="(c) Ozan Arkan Can, 2018. Demonstration of the REINFORCE algorithm on the continuous action space."
    @add_arg_table s begin
        ("--env_id"; default="MountainCarContinuous-v0"; help="environment name")
        ("--episodes"; arg_type=Int; default=100; help="number of episodes")
        ("--gamma"; arg_type=Float64; default=0.9; help="discount factor")
        ("--threshold"; arg_type=Int; default=1000; help="stop the episode even it is not terminal after number of steps exceeds the threshold")
        ("--lr"; arg_type=Float64; default=0.001; help="learning rate")
        ("--render"; help = "render the environment"; action = :store_true)
        ("--hidden"; arg_type=Int; default=120; help="hidden units")
        ("--usegpu"; action=:store_true; help="use GPU or not")
    end

    Knet.seed!(12345)
    isa(args, AbstractString) && (args=split(args))
    if in("--help", args) || in("-h", args)
        ArgParse.show_help(s; exit_when_done=false)
        return
    end

    o = parse_args(args, s)
    o["atype"] = !o["usegpu"] ? Array{Float32} : KnetArray{Float32}

    env = GymEnv(o["env_id"])
    seed!(env, 12345)

    INPUT = env.observation_space.shape[1]
    OUTPUT = env.action_space.shape[1]

    w = init_params(INPUT, o["hidden"], OUTPUT, o["atype"])
    opts = Dict()
    for k in keys(w)
        opts[k] = Rmsprop(lr=o["lr"])
    end

    for i=1:o["episodes"]
        total = train!(w, opts, env, o)
        println("episode $i , total rewards: $total")
    end
end

PROGRAM_FILE=="reinforce_continuous.jl" && main(ARGS)

end


text:
##################################################
# Parsing of the command-line arguments
##################################################
tabargs = ArgParseSettings()
@add_arg_table! tabargs begin
    "--parallel"
    help = "Parallel computation: yes/no"
    arg_type = String
    default = "yes"
    "--eps"
    help = "Disk temperature"
    arg_type = Float64
    default = 0.1
    "--nbK"
    help = "Sampling number for integration"
    arg_type = Int64
    default = 5000
    "--N"
    help = "Number of basis elements"
    arg_type = Int64
    default = 170
    "--m"
    help = "Azimutal wavenumber"
    arg_type = Int64
    default = 2
end
parsed_args = parse_args(tabargs)

max_stars_repo_name: KerwannTEP/LiRA

max_issues_repo_name: KerwannTEP/LiRA

max_forks_repo_name: KerwannTEP/LiRA


text:
using AssociatedLegendrePolynomials
# input lat and long in degrees

function Ylm(l, m, lat, long)
    return λlm(l, m, sind(lat)) * cosd(m * long)
end

# V and W for vector spherical harmonics
function Vlm(l, m, lat)
    return 0.5 * (
        Plm(l, m + 1, sind(lat)) -
        (l + m) * (l - m + 1) * Plm(l, m - 1, sind(lat))
    ) / sqrt(l * (l + 1))
end

function Wlm(l, m, lat)
    return 0.5 * (
        Plm(l - 1, m + 1, sind(lat)) +
        (l + m) * (l + m - 1) * Plm(l - 1, m - 1, sind(lat))
    ) / sqrt(l * (l + 1))
end

function Blm(l, m, lat, long)
    C = 1 / sqrt(2 * pi)
    E = exp((im * m) * deg2rad(long))
    W = Wlm(l, m, lat)
    V = Vlm(l, m, lat)
    return (real(C * im * W * E), real(C * V * E))
end

function Clm(l, m, lat, long)
    C = 1 / sqrt(2 * pi)
    E = exp((im * m) * deg2rad(long))
    W = Wlm(l, m, lat)
    V = Vlm(l, m, lat)
    return (real(-C * V * E), real(C * im * W * E))
end

function VSH(l, m, lat, long)
    C = 1 / sqrt(2 * pi)
    E = exp((im * m) * deg2rad(long))
    W = Wlm(l, m, lat)
    V = Vlm(l, m, lat)
    return (real(C * (im * W - V) * E), real(C * (V + im * W) * E))
end


text:
using JLSponge
using BenchmarkTools
function main()
    bs = ByteStream(3)
    JLSponge.write!(bs, "cb")
    for _ in 1:999999
        JLSponge.write!(bs, "abc")
        JLSponge.read!(bs, 1)
        JLSponge.write!(bs, "bca")
        JLSponge.read!(bs, 1)
        JLSponge.write!(bs, "cab")
        JLSponge.read!(bs, 3) == "abc"
        JLSponge.write!(bs, "bc")
    end
    return nothing
end

@time main()

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
module CDLagrangianFESpacesTests

using Test
using Gridap.Geometry
using Gridap.ReferenceFEs
using Gridap.FESpaces
using Gridap.CellData
using Gridap.TensorValues
using Gridap.Fields

domain =(0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

# order = 1
orders = (2,1)
order = 2*max(orders...)

Ω = Triangulation(model)
dΩ = Measure(Ω,order)

T = VectorValue{2,Float64}

V = FESpace(model,ReferenceFE(lagrangian,T,orders),conformity=CDConformity((DISC,CONT)))
test_single_field_fe_space(V)
u(x) = x
U = TrialFESpace(V,u)
uh = interpolate(u,U)
e = u - uh
el2 = sqrt(sum(∫(e⋅e)*dΩ))
@test el2 < 1.0e-10

reffe = LagrangianRefFE(T,QUAD,2)
V = FESpace(model,ReferenceFE(lagrangian,T,2),conformity=CDConformity((CONT,DISC)))
U = TrialFESpace(V,u)
uh = interpolate(u,U)
e = u - uh
el2 = sqrt(sum(∫(e⋅e)*dΩ))
@test el2 < 1.0e-10

V = FESpace(model,ReferenceFE(lagrangian,T,(2,0)),conformity=CDConformity((CONT,DISC)))
u(x) = VectorValue(x[1],0.0)
U = TrialFESpace(V,u)
uh = interpolate(u,U)
e = u - uh
el2 = sqrt(sum(∫(e⋅e)*dΩ))
@test el2 < 1.0e-10

V = FESpace(model,ReferenceFE(lagrangian,T,(2,0)),conformity=CDConformity((DISC,DISC)))
u(x) = VectorValue(x[1],0.0)
U = TrialFESpace(V,u)
uh = interpolate(u,U)
e = u - uh
el2 = sqrt(sum(∫(e⋅e)*dΩ))
@test el2 < 1.0e-10

# # Skeleton triangulation
# face_own_dofs = get_face_own_dofs(reffe)
# strian = SkeletonTriangulation(model,reffe,face_own_dofs)
# ns = get_normal_vector(strian)
# writevtk(strian,"strian",cellfields=["normal"=>ns])
#
# # Random function for visualization purposes
# model = CartesianDiscreteModel((0,1,0,1),(10,5))
# V = FESpace(model=model,reffe=reffe,conformity=true)
# trian = Triangulation(model)
# vh = FEFunction(V,rand(num_free_dofs(V)))
# writevtk(trian,"trian",nsubcells=20,cellfields=["vh"=>vh])

#using Gridap.Visualization
#
#writevtk(trian,"trian",nsubcells=10,cellfields=["uh"=>uh])

end # module


text:
using UnsteadyFlowSolvers

alphadef = ConstDef(13. *pi/180)

hdef = ConstDef(0.)

lambda = 0.5
k = 0.1

udef = SinDef(1., lambda, k, 0.)

full_kinem = KinemDef(alphadef, hdef, udef)

pvt = 0.25

geometry = "FlatPlate"

lespcrit = [0.25;]

surf = TwoDSurf(geometry, pvt, full_kinem, lespcrit)

curfield = TwoDFlowField()

dtstar = 0.015 #UNSflow.find_tstep(hdef)

t_tot = 2. *pi/udef.k

nsteps = Int(round(t_tot/dtstar))+1

println("nsteps", nsteps)

startflag = 0

writeflag = 1

writeInterval = t_tot/10.

#delvort = delSpalart(500, 12, 1e-5)
delvort = UNSflow.delNone()

mat, surf, curfield = ldvmLin(surf, curfield, nsteps, dtstar,startflag, writeflag, writeInterval, delvort)

UNSflow.makeVortPlots2D()

UNSflow.makeForcePlots2D()

#UNSflow.makeKinemClVortPlots2D()

UNSflow.cleanWrite()

max_stars_repo_name: MaximeUofG/UNSflowTutorials

max_issues_repo_name: MaximeUofG/UNSflowTutorials

max_forks_repo_name: MaximeUofG/UNSflowTutorials


text:
"""
    LinearTransform(A::AbstractMatrix)

Linear transformation of the input realised by the matrix `A`.

The second dimension of `A` must match the number of features of the target.

# Examples

```jldoctest
julia> A = rand(10, 5); t = LinearTransform(A); X = rand(5, 100);

julia> map(t, ColVecs(X)) == ColVecs(A * X)
true
```
"""
struct LinearTransform{T<:AbstractMatrix{<:Real}} <: Transform
    A::T
end

@functor LinearTransform

function set!(t::LinearTransform{<:AbstractMatrix{T}}, A::AbstractMatrix{T}) where {T<:Real}
    size(t.A) == size(A) || error(
        "size of the given matrix ",
        size(A),
        " and of the transformation matrix ",
        size(t.A),
        " are not the same",
    )
    return t.A .= A
end

(t::LinearTransform)(x::Real) = vec(t.A * x)
(t::LinearTransform)(x::AbstractVector{<:Real}) = t.A * x

_map(t::LinearTransform, x::AbstractVector{<:Real}) = ColVecs(t.A * collect(x'))
_map(t::LinearTransform, x::ColVecs) = ColVecs(t.A * x.X)
_map(t::LinearTransform, x::RowVecs) = RowVecs(x.X * t.A')

function Base.show(io::IO, t::LinearTransform)
    return print(io::IO, "Linear transform (size(A) = ", size(t.A), ")")
end


text:
function get_default_law(::Type{Propagator}, arg0::Frames)
    return jcall(Propagator, "getDefaultLaw", AttitudeProvider, (Frames,), arg0)
end


text:
using ArgParse
function parse_commandline()
  s = ArgParseSettings()

  @add_arg_table! s begin
    "--spreadNH_start"
      help = "spreadNH_start"
      arg_type = Float64
      default = 3.0
    "--spreadNH_stop"
      help = "spreadNH_stop"
      arg_type = Float64
      default = 3.0
    "--spreadNH_step"
      help = "spreadNH_step"
      arg_type = Float64
      default = 1.0
    "--inflation_start"
      help = "inflation_start"
      arg_type = Float64
      default = 5.0
    "--inflation_stop"
      help = "inflation_stop"
      arg_type = Float64
      default = 5.0
    "--inflation_step"
      help = "inflation_step"
      arg_type = Float64
      default = 1.0
    "--useMsgLikelihoods"
      help = "useMsgLikelihoods"
      arg_type = Bool
      default = true
    "--anneal"
      help = "anneal the solve value according to param settings iterations"
      action = :store_true
    "fgtargz"
      help = "fgtargz, but leave out the extension .tar.gz"
      required = true
  end

  return parse_args(s)
end

using RoME
using Distributed
addprocs(10)
using RoME
@everywhere using RoME

function main()
  pargs = parse_commandline()
  println("Parsed args:")
  for (arg,val) in pargs
      println("  $arg  =>  $val")
  end

  fname	= replace(pargs["fgtargz"],".tar.gz"=>"")
  fg = loadDFG(fname);

  getSolverParams(fg).useMsgLikelihoods = pargs["useMsgLikelihoods"]

  SNH = pargs["spreadNH_start"]:pargs["spreadNH_step"]:pargs["spreadNH_stop"]
  IFL = pargs["inflation_start"]:pargs["inflation_step"]:pargs["inflation_stop"]

  if pargs["anneal"]
    for i in 1:length(SNH)
      getSolverParams(fg).spreadNH = SNH[i]
      getSolverParams(fg).inflation = IFL[i]
      # solve
      solveTree!(fg, storeOld=true);
    end
    # store the final result
    snhn = replace(string(SNH),'.'=>'_')
    ifln = replace(string(IFL),'.'=>'_')
    saveDFG(pargs["fgtargz"]*"S$snhn"*"I$ifln", fg)
  else
    for snh in SNH, ifl in IFL
      getSolverParams(fg).spreadNH = snh
      getSolverParams(fg).inflation = ifl
      fg_ = deepcopy(fg)
      solveTree!(fg_, storeOld=true);
      # save the result
      snhn = replace("$snh",'.'=>'_')
      ifln = replace("$ifl",'.'=>'_')
      saveDFG(pargs["fgtargz"]*"S$snhn"*"I$ifln", fg_)
    end
  end

  @info "done"
end

main()

#


text:
# This file was generated, do not modify it. # hide
RandomForestRegressor = @load RandomForestRegressor pkg=DecisionTree
rfr = RandomForestRegressor()

rfr_m = machine(rfr, X, y);

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
#unique!(push!(LOAD_PATH, "~/Documents/UGent/PhD/Code_cMPS/julia1.0/CMPSKit.jl/src/"))
unique!(push!(LOAD_PATH, joinpath(pwd(), "src")))

using Revise
using CMPSKit
using KrylovKit
using OptimKit
using LinearAlgebra
using JLD2
using TensorOperations
using Plots

D = 4
k = 1.
μ = 2.3
c = 1.5

leftgauge = false
if leftgauge
    #Tensors product ansatz
    Id = 1*Matrix(I,D,D)
    KL = Constant(randn(D^2,D^2))
    KL = 0.5*(KL-KL')
    R1 = Constant(kron(Id,randn(D,D)))
    R2 = Constant(kron(randn(D,D),Id))
    RLs = (R1,R2)
    QL = KL
    for R in RLs
        mul!(QL, R', R, -1/2, 1)
    end
    #Put them in cMPS form
    Ψ = InfiniteCMPS(QL, (R1,R2); gauge = :left)
else
    #Tensors product ansatz
    Id = 1*Matrix(I,D,D)
    Q = Constant(randn(D^2,D^2))
    R1 = Constant(kron(Id,randn(D,D)))
    R2 = Constant(kron(randn(D,D),Id))
    #Put them in cMPS form
    Ψ = InfiniteCMPS(Q, (R1,R2))
end

alg1 = LBFGS(; verbosity = 20, maxiter = 1000000, gradtol = 1e-3);
linalg = GMRES(krylovdim = 100; tol = 1e-5)
h = k * (∂ψ[1]'*∂ψ[1] + ∂ψ[2]'*∂ψ[2]) - μ * (ψ[1]'*ψ[1] + ψ[2]'*ψ[2]) + c * ((ψ[1]')^2*ψ[1]^2 + (ψ[2]')^2*ψ[2]^2)
H = ∫(h, (-Inf,+Inf))

if leftgauge
    Ψ, ρR, E, e, normgrad, numfg, history = groundstate4_unconstrained(H, Ψ; optalg = alg1, linalg = linalg)
    #αs,fs, dfs1, dfs2 = groundstate4(H, Ψ; optalg = alg1, linalg = GMRES(; tol = 1e-5))
    #αs = (αs[1:end-1] + αs[2:end])/2
    #push!(αs,0.1)
    #display(plot(αs,[dfs1,dfs2]))
    #gui()
else
    Ψ, ρL, ρR, E, e, normgrad, numfg, history = groundstate3(H, Ψ; optalg = alg1, linalg = linalg)
    #αs,fs, dfs1, dfs2 = groundstate3(H, Ψ; optalg = alg1, linalg = GMRES(; tol = 1e-5))
    #αs = (αs[1:end-1] + αs[2:end])/2
    #push!(αs,0.1)
end

@show Q = Ψ.Q
@show R1 = Ψ.Rs[1]
@show R2 = Ψ.Rs[2]
@show R1[]*R2[] - R2[]*R1[]
@show expval(ψ[1]*ψ[2] - ψ[2]*ψ[1],Ψ)[]


text:
# utility/z-test.jl

using Test

@test downsample(:test)
@test map_many(:test)
@test mask_test()
@test ndgrid(:test)
@test interp1(:test)


text:
import Term: Panel, RenderableText

println("Panel no content")
@time Panel()
@time Panel(; fit = true)
@time Panel(; width = 22, height = 5)

println("short text panels")
txt = "This is my first panel!"
@time Panel(txt)
@time Panel(txt; justify = :right)
@time Panel(txt; fit = true)

println("long text panels")
txt = ","^500
@time Panel(txt)
@time Panel(txt; justify = :right)
@time Panel(txt; fit = true)

println("nested panels")
p1 = Panel(; width = 10, height = 5, style = "red")
@time Panel(p1)
@time Panel(p1; fit = true)

# TODO visually inspect panel with different texts/renderables/padding combinations
# TODO improve Panel docstring and functions order


text:
"""
	update(phd::PHDFilter, b0::GaussianMixture, Z::Vector{AbstractVector},
		T::Real, U::Real, J_max::Integer)

Perform an update step using a GMPHD Filter.

Arguments:
- `phd::PHDFilter` PHD filter to step through. [PHD Filter]
- `b0::GaussianMixture` Prior state distribution. [Gaussian Mixture]
- `Z::Matrix` Matrix of measurements. [Measurements]
- `T::Real` Truncation threshold.  Drop distributions with weight less than T
- `U::Real` Merging threshold.  Merge if (μ_1-μ_2)^T * Σ^-1 * (μ_1-μ_2) < U
- `J_max::Integer` Maximum number of features

Returns:
- `bn::GaussianMixture` Describe first return value. [Gaussian Mixture]

References:
1. Vo, B. N., & Ma, W. K. (2006). The Gaussian mixture probability hypothesis
density filter. IEEE Transactions on signal processing, 54(11), 4091-4104.
"""
function update(phd::PHDFilter, b0::GaussianMixture,
	Z::Vector{<:AbstractVector{<:Number}}, T::Real, U::Real, J_max::Integer)

	bp = predict(phd, b0)
	bm = measure(phd, bp, Z)
    bn = prune(bm, T, U, J_max)
    return bn
end

"""
	predict(phd::PHDFilter, b0::GaussianMixture)

Make a prediction on next state based on PHD dynamics

Arguments:
- `phd::PHDFilter` PHD filter to step through. [PHD Filter]
- `b0::GaussianMixture` Prior state distribution. [Gaussian Mixture]

Returns:
- `bp::GaussianMixture` Predicted next state distribution. [Gaussian Mixture]

References:
1. Vo, B. N., & Ma, W. K. (2006). The Gaussian mixture probability hypothesis
density filter. IEEE Transactions on signal processing, 54(11), 4091-4104.
"""
function predict(phd::PHDFilter, b0::GaussianMixture)

    J_k = b0.N
    n = length(b0.μ[1])
    J_γ = phd.γ.N
    J_β = phd.spawn.β.N
    J_dyn = length(phd.dyn)
    tot = J_γ + J_β*J_k + J_k*J_dyn
    w_p = [0.0 for i in 1:tot]
    μ_p = [zeros(n) for i in 1:tot]
    Σ_p = [zeros(n,n) for i in 1:tot]

    # 1. Prediction for birth targets
    i = 0
    for j in 1:J_γ
        i += 1
        w_p[i] = phd.γ.w[j]
        μ_p[i] = phd.γ.μ[j]
        Σ_p[i] = phd.γ.Σ[j]
    end

    for j in 1:J_β, ℓ in 1:J_k
        i += 1
        spawn_dyn = phd.spawn.dyn[j]
        w_p[i] = b0.w[ℓ]*phd.spawn.β.w[j]
        μ_p[i] =  spawn_dyn.d + spawn_dyn.A*b0.μ[ℓ]
        Σ_p[i] = spawn_dyn.Q + spawn_dyn.A*b0.Σ[ℓ]*spawn_dyn.A'
    end

    # 2. Prediction for existing targets
    for j in 1:J_k
        for k in 1:J_dyn
            i += 1
            w_p[i] = phd.Ps*b0.w[j]
            μ_p[i] = phd.dyn[k].A*b0.μ[j]
            Σ_p[i] = phd.dyn[k].Q + phd.dyn[k].A*b0.Σ[j]*phd.dyn[k].A'
        end
    end
    @assert i == tot "Poor Indexing"
    bp = GaussianMixture(w_p, μ_p, Σ_p)
	return bp
end

"""
	measure(phd::PHDFilter, bp::GaussianMixture, Z::Vector{AbstractVector})

Perform a measurement update on a predicted PHD next state.

Arguments:
- `phd::PHDFilter` PHD filter to step through. [PHD Filter]
- `bp::GaussianMixture` Prior state distribution. [Gaussian Mixture]
- `Z::Vector{AbstractVector}` Array of measurements. [Measurements]

Returns:
- `bm::GaussianMixture` Describe first return value. [Gaussian Mixture]

References:
1. Vo, B. N., & Ma, W. K. (2006). The Gaussian mixture probability hypothesis
density filter. IEEE Transactions on signal processing, 54(11), 4091-4104.
"""
function measure(phd::PHDFilter, bp::GaussianMixture,
				Z::Vector{<:AbstractVector{<: Real}})

	# 3. Construction of PHD update components

    J_p = bp.N

    η_n = []
    S_n = []
    K_n = []
    P_n = []
    for j in 1:J_p
        push!(η_n, phd.meas.C*bp.μ[j])

        S_add = phd.meas.R + phd.meas.C*bp.Σ[j]*phd.meas.C'
        push!(S_n, S_add)

        K_add = bp.Σ[j]*phd.meas.C'*inv(S_add)
        push!(K_n, K_add)

        P_add = (1.0I - K_add*phd.meas.C)*bp.Σ[j]
        push!(P_n, P_add)
    end

    # 4. Update
    w_n = similar(bp.w,0)
    μ_n = similar(bp.μ,0)
    Σ_n = similar(bp.Σ,0)
    for j in 1:J_p
        push!(w_n,(1-phd.Pd)*bp.w[j])
        push!(μ_n, bp.μ[j])
        push!(Σ_n, bp.Σ[j])
    end

    for z in Z
        w_ntemp = []
        for j in 1:J_p
            push!(w_ntemp, phd.Pd*bp.w[j]*MvNormalPDF(z,η_n[j],S_n[j]))
            push!(μ_n, bp.μ[j] + K_n[j]*(z-η_n[j]))
            push!(Σ_n, P_n[j])
        end
        append!(w_n, w_ntemp ./ (phd.κ(z) + sum(w_ntemp)))
    end

    xm = GaussianMixture(w_n, μ_n, Σ_n)
    return xm
end

"""
    prune(b::GaussianMixture, T::Real, U::Real, J_max::Integer)

Prune the posterior Gaussian-Mixture next PHD state

Arguments:
- `b::GaussianMixture` Posterior state distribution. [Gaussian Mixture]
- `T::Real` Truncation threshold.  Drop distributions with weight less than T
- `U::Real` Merging threshold.  Merge if (μ_1-μ_2)^T * Σ^-1 * (μ_1-μ_2) < U
- `J_max::Integer` Maximum number of features
"""
function prune(b::GaussianMixture, T::Real, U::Real, J_max::Integer)
	l = 0
	J_k = b.N
	n = length(b.μ[1])

	# prune indices based on truncation threshold
	I = collect(1:J_k)
	I = I[b.w.>T]

	# setup output vectors
	w_new = zeros(J_k)
	μ_new = [zeros(n) for i=1:J_k]
	Σ_new = [zeros(n,n) for i=1:J_k]

	# perform pruning
	while length(I) != 0
		l = l+1
		j = I[argmax(b.w[I])]

		# Merge close features
		L = []
		for i in I
			delta = b.μ[i]-b.μ[j]
			if dot(delta, inv(b.Σ[i])*delta) < U
				push!(L,i)
			end
		end

		# determine merged parameters
		w_tilde = sum(b.w[L])
		μ_tilde = 1/w_tilde * sum(b.w[L].*b.μ[L])
		Σ_tilde = zeros(n,n)
		for i in L
			tmp = b.Σ[i] + (μ_tilde-b.μ[i])*(μ_tilde-b.μ[i])'
			Σ_tilde = Σ_tilde + b.w[i]*tmp
		end
		Σ_tilde = 1/w_tilde * Σ_tilde

		# store values
		w_new[l] = w_tilde
		μ_new[l] = μ_tilde
		Σ_new[l] = Σ_tilde

		# prune merged features
		I = setdiff(I,L)
	end

	if l > J_max 	# only keep the J_max with highest weights
		idxs = sort!([1:length(w_new);], by=i->(w_new[i]), rev=true)
		idxs = idxs[1:J_max]
		w_new = w_new[idxs]
		μ_new = μ_new[idxs]
		Σ_new = Σ_new[idxs]
	else		# drop zeros at end
		w_new = w_new[1:l]
		μ_new = μ_new[1:l]
		Σ_new = Σ_new[1:l]
	end

	bn = GaussianMixture(w_new,μ_new,Σ_new)
	return bn
end

### Utilities

"""
	MvNormalPDF(x::AbstractVector, μ::AbstractVector, Σ::AbstractMatrix)

Evaluate a multivariate normal pdf function parameterized by μ
and Σ at x

Arguments:
- `x::AbstractVector` Point to evaluate pdf
- `μ::AbstractVector` Mean of Multivariate Normal Gaussian
- `Σ::AbstractMatrix` Covariance Multivariate Normal Gaussian

Returns:
- `xp::AbstractVector` evaluation of multivariate normal pdf
"""
MvNormalPDF(x::AbstractVector, μ::AbstractVector, Σ::AbstractMatrix) =
	(2*pi)^(-length(μ)/2)*det(Σ)^(-1/2)*exp(-1/2*(x-μ)'*inv(Σ)*(x-μ))


text:
module Tables

import Tables as TablesPkg

import Term: do_by_line, term_theme, fillin, truncate
import MyterialColors: orange

import ..Renderables: AbstractRenderable, RenderableText
import ..Layout: cvstack, hstack, vstack, pad, vLine, vertical_pad
import ..Measures: Measure, width, height
import ..Style: apply_style
import ..Segments: Segment
import ..Tprint: tprintln
using ..Boxes

export Table

include("_tables.jl")

"""
    Table

A Table renderable.

# Examples
```julia

t = 1:3
data = hcat(t, rand(Int8, length(t)))
Table(data)

┌───────────┬───────────┐
│  Column1  │  Column2  │
├───────────┼───────────┤
│     1     │    -95    │
├───────────┼───────────┤
│     2     │    -85    │
├───────────┼───────────┤
│     3     │    115    │
└───────────┴───────────┘
```
"""
mutable struct Table <: AbstractRenderable
    segments::Vector{Segment}
    measure::Measure
end

"""
    Table(
        tb::TablesPkg.AbstractColumns;
        box::Symbol = :SQUARE,
        style::String = "default",
        hpad::Union{Vector, Int}=2,
        vpad::Union{Vector, Int}=0,
        vertical_justify::Symbol=:center,

        show_header::Bool = true,
        header::Union{Nothing, Vector, Tuple}=nothing,
        header_style::Union{String, Vector, Tuple} = "default",
        header_justify::Union{Nothing, Symbol, Vector, Tuple} = nothing,

        columns_style::Union{String, Vector, Tuple} = "default",
        columns_justify::Union{Symbol, Vector, Tuple} = :center,

        footer::Union{Function, Nothing, Vector, Tuple}=nothing,
        footer_style::Union{String, Vector, Tuple} = "default",
        footer_justify::Union{Nothing, Symbol, Vector, Tuple} = :center
    )

Generic constructo for a Table renderable.
"""
function Table(
    tb::TablesPkg.AbstractColumns;
    box::Symbol = :SQUARE,
    style::String = "default",
    hpad::Union{Vector,Int} = 2,
    vpad::Union{Vector,Int} = 0,
    vertical_justify::Symbol = :center,
    show_header::Bool = true,
    header::Union{Nothing,Vector,Tuple} = nothing,
    header_style::Union{String,Vector,Tuple} = "default",
    header_justify::Union{Nothing,Symbol,Vector,Tuple} = nothing,
    columns_style::Union{String,Vector,Tuple} = "default",
    columns_justify::Union{Symbol,Vector,Tuple} = :center,
    columns_widths::Union{Nothing,Int,Vector} = nothing,
    footer::Union{Function,Nothing,Vector,Tuple} = nothing,
    footer_style::Union{String,Vector,Tuple} = "default",
    footer_justify::Union{Nothing,Symbol,Vector,Tuple} = :center,
    compact::Bool = false,
)

    # prepare some variables
    header_justify = isnothing(header_justify) ? columns_justify : header_justify
    box = eval(box)

    # table info
    rows = TablesPkg.rows(tb)
    sch = TablesPkg.schema(rows)
    N_cols = length(sch.names)
    N_rows = length(rows) + 2

    # make sure arguemnts combination is valud
    valid = assert_table_arguments(
        N_cols,
        N_rows,
        show_header,
        header,
        header_style,
        header_justify,
        columns_style,
        columns_justify,
        columns_widths,
        footer,
        footer_style,
        footer_justify,
        hpad,
        vpad,
    )
    valid || return nothing

    # columns style
    columns_style, columns_justify, hpad =
        expand.([columns_style, columns_justify, hpad], N_cols)
    vpad = expand(vpad, N_rows)

    # headers and headers style
    if show_header
        header = isnothing(header) ? string.(sch.names) : header
        header_style, header_justify = expand.([header_style, header_justify], N_cols)
    end

    # get footer (if it's a function)
    if footer isa Function
        try
            footer_entries = footer.(map(c -> tb[c], sch.names))
            footer = string(footer) * ": " .* string.(footer_entries)
        catch
            @warn "Could not apply function $footer to table - types mismatch?"
            footer = repeat(["couldn't apply"], N_cols)
        end
    end

    # get the max-width of each column
    widths = calc_columns_widths(
        N_cols,
        N_rows,
        columns_widths,
        show_header,
        header,
        tb,
        sch,
        footer,
        hpad,
    )

    # get the table values as vectors of strings
    rows_values = []
    for row in rows
        _row::Vector = []
        TablesPkg.eachcolumn(sch, row) do val, i, nm
            push!(_row, val isa AbstractRenderable ? val : string(val))
        end
        push!(rows_values, _row)
    end

    # get the height of each row
    heights = rows_heights(N_rows, show_header, header, rows_values, footer, vpad)
    # @info "sizes" widths heights

    # create table lines
    nrows = length(rows_values)
    lines::Vector{String} = []
    show_header && push!(
        lines,
        table_row(
            make_row_cells(
                header,
                header_style,
                header_justify,
                widths,
                hpad,
                heights[1],
                vertical_justify,
            ),
            widths,
            box,
            :top,
            :head,
            :head_row,
            style,
            heights[1],
        ),
    )

    for (l, row) in enumerate(rows_values)
        I = l + 1
        row = make_row_cells(
            row,
            columns_style,
            columns_justify,
            widths,
            hpad,
            heights[I],
            vertical_justify,
        )
        if l == 1
            bottom = if nrows < 2
                :bottom
            elseif nrows > 2
                :row
            else
                :foot_row
            end
            push!(
                lines,
                table_row(
                    row,
                    widths,
                    box,
                    show_header ? nothing : :top,
                    :mid,
                    bottom,
                    style,
                    heights[I];
                    compact = show_header ? false : compact,
                ),
            )
        elseif l == nrows
            push!(
                lines,
                table_row(
                    row,
                    widths,
                    box,
                    nothing,
                    :mid,
                    isnothing(footer) ? :bottom : :foot_row,
                    style,
                    heights[I],
                ),
            )
        else
            push!(
                lines,
                table_row(
                    row,
                    widths,
                    box,
                    nothing,
                    :mid,
                    :row,
                    style,
                    heights[I];
                    compact = compact,
                ),
            )
        end
    end

    if !isnothing(footer)
        # get footer style
        footer_justify = isnothing(footer_justify) ? columns_justify : footer_justify
        footer_style, footer_justify = expand.([footer_style, footer_justify], N_cols)

        push!(
            lines,
            table_row(
                make_row_cells(
                    footer,
                    footer_style,
                    footer_justify,
                    widths,
                    hpad,
                    heights[end],
                    vertical_justify,
                ),
                widths,
                box,
                nothing,
                :foot,
                :bottom,
                style,
                heights[end],
            ),
        )
    end

    lines = split(fillin(join(lines, "\n")), "\n")

    segments = Segment.(lines)
    return Table(segments, Measure(segments))
end

"""
    Table(data::Union{AbstractVector, AbstractMatrix}; kwargs...)

Construct `Table` from `Vector` and `Matrix`
"""
function Table(data::Union{AbstractVector,AbstractMatrix}; kwargs...)
    return Table(TablesPkg.table(data); kwargs...)
end

"""
Table(data::AbstractDict; kwargs...)

Construct `Table` from a `Dict`.
The Dict's keys make up the table header if none is assigned.
"""
function Table(data::AbstractDict; kwargs...)
    kwargs = Dict(kwargs...)
    header = pop!(kwargs, :header, string.(collect(keys(data))))

    data = hcat(values(data)...)
    return Table(data; header = header, kwargs...)
end

"""
    table_row(cells, widths, box, top_level, mid_level, bottom_level, box_style, row_height)

Create a single row of a `Table` renderable.

Each row is made up of an upper, mid and bottom section. The
upper and bottom sections are created by rows of the `Box` type.
The mid section is made up of the row's cells and dividing lines.
"""
function table_row(
    cells,
    widths,
    box,
    top_level,
    mid_level,
    bottom_level,
    box_style,
    row_height;
    compact = false,
)
    # get box characters
    mid_level = getfield(box, mid_level)

    l = vLine(row_height; style = box_style, char = mid_level.left)
    m = vLine(row_height; style = box_style, char = mid_level.vertical)
    r = vLine(row_height; style = box_style, char = mid_level.right)

    if row_height == 1
        l, m, r = string(l), string(m), string(r)
    end

    # create row
    mid = if length(widths) > 1
        l * foldl((a, b) -> a * m * b, cells[1:end]) * r
    else
        l * cells[1] * r
    end
    bottom = apply_style(get_row(box, widths, bottom_level), box_style)

    return if !isnothing(top_level)
        top = apply_style(get_row(box, widths, top_level), box_style)
        string(vstack(top, mid, bottom))
    else
        compact ? string(mid) : string(vstack(mid, bottom))
    end
end

"""
    cell(x::AbstractString, w::Int, h::Int, justify::Symbol, style::String, vertical_justify::Symbol)

Create a Table's row's cell from a string - apply styling and vertical/horizontal justification.
"""
function cell(
    x::AbstractString,
    w::Int,
    hor_pad::Int,
    h::Int,
    justify::Symbol,
    style::String,
    vertical_justify::Symbol,
)
    return vertical_pad(
        do_by_line(
            y -> apply_style(" " * pad(y, w - 2, justify) * " ", style),
            truncate(x, w - hor_pad),
        ),
        h,
        vertical_justify,
    )
end

"""
    cell(x::AbstractString, w::Int, h::Int, justify::Symbol, style::String, vertical_justify::Symbol)

Create a Table's row's cell from a renderable - apply styling and vertical/horizontal justification.
"""
function cell(
    x::AbstractRenderable,
    w::Int,
    hor_pad::Int,
    h::Int,
    justify::Symbol,
    style::String,
    vertical_justify::Symbol,
)
    return vertical_pad(do_by_line(y -> pad(y, w, justify), string(x)), h, vertical_justify)
end

"""
    make_row_cells(
        entries::Union{Tuple, Vector},
        style::Vector{String},
        justify::Vector{Symbol},
        widths::Vector{Int},
        height::Int,
        vertical_justify::Symbol,
    )

Create a row's cell from a vector of 'entries' (renderables or strings).
"""
function make_row_cells(
    entries::Union{Tuple,Vector},
    style::Vector{String},
    justify::Vector{Symbol},
    widths::Vector{Int},
    hor_pad::Vector{Int},
    height::Int,
    vertical_justify::Symbol,
)
    N = length(entries)
    cells = map(
        i -> cell(
            entries[i],
            widths[i],
            hor_pad[i],
            height,
            justify[i],
            style[i],
            vertical_justify,
        ),
        1:N,
    )
    return cells
end

end


text:
"""
    FreeBridge{T} <: Bridges.Variable.AbstractBridge

Transforms constrained variables in [`MOI.Reals`](@ref) to the difference of
constrained variables in [`MOI.Nonnegatives`](@ref).
"""
struct FreeBridge{T} <: AbstractBridge
    variables::Vector{MOI.VariableIndex}
    constraint::MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}
end
function bridge_constrained_variable(::Type{FreeBridge{T}},
                                     model::MOI.ModelLike,
                                     set::MOI.Reals) where T
    variables, constraint = MOI.add_constrained_variables(
        model, MOI.Nonnegatives(2MOI.dimension(set)))
    return FreeBridge{T}(variables, constraint)
end

function supports_constrained_variable(
    ::Type{<:FreeBridge}, ::Type{MOI.Reals})
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:FreeBridge})
    return [(MOI.Nonnegatives,)]
end
function MOIB.added_constraint_types(::Type{FreeBridge{T}}) where T
    return Tuple{DataType, DataType}[]
end

# Attributes, Bridge acting as a model
function MOI.get(bridge::FreeBridge, ::MOI.NumberOfVariables)
    return length(bridge.variables)
end
function MOI.get(bridge::FreeBridge, ::MOI.ListOfVariableIndices)
    return vcat(bridge.variables)
end
function MOI.get(bridge::FreeBridge,
                 ::MOI.NumberOfConstraints{MOI.VectorOfVariables,
                                           MOI.Nonnegatives})
    return 1
end
function MOI.get(bridge::FreeBridge,
                 ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,
                                               MOI.Nonnegatives})
    return [bridge.constraint]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::FreeBridge)
    MOI.delete(model, bridge.variables)
end

function MOI.delete(model::MOI.ModelLike, bridge::FreeBridge, i::IndexInVector)
    n = div(length(bridge.variables), 2)
    MOI.delete(model, bridge.variables[i.value])
    MOI.delete(model, bridge.variables[n + i.value])
    deleteat!(bridge.variables, [i.value, n + i.value])
end

# Attributes, Bridge acting as a constraint

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal,
                 bridge::FreeBridge{T}) where T
    n = div(length(bridge.variables), 2)
    primal = MOI.get(model, attr, bridge.constraint)
    return primal[1:n] - primal[n .+ (1:n)]
end
# The transformation is x_free = [I -I] * x
# so the transformation of the dual is
# y = [I; -I] * y_free
# that is
# y[1:n] = -y[n .+ (1:n)] = y_free
# We can therefore compute `y_free` from either of them, let's take `y[1:n]`.
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 bridge::FreeBridge{T}) where T
    n = div(length(bridge.variables), 2)
    return MOI.get(model, attr, bridge.constraint)[1:n]
end

function MOI.get(model::MOI.ModelLike,
                 attr::Union{MOI.VariablePrimal, MOI.VariablePrimalStart},
                 bridge::FreeBridge{T}, i::IndexInVector) where T
    n = div(length(bridge.variables), 2)
    return MOI.get(model, attr, bridge.variables[i.value]) -
        MOI.get(model, attr, bridge.variables[n + i.value])
end

function MOIB.bridged_function(bridge::FreeBridge{T}, i::IndexInVector) where T
    n = div(length(bridge.variables), 2)
    return MOIU.operate(-, T, MOI.SingleVariable(bridge.variables[i.value]),
                        MOI.SingleVariable(bridge.variables[n + i.value]))
end
# x_free has been replaced by x[i] - x[n + i].
# To undo it we replace x[i] by x_free and x[n + i] by 0.
function unbridged_map(bridge::FreeBridge{T}, vi::MOI.VariableIndex,
                       i::IndexInVector) where T
    sv = MOI.SingleVariable(vi)
    # `unbridged_map` is required to return a `MOI.ScalarAffineFunction`.
    func = convert(MOI.ScalarAffineFunction{T}, sv)
    n = div(length(bridge.variables), 2)
    return bridge.variables[i.value] => func,
        bridge.variables[n + i.value] => zero(MOI.ScalarAffineFunction{T})
end
function MOI.supports(model::MOI.ModelLike, attr::MOI.VariablePrimalStart,
                      ::Type{<:FreeBridge})
    return MOI.supports(model, attr, MOI.VariableIndex)
end
function MOI.set(model::MOI.ModelLike, attr::MOI.VariablePrimalStart,
                 bridge::FreeBridge, value, i::IndexInVector)
    if value < 0
        nonneg = zero(value)
        nonpos = -value
    else
        nonneg = value
        nonpos = zero(value)
    end
    n = div(length(bridge.variables), 2)
    MOI.set(model, attr, bridge.variables[i.value], nonneg)
    MOI.set(model, attr, bridge.variables[n + i.value], nonpos)
end


text:
"""
    grdview(cmd0::String="", arg1=nothing, arg2=nothing, arg3=nothing; kwargs...)

Reads a 2-D grid and produces a 3-D perspective plot by drawing a mesh, painting a
colored/grayshaded surface made up of polygons, or by scanline conversion of these polygons
to a raster image.

Full option list at [`grdview`]($(GMTdoc)grdview.html)

- $(GMT.opt_J)
- $(GMT.opt_Jz)
- $(GMT.opt_R)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **G** | **drapefile** :: [Type => Str | GMTgrid | a Tuple with 3 GMTgrid types]

    Drape the image in drapefile on top of the relief provided by relief_file.
    ($(GMTdoc)grdview.html#g)
- **I** | **shade** | **shading** | **intensity** :: [Type => Str | GMTgrid]		``Arg = GMTgrid | filename``

    Gives the name of a grid file or GMTgrid with intensities in the (-1,+1) range,
    or a grdgradient shading flags.
    ($(GMTdoc)grdview.html#i)
- **N** | **plane** :: [Type => Str | Int]		``Arg = (level [,fill])``

    Draws a plane at this z-level.
    ($(GMTdoc)grdview.html#n)
- $(GMT.opt_P)
- **Q** | **surftype** | **surf** :: [Type => Str | Int] ``Arg = mesh=Bool, surface=Bool, image=Bool, wterfall=(:rows|cols,[fill])``

    Specify **m** for mesh plot, **s** for surface, **i** for image.
    ($(GMTdoc)grdview.html#q)
- **S** | **smoothfactor** :: [Type => Number]

    Used to resample the contour lines at roughly every (gridbox_size/smoothfactor) interval..
    ($(GMTdoc)grdview.html#s)
- **T** | **tiles** | **no_interp** :: [Type => Str | NT]	``Arg = (skip|skip_nan=Bool, outlines=Bool|pen)``

    Plot image without any interpolation.
    ($(GMTdoc)grdview.html#t)
- **W** | **pens** | **pen** :: [Type => Str]	``Arg = (contour=Bool|pen, mesh=Bool|pen, facade=Bool|pen)``

    Draw contour, mesh or facade. Append pen attributes.
    ($(GMTdoc)grdview.html#w)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_f)
- $(GMT.opt_n)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function grdview(cmd0::String="", arg1=nothing; first=true, kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("grdview", cmd0, arg1)
	arg2 = nothing;	arg3 = nothing;	arg4 = nothing;	arg5 = nothing;

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode
	common_insert_R!(d, O, cmd0, arg1)			# Set -R in 'd' out of grid/images (with coords) if limits was not used

	has_opt_B = (find_in_dict(d, [:B :frame :axis :axes], false)[1] !== nothing)
	cmd, opt_B, _, opt_R = parse_BJR(d, "", "grdview", O, " -JX" * split(def_fig_size, '/')[1] * "/0")
	(!has_opt_B && isa(arg1, GMTimage) && (isimgsize(arg1) || CTRL.limits[1:4] == zeros(4)) && opt_B == def_fig_axes) &&
		(cmd = replace(cmd, opt_B => ""))	# Dont plot axes for plain images if that was not required

	cmd, = parse_common_opts(d, cmd, [:UVXY :c :f :n :p :t :params], first)
	cmd  = add_opt(d, cmd, 'S', [:S :smooth])
	if ((val = find_in_dict(d, [:N :plane])[1]) !== nothing)
		cmd *= " -N" * parse_arg_and_pen(val, "+g", false)
	end
	cmd = add_opt(d, cmd, 'Q', [:Q :surf :surftype],
				  (mesh=("m", add_opt_fill), surface="_s", surf="_s", img=("i",arg2str), image="i", nan_alpha="_c", monochrome="_+m", waterfall=(rows="my", cols="mx", fill=add_opt_fill)))
	cmd = add_opt(d, cmd, 'W', [:W :pens :pen], (contour=("c", add_opt_pen),
	              mesh=("m", add_opt_pen), facade=("f", add_opt_pen)) )
	cmd = add_opt(d, cmd, 'T', [:T :no_interp :tiles], (skip="_+s", skip_nan="_+s", outlines=("+o", add_opt_pen)) )
	(!occursin(" -T", cmd)) ? cmd = parse_JZ(d, cmd)[1] : del_from_dict(d, [:JZ])	# Means, even if we had one, ignore silently
	cmd = add_opt(d, cmd, "%", [:layout :mem_layout], nothing)

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, arg1)		# Find how data was transmitted

	(isa(arg1, Array{<:Real})) && (arg1 = mat2grid(arg1))

	cmd, N_used, arg1, arg2, arg3 = common_get_R_cpt(d, cmd0, cmd, opt_R, got_fname, arg1, arg2, arg3, "grdview")
	cmd, arg1, arg2, arg3, arg4 = common_shade(d, cmd, arg1, arg2, arg3, arg4, "grdview")
	cmd, arg1, arg2, arg3, arg4, arg5 = parse_G_grdview(d, [:G :drapefile], cmd, arg1, arg2, arg3, arg4, arg5)

	_cmd, K = finish_PS_nested(d, ["grdview " * cmd], K)
    return finish_PS_module(d, _cmd, "", K, O, true, arg1, arg2, arg3, arg4, arg5)
end

# ---------------------------------------------------------------------------------------------------
function parse_G_grdview(d::Dict, symbs::Array{<:Symbol}, cmd::String, arg1, arg2, arg3, arg4, arg5)
	(show_kwargs[1]) && return print_kwarg_opts(symbs, "GMTgrid | Tuple | String"), arg1, arg2, arg3, arg4, arg5
	if ((val = find_in_dict(d, [:G :drapefile])[1]) !== nothing)
		if (isa(val, String))				# Uff, simple. Either a file name or a -A type modifier
			cmd *= " -G" * val
		elseif (isa(val, GMTgrid))			# A single drape grid (arg1-3 may be used already)
			cmd, N_used = put_in_slot(cmd, val, 'G', [arg1, arg2, arg3, arg4])
			if     (N_used == 1)  arg1 = val
			elseif (N_used == 2)  arg2 = val
			elseif (N_used == 3)  arg3 = val
			elseif (N_used == 4)  arg4 = val
			end
		elseif (isa(val, Tuple) && length(val) == 3)
			cmd, N_used = put_in_slot(cmd, val[1], 'G', [arg1, arg2, arg3, arg4, arg5])
			cmd *= " -G -G"		# Because the above only set one -G and we need 3
			if     (N_used == 1)  arg1 = val[1];	arg2 = val[2];		arg3 = val[3]
			elseif (N_used == 2)  arg2 = val[1];	arg3 = val[2];		arg4 = val[3]
			elseif (N_used == 3)  arg3 = val[1];	arg4 = val[2];		arg5 = val[3]
			end
		else
			error("Wrong way of setting the drape (G) option.")
		end
	end
	return cmd, arg1, arg2, arg3, arg4, arg5
end

# ---------------------------------------------------------------------------------------------------
grdview!(cmd0::String="", arg1=nothing; first=false, kw...) = grdview(cmd0, arg1; first=first, kw...)
grdview(arg1; first=true, kw...) = grdview("", arg1; first=first, kw...)
grdview!(arg1; first=false, kw...) = grdview("", arg1; first=first, kw...)


text:
n = rand(1: 10)
@test matrixdepot("moler", n) == matrixdepot("moler", n, -1)

M = matrixdepot("moler", n)
for i = 1:n, j = 1:n
    i != j ? (@test_approx_eq M[i,j] min(i,j) - 2) : (@test_approx_eq M[i,i] i)
end

println("'moler' passed test...")


text:
function IsotropicRadiationClassicalConvention(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return IsotropicRadiationClassicalConvention((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function get_radiation_parameters_drivers(obj::IsotropicRadiationClassicalConvention)
    return jcall(obj, "getRadiationParametersDrivers", Vector{ParameterDriver}, ())
end

function radiation_pressure_acceleration(obj::IsotropicRadiationClassicalConvention, arg0::AbsoluteDate, arg1::Frame, arg2::Vector3D, arg3::Rotation, arg4::jdouble, arg5::Vector3D, arg6::Vector{jdouble})
    return jcall(obj, "radiationPressureAcceleration", Vector3D, (AbsoluteDate, Frame, Vector3D, Rotation, jdouble, Vector3D, Vector{jdouble}), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function radiation_pressure_acceleration(obj::IsotropicRadiationClassicalConvention, arg0::FieldAbsoluteDate, arg1::Frame, arg2::FieldVector3D, arg3::FieldRotation, arg4::RealFieldElement, arg5::FieldVector3D, arg6::Vector{RealFieldElement})
    return jcall(obj, "radiationPressureAcceleration", FieldVector3D, (FieldAbsoluteDate, Frame, FieldVector3D, FieldRotation, RealFieldElement, FieldVector3D, Vector{RealFieldElement}), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end


text:
"""
`ReMat` - a scalar random effects matrix

The matrix consists of the grouping factor, `f`, and the transposed dense model
matrix `z`.  The length of `f` must be equal to the number of columns of `z`
"""
abstract AbstractReMat

immutable ReMat <: AbstractReMat
    f::PooledDataVector
    z::Vector
    function ReMat(p::PooledDataVector,z::Vector)
        length(p) == length(z) || throw(DimensionMismatch())
        new(p,z)
    end
end

immutable VectorReMat <: AbstractReMat
    f::PooledDataVector                 # grouping factor
    z::Matrix
    function VectorReMat(p::PooledDataVector,z::Matrix)
        length(p) == size(z,2) || throw(DimensionMismatch())
        size(z,1) > 1 || error("use ReMat instead")
        new(p,z)
    end
end

ReMat(p::PooledDataVector) = ReMat(p,ones(length(p)))

ReMat{T<:Integer}(v::Vector{T}) = ReMat(compact(pool(v)))

Base.eltype(A::AbstractReMat) = eltype(A.z)

vsize(A::AbstractReMat) = isa(A,ReMat) ? 1 : size(A.z,1)

Base.size(A::AbstractReMat) = (length(A.f),vsize(A)*length(A.f.pool))

Base.size(A::AbstractReMat,i::Integer) =
    i < 1 ? throw(BoundsError()) :
    i == 1 ? length(A.f) :
    i == 2 ? vsize(A)*length(A.f.pool) : 1

==(A::AbstractReMat,B::AbstractReMat) = (A.f == B.f) && (A.z == B.z)

# FIXME add a tA boolean argument to combine the code for A_mul_B! and Ac_mul_B!
function Base.A_mul_B!{T}(α::Real,A::AbstractReMat,B::StridedVecOrMat{T},β::Real,R::StridedVecOrMat{T})
    n,q = size(A)
    k = size(B,2)
    size(R,1) == n && size(B,1) == q && size(R,2) == k || throw(DimensionMismatch())
    if β ≠ 1
        scale!(β,R)
    end
    rr = A.f.refs
    zz = A.z
    if isa(A,ReMat)
        for j in 1:k, i in 1:n
            R[i,j] += α * zz[i] * B[rr[i],j]
        end
    else
        l = size(zz,1)
        Bt = reshape(B,(l,div(q,l),k))
        for j in 1:k, i in 1:n
            R[i,j] += α * dot(sub(Bt,:,Int(rr[i]),j),sub(zz,:,i))
        end
    end
    R
end

Base.A_mul_B!{T}(A::AbstractReMat,B::StridedVecOrMat{T},R::StridedVecOrMat{T}) = A_mul_B!(one(T),A,B,zero(T),R)

function Base.Ac_mul_B!{T}(α::Real,A::AbstractReMat,B::StridedVecOrMat{T},β::Real,R::StridedVecOrMat{T})
    n,q = size(A)
    k = size(B,2)
    size(R,1) == q && size(B,1) == n && size(R,2) == k || throw(DimensionMismatch())
    scale!(β,R)
    rr = A.f.refs
    zz = A.z
    if isa(A,ReMat)
        for j in 1:k, i in 1:n
            R[rr[i],j] += α * zz[i] * B[i,j]
        end
    else
        l = size(zz,1)
        rt = reshape(R,(l,div(q,l),k))
        for j in 1:k, i in 1:n
            Base.axpy!(α*B[i,j],sub(zz,:,i),sub(rt,:,Int(rr[i]),j))
        end
    end
    R
end

Base.Ac_mul_B!{T}(R::StridedVecOrMat{T},A::AbstractReMat,B::StridedVecOrMat{T}) = Ac_mul_B!(one(T),A,B,zero(T),R)
function Base.Ac_mul_B(A::AbstractReMat,B::DenseVecOrMat)
    k = size(A,2)
    Ac_mul_B!(Array(eltype(B), isa(B,Vector) ? (k,) : (k, size(B,2))), A, B)
end

function Base.Ac_mul_B(A::ReMat, B::ReMat)
    Az = A.z
    Ar = A.f.refs
    if is(A,B)
        v = zeros(eltype(Az),length(A.f.pool))
        for i in eachindex(Ar)
            v[Ar[i]] += abs2(Az[i])
        end
        return Diagonal(v)
    end
    densify(sparse(convert(Vector{Int32},Ar),convert(Vector{Int32},B.f.refs),Az .* B.z))
end

function Base.Ac_mul_B(A::VectorReMat,B::VectorReMat)
    Az = A.z
    Ar = convert(Vector{Int},A.f.refs)
    if is(A,B)
        l,n = size(Az)
        T = eltype(Az)
        np = length(A.f.pool)
        a = zeros(T,(l,l,np))
        for i in eachindex(Ar)
            Base.LinAlg.BLAS.syr!('L',one(T),sub(Az,:,i),sub(a,:,:,Ar[i]))
        end
        for k in 1:np
            Base.LinAlg.copytri!(sub(a,:,:,k),'L')
        end
        return HBlkDiag(a)
    end
    Bz = B.z
    Br = convert(Vector{Int},B.f.refs)
    (m = length(Ar)) == length(Br) || throw(DimensionMismatch())
    sparse(Ar,Br,[sub(Az,:,i)*sub(Bz,:,i)' for i in 1:m])
end

function Base.Ac_mul_B!{T}(R::DenseVecOrMat{T},A::DenseVecOrMat{T},B::AbstractReMat)
    m = size(A,1)
    n = size(A,2)
    p,q = size(B)
    m == p && size(R,1) == n && size(R,2) == q || throw(DimensionMismatch(""))
    fill!(R,zero(T))
    rr = B.f.refs
    zz = B.z
    if isa(B,ReMat)
        for j in 1:n, i in 1:m
            R[j,rr[i]] += A[i,j] * zz[i]
        end
    else
        l = size(zz,1)
        for j in 1:n, i in 1:m
            Base.LinAlg.axpy!(A[i,j],sub(zz,:,i),sub(R,j,(rr[i]-1)*l + (1:l)))
        end
    end
    R
end

Base.Ac_mul_B(A::DenseVecOrMat,B::AbstractReMat) =
    Ac_mul_B!(Array(eltype(A),(size(A,2),size(B,2))),A,B)


text:
using Documenter
using DeepUncertainty

DocMeta.setdocmeta!(
    DeepUncertainty,
    :DocTestSetup,
    :(using DeepUncertainty);
    recursive = true,
)

makedocs(
    sitename = "DeepUncertainty.jl",
    modules = [DeepUncertainty],
    repo = "https://github.com/aced-differentiate/DeepUncertainty.jl/blob/{commit}{path}#{line}",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        assets = String["assets/flux.css"],
    ),
    pages = [
        "Home" => "index.md",
        "Trainable Distributions" => "initializers.md",
        "MC Dropout" => "layers/mclayers.md",
        "BatchEnsemble" => "layers/batchensemble.md",
        "Bayesian BatchEnsemble" => "layers/bayesianbe.md",
        "Variational Inference" => "layers/variational.md",
    ],
)

deploydocs(
    repo = "github.com/aced-differentiate/DeepUncertainty.jl",
    target = "build",
    branch = "gh-pages",
    devbranch = "main",
    push_preview = true,
)


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

using Markdown
using InteractiveUtils

# ╔═╡ f4b7d3df-213c-46a6-8156-a83fcf6ffb34
begin

	# setup paths -
	_PATH_TO_ROOT = pwd()
	_PATH_TO_SRC = joinpath(_PATH_TO_ROOT, "src")
	_PATH_TO_DATA = joinpath(_PATH_TO_ROOT, "data")
	_PATH_TO_CONFIG = joinpath(_PATH_TO_ROOT, "config")

	# load required packages -
	import Pkg

	# If you do NOT have a Project.toml file -
	# Pkg.activate(mktempdir())
	# Pkg.add([
	# 	Pkg.PackageSpec(name="DataFrames"),
	# 	Pkg.PackageSpec(name="DelimitedFiles"),
	# 	Pkg.PackageSpec(name="CSV"),
	# 	Pkg.PackageSpec(name="Plots"),
	# 	Pkg.PackageSpec(name="TOML"),
	# ])

	# If you have already activated the project -
	Pkg.activate(_PATH_TO_ROOT);
	Pkg.instantiate()
	Pkg.update()

	using DataFrames
	using DelimitedFiles
	using CSV
	using Plots
	using TOML
	using Optim
end

# ╔═╡ 4b710f82-b1cc-11eb-2d14-11b6421a813f
html"""<style>
main {
    max-width: 1200px;
  	width: 70%;
	margin: auto;
}
"""

# ╔═╡ 0df967d7-ca1a-4a19-9c40-7f3a94764504
md""" #### Setup the Julia environment
"""

# ╔═╡ e0a279c1-7677-4c5b-a6e1-ee450b27020e
md""" ### Solution: Enzyme Activity Model

Derive and test an expression for the allosteric regulation of enzyme activity.
Allosteric regulation is a fast mechanism that regulates the catalytic activity of metabolic enzymes.
In this type of regulation, the ability of an enzyme to catalyze a reaction depends upon the concentration of
effector molecules which do not directly participate in the chemical reaction.

We will be developing a model for Phosphofructokinase (PFK) with and without the effector 3$^{\prime}$-5$^{\prime}$-AMP.

"""

# ╔═╡ aa6bba31-0d33-4476-acb8-707128e63a72
md"""##### a) Describe the Model

Let the model take the form $\hat{r}_{j} = r_{j}v\left(...\right)_{j}$
where $\hat{r}_{j}$, the overall rate of the PFK reaction ($\mu$M h$^{-1}$),
is the product of a kinetic limit $r_{j}$ ($\mu$M h$^{-1}$), i.e., the maximum rate of conversion,
and a control variable $0\leq v\left(...\right)_{j}\leq 1$ (dimensionless) that describes the influence
of effector molecules. Since we have only a single enzyme let j = 1.

The model for the kinetic limit was given in the problem:

$r_{1} = k_{cat}E_{1}\left(\frac{F6P}{K_{F6P}+F6P}\right)\left(\frac{ATP}{K_{ATP}+ATP}\right)$

along with all the parameters associated with this model, and the concentration of enzyme ($E_{1}$), F6P and ATP. Thus, we can calculate the kinetic limit. However, we need to postulate a form for the v-varible, and estimate the parameters that appear in the v-variable expression. We'll do this from the data.

###### v-model
In the problem we proposed a three state model for $v\left(...\right)_{j}$: State 0: no effector+no substrate, State 1: no effector+substrate and State 2: effector+substrate. Following the procedure discussed in class, given this three state model, the v-variable takes the form:

$v_{1} = \frac{W_{1}+W_{2}f_{2}}{1+W_{1}+W_{2}f_{2}}$

where $W_{i}$ denotes the weight of state i (dimensionless), and $f_{2}$ denote the fraction of bound effector given by: $f_{2} = (x/K_{2})^{n_{2}}/(1+(x/K_{2})^{n_{2}})$. The quanity x denotes the concentration of the effector 3$^{\prime}$-5$^{\prime}$-AMP, $K_{2}$ denotes a binding constant (units of concentration) and $n_{2}$ denotes a binding cooperativity parameter (dimensionless). State 0 has no activity, while State 1 (with weight $W_{1}$) has basal enzyme activity, and State 2 describes the enhanced activity following from effector binding (with weight $W_{2}f_{2}$).

 - The probability of state 0 is $p_{0}$ = $1/Z$ where $Z$=$1+W_{1}+W_{2}f_{2}$.
 - The probability of state 1 is $p_{1}$ = $W_{1}/Z$ where $Z$=$1+W_{1}+W_{2}f_{2}$.
 - The probability of state 2 is $p_{2}$ = $W_{2}f_{2}/Z$ where $Z$=$1+W_{1}+W_{2}f_{2}$.

"""

# ╔═╡ c4077a00-bef2-4003-ba51-8596550335b2
md""" ##### b) How do we get the parameters in the activity model?

We can estimate the parameters in the activity model from the experimental data given in the problem. This data set measured the rate of PFK at the same concentration(s) of ATP and F6P but different concentrations of 3$^{\prime}$-5$^{\prime}$-AMP. Thus, the kinetic limit was constant in this experiment.
Thus, since the kinetic limit is constant, we can estimate a value of the v-variable as a function of 3$^{\prime}$-5$^{\prime}$-AMP:

$v_{1} \simeq \frac{\hat{r}_{1}^{*}}{r_{1}}$

where $\hat{r}_{1}^{*}$ denotes the _measured value_ of the overall rate.

Lets look at the ends of the experimental data. Let $\alpha$ denote the ratio of rates when there is no effector, and let $\beta$ denote the ratio of rates when there is saturating levels of effector.

###### Case: no effector
When there is no effector, the v-variable model becomes:

$\alpha \simeq \frac{W_{1}}{1+W_{1}}$

which we can solve for $W_{1}$:

$W_{1} \simeq \frac{\alpha}{1-\alpha}$

###### Case: saturating effector

When the effector levels are saturating we expect $f_{2}\simeq{1}$ which gives the expression:

$\beta\simeq\frac{W_{1}+W_{2}}{1+W_{1}+W_{2}}$

which can be solved for $W_{2}$:

$W_{2} \simeq \left(1-\beta\right)^{-1}\Bigl\{\beta\left(1+W_{1}\right)-W_{1}\Bigr\}$

"""

# ╔═╡ 50886df5-d361-4646-898f-795a064c6001
md""" ###### Estimate the configuration weights
"""

# ╔═╡ a3264531-d9c8-42cf-b671-c6f46980a5a4
md""" ###### Estimate the binding and the cooperativity parameters

To estimate the remaining model parameters, we'll use an optimization approach e.g., Solver in xls or the [Optim](https://github.com/JuliaNLSolvers/Optim.jl) package in Julia. In this case, we want to minimize the squared difference between the model predicted v-variable and the v-variable estimated from the experimental data where the decision variables are $K_{2}$ and $n_{2}$:

$\min_{K_{2},n_{2}} \sum_{m=1}^{M}\left(v_{m,1} - v^{*}_{m,1}\right)^2$

where the summation is over the experimental data (the different values of the 3$^{\prime}$-5$^{\prime}$-AMP concentrations), $v_{m,1}$ denotes the model v-variable and $v^{*}_{m,1}$ denotes the experimentally estimated v-variable.

Note: these estimates will __may not be unique__ because the problem is underspecified. Thus, your values may be different.

"""

# ╔═╡ dc395064-55a7-42fe-8fb0-880692e7d81a
md""" ###### c) Plot model versus data. Does it fit?

Now that have the model parameters, let's plot the model v-variable versus the estimated v-variable and see if the model describes the data. In real life, we would do a statistical test on the residual (difference between the model and the data) and/or a cross-validation procedure (partition the data into test and validation sets) and calculate a training and prediction error. However, in this case let's just look - does the model generally describe the trend in the data?

"""

# ╔═╡ b1c395bd-6fa7-4e07-a339-894675791d91
md""" ##### Helper Functions Below Here
"""

# ╔═╡ 169f11a1-34dd-45a9-8f9a-6b18f4957998
function compute_activity_variable(effector_array::Array{Float64,1}, 	 parameter_dictionary::Dict{String,Any})::Array{Float64,1}

    # initialize -
    v_variable_array = Array{Float64,1}()

    try

        # process each effector -
        for effector_value in effector_array

            # compute -
            value = compute_activity_variable(effector_value, parameter_dictionary)

            # cache -
            push!(v_variable_array, value)
        end

        # return -
        return v_variable_array
    catch error
        throw(error)
    end
end

# ╔═╡ de9f574f-fc03-442c-90bd-810f44673466
function compute_activity_variable(effector::Float64, parameter_dictionary::Dict{String,Any})::Float64

    try

        # grab the parameter values from the dictionary -
        W1 = parameter_dictionary["W1"]
        W2 = parameter_dictionary["W2"]
        KD = parameter_dictionary["KD"]
        n = parameter_dictionary["n"]

        # compute the f -
        f = ((effector / KD)^n) / (1 + (effector / KD)^n)

        # compute the v-variable -
        v_variable = (W1 + W2 * f) / (1 + W1 + W2 * f)

        # return -
        return v_variable
    catch error
        throw(error)
    end
end

# ╔═╡ 820c02ef-da91-482d-8bab-36de9b524fa1
function _obj_function(p,parameter_dictionary::Dict{String,Any}, v_variable_data::Array{Float64,1},
    experimental_data_array::Array{Float64,2})

    parameter_dictionary["KD"] = p[1]
    parameter_dictionary["n"] = p[2]

    # compute the model v-variable -
    model_v_variable = compute_activity_variable(experimental_data_array[:,1], parameter_dictionary)

    # compute the residual between the data and model estimated v-variable -
    error = sum((v_variable_data .- model_v_variable).^2)

    # return -
    return error
end

# ╔═╡ b2fe595d-e90e-4041-bba4-8a1558bb83e2
function estimate_parameters(parameter_dictionary::Dict{String,Any}, v_variable_data::Array{Float64,1},
    experimental_data_array::Array{Float64,2})

    try

        # get unkown parameter guess from config -
        KD = parameter_dictionary["KD"]
        n = parameter_dictionary["n"]

        # setup the calculation -
        xinitial = [KD n]
        OF(p) = _obj_function(p,parameter_dictionary, v_variable_data, experimental_data_array)

        # call the optimizer -
        opt_result = optimize(OF,xinitial,BFGS())

        # return -
        return Optim.minimizer(opt_result)
    catch error
        throw(error)
    end
end

# ╔═╡ aa2267d4-ae47-4532-9fa0-460498b40fea
function load_experimental_data_set(path_to_data_file::String)::Array{Float64,2}

    try

        # load the experimental data -
        df = readdlm(path_to_data_file)

        # return -
        return df
    catch error
        throw(error)
    end
end

# ╔═╡ b25aeabb-2be1-48ef-bd72-3680ae8e9407
function load_configuration_dictionary(path_to_configuration_file)::Dict{String,Any}

    try

        # load the config TOML file -
        configuration_dictionary = TOML.parsefile(path_to_configuration_file)

        # return -
        return configuration_dictionary
    catch error
        throw(error)
    end
end

# ╔═╡ 5f846f30-f20f-4153-b87f-71ab429019cb
function compute_kinetic_limit(parameter_dictionary::Dict{String,Any})::Float64

    try

        # grab the parameter values from the dictionary -
        E1 = parameter_dictionary["PFK_concentration"]
        ATP = parameter_dictionary["ATP_concentration"]
        F6P = parameter_dictionary["F6P_concentration"]
        kcat = parameter_dictionary["kcat"]
        K_ATP = parameter_dictionary["K_ATP"]
        K_F6P = parameter_dictionary["K_F6P"]

        # compute the kinetic limit -
        kinetic_limit = kcat * E1 * (ATP / (K_ATP + ATP)) * (F6P / (K_F6P + F6P))

        # return -
        return kinetic_limit
    catch error
        throw(error)
    end
end

# ╔═╡ 5bab3b76-be8a-44c9-8de3-ac35992bda6a
begin

	# load the experimental_data file -
	path_to_experimental_data_file = joinpath(_PATH_TO_DATA, "Data-3-5-AMP.txt")
	experimental_data_array = load_experimental_data_set(path_to_experimental_data_file)

	# load my parameters -> configuration_dictionary contains, two sub-dictionaries
	path_to_configuration_file = joinpath(_PATH_TO_CONFIG, "Configuration.toml")
	configuration_dictionary = load_configuration_dictionary(path_to_configuration_file)
	v_variable_parameter_dictionary = configuration_dictionary["vvariableparameters"]
	kinetic_limit_parameter_dictionary = configuration_dictionary["kineticparameters"]

	# lets compute the kinetic limit -
	kinetic_limit_value = compute_kinetic_limit(kinetic_limit_parameter_dictionary)

	# ok, so let's compute a data driven v-variable -
	data_driven_v_variable = Array{Float64,1}()
	(number_of_effector_values, number_of_cols) = size(experimental_data_array)
	for row_index = 1:number_of_effector_values
    	value = experimental_data_array[row_index,2]
    	scaled_value = value / kinetic_limit_value
    	push!(data_driven_v_variable, scaled_value)
	end

	# return -
	nothing
end

# ╔═╡ 2bb3a988-e054-4bfb-a306-cba7adb278f6
begin

	# compute W1 and W2 estimates -
	# W1 -
	⍺ = data_driven_v_variable[1]
	W1 = (⍺)/(1-⍺)

	# W2 -
	β = data_driven_v_variable[end]
	W2 = ((1-β)^-1)*(β*(1+W1)-W1)

	# return -
	nothing
end

# ╔═╡ 11df9c9f-e58f-450a-b630-3010db01391b
(W1,W2)

# ╔═╡ 8d8f4af4-8519-4a83-9027-bbf63b56e038
# estimate parameters -
parameter_array = estimate_parameters(v_variable_parameter_dictionary, data_driven_v_variable, experimental_data_array)

# ╔═╡ 81e97235-a5d6-4404-a701-c089f76aa791
begin

	# note: I put the values for K and n into the config file -

	# let's make an array of effector values -
	effector_concentration_array = range(0,stop=1.0, length=100) |> collect

	# compute the model v-variable -
	model_v_variable = compute_activity_variable(effector_concentration_array, v_variable_parameter_dictionary)

	# plot the data -
	plot(effector_concentration_array,model_v_variable,label="model",legend=:bottomright)
	scatter!(experimental_data_array[:,1],data_driven_v_variable,label="mean data")
	xlabel!("Effector 3-5-AMP (mM)", fontsize=16)
	ylabel!("Estimated v-variable (AU)", fontsize=16)

end

# ╔═╡ Cell order:
# ╟─4b710f82-b1cc-11eb-2d14-11b6421a813f
# ╟─0df967d7-ca1a-4a19-9c40-7f3a94764504
# ╠═f4b7d3df-213c-46a6-8156-a83fcf6ffb34
# ╟─e0a279c1-7677-4c5b-a6e1-ee450b27020e
# ╟─aa6bba31-0d33-4476-acb8-707128e63a72
# ╟─c4077a00-bef2-4003-ba51-8596550335b2
# ╟─50886df5-d361-4646-898f-795a064c6001
# ╠═5bab3b76-be8a-44c9-8de3-ac35992bda6a
# ╠═2bb3a988-e054-4bfb-a306-cba7adb278f6
# ╠═11df9c9f-e58f-450a-b630-3010db01391b
# ╟─a3264531-d9c8-42cf-b671-c6f46980a5a4
# ╠═8d8f4af4-8519-4a83-9027-bbf63b56e038
# ╟─dc395064-55a7-42fe-8fb0-880692e7d81a
# ╠═81e97235-a5d6-4404-a701-c089f76aa791
# ╟─b1c395bd-6fa7-4e07-a339-894675791d91
# ╠═820c02ef-da91-482d-8bab-36de9b524fa1
# ╠═169f11a1-34dd-45a9-8f9a-6b18f4957998
# ╠═de9f574f-fc03-442c-90bd-810f44673466
# ╠═b2fe595d-e90e-4041-bba4-8a1558bb83e2
# ╠═aa2267d4-ae47-4532-9fa0-460498b40fea
# ╠═b25aeabb-2be1-48ef-bd72-3680ae8e9407
# ╠═5f846f30-f20f-4153-b87f-71ab429019cb