text:
# sqrt:
# tests special values from csqrt man page
# as well as conj(sqrt(z)) = sqrt(conj(z))
@test isequal(sqrt(complex( 0.0, 0.0)), complex( 0.0, 0.0))
@test isequal(sqrt(complex( 0.0,-0.0)), complex( 0.0,-0.0))
@test isequal(sqrt(complex( 0.0, Inf)), complex( Inf, Inf))
@test isequal(sqrt(complex( 0.0,-Inf)), complex( Inf,-Inf))
@test isequal(sqrt(complex( 0.0, NaN)), complex( NaN, NaN))
@test isequal(sqrt(complex(-0.0, 0.0)), complex( 0.0, 0.0))
@test isequal(sqrt(complex(-0.0,-0.0)), complex( 0.0,-0.0))

@test isequal(sqrt(complex( 5.0, 0.0)), complex(sqrt(5.0), 0.0))
@test isequal(sqrt(complex( 5.0,-0.0)), complex(sqrt(5.0),-0.0))
@test isequal(sqrt(complex(-5.0, 0.0)), complex( 0.0, sqrt(5.0)))
@test isequal(sqrt(complex(-5.0,-0.0)), complex( 0.0,-sqrt(5.0)))

@test isequal(sqrt(complex( Inf, 0.0)), complex( Inf, 0.0))
@test isequal(sqrt(complex( Inf,-0.0)), complex( Inf,-0.0))
@test isequal(sqrt(complex( Inf, 5.0)), complex( Inf, 0.0))
@test isequal(sqrt(complex( Inf,-5.0)), complex( Inf,-0.0))
@test isequal(sqrt(complex( Inf, Inf)), complex( Inf, Inf))
@test isequal(sqrt(complex( Inf,-Inf)), complex( Inf,-Inf))
@test isequal(sqrt(complex( Inf, NaN)), complex( Inf, NaN))

@test isequal(sqrt(complex(-Inf, 0.0)), complex( 0.0, Inf))
@test isequal(sqrt(complex(-Inf,-0.0)), complex( 0.0,-Inf))
@test isequal(sqrt(complex(-Inf, 5.0)), complex( 0.0, Inf))
@test isequal(sqrt(complex(-Inf,-5.0)), complex( 0.0,-Inf))
@test isequal(sqrt(complex(-Inf, Inf)), complex( Inf, Inf))
@test isequal(sqrt(complex(-Inf,-Inf)), complex( Inf,-Inf))
@test isequal(sqrt(complex(-Inf, NaN)), complex( NaN, Inf))

@test isequal(sqrt(complex( NaN, 0.0)), complex( NaN, NaN))
@test isequal(sqrt(complex( NaN, 0.0)), complex( NaN, NaN))
@test isequal(sqrt(complex( NaN, Inf)), complex( Inf, Inf))
@test isequal(sqrt(complex( NaN,-Inf)), complex( Inf,-Inf))

# log:
#  log(conj(z)) = conj(log(z))

@test isequal(log(complex( 0.0, 0.0)), complex(-Inf, 0.0))
@test isequal(log(complex( 0.0,-0.0)), complex(-Inf,-0.0))
@test isequal(log(complex( 0.0, 1.0)), complex( 0.0, pi/2))
@test isequal(log(complex( 0.0,-1.0)), complex( 0.0,-pi/2))
@test isequal(log(complex( 0.0, Inf)), complex( Inf, pi/2))
@test isequal(log(complex( 0.0,-Inf)), complex( Inf,-pi/2))
@test isequal(log(complex( 0.0, NaN)), complex( NaN, NaN))
@test isequal(log(complex(-0.0, 0.0)), complex(-Inf, pi))
@test isequal(log(complex(-0.0,-0.0)), complex(-Inf,-pi))

@test isequal(log(complex( 5.0, 0.0)),complex(log(5.0), 0.0))
@test isequal(log(complex( 5.0,-0.0)),complex(log(5.0),-0.0))

@test isequal(log(complex( Inf, 5.0)), complex( Inf, 0.0))
@test isequal(log(complex( Inf,-5.0)), complex( Inf,-0.0))
@test isequal(log(complex( Inf, Inf)), complex( Inf, pi/4))
@test isequal(log(complex( Inf,-Inf)), complex( Inf,-pi/4))
@test isequal(log(complex( Inf, NaN)), complex( Inf, NaN))

@test isequal(log(complex(-Inf, 5.0)), complex( Inf, pi))
@test isequal(log(complex(-Inf,-5.0)), complex( Inf,-pi))
@test isequal(log(complex(-Inf, Inf)), complex( Inf, 3*pi/4))
@test isequal(log(complex(-Inf,-Inf)), complex( Inf,-3*pi/4))
@test isequal(log(complex(-Inf, NaN)), complex( Inf, NaN))

@test isequal(log(complex( NaN, 0.0)), complex( NaN, NaN))
@test isequal(log(complex( NaN, Inf)), complex( Inf, NaN))
@test isequal(log(complex( NaN,-Inf)), complex( Inf, NaN))
@test isequal(log(complex( NaN, NaN)), complex( NaN, NaN))

# exp:
#  exp(conj(z)) = conj(exp(z))

@test isequal(exp(complex( 0.0, 0.0)), complex(1.0, 0.0))
@test isequal(exp(complex( 0.0,-0.0)), complex(1.0,-0.0))
@test isequal(exp(complex( 0.0, Inf)), complex(NaN, NaN))
@test isequal(exp(complex( 0.0,-Inf)), complex(NaN, NaN))
@test isequal(exp(complex( 0.0, NaN)), complex(NaN, NaN))

@test isequal(exp(complex(-0.0, 0.0)), complex(1.0, 0.0))
@test isequal(exp(complex(-0.0,-0.0)), complex(1.0,-0.0))

@test isequal(exp(complex( 5.0, Inf)), complex(NaN, NaN))

@test isequal(exp(complex( Inf, 0.0)), complex(Inf, 0.0))
@test isequal(exp(complex( Inf,-0.0)), complex(Inf,-0.0))
@test isequal(exp(complex( Inf, 5.0)), complex(cos(5.0)*Inf,sin(5.0)* Inf))
@test isequal(exp(complex( Inf,-5.0)), complex(cos(5.0)*Inf,sin(5.0)*-Inf))
@test isequal(exp(complex( Inf, NaN)), complex(-Inf, NaN))
@test isequal(exp(complex( Inf, Inf)), complex(-Inf, NaN))
@test isequal(exp(complex( Inf,-Inf)), complex(-Inf, NaN))

@test isequal(exp(complex(-Inf, 0.0)), complex(0.0, 0.0))
@test isequal(exp(complex(-Inf,-0.0)), complex(0.0,-0.0))
@test isequal(exp(complex(-Inf, 5.0)), complex(cos(5.0)*0.0,sin(5.0)* 0.0))
@test isequal(exp(complex(-Inf,-5.0)), complex(cos(5.0)*0.0,sin(5.0)*-0.0))
@test isequal(exp(complex(-Inf, Inf)), complex(-0.0, 0.0))
@test isequal(exp(complex(-Inf,-Inf)), complex(-0.0,-0.0))
@test isequal(exp(complex(-Inf, NaN)), complex(-0.0, 0.0))

@test isequal(exp(complex( NaN, 0.0)), complex( NaN, 0.0))
@test isequal(exp(complex( NaN,-0.0)), complex( NaN,-0.0))
@test isequal(exp(complex( NaN, 5.0)), complex( NaN, NaN))
@test isequal(exp(complex( NaN, NaN)), complex( NaN, NaN))

# ^ (cpow)
#  equivalent to exp(y*log(x))
#    except for 0^0?
#  conj(x)^conj(y) = conj(x^y)
@test isequal(complex( 0.0, 0.0)^complex( 0.0, 0.0), complex(1.0, 0.0))
@test isequal(complex( 0.0, 0.0)^complex( 0.0,-0.0), complex(1.0, 0.0))
@test isequal(complex( 0.0, 0.0)^complex(-0.0, 0.0), complex(1.0,-0.0))
@test isequal(complex( 0.0, 0.0)^complex(-0.0,-0.0), complex(1.0,-0.0))

@test isequal(complex( 0.0,-0.0)^complex( 0.0, 0.0), complex(1.0,-0.0))
@test isequal(complex( 0.0,-0.0)^complex( 0.0,-0.0), complex(1.0,-0.0))
@test isequal(complex( 0.0,-0.0)^complex(-0.0, 0.0), complex(1.0, 0.0))
@test isequal(complex( 0.0,-0.0)^complex(-0.0,-0.0), complex(1.0, 0.0))

@test isequal(complex(-0.0, 0.0)^complex( 0.0, 0.0), complex(1.0, 0.0))
@test isequal(complex(-0.0, 0.0)^complex( 0.0,-0.0), complex(1.0, 0.0))
@test isequal(complex(-0.0, 0.0)^complex(-0.0, 0.0), complex(1.0,-0.0))
@test isequal(complex(-0.0, 0.0)^complex(-0.0,-0.0), complex(1.0,-0.0))

@test isequal(complex(-0.0,-0.0)^complex( 0.0, 0.0), complex(1.0,-0.0))
@test isequal(complex(-0.0,-0.0)^complex( 0.0,-0.0), complex(1.0,-0.0))
@test isequal(complex(-0.0,-0.0)^complex(-0.0, 0.0), complex(1.0, 0.0))
@test isequal(complex(-0.0,-0.0)^complex(-0.0,-0.0), complex(1.0, 0.0))

# sinh: has properties
#  sinh(conj(z)) = conj(sinh(z))
#  sinh(-z) = -sinh(z)

@test isequal(sinh(complex( 0.0, 0.0)), complex( 0.0, 0.0))
@test isequal(sinh(complex( 0.0,-0.0)), complex( 0.0,-0.0))
@test isequal(sinh(complex( 0.0, Inf)), complex( 0.0, NaN))
@test isequal(sinh(complex( 0.0,-Inf)), complex( 0.0, NaN))
@test isequal(sinh(complex( 0.0, NaN)), complex( 0.0, NaN))

@test isequal(sinh(complex(-0.0, 0.0)), complex(-0.0, 0.0))
@test isequal(sinh(complex(-0.0,-0.0)), complex(-0.0,-0.0))
@test isequal(sinh(complex(-0.0, Inf)), complex(-0.0, NaN))
@test isequal(sinh(complex(-0.0,-Inf)), complex(-0.0, NaN))
@test isequal(sinh(complex(-0.0, NaN)), complex(-0.0, NaN))

@test isequal(sinh(complex( 5.0, Inf)), complex( NaN, NaN))
@test isequal(sinh(complex( 5.0, NaN)), complex( NaN, NaN))

@test isequal(sinh(complex( Inf, 0.0)), complex( Inf, 0.0))
@test isequal(sinh(complex( Inf,-0.0)), complex( Inf,-0.0))
@test isequal(sinh(complex( Inf, 5.0)), complex(cos(5.0)*Inf,sin(5.0)*Inf))
@test isequal(sinh(complex( Inf, Inf)), complex( Inf, NaN))
@test isequal(sinh(complex( Inf,-Inf)), complex( Inf, NaN))
@test isequal(sinh(complex( Inf, NaN)), complex( Inf, NaN))

@test isequal(sinh(complex(-Inf, 0.0)), complex(-Inf, 0.0))
@test isequal(sinh(complex(-Inf,-0.0)), complex(-Inf,-0.0))
@test isequal(sinh(complex(-Inf, 5.0)), complex(cos(5.0)*-Inf,sin(5.0)*Inf))
@test isequal(sinh(complex(-Inf, Inf)), complex(-Inf, NaN))
@test isequal(sinh(complex(-Inf,-Inf)), complex(-Inf, NaN))
@test isequal(sinh(complex(-Inf, NaN)), complex(-Inf, NaN))

@test isequal(sinh(complex( NaN, 0.0)), complex( NaN, 0.0))
@test isequal(sinh(complex( NaN,-0.0)), complex( NaN,-0.0))
@test isequal(sinh(complex( NaN, 5.0)), complex( NaN, NaN))
@test isequal(sinh(complex( NaN, NaN)), complex( NaN, NaN))

# sin: defined in terms of sinh
#  sin(z) = -i*sinh(i*z)
#  i.e. if sinh(a+ib) = x+iy
#    then  sin(b-ia) = y-ix
#  sin(conj(z)) = conj(sin(z))
#  sin(-z) = -sin(z)

@test isequal(sin(complex( 0.0,-0.0)),complex( 0.0,-0.0))
@test isequal(sin(complex( 0.0, 0.0)),complex( 0.0, 0.0))
@test isequal(sin(complex( 0.0, Inf)),complex( 0.0, Inf))
@test isequal(sin(complex( 0.0,-Inf)),complex( 0.0,-Inf))
@test isequal(sin(complex( 0.0, NaN)),complex( 0.0, NaN))

@test isequal(sin(complex(-0.0, 0.0)),complex(-0.0, 0.0))
@test isequal(sin(complex(-0.0,-0.0)),complex(-0.0,-0.0))
@test isequal(sin(complex(-0.0, Inf)),complex(-0.0, Inf))
@test isequal(sin(complex(-0.0,-Inf)),complex(-0.0,-Inf))
@test isequal(sin(complex(-0.0, NaN)),complex(-0.0, NaN))

@test isequal(sin(complex( 5.0, Inf)),complex(sin(5.0)*Inf,cos(5.0)* Inf))
@test isequal(sin(complex( 5.0,-Inf)),complex(sin(5.0)*Inf,cos(5.0)*-Inf))
@test isequal(sin(complex( 5.0, NaN)),complex( NaN, NaN))

@test isequal(sin(complex( Inf, 0.0)), complex( NaN, 0.0))
@test isequal(sin(complex( Inf,-0.0)), complex( NaN,-0.0))
@test isequal(sin(complex( Inf, 5.0)), complex( NaN, NaN))
@test isequal(sin(complex( Inf, Inf)), complex( Inf, NaN))
@test isequal(sin(complex( Inf,-Inf)), complex( Inf, NaN))

@test isequal(sin(complex(-Inf, 0.0)), complex( NaN, 0.0))
@test isequal(sin(complex(-Inf,-0.0)), complex( NaN,-0.0))
@test isequal(sin(complex(-Inf, 5.0)), complex( NaN, NaN))
@test isequal(sin(complex(-Inf, Inf)), complex(-Inf, NaN))
@test isequal(sin(complex(-Inf,-Inf)), complex(-Inf, NaN))

@test isequal(sin(complex( NaN, 0.0)),complex( NaN, 0.0))
@test isequal(sin(complex( NaN,-0.0)),complex( NaN,-0.0))
@test isequal(sin(complex( NaN, 5.0)),complex( NaN, NaN))
@test isequal(sin(complex( NaN, Inf)),complex( NaN, Inf))
@test isequal(sin(complex( NaN,-Inf)),complex( NaN,-Inf))
@test isequal(sin(complex( NaN, NaN)),complex( NaN, NaN))

# cosh: has properties
#  cosh(conj(z)) = conj(cosh(z))
#  coshh(-z) = cosh(z)

@test isequal(cosh(complex( 0.0, 0.0)), complex( 1.0, 0.0))
@test isequal(cosh(complex( 0.0,-0.0)), complex( 1.0,-0.0))
@test isequal(cosh(complex( 0.0, Inf)), complex( NaN, 0.0))
@test isequal(cosh(complex( 0.0,-Inf)), complex( NaN,-0.0))
@test isequal(cosh(complex( 0.0, NaN)), complex( NaN, 0.0))

@test isequal(cosh(complex(-0.0,-0.0)), complex( 1.0, 0.0))
@test isequal(cosh(complex(-0.0, 0.0)), complex( 1.0,-0.0))
@test isequal(cosh(complex(-0.0, NaN)), complex( NaN, 0.0))

@test isequal(cosh(complex( 5.0, Inf)), complex( NaN, NaN))
@test isequal(cosh(complex( 5.0, NaN)), complex( NaN, NaN))

@test isequal(cosh(complex( Inf, 0.0)), complex( Inf, 0.0))
@test isequal(cosh(complex( Inf,-0.0)), complex( Inf,-0.0))
@test isequal(cosh(complex( Inf, 5.0)), complex(cos(5.0)*Inf,sin(5.0)* Inf))
@test isequal(cosh(complex( Inf,-5.0)), complex(cos(5.0)*Inf,sin(5.0)*-Inf))
@test isequal(cosh(complex( Inf, Inf)), complex( Inf, NaN))
@test isequal(cosh(complex( Inf,-Inf)), complex( Inf, NaN))
@test isequal(cosh(complex( Inf, NaN)), complex( Inf, NaN))

@test isequal(cosh(complex(-Inf,-0.0)), complex( Inf, 0.0))
@test isequal(cosh(complex(-Inf, 0.0)), complex( Inf,-0.0))
@test isequal(cosh(complex(-Inf,-5.0)), complex(cos(5.0)*Inf,sin(5.0)* Inf))
@test isequal(cosh(complex(-Inf, 5.0)), complex(cos(5.0)*Inf,sin(5.0)*-Inf))
@test isequal(cosh(complex(-Inf, Inf)), complex( Inf, NaN))
@test isequal(cosh(complex(-Inf,-Inf)), complex( Inf, NaN))
@test isequal(cosh(complex(-Inf, NaN)), complex( Inf, NaN))

@test isequal(cosh(complex( NaN, 0.0)), complex( NaN, 0.0))
@test isequal(cosh(complex( NaN, 5.0)), complex( NaN, NaN))
@test isequal(cosh(complex( NaN,-0.0)), complex( NaN,-0.0))
@test isequal(cosh(complex( NaN,-5.0)), complex( NaN, NaN))
@test isequal(cosh(complex( NaN, NaN)), complex( NaN, NaN))

# cos
#  cos(z) = cosh(iz)
#   i.e cos(b-ia) = cosh(a+ib)
#   and cos(b+ia) = cosh(a-ib)
#  cos(conj(z)) = conj(cos(z))
#  cos(-z) = cos(z)

@test isequal(cos(complex( 0.0, 0.0)), complex( 1.0,-0.0))
@test isequal(cos(complex( 0.0,-0.0)), complex( 1.0, 0.0))
@test isequal(cos(complex( 0.0,-Inf)), complex( Inf, 0.0))
@test isequal(cos(complex( 0.0, Inf)), complex( Inf,-0.0))
@test isequal(cos(complex( 0.0, NaN)), complex( NaN, 0.0))
@test isequal(cos(complex(-0.0, 0.0)), complex( 1.0, 0.0))
@test isequal(cos(complex(-0.0,-0.0)), complex( 1.0,-0.0))
@test isequal(cos(complex(-0.0,-Inf)), complex( Inf,-0.0))
@test isequal(cos(complex(-0.0, Inf)), complex( Inf, 0.0))
@test isequal(cos(complex(-0.0, NaN)), complex( NaN,-0.0))

@test isequal(cos(complex( 5.0,-Inf)), complex(cos(5.0)*Inf,sin(5.0)* Inf))
@test isequal(cos(complex( 5.0, Inf)), complex(cos(5.0)*Inf,sin(5.0)*-Inf))
@test isequal(cos(complex( 5.0, NaN)), complex( NaN, NaN))
@test isequal(cos(complex(-5.0,-Inf)), complex(cos(5.0)*Inf,sin(5.0)*-Inf))
@test isequal(cos(complex(-5.0, Inf)), complex(cos(5.0)*Inf,sin(5.0)* Inf))
@test isequal(cos(complex(-5.0, NaN)), complex( NaN, NaN))

@test isequal(cos(complex( Inf, 0.0)), complex( NaN,-0.0))
@test isequal(cos(complex( Inf, 5.0)), complex( NaN, NaN))
@test isequal(cos(complex( Inf, Inf)), complex( Inf, NaN))
@test isequal(cos(complex( Inf,-Inf)), complex( Inf, NaN))
@test isequal(cos(complex(-Inf, 0.0)), complex( NaN, 0.0))
@test isequal(cos(complex(-Inf, Inf)), complex( Inf, NaN))
@test isequal(cos(complex(-Inf,-Inf)), complex( Inf, NaN))

@test isequal(cos(complex( NaN, 0.0)), complex( NaN, 0.0))
@test isequal(cos(complex( NaN,-0.0)), complex( NaN, 0.0))
@test isequal(cos(complex( NaN, 5.0)), complex( NaN, NaN))
@test isequal(cos(complex( NaN, Inf)), complex( Inf, NaN))
@test isequal(cos(complex( NaN,-Inf)), complex( Inf, NaN))
@test isequal(cos(complex( NaN, NaN)), complex( NaN, NaN))

# tanh
#  tanh(conj(z)) = conj(tanh(z))
#  tanh(-z) = -tanh(z)
@test isequal(tanh(complex( 0.0, 0.0)),complex(0.0,0.0))
@test isequal(tanh(complex( 0.0,-0.0)),complex(0.0,-0.0))
@test_throws   tanh(complex( 0.0, Inf)) #complex(NaN,0.0)
@test_throws   tanh(complex( 0.0,-Inf)) #complex(NaN,-0.0)
@test isequal(tanh(complex( 0.0, NaN)),complex(NaN,NaN))

@test isequal(tanh(complex(-0.0, 0.0)),complex(-0.0,0.0))
@test isequal(tanh(complex(-0.0,-0.0)),complex(-0.0,-0.0))

@test_throws   tanh(complex( 5.0, Inf)) #complex(NaN,NaN)
@test isequal(tanh(complex( 5.0, NaN)),complex(NaN,NaN))

@test isequal(tanh(complex( Inf, 0.0)),complex(1.0, 0.0))
@test isequal(tanh(complex( Inf,-0.0)),complex(1.0,-0.0))
@test isequal(tanh(complex( Inf, 5.0)),complex(1.0,sin(2*5.0)* 0.0))
@test isequal(tanh(complex( Inf,-5.0)),complex(1.0,sin(2*5.0)*-0.0))
@test isequal(tanh(complex( Inf, Inf)),complex(1.0, 0.0))
@test isequal(tanh(complex( Inf,-Inf)),complex(1.0,-0.0))
@test isequal(tanh(complex( Inf, NaN)),complex(1.0, 0.0))

@test isequal(tanh(complex(-Inf, 0.0)),complex(-1.0, 0.0))
@test isequal(tanh(complex(-Inf,-0.0)),complex(-1.0,-0.0))
@test isequal(tanh(complex(-Inf, 5.0)),complex(-1.0,sin(2*5.0)* 0.0))
@test isequal(tanh(complex(-Inf,-5.0)),complex(-1.0,sin(2*5.0)*-0.0))
@test isequal(tanh(complex(-Inf, Inf)),complex(-1.0, 0.0))
@test isequal(tanh(complex(-Inf,-Inf)),complex(-1.0,-0.0))
@test isequal(tanh(complex(-Inf, NaN)),complex(-1.0, 0.0))

@test isequal(tanh(complex( NaN, 0.0)),complex(NaN, 0.0))
@test isequal(tanh(complex( NaN,-0.0)),complex(NaN,-0.0))
@test isequal(tanh(complex( NaN, 5.0)),complex(NaN, NaN))
@test isequal(tanh(complex( NaN,-5.0)),complex(NaN, NaN))
@test isequal(tanh(complex( NaN, NaN)),complex(NaN, NaN))

# tan
#  tan(z) = -i tanh(iz)

@test isequal(tan(complex( 0.0, Inf)),complex( 0.0, 1.0))
@test isequal(tan(complex( 0.0,-Inf)),complex( 0.0,-1.0))
@test isequal(tan(complex( 0.0, NaN)),complex( 0.0, NaN))
@test isequal(tan(complex(-0.0,-Inf)),complex(-0.0,-1.0))
@test isequal(tan(complex(-0.0, Inf)),complex(-0.0, 1.0))
@test isequal(tan(complex(-0.0, NaN)),complex(-0.0, NaN))

@test isequal(tan(complex( 5.0, Inf)),complex(sin(2*5.0)* 0.0, 1.0))
@test isequal(tan(complex( 5.0,-Inf)),complex(sin(2*5.0)* 0.0,-1.0))
@test isequal(tan(complex( 5.0, NaN)),complex( NaN, NaN))
@test isequal(tan(complex(-5.0, Inf)),complex(sin(2*5.0)*-0.0, 1.0))
@test isequal(tan(complex(-5.0,-Inf)),complex(sin(2*5.0)*-0.0,-1.0))
@test isequal(tan(complex(-5.0, NaN)),complex( NaN, NaN))

@test_throws   tan(complex( Inf, 5.0)) #complex(NaN,NaN)
@test isequal(tan(complex( Inf, Inf)),complex( 0.0, 1.0))
@test isequal(tan(complex( Inf,-Inf)),complex (0.0,-1.0))
@test isequal(tan(complex(-Inf, Inf)),complex(-0.0, 1.0))
@test isequal(tan(complex(-Inf,-Inf)),complex(-0.0,-1.0))

@test isequal(tan(complex( NaN, 5.0)),complex( NaN, NaN))
@test isequal(tan(complex( NaN, Inf)),complex( 0.0, 1.0))
@test isequal(tan(complex( NaN,-Inf)),complex( 0.0,-1.0))
@test isequal(tan(complex( NaN, NaN)),complex( NaN, NaN))

# acosh
#  acosh(conj(z)) = conj(acosh(z))

@test isequal(acosh(complex( 0.0, 0.0)), complex( 0.0, pi/2))
@test isequal(acosh(complex( 0.0,-0.0)), complex( 0.0,-pi/2))
@test isequal(acosh(complex( 0.0, Inf)), complex( Inf, pi/2))
@test isequal(acosh(complex( 0.0,-Inf)), complex( Inf,-pi/2))
@test isequal(acosh(complex(-0.0, 0.0)), complex( 0.0, pi/2))
@test isequal(acosh(complex(-0.0,-0.0)), complex( 0.0,-pi/2))
@test isequal(acosh(complex( 5.0, Inf)), complex( Inf, pi/2))
@test isequal(acosh(complex( 5.0,-Inf)), complex( Inf,-pi/2))
@test isequal(acosh(complex( 5.0, NaN)), complex( NaN, NaN))

@test isequal(acosh(complex( Inf, 0.0)), complex( Inf, 0.0))
@test isequal(acosh(complex( Inf,-0.0)), complex( Inf,-0.0))
@test isequal(acosh(complex( Inf, 5.0)), complex( Inf, 0.0))
@test isequal(acosh(complex( Inf,-5.0)), complex( Inf,-0.0))
@test isequal(acosh(complex( Inf, Inf)), complex( Inf, pi/4))
@test isequal(acosh(complex( Inf,-Inf)), complex( Inf,-pi/4))
@test isequal(acosh(complex( Inf, NaN)), complex( Inf, NaN))

@test isequal(acosh(complex(-Inf, 0.0)), complex( Inf, pi))
@test isequal(acosh(complex(-Inf,-0.0)), complex( Inf,-pi))
@test isequal(acosh(complex(-Inf, 5.0)), complex( Inf, pi))
@test isequal(acosh(complex(-Inf,-5.0)), complex( Inf,-pi))
@test isequal(acosh(complex(-Inf, Inf)), complex( Inf, 3*pi/4))
@test isequal(acosh(complex(-Inf,-Inf)), complex( Inf,-3*pi/4))
@test isequal(acosh(complex(-Inf, NaN)), complex( Inf, NaN))

@test isequal(acosh(complex( NaN, Inf)), complex( Inf, NaN))
@test isequal(acosh(complex( NaN,-Inf)), complex( Inf, NaN))
@test isequal(acosh(complex( NaN, NaN)), complex( NaN, NaN))

#
## acos
##  acos(conj(z)) = conj(acos(z))

@test isequal(acos(complex( 0.0, 0.0)),complex(pi/2,-0.0))
@test isequal(acos(complex( 0.0,-0.0)),complex(pi/2, 0.0))
@test isequal(acos(complex( 0.0, Inf)),complex(pi/2,-Inf))
@test isequal(acos(complex( 0.0,-Inf)),complex(pi/2, Inf))
@test isequal(acos(complex( 0.0, NaN)),complex(pi/2, NaN))

@test isequal(acos(complex(-0.0, 0.0)),complex(pi/2,-0.0))
@test isequal(acos(complex(-0.0,-0.0)),complex(pi/2, 0.0))
@test isequal(acos(complex(-0.0, NaN)),complex(pi/2, NaN))

@test isequal(acos(complex( 5.0, Inf)),complex(pi/2,-Inf))
@test isequal(acos(complex( 5.0,-Inf)),complex(pi/2, Inf))
@test isequal(acos(complex( 5.0, NaN)),complex( NaN, NaN))

@test isequal(acos(complex( Inf, 0.0)),complex( 0.0,-Inf))
@test isequal(acos(complex( Inf,-0.0)),complex( 0.0, Inf))
@test isequal(acos(complex( Inf, 5.0)),complex( 0.0,-Inf))
@test isequal(acos(complex( Inf,-5.0)),complex( 0.0, Inf))
@test isequal(acos(complex( Inf, Inf)),complex(pi/4,-Inf))
@test isequal(acos(complex( Inf,-Inf)),complex(pi/4, Inf))
@test isequal(acos(complex( Inf, NaN)),complex( NaN, Inf))

@test isequal(acos(complex(-Inf, 0.0)),complex(pi,-Inf))
@test isequal(acos(complex(-Inf,-0.0)),complex(pi, Inf))
@test isequal(acos(complex(-Inf, 5.0)),complex(pi,-Inf))
@test isequal(acos(complex(-Inf,-5.0)),complex(pi, Inf))
@test isequal(acos(complex(-Inf, Inf)),complex(3*pi/4,-Inf))
@test isequal(acos(complex(-Inf,-Inf)),complex(3*pi/4, Inf))
@test isequal(acos(complex(-Inf, NaN)),complex( NaN, Inf))

@test isequal(acos(complex( NaN, 0.0)),complex( NaN, NaN))
@test isequal(acos(complex( NaN, 5.0)),complex( NaN, NaN))
@test isequal(acos(complex( NaN, Inf)),complex( NaN,-Inf))
@test isequal(acos(complex( NaN,-Inf)),complex( NaN, Inf))
@test isequal(acos(complex( NaN, NaN)),complex( NaN, NaN))

## asinh
##  asinh(conj(z)) = conj(asinh(z))
##  asinh(-z) = -asinh(z)
@test isequal(asinh(complex( 0.0, 0.0)),complex( 0.0, 0.0))
@test isequal(asinh(complex( 0.0,-0.0)),complex( 0.0,-0.0))
@test isequal(asinh(complex( 0.0, Inf)),complex( Inf, pi/2))
@test isequal(asinh(complex( 0.0,-Inf)),complex( Inf,-pi/2))
@test isequal(asinh(complex( 0.0, NaN)),complex( NaN, NaN))

@test isequal(asinh(complex(-0.0, 0.0)),complex(-0.0, 0.0))
@test isequal(asinh(complex(-0.0,-0.0)),complex(-0.0,-0.0))
@test isequal(asinh(complex(-0.0, Inf)),complex(-Inf, pi/2))
@test isequal(asinh(complex(-0.0,-Inf)),complex(-Inf,-pi/2))

@test isequal(asinh(complex( 5.0, Inf)),complex( Inf, pi/2))
@test isequal(asinh(complex( 5.0,-Inf)),complex( Inf,-pi/2))
@test isequal(asinh(complex( 5.0, NaN)),complex( NaN, NaN))
@test isequal(asinh(complex(-5.0, Inf)),complex(-Inf, pi/2))
@test isequal(asinh(complex(-5.0,-Inf)),complex(-Inf,-pi/2))

@test isequal(asinh(complex( Inf, Inf)),complex( Inf, pi/4))
@test isequal(asinh(complex( Inf,-Inf)),complex( Inf,-pi/4))
@test isequal(asinh(complex( Inf, NaN)),complex( Inf, NaN))
@test isequal(asinh(complex(-Inf, Inf)),complex(-Inf, pi/4))
@test isequal(asinh(complex(-Inf,-Inf)),complex(-Inf,-pi/4))
@test isequal(asinh(complex(-Inf, NaN)),complex(-Inf, NaN))

@test isequal(asinh(complex( NaN, 0.0)),complex( NaN, 0.0))
@test isequal(asinh(complex( NaN,-0.0)),complex( NaN,-0.0))
@test isequal(asinh(complex( NaN, 5.0)),complex( NaN, NaN))
@test isequal(asinh(complex( NaN, Inf)),complex( Inf, NaN))
@test isequal(asinh(complex( NaN,-Inf)),complex( Inf, NaN))
@test isequal(asinh(complex( NaN, NaN)),complex( NaN, NaN))

# asin
#  asin(z) = -i*asinh(iz)

@test isequal(asin(complex( 0.0, 0.0)),complex( 0.0, 0.0))
@test isequal(asin(complex( 0.0,-0.0)),complex( 0.0,-0.0))
@test isequal(asin(complex(-0.0, 0.0)),complex(-0.0, 0.0))
@test isequal(asin(complex( 0.0, NaN)),complex( 0.0, NaN))
@test isequal(asin(complex(-0.0,-0.0)),complex(-0.0,-0.0))
@test isequal(asin(complex(-0.0, NaN)),complex(-0.0, NaN))
@test isequal(asin(complex( 5.0, NaN)),complex( NaN, NaN))

@test isequal(asin(complex( Inf, 0.0)),complex( pi/2, Inf))
@test isequal(asin(complex( Inf,-0.0)),complex( pi/2,-Inf))
@test isequal(asin(complex( Inf, 5.0)),complex( pi/2, Inf))
@test isequal(asin(complex( Inf,-5.0)),complex( pi/2,-Inf))
@test isequal(asin(complex( Inf, Inf)),complex( pi/4, Inf))
@test isequal(asin(complex( Inf,-Inf)),complex( pi/4,-Inf))
@test isequal(asin(complex( Inf, NaN)),complex( NaN, Inf))

@test isequal(asin(complex(-Inf, 0.0)),complex(-pi/2, Inf))
@test isequal(asin(complex(-Inf,-0.0)),complex(-pi/2,-Inf))
@test isequal(asin(complex(-Inf, 5.0)),complex(-pi/2, Inf))
@test isequal(asin(complex(-Inf,-5.0)),complex(-pi/2,-Inf))
@test isequal(asin(complex(-Inf, Inf)),complex(-pi/4, Inf))
@test isequal(asin(complex(-Inf,-Inf)),complex(-pi/4,-Inf))
@test isequal(asin(complex(-Inf, NaN)),complex( NaN, Inf))

@test isequal(asin(complex( NaN, 0.0)),complex( NaN, NaN))
@test isequal(asin(complex( NaN, 5.0)),complex( NaN, NaN))
@test isequal(asin(complex( NaN, Inf)),complex( NaN, Inf))
@test isequal(asin(complex( NaN,-Inf)),complex( NaN,-Inf))
@test isequal(asin(complex( NaN, NaN)),complex( NaN, NaN))

# atanh
#  atanh(conj(z)) = conj(atanh(z))
#  atanh(-z) = -atanh(z)

@test isequal(atanh(complex( 0.0, 0.0)),complex( 0.0, 0.0))
@test isequal(atanh(complex( 0.0,-0.0)),complex( 0.0,-0.0))
@test isequal(atanh(complex( 0.0, NaN)),complex( 0.0, NaN))
@test isequal(atanh(complex( 0.0, Inf)),complex( 0.0, pi/2))
@test isequal(atanh(complex( 0.0,-Inf)),complex( 0.0,-pi/2))

@test isequal(atanh(complex(-0.0, NaN)),complex(-0.0, NaN))
@test isequal(atanh(complex(-0.0, 0.0)),complex(-0.0, 0.0))
@test isequal(atanh(complex(-0.0,-0.0)),complex(-0.0,-0.0))
@test isequal(atanh(complex(-0.0, Inf)),complex(-0.0, pi/2))
@test isequal(atanh(complex(-0.0,-Inf)),complex(-0.0,-pi/2))

@test isequal(atanh(complex( 1.0, 0.0)),complex( Inf, 0.0))
@test isequal(atanh(complex(-1.0, 0.0)),complex(-Inf, 0.0))
@test isequal(atanh(complex( 5.0, Inf)),complex( 0.0, pi/2))
@test isequal(atanh(complex( 5.0,-Inf)),complex( 0.0,-pi/2))
@test isequal(atanh(complex( 5.0, NaN)),complex( NaN, NaN))
@test isequal(atanh(complex(-5.0, Inf)),complex(-0.0, pi/2))
@test isequal(atanh(complex(-5.0,-Inf)),complex(-0.0,-pi/2))
@test isequal(atanh(complex(-5.0, NaN)),complex( NaN, NaN))

@test isequal(atanh(complex( Inf, 0.0)),complex(0.0, pi/2))
@test isequal(atanh(complex( Inf,-0.0)),complex(0.0,-pi/2))
@test isequal(atanh(complex( Inf, 5.0)),complex(0.0, pi/2))
@test isequal(atanh(complex( Inf,-5.0)),complex(0.0,-pi/2))
@test isequal(atanh(complex( Inf, Inf)),complex(0.0, pi/2))
@test isequal(atanh(complex( Inf,-Inf)),complex(0.0,-pi/2))
@test isequal(atanh(complex( Inf, NaN)),complex(0.0, NaN))

@test isequal(atanh(complex(-Inf, 0.0)),complex(-0.0, pi/2))
@test isequal(atanh(complex(-Inf,-0.0)),complex(-0.0,-pi/2))
@test isequal(atanh(complex(-Inf, 5.0)),complex(-0.0, pi/2))
@test isequal(atanh(complex(-Inf,-5.0)),complex(-0.0,-pi/2))
@test isequal(atanh(complex(-Inf, Inf)),complex(-0.0, pi/2))
@test isequal(atanh(complex(-Inf,-Inf)),complex(-0.0,-pi/2))
@test isequal(atanh(complex(-Inf, NaN)),complex(-0.0, NaN))

@test isequal(atanh(complex( NaN, 0.0)),complex( NaN, NaN))
@test isequal(atanh(complex( NaN,-0.0)),complex( NaN, NaN))
@test isequal(atanh(complex( NaN, 5.0)),complex( NaN, NaN))
@test isequal(atanh(complex( NaN,-5.0)),complex( NaN, NaN))
@test isequal(atanh(complex( NaN, Inf)),complex( 0.0, pi/2))
@test isequal(atanh(complex( NaN,-Inf)),complex( 0.0,-pi/2))
@test isequal(atanh(complex( NaN, NaN)),complex( NaN, NaN))

# atan
#  atan(z) = -i*atanh(iz)

@test isequal(atan(complex( 0.0, 0.0)),complex( 0.0, 0.0))
@test isequal(atan(complex( 0.0,-0.0)),complex( 0.0,-0.0))
@test isequal(atan(complex( 0.0, 1.0)),complex( 0.0, Inf))
@test isequal(atan(complex( 0.0, Inf)),complex( pi/2, 0.0))
@test isequal(atan(complex( 0.0,-Inf)),complex( pi/2,-0.0))
@test isequal(atan(complex( 0.0, NaN)),complex( NaN, NaN))

@test isequal(atan(complex(-0.0, 0.0)),complex(-0.0, 0.0))
@test isequal(atan(complex(-0.0,-0.0)),complex(-0.0,-0.0))
@test isequal(atan(complex(-0.0, Inf)),complex(-pi/2, 0.0))
@test isequal(atan(complex(-0.0,-Inf)),complex(-pi/2,-0.0))
@test isequal(atan(complex(-0.0, NaN)),complex( NaN, NaN))

@test isequal(atan(complex( 5.0, Inf)),complex( pi/2, 0.0))
@test isequal(atan(complex( 5.0,-Inf)),complex( pi/2,-0.0))
@test isequal(atan(complex( 5.0, NaN)),complex( NaN, NaN))

@test isequal(atan(complex(-5.0, Inf)),complex(-pi/2, 0.0))
@test isequal(atan(complex(-5.0,-Inf)),complex(-pi/2,-0.0))
@test isequal(atan(complex(-5.0, NaN)),complex( NaN, NaN))

@test isequal(atan(complex( Inf, 0.0)),complex( pi/2, 0.0))
@test isequal(atan(complex( Inf,-0.0)),complex( pi/2,-0.0))
@test isequal(atan(complex( Inf, 5.0)),complex( pi/2, 0.0))
@test isequal(atan(complex( Inf,-5.0)),complex( pi/2,-0.0))
@test isequal(atan(complex( Inf, Inf)),complex( pi/2, 0.0))
@test isequal(atan(complex( Inf,-Inf)),complex( pi/2,-0.0))
@test isequal(atan(complex( Inf, NaN)),complex( pi/2, 0.0))

@test isequal(atan(complex(-Inf, 0.0)),complex(-pi/2, 0.0))
@test isequal(atan(complex(-Inf,-0.0)),complex(-pi/2,-0.0))
@test isequal(atan(complex(-Inf, 5.0)),complex(-pi/2, 0.0))
@test isequal(atan(complex(-Inf,-5.0)),complex(-pi/2,-0.0))
@test isequal(atan(complex(-Inf, Inf)),complex(-pi/2, 0.0))
@test isequal(atan(complex(-Inf,-Inf)),complex(-pi/2,-0.0))
@test isequal(atan(complex(-Inf, NaN)),complex(-pi/2, 0.0))

@test isequal(atan(complex( NaN, 0.0)),complex( NaN, 0.0))
@test isequal(atan(complex( NaN,-0.0)),complex( NaN,-0.0))
@test isequal(atan(complex( NaN, 5.0)),complex( NaN, NaN))
@test isequal(atan(complex( NaN,-5.0)),complex( NaN, NaN))
@test isequal(atan(complex( NaN, Inf)),complex( NaN, 0.0))
@test isequal(atan(complex( NaN,-Inf)),complex( NaN,-0.0))
@test isequal(atan(complex( NaN, NaN)),complex( NaN, NaN))

# lexcmp
@test lexcmp(1.0-1.0im, 1.0+0.0im) == -1
@test lexcmp(0.0+0.0im, 0.0+0.0im) == 0
@test lexcmp(1.0-1.0im, 0.0+0.0im) == 1

# misc.

@test complex(1//2,1//3)^2 === complex(5//36, 1//3)
@test complex(2,2)^2 === complex(0,8)
@test_throws complex(2,2)^(-2)
@test complex(2.0,2.0)^(-2) === complex(0.0, -0.125)

# robust division of Float64
# hard complex divisions from Fig 6 of arxiv.1210.4539
z7 = Complex{Float64}(3.898125604559113300e289, 8.174961907852353577e295)
z9 = Complex{Float64}(0.001953125, -0.001953125)
z10 = Complex{Float64}( 1.02951151789360578e-84, 6.97145987515076231e-220)
harddivs = ((1.0+im*1.0, 1.0+im*2^1023.0, 2^-1023.0-im*2^-1023.0), #1
      (1.0+im*1.0, 2^-1023.0+im*2^-1023.0, 2^1023.0+im*0.0), #2
      (2^1023.0+im*2^-1023.0, 2^677.0+im*2^-677.0, 2^346.0-im*2^-1008.0), #3
      (2^1023.0+im*2^1023.0, 1.0+im*1.0, 2^1023.0+im*0.0), #4
      (2^1020.0+im*2^-844., 2^656.0+im*2^-780.0, 2^364.0-im*2^-1072.0), #5
      (2^-71.0+im*2^1021., 2^1001.0+im*2^-323.0, 2^-1072.0+im*2^20.0), #6
      (2^-347.0+im*2^-54., 2^-1037.0+im*2^-1058.0, z7), #7
      (2^-1074.0+im*2^-1074., 2^-1073.0+im*2^-1074., 0.6+im*0.2), #8
      (2^1015.0+im*2^-989., 2^1023.0+im*2^1023.0, z9), #9
      (2^-622.0+im*2^-1071., 2^-343.0+im*2^-798.0, z10)#10
      )

# calculate "accurate bits" in range 0:53 by algorithm given in arxiv.1210.4539
function sb_accuracy(x,expected)
  min(logacc(real(x),real(expected)),
    logacc(imag(x),imag(expected)) )
end
relacc(x,expected) = abs(x-expected)/abs(expected)
function logacc(x::Float64,expected::Float64)
  x == expected && (return 53)
  expected == 0 && (return 0)
  (x == Inf || x == -Inf) && (return 0)
  isnan(x) && (return 0)
  ra = relacc(BigFloat(x),BigFloat(expected))
  max(ifloor(-log2(ra)),0)
end
# the robust division algorithm should have 53 or 52
# bits accuracy for each of the hard divisions
@test 52 <= minimum([sb_accuracy(h[1]/h[2],h[3]) for h in harddivs])

# division of non-Float64
function cdiv_test(a,b)
  c=convert(Complex{Float64},a)/convert(Complex{Float64},b)
  50 <= sb_accuracy(c,convert(Complex{Float64},a/b))
end
@test cdiv_test(complex(1//2, 3//4), complex(17//13, 4//5))
@test cdiv_test(complex(1,2), complex(8997,2432))

# inv
@test inv(1e300+0im) == 1e-300 - 0.0im

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Fixes problem where members of leaf nodes are not included in members log.
This tool is used to update old results to new format using the tree
saved in save.jld
"""
module RegenMembers

export regen_members, regen_members_recursive, DecisionTrees

using Compat
import Compat.ASCIIString

using GrammarExpts
using RLESUtils, TreeIterators
import TreeIterators.get_children
using JLD, DataFrames

#################
#recreate for loading JLD, real got updated and breaks JLD
module DecisionTrees

export DTNode, DecisionTree

type DTNode{T1,T2}
    depth::Int64
    members::Vector{Int64} #indices into data starting at 1
    split_rule::Any #object used in callback for split rule
    children::Dict{T1,DTNode{T1,T2}} #key=split_rule predicts, value=child node, T1=predict label type
    label::T2 #T2=label type
    confidence::Float64
end
type DecisionTree{T1,T2}
    root::DTNode{T1,T2}
end
end #module
######################

using .DecisionTrees

get_children(node::DTNode) = collect(values(node.children))

function regen_members(jldfile::AbstractString, memberfile::AbstractString)
    dtree = load(jldfile, "dtree")
    D = DataFrame([Int64,ASCIIString], [:decision_id,:members], 0)
    node_id = 1
    for node in tree_iter(dtree.root)
        push!(D, [node_id, string(node.members)])
        node_id += 1
    end
    writetable(memberfile, D)
    D
end

function regen_members_recursive(topdir::AbstractString, jldfile::AbstractString,
    memberfile::AbstractString)
    for (root, dirs, files) in walkdir(topdir)
        fin = joinpath(root, jldfile)
        if isfile(fin)
            fout = joinpath(root, memberfile)
            regen_members(fin, fout)
        end
    end
end

end #module

ext: jl


text:
############################
#                          #
# Generate Figures         #
#                          #
############################

### Figure 5
# Using user-provided results
"""
    createFig5(noint, notest, untarget)
    createFig5()

Replicate Figure 5 in Chari, Kirpalani, and Phelan (2021).
If no arguments are provided, then stored results obtained with default parameters are used.

# Arguments:
- `noint::Dict`: Output from [`nopolicy`](@ref).
- `notest::Dict`: Output from [`withpolicy(::String)`](@ref) with policy `notest`.
- `untarget::Dict`: Output from [`withpolicy(::String)`](@ref) with policy `untargettest`.

"""
function createFig5(noint::Dict, notest::Dict, untarget::Dict)

    println("Replicating Figure 5...")

    # Create variables
    TT = noint["TT"];
    IN_1 = 1 .- noint["ST"] .- noint["RT"];
    IN_2 = 1 .- notest["ST"] .- notest["RT"];
    IN_3 = 1 .- untarget["ST"] .- untarget["RT"];
    death_1 = cumsum(noint["Par"].gamma .* noint["Par"].delta .* IN_1);
    death_2 = cumsum(notest["Par"].gamma .* notest["Par"].delta .* IN_2);
    death_3 = cumsum(untarget["Par"].gamma .* untarget["Par"].delta .* IN_3);
    Yss = 1.01;
    Cdev_1 = ((noint["YT"]-noint["C1T"]).-Yss)./Yss.*100;
    Cdev_2 = ((notest["YT"]-notest["C1T"]).-Yss)./Yss.*100;
    Cdev_3 = ((untarget["YT"]-untarget["C1T"]).-Yss)./Yss.*100;

    # Figure 5
    args_first = (titlefontsize=8, guidefontsize=6, linecolor=["black" "green" "red"], legendfontsize=6, label=["no policy" "no testing" "untargeted"]);
    args_rest  = (titlefontsize=8, guidefontsize=6, linecolor=["black" "green" "red"], legend=false);
    p1 = plot(1:TT, [IN_1.*100 IN_2.*100 IN_3.*100], title="Infected, I_t", ylabel="%"; args_first...)
    p2 = plot(1:TT, [noint["ST"].*100 notest["ST"].*100 untarget["ST"].*100], title="Susceptible, S_t", ylabel="%"; args_rest...)
    p3 = plot(1:TT, [death_1.*100 death_2.*100 death_3.*100], title="Cumulative Deaths", ylabel="%"; args_rest...)
    p4 = plot(1:TT, [noint["Rknot"] notest["Rknot"] untarget["Rknot"]], title="R_0"; args_rest...)
    p5 = plot(1:TT, [noint["IT"].*100 notest["IT"].*100 untarget["IT"].*100], title="Unknown Infected, ~I_t", ylabel="%"; args_rest...)
    p6 = plot(1:TT, [Cdev_1 Cdev_2 Cdev_3], title="Agg. Consumption, C_t", ylabel="% dev. from steady state"; args_rest...)
    p7 = plot(1:TT, [noint["masstest"].*100 notest["masstest"].*100 untarget["masstest"].*100], title="Mass tested", ylabel="%"; args_rest...)
    p8 = plot(1:TT, [noint["mctest"] notest["mctest"] untarget["mctest"]], title="Marginal testing cost", ylabel="Dollars"; args_rest...)
    p = plot(p1, p2, p3, p4, p5, p6, p7, p8, layout = (2, 4));
    display(p)

end

# Using stored results
function createFig5()

    println("Replicating Figure 5...")

    # Load stored data
    noint = load("src/data/calib1_nointervention.jld2");
    notest = load("src/data/calib1_notest_small.jld2");
    untarget = load("src/data/calib1_untarget_small.jld2");

    # Create variables
    TT = noint["TT"];
    IN_1 = 1 .- noint["ST"] .- noint["RT"];
    IN_2 = 1 .- notest["ST"] .- notest["RT"];
    IN_3 = 1 .- untarget["ST"] .- untarget["RT"];
    death_1 = cumsum(noint["Par"].gamma .* noint["Par"].delta .* IN_1);
    death_2 = cumsum(notest["Par"].gamma .* notest["Par"].delta .* IN_2);
    death_3 = cumsum(untarget["Par"].gamma .* untarget["Par"].delta .* IN_3);
    Yss = 1.01;
    Cdev_1 = ((noint["YT"]-noint["C1T"]).-Yss)./Yss.*100;
    Cdev_2 = ((notest["YT"]-notest["C1T"]).-Yss)./Yss.*100;
    Cdev_3 = ((untarget["YT"]-untarget["C1T"]).-Yss)./Yss.*100;

    # Figure 5
    args_first = (titlefontsize=8, guidefontsize=6, linecolor=["black" "green" "red"], legendfontsize=6, label=["no policy" "no testing" "untargeted"]);
    args_rest  = (titlefontsize=8, guidefontsize=6, linecolor=["black" "green" "red"], legend=false);
    p1 = plot(1:TT, [IN_1.*100 IN_2.*100 IN_3.*100], title="Infected, I_t", ylabel="%"; args_first...)
    p2 = plot(1:TT, [noint["ST"].*100 notest["ST"].*100 untarget["ST"].*100], title="Susceptible, S_t", ylabel="%"; args_rest...)
    p3 = plot(1:TT, [death_1.*100 death_2.*100 death_3.*100], title="Cumulative Deaths", ylabel="%"; args_rest...)
    p4 = plot(1:TT, [noint["Rknot"] notest["Rknot"] untarget["Rknot"]], title="R_0"; args_rest...)
    p5 = plot(1:TT, [noint["IT"].*100 notest["IT"].*100 untarget["IT"].*100], title="Unknown Infected, ~I_t", ylabel="%"; args_rest...)
    p6 = plot(1:TT, [Cdev_1 Cdev_2 Cdev_3], title="Agg. Consumption, C_t", ylabel="% dev. from steady state"; args_rest...)
    p7 = plot(1:TT, [noint["masstest"].*100 notest["masstest"].*100 untarget["masstest"].*100], title="Mass tested", ylabel="%"; args_rest...)
    p8 = plot(1:TT, [noint["mctest"] notest["mctest"] untarget["mctest"]], title="Marginal testing cost", ylabel="Dollars"; args_rest...)
    p = plot(p1, p2, p3, p4, p5, p6, p7, p8, layout = (2, 4));
    display(p)

end

##########################################################################################
### Figure 6
# Using user-provided results
"""
    createFig6(notest, untarget, target)
    createFig6()

Analogous to [`createFig5`](@ref).
"""
function createFig6(notest::Dict, untarget::Dict, target::Dict)

    println("Replicating Figure 6...")

    # Create variables
    TT = notest["TT"];
    IN_2 = 1 .- notest["ST"] .- notest["RT"];
    IN_3 = 1 .- untarget["ST"] .- untarget["RT"];
    IN_4 = 1 .- target["ST"] .- target["RT"];
    death_2 = cumsum(notest["Par"].gamma .* notest["Par"].delta .* IN_2);
    death_3 = cumsum(untarget["Par"].gamma .* untarget["Par"].delta .* IN_3);
    death_4 = cumsum(target["Par"].gamma .* target["Par"].delta .* IN_4);
    Yss = 1.01;
    Cdev_2 = ((notest["YT"]-notest["C1T"]).-Yss)./Yss.*100;
    Cdev_3 = ((untarget["YT"]-untarget["C1T"]).-Yss)./Yss.*100;
    Cdev_4 = ((target["YT"]-target["C1T"]).-Yss)./Yss.*100;

    # Figure 6
    args_first = (titlefontsize=8, guidefontsize=6, linecolor=["green" "red" "blue"], legendfontsize=6, label=["no testing" "untargeted" "targeted"]);
    args_rest  = (titlefontsize=8, guidefontsize=6, linecolor=["green" "red" "blue"], legend=false);
    p1 = plot(1:TT, [IN_2.*100 IN_3.*100 IN_4.*100], title="Infected, I_t", ylabel="%"; args_first...)
    p2 = plot(1:TT, [notest["ST"].*100 untarget["ST"].*100 target["ST"].*100], title="Susceptible, S_t", ylabel="%"; args_rest...)
    p3 = plot(1:TT, [death_2.*100 death_3.*100 death_4.*100], title="Cumulative Deaths", ylabel="%"; args_rest...)
    p4 = plot(1:TT, [notest["Rknot"] untarget["Rknot"] target["Rknot"]], title="R_0"; args_rest...)
    p5 = plot(1:TT, [notest["IT"].*100 untarget["IT"].*100 target["IT"].*100], title="Unknown Infected, ~I_t", ylabel="%"; args_rest...)
    p6 = plot(1:TT, [Cdev_2 Cdev_3 Cdev_4], title="Agg. Consumption, C_t", ylabel="% dev. from steady state"; args_rest...)
    p7 = plot(1:TT, [notest["masstest"].*100 untarget["masstest"].*100 target["masstest"].*100], title="Mass tested", ylabel="%"; args_rest...)
    p8 = plot(1:TT, [notest["mctest"] untarget["mctest"] target["mctest"]], title="Marginal testing cost", ylabel="Dollars"; args_rest...)
    p = plot(p1, p2, p3, p4, p5, p6, p7, p8, layout = (2, 4));
    display(p)

end

# Using stored results
function createFig6()

    println("Replicating Figure 6...")

    # Load stored data
    notest = load("src/data/calib1_notest_small.jld2");
    untarget = load("src/data/calib1_untarget_small.jld2");
    target = load("src/data/calib1_target_small.jld2");

    # Create variables
    TT = notest["TT"];
    IN_2 = 1 .- notest["ST"] .- notest["RT"];
    IN_3 = 1 .- untarget["ST"] .- untarget["RT"];
    IN_4 = 1 .- target["ST"] .- target["RT"];
    death_2 = cumsum(notest["Par"].gamma .* notest["Par"].delta .* IN_2);
    death_3 = cumsum(untarget["Par"].gamma .* untarget["Par"].delta .* IN_3);
    death_4 = cumsum(target["Par"].gamma .* target["Par"].delta .* IN_4);
    Yss = 1.01;
    Cdev_2 = ((notest["YT"]-notest["C1T"]).-Yss)./Yss.*100;
    Cdev_3 = ((untarget["YT"]-untarget["C1T"]).-Yss)./Yss.*100;
    Cdev_4 = ((target["YT"]-target["C1T"]).-Yss)./Yss.*100;

    # Figure 6
    args_first = (titlefontsize=8, guidefontsize=6, linecolor=["green" "red" "blue"], legendfontsize=6, label=["no testing" "untargeted" "targeted"]);
    args_rest  = (titlefontsize=8, guidefontsize=6, linecolor=["green" "red" "blue"], legend=false);
    p1 = plot(1:TT, [IN_2.*100 IN_3.*100 IN_4.*100], title="Infected, I_t", ylabel="%"; args_first...)
    p2 = plot(1:TT, [notest["ST"].*100 untarget["ST"].*100 target["ST"].*100], title="Susceptible, S_t", ylabel="%"; args_rest...)
    p3 = plot(1:TT, [death_2.*100 death_3.*100 death_4.*100], title="Cumulative Deaths", ylabel="%"; args_rest...)
    p4 = plot(1:TT, [notest["Rknot"] untarget["Rknot"] target["Rknot"]], title="R_0"; args_rest...)
    p5 = plot(1:TT, [notest["IT"].*100 untarget["IT"].*100 target["IT"].*100], title="Unknown Infected, ~I_t", ylabel="%"; args_rest...)
    p6 = plot(1:TT, [Cdev_2 Cdev_3 Cdev_4], title="Agg. Consumption, C_t", ylabel="% dev. from steady state"; args_rest...)
    p7 = plot(1:TT, [notest["masstest"].*100 untarget["masstest"].*100 target["masstest"].*100], title="Mass tested", ylabel="%"; args_rest...)
    p8 = plot(1:TT, [notest["mctest"] untarget["mctest"] target["mctest"]], title="Marginal testing cost", ylabel="Dollars"; args_rest...)
    p = plot(p1, p2, p3, p4, p5, p6, p7, p8, layout = (2, 4));
    display(p)

end

##########################################################################################
### Figure 7
# Using user-provided results
"""
    createFig7(untarget, target, isolate)
    createFig7()

Analogous to [`createFig5`](@ref).
"""
function createFig7(untarget::Dict, target::Dict, isolate::Dict)

    println("Replicating Figure 7...")

    # Create variables
    TT = untarget["TT"];
    IN_3 = 1 .- untarget["ST"] .- untarget["RT"];
    IN_4 = 1 .- target["ST"] .- target["RT"];
    IN_5 = 1 .- isolate["ST"] .- isolate["RT"];
    death_3 = cumsum(untarget["Par"].gamma .* untarget["Par"].delta .* IN_3);
    death_4 = cumsum(target["Par"].gamma .* target["Par"].delta .* IN_4);
    death_5 = cumsum(isolate["Par"].gamma .* isolate["Par"].delta .* IN_5);
    Yss = 1.01;
    Cdev_3 = ((untarget["YT"]-untarget["C1T"]).-Yss)./Yss.*100;
    Cdev_4 = ((target["YT"]-target["C1T"]).-Yss)./Yss.*100;
    Cdev_5 = ((isolate["YT"]-isolate["C1T"]).-Yss)./Yss.*100;

    # Figure 7
    args_first = (titlefontsize=8, guidefontsize=6, linecolor=["red" "blue" "purple"], legendfontsize=6, label=["untargeted" "targeted" "isolate"]);
    args_rest  = (titlefontsize=8, guidefontsize=6, linecolor=["red" "blue" "purple"], legend=false);
    p1 = plot(1:TT, [IN_3.*100 IN_4.*100 IN_5.*100], title="Infected, I_t", ylabel="%"; args_first...)
    p2 = plot(1:TT, [untarget["ST"].*100 target["ST"].*100 isolate["ST"].*100], title="Susceptible, S_t", ylabel="%"; args_rest...)
    p3 = plot(1:TT, [death_3.*100 death_4.*100 death_5.*100], title="Cumulative Deaths", ylabel="%"; args_rest...)
    p4 = plot(1:TT, [untarget["Rknot"] target["Rknot"] isolate["Rknot"]], title="R_0"; args_rest...)
    p5 = plot(1:TT, [untarget["IT"].*100 target["IT"].*100 isolate["IT"].*100], title="Unknown Infected, ~I_t", ylabel="%"; args_rest...)
    p6 = plot(1:TT, [Cdev_3 Cdev_4 Cdev_5], title="Agg. Consumption, C_t", ylabel="% dev. from steady state"; args_rest...)
    p7 = plot(1:TT, [untarget["masstest"].*100 target["masstest"].*100 isolate["masstest"].*100], title="Mass tested", ylabel="%"; args_rest...)
    p8 = plot(1:TT, [untarget["mctest"] target["mctest"] isolate["mctest"] ], title="Marginal testing cost", ylabel="Dollars"; args_rest...)
    p = plot(p1, p2, p3, p4, p5, p6, p7, p8, layout = (2, 4));
    display(p)

end

# Using stored results
function createFig7()

    println("Replicating Figure 7...")

    # Load stored data
    untarget = load("src/data/calib1_untarget_small.jld2");
    target = load("src/data/calib1_target_small.jld2");
    isolate = load("src/data/calib1_isolate_small.jld2");

    # Create variables
    TT = untarget["TT"];
    IN_3 = 1 .- untarget["ST"] .- untarget["RT"];
    IN_4 = 1 .- target["ST"] .- target["RT"];
    IN_5 = 1 .- isolate["ST"] .- isolate["RT"];
    death_3 = cumsum(untarget["Par"].gamma .* untarget["Par"].delta .* IN_3);
    death_4 = cumsum(target["Par"].gamma .* target["Par"].delta .* IN_4);
    death_5 = cumsum(isolate["Par"].gamma .* isolate["Par"].delta .* IN_5);
    Yss = 1.01;
    Cdev_3 = ((untarget["YT"]-untarget["C1T"]).-Yss)./Yss.*100;
    Cdev_4 = ((target["YT"]-target["C1T"]).-Yss)./Yss.*100;
    Cdev_5 = ((isolate["YT"]-isolate["C1T"]).-Yss)./Yss.*100;

    # Figure 7
    args_first = (titlefontsize=8, guidefontsize=6, linecolor=["red" "blue" "purple"], legendfontsize=6, label=["untargeted" "targeted" "isolate"]);
    args_rest  = (titlefontsize=8, guidefontsize=6, linecolor=["red" "blue" "purple"], legend=false);
    p1 = plot(1:TT, [IN_3.*100 IN_4.*100 IN_5.*100], title="Infected, I_t", ylabel="%"; args_first...)
    p2 = plot(1:TT, [untarget["ST"].*100 target["ST"].*100 isolate["ST"].*100], title="Susceptible, S_t", ylabel="%"; args_rest...)
    p3 = plot(1:TT, [death_3.*100 death_4.*100 death_5.*100], title="Cumulative Deaths", ylabel="%"; args_rest...)
    p4 = plot(1:TT, [untarget["Rknot"] target["Rknot"] isolate["Rknot"]], title="R_0"; args_rest...)
    p5 = plot(1:TT, [untarget["IT"].*100 target["IT"].*100 isolate["IT"].*100], title="Unknown Infected, ~I_t", ylabel="%"; args_rest...)
    p6 = plot(1:TT, [Cdev_3 Cdev_4 Cdev_5], title="Agg. Consumption, C_t", ylabel="% dev. from steady state"; args_rest...)
    p7 = plot(1:TT, [untarget["masstest"].*100 target["masstest"].*100 isolate["masstest"].*100], title="Mass tested", ylabel="%"; args_rest...)
    p8 = plot(1:TT, [untarget["mctest"] target["mctest"] isolate["mctest"] ], title="Marginal testing cost", ylabel="Dollars"; args_rest...)
    p = plot(p1, p2, p3, p4, p5, p6, p7, p8, layout = (2, 4));
    display(p)

end


text:
using ChainRulesCore

function foo(x)
    a = sin(x)
    b = 0.2 + a
    c = asin(b)
    return c
end

# Define rules (alternatively get them for free via `using ChainRules`)
@scalar_rule(sin(x), cos(x))
@scalar_rule(+(x, y), (1.0, 1.0))
@scalar_rule(asin(x), inv(sqrt(1 - x^2)))

#### Find dfoo/dx via rrules
#### First the forward pass, gathering up the pullbacks
x = 3;
a, a_pullback = rrule(sin, x); # a_pullback(x) = (NoTangent(), cos(x) * x)
b, b_pullback = rrule(+, 0.2, a);
c, c_pullback = rrule(asin, b)

#### Then the backward pass calculating gradients
c̄ = 1;                    # ∂c/∂c
_, b̄ = c_pullback(c̄);     # ∂c/∂b = ∂c/∂b ⋅ ∂c/∂c
println(b̄)
println(inv(sqrt(1 - b^2)))
_, ā2, ā = b_pullback(b̄);  # ∂c/∂a = ∂c/∂b ⋅ ∂b/∂a
println((ā2, ā))
println(b̄ .* (1.0, 1.0))
_, x̄ = a_pullback(ā);     # ∂c/∂x = ∂c/∂a ⋅ ∂a/∂x
println(x̄)                # ∂c/∂x = ∂foo/∂x
# ∂c/∂x = ∂c/∂b ⋅ ∂b/∂a ⋅ ∂a/∂x
#### Find dfoo/dx via frules
x = 3;
ẋ = 1;              # ∂x/∂x
nofields = ZeroTangent();  # ∂self/∂self

a, ȧ = frule((nofields, ẋ), sin, x);                    # ∂a/∂x = ∂a/∂x ⋅ ∂x/∂x
println(ȧ)
println(cos(x))
b, ḃ = frule((nofields, ZeroTangent(), ȧ), +, 0.2, a);  # ∂b/∂x = ∂b/∂a ⋅ ∂a/∂x

c, ċ = frule((nofields, ḃ), asin, b);                   # ∂c/∂x = ∂c/∂b ⋅ ∂b/∂x
ċ                                                       # ∂c/∂x = ∂foo/∂x
# ∂c/∂x = ∂c/∂b ⋅ ∂b/∂a ⋅ ∂a/∂x
#### Find dfoo/dx via FiniteDifferences.jl
using FiniteDifferences
central_fdm(5, 1)(foo, x)
# output
-1.0531613736418257

#### Find dfoo/dx via ForwardDiff.jl
using ForwardDiff
ForwardDiff.derivative(foo, x)
# output
-1.0531613736418153

#### Find dfoo/dx via Zygote.jl
using Zygote
Zygote.gradient(foo, x)

max_stars_repo_name: sandreza/GradientLearning

max_issues_repo_name: sandreza/GradientLearning

max_forks_repo_name: sandreza/GradientLearning


text:
using Downloads: download
using BenchmarkTools

function create_template_if_missing(filename)
    isfile(filename) && return
    write(filename, """function part1(input)

                       end

                       function part2(input)

                       end""")
end

args = copy(ARGS)
benchmark_arg = filter(startswith("--benchmark"), args)
filter!(!startswith("--benchmark"), args)
length(args) >= 2 || error("Usage: julia aoc.jl YEAR DAY [TEST] [EXTRA_ARGS] [--benchmark=SUFFIX]")
year, day = args[1:2]
test = length(args) > 2
input_file = joinpath(@__DIR__, year, "input", "day$(day)")
if !isfile(input_file)
    mkpath(dirname(input_file))
    println("Downloading input for $year day $day.")
    session_cookie_file = joinpath(@__DIR__, "session_cookie")
    if isfile(session_cookie_file)
        session_cookie = read(session_cookie_file, String)
        download("https://adventofcode.com/$(year)/day/$(day)/input",
                 input_file,
                 headers = Dict("cookie" => "session=$(session_cookie)"))
    else
        println("No session cookie available. Download the input file manually and save it to $(input_file)")
    end
end
code_filename = "day$(day).jl"
create_template_if_missing(joinpath(@__DIR__, year, code_filename))
if !isempty(benchmark_arg) && occursin("=", only(benchmark_arg))
    suffix = last(split(only(benchmark_arg), "="))
    code_filename = "day$(day)$(suffix).jl"
end
include(joinpath(@__DIR__, year, code_filename))
if test
    input_file *= args[3]
end
for part in (part1, part2)
    part == part2 && println()
    print("Part ", last(string(part)), ": ")
    test && print(args[3])
    println()
    data = read(input_file)
    println(part(IOBuffer(data), args[4:end]...))
    if !isempty(benchmark_arg)
        length(args) < 4 || error("EXTRA_ARGS are not supported for benchmarking (due to interference with measurement).")
        @btime $(part)(input) setup = (input = IOBuffer($data)) evals = 1
    end
end


text:
## Exercise 8-2
## Modify the program to fix this error.
prefixes = "JKLMNOPQ"
suffix = "ack"

println("Ans.")

for letter in prefixes
    if letter == 'O' || letter == 'Q'
        letter *= 'u'
    end
    println(letter * suffix)
end

println("End.")


text:
using LinearAlgebra
#
# This script should return a mxn Matrix{Float64}.
# Example:
# return [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0];
#

m = 100
n = m

A = rand(m, n);

println("Julia: sum(A)=", sum(A));
return A

max_stars_repo_name: yundantianchang/dynload-julia

max_issues_repo_name: yundantianchang/dynload-julia

max_forks_repo_name: yundantianchang/dynload-julia


text:
"""
    addColOverflow!(df::DataFrame, ColDuration::String, NewCol::String, ColNum::Int64)

Adds a column of 0's and 1's indicating the absence or presence of overflow using overflows' duration.

Data provided must not contain any missing values for overflows' duration.

# Arguments

- `df::DataFrame`: The dataframe containing the data.
- `ColDuration::String`: The name of the column of `df` containing the duration of overflows.
- `NewCol::String`: The name of the new column of `df` that will contain the 0's and 1's.
- `ColNum::Int64`: The column number of `NewCol`.
"""
function addColOverflow!(df::DataFrame, ColDuration::String, NewCol::String, ColNum::Int64)
    insertcols!(df, ColNum, Symbol(NewCol) => Int.(df[:, Symbol(ColDuration)] .> 0.0))
end

"""
    countOverflow(df::DataFrame, ColOverflow::String, Print::Bool=true)

Returns the number of non-overflows, of overflows, and the percentage of overflows.

# Arguments

- `df::DataFrame`: The dataframe containing the data.
- `ColOverflow::String`: The name of the column of `df` containing the binary vector of overflows.
- `Print::Bool=true`: Whether the number of non-overflows, overflows and percentage of overflows are printed or not. Default = true.
"""
function countOverflow(df::DataFrame, ColOverflow::String, Print::Bool=true)
    NO = countmap(df[:, Symbol(ColOverflow)])[0]
    O = countmap(df[:, Symbol(ColOverflow)])[1]
    Per = O/(NO + O)*100

    if Print
        println("The number of non-overflows equals $(NO), the number of overflows equals $O, the percentage of overflows equals around $(round(Per, digits=2)) %.")
    end

    return NO, O, Per
end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule protoc_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("protoc")
JLLWrappers.@generate_main_file("protoc", UUID("c7845625-083e-5bbe-8504-b32d602b7110"))
end  # module protoc_jll


text:
export Analytical, CG,
        Newton, NewtonCG,
        LBFGS,
        ProxGrad, FISTA, ISTA,
        IWLSCG,
        ADMM, FADMM # NOTE: these do not work currently

# =====
# TODO
# * all - pick linesearch
# * NewtonCG number of inner iter
# * FISTA field to enforce descent
# ====

abstract type Solver end

# ===================== analytical.jl

@with_kw struct Analytical <: Solver
    iterative::Bool = false
    max_inner::Int  = 200
end

CG() = Analytical(; iterative=true)

# ===================== newton.jl

struct Newton <: Solver end

struct NewtonCG <: Solver end

struct LBFGS <: Solver end

# struct BFGS <: Solver end

# ===================== pgrad.jl

@with_kw struct ProxGrad <: Solver
    accel::Bool    = false # use Nesterov style acceleration (see also FISTA)
    max_iter::Int  = 1000  # max number of overall iterations
    tol::Float64   = 1e-4  # tolerance over relative change of θ i.e. norm(θ-θ_)/norm(θ)
    max_inner::Int = 100   # β^max_inner should be > 1e-10
    beta::Float64  = 0.8   # in (0, 1); shrinkage in the backtracking step
end

FISTA(; kwa...) = ProxGrad(;accel = true, kwa...)
ISTA(; kwa...)  = ProxGrad(;accel = false, kwa...)

# ===================== iwls.jl

@with_kw struct IWLSCG <: Solver
    max_iter::Int      = 100
    max_inner::Int     = 200
    tol::Float64       = 1e-4
    damping::Float64   = 1.0   # should be between 0 and 1, 1 = trust iterates
    threshold::Float64 = 1e-6  # threshold for the residuals used for instance in quantile reg
end

# ===================== admm.jl

# @with_kw struct ADMM <: Solver
#     max_iter::Int  = 100
#     tol::Float64   = 1e-4
#     alpha::Float64 = 1.5  # over-relaxation (recommended between 1.5 and 1.8)
#     rho::Float64   = 1.0  # Lagrangian parameter (should be decreased if poor condition)
# end
#
# @with_kw struct FADMM <: Solver
#     max_iter::Int  = 100
#     tol::Float64   = 1e-4
#     eta::Float64   = 0.999 # restart parameter
#     rho::Float64   = 1.0   # Lagrangian parameter (should be decreased if poor condition)
#     tau::Float64   = 2.0   # Increase / decrease of ρ should be > 1
#     mu::Float64    = 10.0  #
# end


text:
function c_R_to_G(Ns::Array{Int,1}, fR)
  out = zeros(Complex128,Ns[1],Ns[2],Ns[3])
  ccall( (:fftw_fw_fft3d, "fft3d"), Void,
    (Ptr{Complex128}, Ptr{Complex128}, Int64, Int64, Int64), fR, out,
    Ns[3],Ns[2],Ns[1] )
  return out
end

function c_G_to_R(Ns::Array{Int,1}, fG)
  out = zeros(Complex128,Ns[1],Ns[2],Ns[3])
  ccall( (:fftw_inv_fft3d, "fft3d"), Void,
    (Ptr{Complex128}, Ptr{Complex128}, Int64, Int64, Int64), fG, out,
    Ns[3],Ns[2],Ns[1] )
  return out
end

# generic, one-column version
function R_to_G( Ns::Array{Int,1}, fR )
  out = reshape( fft( reshape(fR,Ns[1],Ns[2],Ns[3]) ),
                   size(fR) )
end

function G_to_R( Ns::Array{Int,1}, fG )
  out = reshape( ifft( reshape(fG,Ns[1],Ns[2],Ns[3]) ),
                 size(fG) )
  return out
end

function test_main()
  Ns = [150,150,150]
  in1 = rand(Complex128,Ns[1],Ns[2],Ns[3])

  @time out1c = c_R_to_G(Ns,in1)
  @time in2c = c_G_to_R(Ns,out1c)

  @time out1 = R_to_G(Ns,in1)
  @time in2 = G_to_R(Ns,out1)

  println("diff out1 = ", sum( abs(out1c - out1) ))
  println("diff out1 = ", sum( abs(in2c - in2) ))
end

test_main()


text:
@testset "OpenCL.Program" begin

    test_source = "
    __kernel void sum(__global const float *a,
                      __global const float *b,
                      __global float *c)
    {
      uint gid = get_global_id(0);
      c[gid] = a[gid] + b[gid];
    }
    "

    function create_test_program()
        ctx = cl.create_some_context()
        cl.Program(ctx, source=test_source)
    end

    @testset "OpenCL.Program source constructor" begin
        for device in cl.devices()
            ctx = cl.Context(device)
            prg = cl.Program(ctx, source=test_source)
            @test prg != nothing
        end
    end
    @testset "OpenCL.Program info" begin
        for device in cl.devices()
            ctx = cl.Context(device)
            prg = cl.Program(ctx, source=test_source)

            @test prg[:context] == ctx

            @test typeof(prg[:devices]) == Vector{cl.Device}
            @test length(prg[:devices]) > 0
            @test device in prg[:devices]

            @test typeof(prg[:source]) == String
            @test prg[:source] == test_source

            @test prg[:reference_count] > 0
            @test isempty(strip(prg[:build_log][device]))
         end
    end

    @testset "OpenCL.Program build" begin
        for device in cl.devices()
            ctx = cl.Context(device)
            prg = cl.Program(ctx, source=test_source)
            @test cl.build!(prg) != nothing

            # BUILD_SUCCESS undefined in POCL implementation..
            if device[:platform][:name] == "Portable Computing Language"
                @warn("Skipping OpenCL.Program build for Portable Computing Language Platform")
                continue
            end
            @test prg[:build_status][device] == cl.CL_BUILD_SUCCESS

            # test build by methods chaining
            @test prg[:build_status][device] == cl.CL_BUILD_SUCCESS
            if device[:platform][:name] != "Intel(R) OpenCL"
                # The intel CPU driver is very verbose on Linux and output
                # compilation status even without any warnings
                @test isempty(strip(prg[:build_log][device]))
            end
        end
    end

    @testset "OpenCL.Program source code" begin
        for device in cl.devices()
           ctx = cl.Context(device)
           prg = cl.Program(ctx, source=test_source)
           @test prg[:source] == test_source
        end
    end

    @testset "OpenCL.Program binaries" begin
        for device in cl.devices()
            ctx = cl.Context(device)
            prg = cl.Program(ctx, source=test_source) |> cl.build!

            @test device in collect(keys(prg[:binaries]))
            binaries = prg[:binaries]
            @test device in collect(keys(binaries))
            @test binaries[device] != nothing
            @test length(binaries[device]) > 0
            prg2 = cl.Program(ctx, binaries=binaries)
            @test prg2[:binaries] == binaries
            try
                prg2[:source]
                error("should not happen")
            catch err
                @test isa(err, cl.CLError)
                @test err.code == -45
                @test err.desc == :CL_INVALID_PROGRAM_EXECUTABLE
            end
        end
    end
end


text:
# symmetrization

act(mon :: SoSMonom, perm :: Dict{Symbol,Symbol}) =
    [ get(perm,k,k) => v for (k,v) in mon ] :: SoSMonom

# XXX currently unused
act(poly :: SoSPoly, perm :: Dict{Symbol,Symbol}) =
    [ act(k,perm) => v for (k,v) in poly ] :: SoSPoly

# maps monomials to their orbits; returns maps in both directions
# note: this is roughly order-preserving, at least to the extent that the first
# monomial given will lie in orbit 1, which is relied on in sossolve().
function monom_orbits(monoms :: Array{SoSMonom}, genperms :: Array{Dict{Symbol,Symbol}})
    ret1 = Dict{SoSMonom,Int64}()
    ret2 = Set{SoSMonom}[]

    orbitidx = 0

    seen = Set{SoSMonom}()
    for b in monoms
        if in(b, seen)
            continue
        end

        orbitidx += 1
        orbit = Set{SoSMonom}()
        q = SoSMonom[]

        # start off with b
        push!(orbit, b)
        push!(q, b)

        # explore the orbit
        while !isempty(q)
            c = shift!(q)
            ret1[c] = orbitidx # here's the actual orbit-mapping assignment

            # act every possible way on the current element
            for perm in genperms
                i = act(c,perm)

                if !in(i,orbit)
                    push!(orbit,i)
                    push!(q,i)
                end
            end
        end

        union!(seen,orbit)
        push!(ret2, orbit)
    end

    ret1, ret2
end

symmetrize!(prog :: Program, genperms :: Array{Dict{Symbol,Symbol}}) = append!(prog.symmetries, genperms)
symmetrize!(prog :: Program, genperm :: Dict{Symbol,Symbol}) = push!(prog.symmetries, genperm)

# symmetry according to dihedral group (rotation and reflection)
function symmetrize_dihedral!(prog :: Program, cycles...)
    n = length(cycles[1])

    symmetrize_cyclic!(prog, cycles...)

    if n >= 2
        reflect = Dict{Symbol,Symbol}()
        for cycle in cycles
            for i in 1:n
                reflect[cycle[i]] = cycle[n-i+1]
            end
        end
        symmetrize!(prog, reflect)
    end
end

# symmetry according to cyclic group (rotation only)
function symmetrize_cyclic!(prog :: Program, cycles...)
    n = length(cycles[1])
    if n < 2 return end

    rotate = Dict{Symbol,Symbol}()
    for cycle in cycles
        for i in 0:(n-1)
            j = (i+1) % n
            rotate[cycle[i+1]] = cycle[j+1]
        end
    end
    symmetrize!(prog, rotate)
end

# symmetry according to hyperoctahedral group
# warning: bit hacks
function symmetrize_hypercube!(prog :: Program, cubes...)
    n = length(cubes[1])
    if (n & (n - 1)) != 0 # n is not a power of 2
        throw(ArgumentException("argument list in symmetrize_hypercube must be a power of 2"))
    end

    # reflection along first coordinate
    if n >= 2
        refl = [ i => (i$1) for i in 0:(n-1) ] :: Dict{Symbol,Symbol}
        symmetrize!(prog, refl)
    end

    # transposition of first two coordinates
    if n >= 4
        trans = Dict{Symbol,Symbol}()
        for cube in cubes
            for i in 0:(n-1)
                bot = i & 3
                j = (bot == 1 || bot == 2) ? i$3 : i
                trans[cube[i+1]] = cube[j+1]
            end
        end
        symmetrize!(prog, trans)
    end

    # cyclic permutation of all coordinates
    if n >= 8
        cycle = Dict{Symbol,Symbol}()
        for cube in cubes
            for i in 0:(n-1)
                j = i<<1
                if ( (j & n) != 0)
                    j $= n
                    j $= 1
                end
                cycle[cube[i+1]] = cube[j+1]
            end
        end
        symmetrize!(prog, cycle)
    end
end

# symmetry according to full symmetric group
function symmetrize_full!(prog :: Program, symsets...)
    n = length(symsets[1])
    if n < 2 return end

    symmetrize!(prog, Dict{Symbol,Symbol}(set[1] => set[2], set[2] => set[1]))

    if n >= 3
        symmetrize_cyclic!(prog, symsets...)
    end
end

max_stars_repo_name: willperry/SumOfSquaresOptimization


text:
import Base: inv
import LinearAlgebra: det, diag, logdet, mul!, lmul!, rmul!, rdiv!
using LinearAlgebra, FillArrays

####### linear algebra  ######
inv(M::IMatrix) = M
det(M::IMatrix) = 1
diag(M::IMatrix{N,T}) where {N,T} = ones(T, N)
logdet(M::IMatrix) = 0
Base.sqrt(x::PermMatrix) = sqrt(Matrix(x))
Base.sqrt(x::IMatrix) = x
Base.exp(x::PermMatrix) = exp(Matrix(x))
Base.exp(x::IMatrix) = ℯ * x

#det(M::PermMatrix) = parity(M.perm)*prod(M.vals)
function inv(M::PermMatrix)
    new_perm = fast_invperm(M.perm)
    return PermMatrix(new_perm, 1.0 ./ M.vals[new_perm])
end

####### multiply ###########
*(A::IMatrix{N}, B::AbstractVector) where {N} = size(A, 2) == size(B, 1) ? B :
    throw(DimensionMismatch("matrix A has dimensions $(size(A)), matrix B has dimensions $((size(B, 1), 1))"))

for MATTYPE in
    [:AbstractMatrix, :StridedMatrix, :Diagonal, :SparseMatrixCSC, :Matrix, :PermMatrix]
    @eval *(A::IMatrix{N}, B::$MATTYPE) where {N} = N == size(B, 1) ? B :
        throw(DimensionMismatch("matrix A has dimensions $(size(A)), matrix B has dimensions $(size(B))"))

    @eval *(A::$MATTYPE, B::IMatrix{N}) where {N} = size(A, 2) == N ? A :
        throw(DimensionMismatch("matrix A has dimensions $(size(A)), matrix B has dimensions $(size(B))"))
end

*(A::Adjoint{T,<:AbstractVector{T}}, B::IMatrix) where {T} = size(A, 2) == size(B, 1) ? A :
    throw(DimensionMismatch("matrix A has dimensions $(size(A)), matrix B has dimensions $(size(B))"))

*(A::IMatrix, B::IMatrix) = size(A, 2) == size(B, 1) ? A :
    throw(DimensionMismatch("matrix A has dimensions $(size(A)), matrix B has dimensions $(size(B))"))

########## Multiplication #############

# NOTE: making them dry?
# to vector
function *(A::PermMatrix{Ta}, X::AbstractVector{Tx}) where {Ta,Tx}
    nX = length(X)
    nX == size(A, 2) || throw(DimensionMismatch())
    v = similar(X, promote_type(Ta, Tx))
    @simd for i = 1:nX
        @inbounds v[i] = A.vals[i] * X[A.perm[i]]
    end
    v
end

function *(X::Adjoint{Tx,<:AbstractVector{Tx}}, A::PermMatrix{Ta}) where {Tx,Ta}
    nX = length(X)
    nX == size(A, 1) || throw(DimensionMismatch())
    v = similar(X, promote_type(Tx, Ta))
    @simd for i = 1:nX
        @inbounds v[A.perm[i]] = A.vals[i] * X[i]
    end
    v
end

# to diagonal
function *(D::Diagonal{Td}, A::PermMatrix{Ta}) where {Td,Ta}
    T = Base.promote_op(*, Td, Ta)
    PermMatrix(A.perm, A.vals .* D.diag)
end

function *(A::PermMatrix{Ta}, D::Diagonal{Td}) where {Td,Ta}
    T = Base.promote_op(*, Td, Ta)
    PermMatrix(A.perm, A.vals .* view(D.diag, A.perm))
end

# to self
function *(A::PermMatrix, B::PermMatrix)
    size(A, 1) == size(B, 1) || throw(DimensionMismatch())
    PermMatrix(B.perm[A.perm], A.vals .* view(B.vals, A.perm))
end

# to matrix
function *(A::PermMatrix, X::AbstractMatrix)
    size(X, 1) == size(A, 2) || throw(DimensionMismatch())
    return @views A.vals .* X[A.perm, :]   # this may be inefficient for sparse CSC matrix.
end

function *(X::AbstractMatrix, A::PermMatrix)
    mX, nX = size(X)
    nX == size(A, 1) || throw(DimensionMismatch())
    return @views (A.vals'.*X)[:, fast_invperm(A.perm)]
end

# NOTE: this is just a temperory fix for v0.7. We should overload mul! in
# the future (when we start to drop v0.6) to enable buildin lazy evaluation.

*(x::Adjoint{<:Any,<:AbstractVector}, D::PermMatrix) = Matrix(x) * D
*(x::Transpose{<:Any,<:AbstractVector}, D::PermMatrix) = Matrix(x) * D
*(A::Adjoint{<:Any,<:AbstractArray}, D::PermMatrix) = Adjoint(adjoint(D) * parent(A))
*(A::Transpose{<:Any,<:AbstractArray}, D::PermMatrix) = Transpose(transpose(D) * parent(A))
*(A::Adjoint{<:Any,<:PermMatrix}, D::PermMatrix) = adjoint(parent(A)) * D
*(A::Transpose{<:Any,<:PermMatrix}, D::PermMatrix) = transpose(parent(A)) * D
*(A::PermMatrix, D::Adjoint{<:Any,<:PermMatrix}) = A * adjoint(parent(D))
*(A::PermMatrix, D::Transpose{<:Any,<:PermMatrix}) = A * transpose(parent(D))

# for MAT in [:AbstractArray, :Matrix, :SparseMatrixCSC, :PermMatrix]
#     @eval begin
#         *(A::Adjoint{<:Any, <:$MAT}, D::PermMatrix) = copy(A) * D
#         *(A::Transpose{<:Any, <:$MAT}, D::PermMatrix) = copy(A) * D
#         *(A::PermMatrix, D::Adjoint{<:Any, <:$MAT}) = A * copy(D)
#         *(A::PermMatrix, D::Transpose{<:Any, <:$MAT}) = A * copy(D)
#     end
# end

############### Transpose, Adjoint for IMatrix ###############
for MAT in
    [:AbstractArray, :AbstractVector, :Matrix, :SparseMatrixCSC, :PermMatrix, :IMatrix]
    @eval *(A::Adjoint{<:Any,<:$MAT}, D::IMatrix) = Adjoint(D * parent(A))
    @eval *(A::Transpose{<:Any,<:$MAT}, D::IMatrix) = Transpose(D * parent(A))
    if MAT != :AbstactVector
        @eval *(A::IMatrix, D::Transpose{<:Any,<:$MAT}) = Transpose(parent(D) * A)
        @eval *(A::IMatrix, D::Adjoint{<:Any,<:$MAT}) = Adjoint(parent(D) * A)
    end
end

# to sparse
function *(A::PermMatrix, X::SparseMatrixCSC)
    nA = size(A, 1)
    mX, nX = size(X)
    mX == nA || throw(DimensionMismatch())
    perm = fast_invperm(A.perm)
    nzval = similar(X.nzval)
    rowval = similar(X.rowval)
    @inbounds for j = 1:nX
        @inbounds for k = X.colptr[j]:X.colptr[j+1]-1
            r = perm[X.rowval[k]]
            nzval[k] = X.nzval[k] * A.vals[r]
            rowval[k] = r
        end
    end
    SparseMatrixCSC(mX, nX, X.colptr, rowval, nzval)
end

function *(X::SparseMatrixCSC, A::PermMatrix)
    nA = size(A, 1)
    mX, nX = size(X)
    nX == nA || throw(DimensionMismatch())
    perm = fast_invperm(A.perm)
    nzval = similar(X.nzval)
    colptr = similar(X.colptr)
    rowval = similar(X.rowval)
    colptr[1] = 1
    z = 1
    @inbounds for j = 1:nA
        pk = perm[j]
        va = A.vals[pk]
        @inbounds @simd for k = X.colptr[pk]:X.colptr[pk+1]-1
            nzval[z] = X.nzval[k] * va
            rowval[z] = X.rowval[k]
            z += 1
        end
        colptr[j+1] = z
    end
    SparseMatrixCSC(mX, nX, colptr, rowval, nzval)
end

############ SparseMatrixCSC and Diagonal #############
function mul!(Y::SparseMatrixCSC, X::SparseMatrixCSC, A::Diagonal{T,Vector{T}}) where {T}
    nA = size(A, 1)
    mX, nX = size(X)
    nX == nA || throw(DimensionMismatch())
    j = 1
    va = A.diag[1]
    @inbounds for k = 1:nnz(X)
        if k == X.colptr[j]
            va = A.diag[j]
            j = j + 1
        end
        Y.nzval[k] *= va
    end
    Y
end

function mul!(Y::SparseMatrixCSC, A::Diagonal{T,Vector{T}}, X::SparseMatrixCSC) where {T}
    nA = size(A, 2)
    mX, nX = size(X)
    mX == nA || throw(DimensionMismatch())
    @inbounds @simd for i = 1:nnz(X)
        Y.nzval[i] *= A.diag[X.rowval[i]]
    end
    Y
end

rmul!(A::SparseMatrixCOO, B::Int) = (A.vs *= B; A)
lmul!(B::Int, A::SparseMatrixCOO) = (A.vs *= B; A)
rdiv!(A::SparseMatrixCOO, B::Int) = (A.vs /= B; A)

*(A::SparseMatrixCOO, B::Int) = rmul!(copy(A), B)
*(B::Int, A::SparseMatrixCOO) = lmul!(B, copy(A))
/(A::SparseMatrixCOO, B::Int) = rdiv!(copy(A), B)

-(ii::IMatrix) = (-1) * ii
-(pm::PermMatrix) = (-1) * pm

for FUNC in [:randn!, :rand!]
    @eval function Random.$FUNC(m::Diagonal)
        $FUNC(m.diag)
        return m
    end

    @eval function Random.$FUNC(m::SparseMatrixCSC)
        $FUNC(m.nzval)
        return m
    end

    @eval function Random.$FUNC(m::PermMatrix)
        $FUNC(m.vals)
        return m
    end
end


text:
using Test, Documenter
using IndividualDisplacements, MeshArrays, OrdinaryDiffEq, OceanStateEstimation
include(joinpath(dirname(pathof(IndividualDisplacements)),"../test/helper_functions.jl"))

IndividualDisplacements.get_ecco_velocity_if_needed();
IndividualDisplacements.get_occa_velocity_if_needed();

@testset "test1" begin
    uvetc,sol=test1_setup()
    @test isapprox(sol[1,end],23.4474; atol=0.01)
    @test isapprox(sol[2,end],8.0896; atol=0.01)
end

@testset "test2" begin
    df,𝑃=test2_periodic_domain()
    @test prod(isapprox.(df[end-35:6:end,:y],12*(0.4:0.04:0.6)))
    @test prod(isapprox.(df[end-5:end,:x],12*(0.4:0.04:0.6).+4.0))
end

@testset "doctests" begin
    doctest(IndividualDisplacements; manual = false)
end


text:
import FFTW: fft, ifft, fft!, ifft!, rfft, irfft
export fft, fft1, fft!, ifft!, rfft, irfft, fft1!, fft2, fft2!, fft3, fft3!, ifft1, ifft2, ifft3, rfft1, rfft2, rfft3
export ifft, irfft2, irfft3, ifft1!, ifft2!, ifft3!
export fft_convolve1, fft_convolve2, fft_convolve3
export set_fft_plan_cache_size

function fft(_in::AFArray{T,1},norm_factor=1.0,odim0::dim_t=0) where {T<:Complex}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t),out,_in.arr,Cdouble(norm_factor),odim0))
    AFArray{T,1}(out[])
end

function fft(_in::AFArray{T,1},norm_factor=1.0,odim0::dim_t=0) where {T<:Real}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t),out,_in.arr,Cdouble(norm_factor),odim0))
    AFArray{Complex{T},1}(out[])
end

function fft1(_in::AFArray{T,N},norm_factor=1.0,odim0::dim_t=0) where {T<:Complex,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t),out,_in.arr,Cdouble(norm_factor),odim0))
    AFArray{T,N}(out[])
end

function fft1(_in::AFArray{T,N},norm_factor=1.0,odim0::dim_t=0) where {T<:Real,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t),out,_in.arr,Cdouble(norm_factor),odim0))
    AFArray{Complex{T},N}(out[])
end

function fft!(_in::AFVector,norm_factor=1.0)
    _error(ccall((:af_fft_inplace,af_lib),af_err,(af_array,Cdouble),_in.arr,Cdouble(norm_factor)))
    _in
end

function fft1!(_in::AFArray,norm_factor=1.0)
    _error(ccall((:af_fft_inplace,af_lib),af_err,(af_array,Cdouble),_in.arr,Cdouble(norm_factor)))
    _in
end

function fft(_in::AFArray{T,2},norm_factor=1.0,odim0::dim_t=0,odim1::dim_t=0) where {T<:Complex}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft2,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),odim0,odim1))
    AFArray{T,2}(out[])
end

function fft(_in::AFArray{T,2},norm_factor=1.0,odim0::dim_t=0,odim1::dim_t=0) where {T<:Real}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft2,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),odim0,odim1))
    AFArray{Complex{T},2}(out[])
end

function fft2(_in::AFArray{T,N},norm_factor=1.0,odim0::dim_t=0,odim1::dim_t=0) where {T<:Complex,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft2,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),odim0,odim1))
    AFArray{T,N}(out[])
end

function fft2(_in::AFArray{T,N},norm_factor=1.0,odim0::dim_t=0,odim1::dim_t=0) where {T<:Real,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft2,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),odim0,odim1))
    AFArray{Complex{T},N}(out[])
end

function fft!(_in::AFMatrix,norm_factor=1.0)
    _error(ccall((:af_fft2_inplace,af_lib),af_err,(af_array,Cdouble),_in.arr,Cdouble(norm_factor)))
    _in
end

function fft2!(_in::AFArray,norm_factor=1.0)
    _error(ccall((:af_fft2_inplace,af_lib),af_err,(af_array,Cdouble),_in.arr,Cdouble(norm_factor)))
    _in
end

function fft(_in::AFArray{T,3},norm_factor=1.0,odim0::dim_t=0,odim1::dim_t=0,odim2::dim_t=0) where {T<:Complex}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft3,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),odim0,odim1,odim2))
    AFArray{T,3}(out[])
end

function fft(_in::AFArray{T,3},norm_factor=1.0,odim0::dim_t=0,odim1::dim_t=0,odim2::dim_t=0) where {T<:Real}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft3,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),odim0,odim1,odim2))
    AFArray{Complex{T},3}(out[])
end

function fft3(_in::AFArray{T,N},norm_factor=1.0,odim0::dim_t=0,odim1::dim_t=0,odim2::dim_t=0) where {T<:Complex,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft3,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),odim0,odim1,odim2))
    AFArray{T,N}(out[])
end

function fft3(_in::AFArray{T,N},norm_factor=1.0,odim0::dim_t=0,odim1::dim_t=0,odim2::dim_t=0) where {T<:Real,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft3,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),odim0,odim1,odim2))
    AFArray{Complex{T},N}(out[])
end

function fft!(_in::AFVolume,norm_factor=1.0)
    _error(ccall((:af_fft3_inplace,af_lib),af_err,(af_array,Cdouble),_in.arr,Cdouble(norm_factor)))
    _in
end

function fft3!(_in::AFArray,norm_factor=1.0)
    _error(ccall((:af_fft3_inplace,af_lib),af_err,(af_array,Cdouble),_in.arr,Cdouble(norm_factor)))
    _in
end

function ifft(_in::AFArray{T,1},norm_factor=-1.0,odim0::dim_t=0) where {T}
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    out = RefValue{af_array}(0)
    _error(ccall((:af_ifft,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t),out,_in.arr,Cdouble(norm_factor),odim0))
    AFArray{T,1}(out[])
end

function ifft1(_in::AFArray{T,N},norm_factor=-1.0,odim0::dim_t=0) where {T,N}
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    out = RefValue{af_array}(0)
    _error(ccall((:af_ifft,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t),out,_in.arr,Cdouble(norm_factor),odim0))
    AFArray{T,N}(out[])
end

function ifft!(_in::AFVector,norm_factor=-1.0)
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    _error(ccall((:af_ifft_inplace,af_lib),af_err,(af_array,Cdouble),_in.arr,Cdouble(norm_factor)))
    _in
end

function ifft1!(_in::AFArray,norm_factor=-1.0)
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    _error(ccall((:af_ifft_inplace,af_lib),af_err,(af_array,Cdouble),_in.arr,Cdouble(norm_factor)))
    _in
end

function ifft(_in::AFArray{T,2},norm_factor=-1.0,odim0::dim_t=0,odim1::dim_t=0) where {T}
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    out = RefValue{af_array}(0)
    _error(ccall((:af_ifft2,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),odim0,odim1))
    AFArray{T,2}(out[])
end

function ifft2(_in::AFArray{T,N},norm_factor=-1.0,odim0::dim_t=0,odim1::dim_t=0) where {T,N}
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    out = RefValue{af_array}(0)
    _error(ccall((:af_ifft2,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),odim0,odim1))
    AFArray{T,N}(out[])
end

function ifft!(_in::AFMatrix,norm_factor=-1.0)
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    _error(ccall((:af_ifft2_inplace,af_lib),af_err,(af_array,Cdouble),_in.arr,Cdouble(norm_factor)))
    _in
end

function ifft2!(_in::AFArray,norm_factor=-1.0)
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    _error(ccall((:af_ifft2_inplace,af_lib),af_err,(af_array,Cdouble),_in.arr,Cdouble(norm_factor)))
    _in
end

function ifft(_in::AFArray{T,3},norm_factor=-1.0,odim0::dim_t=0,odim1::dim_t=0,odim2::dim_t=0) where {T}
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    out = RefValue{af_array}(0)
    _error(ccall((:af_ifft3,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),odim0,odim1,odim2))
    AFArray{T,3}(out[])
end

function ifft3(_in::AFArray{T,N},norm_factor=-1.0,odim0::dim_t=0,odim1::dim_t=0,odim2::dim_t=0) where {T,N}
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    out = RefValue{af_array}(0)
    _error(ccall((:af_ifft3,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),odim0,odim1,odim2))
    AFArray{T,N}(out[])
end

function ifft!(_in::AFVolume,norm_factor=-1.0)
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    _error(ccall((:af_ifft3_inplace,af_lib),af_err,(af_array,Cdouble),_in.arr,Cdouble(norm_factor)))
    _in
end

function ifft3!(_in::AFArray,norm_factor=-1.0)
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    _error(ccall((:af_ifft3_inplace,af_lib),af_err,(af_array,Cdouble),_in.arr,Cdouble(norm_factor)))
    _in
end

function rfft(_in::AFArray{T,1},norm_factor=1.0,pad0::dim_t=0) where {T<:Real}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft_r2c,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t),out,_in.arr,Cdouble(norm_factor),pad0))
    AFArray{Complex{T},1}(out[])
end

function rfft1(_in::AFArray{T,N},norm_factor=1.0,pad0::dim_t=0) where {T<:Real,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft_r2c,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t),out,_in.arr,Cdouble(norm_factor),pad0))
    AFArray{Complex{T},N}(out[])
end

function rfft(_in::AFArray{T,2},norm_factor=1.0,pad0::dim_t=0,pad1::dim_t=0) where {T<:Real}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft2_r2c,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),pad0,pad1))
    AFArray{Complex{T},2}(out[])
end

function rfft2(_in::AFArray{T,N},norm_factor=1.0,pad0::dim_t=0,pad1::dim_t=0) where {T<:Real,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft2_r2c,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),pad0,pad1))
    AFArray{Complex{T},N}(out[])
end

function rfft(_in::AFArray{T,3},norm_factor=1.0,pad0::dim_t=0,pad1::dim_t=0,pad2::dim_t=0) where {T<:Real}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft3_r2c,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),pad0,pad1,pad2))
    AFArray{Complex{T},3}(out[])
end

function rfft3(_in::AFArray{T,N},norm_factor=1.0,pad0::dim_t=0,pad1::dim_t=0,pad2::dim_t=0) where {T<:Real,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft3_r2c,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,dim_t,dim_t,dim_t),out,_in.arr,Cdouble(norm_factor),pad0,pad1,pad2))
    AFArray{Complex{T},N}(out[])
end

function irfft(_in::AFArray{Complex{T},N},norm_factor=-1.0,is_odd::Bool=false) where {T,N}
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft_c2r,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,Bool),out,_in.arr,Cdouble(norm_factor),is_odd))
    AFArray{T,N}(out[])
end

function irfft2(_in::AFArray{Complex{T},N},norm_factor=-1.0,is_odd::Bool=false) where {T,N}
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft2_c2r,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,Bool),out,_in.arr,Cdouble(norm_factor),is_odd))
    AFArray{T,N}(out[])
end

function irfft3(_in::AFArray{Complex{T},N},norm_factor=-1.0,is_odd::Bool=false) where {T,N}
    if norm_factor < 0
        norm_factor = 1/length(_in)
    end
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft3_c2r,af_lib),af_err,(Ptr{af_array},af_array,Cdouble,Bool),out,_in.arr,Cdouble(norm_factor),is_odd))
    AFArray{T,N}(out[])
end

function fft_convolve1(signal::AFArray,filter::AFArray,mode::af_conv_mode)
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft_convolve1,af_lib),af_err,(Ptr{af_array},af_array,af_array,af_conv_mode),out,signal.arr,filter.arr,mode))
    AFArray!(out[])
end

function fft_convolve2(signal::AFArray,filter::AFArray,mode::af_conv_mode)
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft_convolve2,af_lib),af_err,(Ptr{af_array},af_array,af_array,af_conv_mode),out,signal.arr,filter.arr,mode))
    AFArray!(out[])
end

function fft_convolve3(signal::AFArray,filter::AFArray,mode::af_conv_mode)
    out = RefValue{af_array}(0)
    _error(ccall((:af_fft_convolve3,af_lib),af_err,(Ptr{af_array},af_array,af_array,af_conv_mode),out,signal.arr,filter.arr,mode))
    AFArray!(out[])
end

function set_fft_plan_cache_size(cache_size)
    _error(ccall((:af_set_fft_plan_cache_size,af_lib),af_err,(Csize_t,),Csize_t(cache_size)))
end


text:
import LinearAlgebra: dot

function intersect(triangles::Vector{Triangle}, plane::Plane)::Vector{LineSegment}
    filter(l -> l !== nothing, map(t -> intersect(t, plane), triangles))
end

function intersect(triangle::Triangle, plane::Plane)::Union{LineSegment, Nothing}
    intersections = Set(union(intersect(triangle.a, triangle.b, plane),
                              intersect(triangle.b, triangle.c, plane),
                              intersect(triangle.c, triangle.a, plane)))

    length(intersections) == 2 ? LineSegment(intersections...) : nothing
end

"""
    intersect(a::Vertex, b::Vertex, plane::Plane)

Intersect line segment and plane. Two intersections are returned when line
segment is contained in plane.
"""
function intersect(a::Vertex, b::Vertex, plane::Plane)
    intersections = []

    adistance = distance(a, plane)
    bdistance = distance(b, plane)

    if adistance * bdistance > eps(Float32)
        return intersections
    end

    ainplane = isapprox(adistance, 0, atol=eps(Float32))
    binplane = isapprox(bdistance, 0, atol=eps(Float32))

    if ainplane push!(intersections, a) end
    if binplane push!(intersections, b) end

    if ainplane || binplane
        return intersections
    end

    t = adistance / (adistance - bdistance)
    push!(intersections, a + t * (b - a))
end

function distance(point::Vertex, plane::Plane)
    dot(point - plane.point, plane.normal)
end


text:
using Oceananigans: AbstractModel
using Oceananigans.Grids
using Oceananigans.Utils: tupleit

struct KernelComputedField{X, Y, Z, S, A, G, K, C, F, P} <: AbstractField{X, Y, Z, A, G}
                     data :: A
                     grid :: G
                   kernel :: K
      boundary_conditions :: C
    computed_dependencies :: F
               parameters :: P
                   status :: S

    function KernelComputedField{X, Y, Z}(kernel::K, arch, grid;
                                          boundary_conditions = ComputedFieldBoundaryConditions(grid, (X, Y, Z)),
                                          computed_dependencies = (),
                                          parameters::P = nothing,
                                          data = nothing,
                                          recompute_safely = true) where {X, Y, Z, K, P}

        computed_dependencies = tupleit(computed_dependencies)

        if isnothing(data)
            data = new_data(arch, grid, (X, Y, Z))
        end

        # Use FieldStatus if we want to avoid always recomputing
        status = recompute_safely ? nothing : FieldStatus(0.0)

        G = typeof(grid)
        A = typeof(data)
        S = typeof(status)
        F = typeof(computed_dependencies)
        C = typeof(boundary_conditions)

        return new{X, Y, Z, S,
                   A, G, K, C, F, P}(data, grid, kernel, boundary_conditions, computed_dependencies, parameters)
    end
end

"""
    KernelComputedField(X, Y, Z, kernel, model;
                        boundary_conditions = ComputedFieldBoundaryConditions(grid, (X, Y, Z)),
                        computed_dependencies = (),
                        parameters = nothing,
                        data = nothing,
                        recompute_safely = true)

Builds a `KernelComputedField` at `X, Y, Z` computed with `kernel` and `model.architecture` and `model.grid`, with `boundary_conditions`.

`computed_dependencies` are an iterable of `AbstractField`s or other objects on which `compute!` is called prior to launching `kernel`.

`data` is a three-dimensional `OffsetArray` of scratch space where the kernel computation is stored.

If `data=nothing` (the default) then additional memory will be allocated to store the `data` of `KernelComputedField`.

If `isnothing(parameters)`, `kernel` is launched with the function signature

`kernel(data, grid, computed_dependencies...)`

Otherwise, `kernel` is launched with the function signature

`kernel(data, grid, computed_dependencies..., parameters)`

`recompute_safely` (default: `true`) determines whether the `KernelComputedField` is "recomputed" if embedded in the expression
tree of another operation.
    - If `recompute_safely=true`, the `KernelComputedField` is always recomputed.
    - If `recompute_safely=false`, the `KernelComputedField` will not be recomputed if its status is up-to-date.
    - If `data=nothing`, then `recompute_safely` is switched to `false`.

Example
=======

```julia
using KernelAbstractions: @index, @kernel
using Oceananigans.Fields: AveragedField, KernelComputedField, compute!
using Oceananigans.Grids: Center, Face

@inline ψ²(i, j, k, grid, ψ, Ψ) = @inbounds (ψ[i, j, k] - Ψ[i, j, k])^2
@kernel function compute_variance!(var, grid, ϕ, Φ)
    i, j, k = @index(Global, NTuple)

    @inbounds var[i, j, k] = ψ′²(i, j, k, grid, ϕ, Φ)
end

u, v, w = model.velocities

U = AveragedField(u, dims=(1, 2))
V = AveragedField(v, dims=(1, 2))

u′² = KernelComputedField(Face, Center, Center, compute_variance!, model; computed_dependencies=(u, U,))
v′² = KernelComputedField(Center, Face, Center, compute_variance!, model; computed_dependencies=(v, V,))
w′² = KernelComputedField(Center, Center, Face, compute_variance!, model; computed_dependencies=(w, 0,))

compute!(u′²)
compute!(v′²)
compute!(w′²)
```
"""
KernelComputedField(X, Y, Z, kernel, model::AbstractModel; kwargs...) =
    KernelComputedField{X, Y, Z}(kernel, model.architecture, model.grid; kwargs...)

KernelComputedField(X, Y, Z, kernel, arch::AbstractArchitecture, grid::AbstractGrid; kwargs...) =
    KernelComputedField{X, Y, Z}(kernel, arch, grid; kwargs...)

function compute!(kcf::KernelComputedField{X, Y, Z}) where {X, Y, Z}

    for dependency in kcf.computed_dependencies
        compute!(dependency)
    end

    arch = architecture(kcf.data)

    workgroup, worksize = work_layout(kcf.grid,
                                      :xyz,
                                      location=(X, Y, Z),
                                      include_right_boundaries=true)

    compute_kernel! = kcf.kernel(device(arch), workgroup, worksize)

    event = isnothing(kcf.parameters) ?
        compute_kernel!(kcf.data, kcf.grid, kcf.computed_dependencies...) :
        compute_kernel!(kcf.data, kcf.grid, kcf.computed_dependencies..., kcf.parameters)

    wait(device(arch), event)

    fill_halo_regions!(kcf, arch)

    return nothing
end

compute!(field::KernelComputedField{X, Y, Z, <:FieldStatus}, time) where {X, Y, Z} =
    conditional_compute!(field, time)

Adapt.adapt_structure(to, kcf::KernelComputedField) = Adapt.adapt(to, kcf.data)


text:
# ------------------------------------------------------------
# Truncated cosine and its Fourier transform

# Support = [-1/2,1/2]
function tcos(x)
	if 0.0 <= x <= 0.5
		return cos(2*pi*x)
	else
		return 0.0
	end
end

function ftcos(xi)
	if abs(xi) == 1.0
		return 0.25 + 0.0*im
	else
		return im*xi*( 1 + exp(-pi*xi*im) ) / (2*pi*(1 - xi^2))
	end
end

# Overload ftcos to vector input
@vectorize_1arg Float64 ftcos

# ------------------------------------------------------------
# Reconstruct truncated cosine in Haar basis

using GeneralizedSampling

J = 5
M = 2^(J+2)
# Both GeneralizedSampling and Plots (below) have a grid function
xi = GeneralizedSampling.grid(M, 0.5)
f = ftcos(xi)

T = Freq2Wave(xi, "haar", J)
wcoef = T \ f

# ------------------------------------------------------------
# Plot reconstruction

using IntervalWavelets

x, yw = weval( real(wcoef), "haar", 10 )

using Plots

p = plot(x, tcos,
         linecolor=:red,
         linewidth=2,
         label="true"
        )

plot!(p, x, yw,
      linecolor=:black,
      linewidth=2,
      label="approximation"
     )

display(p)


text:
function _degree_centrality(g::AbstractGraph, gtype::Integer; normalize=true)
    n_v = nv(g)
    c = zeros(n_v)
    for v in vertices(g)
        if gtype == 0    # count both in and out degree if appropriate
            deg = is_directed(g) ? outdegree(g, v) + indegree(g, v) : outdegree(g, v)
        elseif gtype == 1    # count only in degree
            deg = indegree(g, v)
        else                 # count only out degree
            deg = outdegree(g, v)
        end
        s = normalize ? (1.0 / (n_v - 1.0)) : 1.0
        c[v] = deg * s
    end
    return c
end

"""
    degree_centrality(g)
    indegree_centrality(g)
    outdegree_centrality(g)

Calculate the [degree centrality](https://en.wikipedia.org/wiki/Centrality#Degree_centrality)
of graph `g`. Return a vector representing the centrality calculated for each node in `g`.

### Optional Arguments
- `normalize=true`: If true, normalize each centrality measure by ``\\frac{1}{|V|-1}``.

# Examples
```jldoctest
julia> using LightGraphs

julia> degree_centrality(star_graph(4))
4-element Array{Float64,1}:
 1.0
 0.3333333333333333
 0.3333333333333333
 0.3333333333333333

julia> degree_centrality(path_graph(3))
3-element Array{Float64,1}:
 0.5
 1.0
 0.5
```
"""
degree_centrality(g::AbstractGraph; all...) = _degree_centrality(g, 0; all...)
indegree_centrality(g::AbstractGraph; all...) = _degree_centrality(g, 1; all...)
outdegree_centrality(g::AbstractGraph; all...) = _degree_centrality(g, 2; all...)


text:
using Test

@testset "./src/Scenario/QuantumBehavior.jl" begin

using BellScenario: QuantumBehavior

@testset "simple separable cases" begin

    test_cases = [
        ([1;0;0;0],[1;1;1;1;1;1;1;1;1]),
        ([0;1;0;0],[1;1;1;0;0;0;0;0;0]),
        ([0;0;1;0],[1;0;0;1;1;0;0;0;0]),
        ([0;0;0;1],[1;0;0;0;0;0;0;0;0])
    ]

    @testset "ψ = $(case[1])" for case in test_cases
        povm = [[1 0;0 0],[0 0;0 1]]

        behavior = QuantumBehavior.bipartite_scenario([povm,povm], [povm,povm], case[1])

        @test behavior == reshape(case[2],(9,1))
    end
end

@testset "prepare_and_measure()" begin
    @test [1 1 0]' == QuantumBehavior.prepare_and_measure(
        [[1 0;0 0],[0 0;0 1]],
        [[[1 0;0 0],[0 0;0 1]]]
    )
    @test [1 1 0.5 0 0.5]' == QuantumBehavior.prepare_and_measure(
        [[1 0;0 0],[0 0;0 1]],
        [[[1 0;0 0],[0 0;0 1]], [[0.5 0.5;0.5 0.5], [0.5 -0.5;-0.5 0.5]]]
    )
end

end


text:
import Multirate
import PyPlot

N𝜙              = 32
tapsPer𝜙        = 10
resampleRate    = 3.14
xƒ              = 0.1
Nx              = 25
x               = cos.(2*pi*xƒ*(0:Nx-1))
# x             = zeros( Nx); x[1] = 1.0
cutoffFreq      = min( 0.45, resampleRate )
transitionWidth = 0.05
hLen            = tapsPer𝜙*N𝜙
h               = Multirate.firdes( hLen, cutoffFreq, DSP.kaiser, samplerate = N𝜙 ) .* N𝜙
myfilter        = Multirate.FIRFilter( h, resampleRate, N𝜙 )
δfilter         = (hLen-1)/(2*N𝜙)
y               = Multirate.filt( myfilter, x )
tx              = (0:length(x)-1)
ty              = (0:length(y)-1)./resampleRate-δfilter

PyPlot.stem( tx, x, linefmt = "r-", markerfmt = "r." )
PyPlot.stem( ty, y, linefmt = "b-", markerfmt = "b." )
PyPlot.show()