text:
@inline function reim(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    return x.re, x.im
end

@inline function csgn(r,i)
    HI(r) > 0 && return  1
    HI(r) < 0 && return -1
    return HI(i) > 0 ? 1 : -1
end

square(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat} = x*x

sqrt(x::ComplexDF64) = ComplexDF64ComplexF128(sqrt, x)

# development from functions.wolfram.com

#=
Julia(evalc(Re(sqrt(r+i*sqrt(-1)))), optimize = true);
t1 = i ^ 2
t2 = r ^ 2
t4 = sqrt(t1 + t2)
t7 = sqrt(2) * sqrt(t4 + r) / 2
Julia(evalc(Im(sqrt(r+i*sqrt(-1)))), optimize = true);
Warning, type signature [complex] for function sqrt is not recognized
Warning, the function names {csgn} are not recognized in the target language
t3 = csgn(-1im * r + i)
t4 = i ^ 2
t5 = r ^ 2
t7 = sqrt(t4 + t5)
t11 = sqrt(2) * sqrt(t7 - r) * t3 / 2
=#
#=
function sqrt(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    fourthroot = sqrt(hypot(rea, ima))
    halfatan = atan(ima, rea) * 0.5
    rea = fourthroot * cos(halfatan)
    ima = fourthroot * sin(halfatan)
    return Complex{DoubleFloat{T}}(rea, ima)
end
=#

#=
Julia(evalc(Re(exp(r+i*sqrt(-1)))), optimize = true);
t1 = exp(r)
t2 = cos(i)
t3 = t2 * t1
Julia(evalc(Im(exp(r+i*sqrt(-1)))), optimize = true);
t1 = exp(r)
t2 = sin(i)
t3 = t2 * t1
=#
function exp(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    erea = exp(rea)
    realpart = erea * cos(ima)
    imagpart = erea * sin(ima)
    return Complex{DoubleFloat{T}}(realpart, imagpart)
end

#=

Julia(evalc(Re(log(r+i*sqrt(-1)))), optimize = true);
t1 = i ^ 2
t2 = r ^ 2
t4 = log(t1 + t2)
t5 = t4 / 2
Julia(evalc(Im(log(r+i*sqrt(-1)))), optimize = true);
t1 = atan2(i, r)
=#
function log(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    realpart = log(square(rea) + square(ima)) * 0.5
    imagpart = atan(ima, rea)
    return Complex{DoubleFloat{T}}(realpart, imagpart)
end

#=
Julia(evalc(Re(sin(r+i*sqrt(-1)))), optimize = true);
t1 = sin(r)
t2 = cosh(i)
t3 = t2 * t1
Julia(evalc(Im(sin(r+i*sqrt(-1)))), optimize = true);
t1 = cos(r)
t2 = sinh(i)
t3 = t2 * t1
=#
function sin(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    rea, ima = sin(rea) * cosh(ima), cos(rea) * sinh(ima)
    return Complex{DoubleFloat{T}}(rea, ima)
end

#=
Julia(evalc(Re(cos(r+i*sqrt(-1)))), optimize = true);
t1 = cos(r)
t2 = cosh(i)
t3 = t2 * t1
Julia(evalc(Im(cos(r+i*sqrt(-1)))), optimize = true);
t1 = sin(r)
t2 = sinh(i)
t4 = -t2 * t1
=#
function cos(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    rea, ima = cos(rea) * cosh(ima), -sin(rea) * sinh(ima)
    return Complex{DoubleFloat{T}}(rea, ima)
end

#=
Julia(evalc(Re(tan(r+i*sqrt(-1)))), optimize = true);
t1 = sin(r)
t2 = cos(r)
t4 = t2 ^ 2
t5 = sinh(i)                        # real(x), t5=0
t6 = t5 ^ 2                         # real(x), t6=0
t9 = 0.1e1 / (t4 + t6) * t2 * t1    # real(x), 1 / t4 * t2 * t1 = t2*t1/t4
Julia(evalc(Im(tan(r+i*sqrt(-1)))), optimize = true);
t1 = sinh(i)
t2 = cosh(i)
t4 = cos(r)
t5 = t4 ^ 2
t6 = t1 ^ 2
t9 = 0.1e1 / (t5 + t6) * t2 * t1
=#
function tan(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = 2*real(x), 2*imag(x)
    den = cos(rea) + cosh(ima)
    rea = sin(rea) / den
    ima = sinh(ima) / den
    return Complex{DoubleFloat{T}}(rea, ima)
end

#=
Julia(evalc(Re(csc(r+i*sqrt(-1)))), optimize = true);
t1 = sin(r)
t2 = cosh(i)
t4 = t1 ^ 2
t5 = sinh(i)
t6 = t5 ^ 2
t9 = 0.1e1 / (t4 + t6) * t2 * t1
Julia(evalc(Im(csc(r+i*sqrt(-1)))), optimize = true);
t1 = cos(r)
t2 = sinh(i)
t4 = sin(r)
t5 = t4 ^ 2
t6 = t2 ^ 2
t10 = -0.1e1 / (t5 + t6) * t2 * t1
=#
function csc(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    den = cos(2*rea) - cosh(2*ima)
    rea, ima = 2*(sin(rea)*cosh(ima)), 2*(cos(rea)*sinh(ima))
    rea = -rea / den
    ima = ima / den
    return Complex{DoubleFloat{T}}(rea, ima)
end

#=
Julia(evalc(Re(sec(r+i*sqrt(-1)))), optimize = true);
t1 = cos(r)
t2 = cosh(i)
t4 = t1 ^ 2
t5 = sinh(i)
t6 = t5 ^ 2
t9 = 0.1e1 / (t4 + t6) * t2 * t1
Julia(evalc(Im(sec(r+i*sqrt(-1)))), optimize = true);
t1 = sin(r)
t2 = sinh(i)
t4 = cos(r)
t5 = t4 ^ 2
t6 = t2 ^ 2
t9 = 0.1e1 / (t5 + t6) * t2 * t1
=#
function sec(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    den = cos(2*rea) + cosh(2*ima)
    rea, ima = 2*(cos(rea)*cosh(ima)), 2*(sin(rea)*sinh(ima))
    rea = rea / den
    ima = ima / den
    return Complex{DoubleFloat{T}}(rea, ima)
end

#=
Julia(evalc(Re(cot(r+i*sqrt(-1)))), optimize = true);
t1 = sin(r)
t2 = cos(r)
t4 = t1 ^ 2
t5 = sinh(i)
t6 = t5 ^ 2
t9 = 0.1e1 / (t4 + t6) * t2 * t1
Julia(evalc(Im(cot(r+i*sqrt(-1)))), optimize = true);
t1 = sinh(i)
t2 = cosh(i)
t4 = sin(r)
t5 = t4 ^ 2
t6 = t1 ^ 2
t10 = -0.1e1 / (t5 + t6) * t2 * t1
=#
function cot(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = 2*real(x), 2*imag(x)
    den = cos(rea) - cosh(ima)
    rea = -sin(rea) / den
    ima = sinh(ima) / den
    return Complex{DoubleFloat{T}}(rea, ima)
end

#=
Julia(evalc(Re(sinh(r+i*sqrt(-1)))), optimize = true);
t1 = sinh(r)
t2 = cos(i)
t3 = t2 * t1
Julia(evalc(Im(sinh(r+i*sqrt(-1)))), optimize = true);
t1 = cosh(r)
t2 = sin(i)
t3 = t2 * t1
=#
function sinh(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    rea, ima = sinh(rea) * cos(ima), cosh(rea) * sin(ima)
    return Complex{DoubleFloat{T}}(rea, ima)
end

#=
Julia(evalc(Re(cosh(r+i*sqrt(-1)))), optimize = true);
t1 = cosh(r)
t2 = cos(i)
t3 = t2 * t1
Julia(evalc(Im(cosh(r+i*sqrt(-1)))), optimize = true);
t1 = sinh(r)
t2 = sin(i)
t3 = t2 * t1
=#
function cosh(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    rea, ima = cosh(rea) * cos(ima), sinh(rea) * sin(ima)
    return Complex{DoubleFloat{T}}(rea, ima)
end

#=
Julia(evalc(Re(tanh(r+i*sqrt(-1)))), optimize = true);
t1 = sinh(r)
t2 = cosh(r)
t4 = t1 ^ 2
t5 = cos(i)
t6 = t5 ^ 2
t9 = 0.1e1 / (t4 + t6) * t2 * t1
Julia(evalc(Im(tanh(r+i*sqrt(-1)))), optimize = true);
t1 = sin(i)
t2 = cos(i)
t4 = sinh(r)
t5 = t4 ^ 2
t6 = t2 ^ 2
t9 = 0.1e1 / (t5 + t6) * t2 * t1
=#
function tanh(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = 2*real(x), 2*imag(x)
    den = cosh(rea) + cos(ima)
    rea = sinh(rea) / den
    ima = sin(ima) / den
    return Complex{DoubleFloat{T}}(rea, ima)
end

#=
Julia(evalc(Re(csch(r+i*sqrt(-1)))), optimize = true);
t1 = sinh(r)
t2 = cos(i)
t4 = t1 ^ 2
t5 = sin(i)
t6 = t5 ^ 2
t9 = 0.1e1 / (t4 + t6) * t2 * t1
Julia(evalc(Im(csch(r+i*sqrt(-1)))), optimize = true);
t1 = cosh(r)
t2 = sin(i)
t4 = sinh(r)
t5 = t4 ^ 2
t6 = t2 ^ 2
=#
function csch(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    den = cos(2*ima) - cosh(2*rea)
    rea, ima = 2*(sinh(rea)*cos(ima)), 2*(cosh(rea)*sin(ima))
    rea = -rea / den
    ima = ima / den
    return Complex{DoubleFloat{T}}(rea, ima)
end

#=
Julia(evalc(Re(sech(r+i*sqrt(-1)))), optimize = true);
t1 = cos(i)
t2 = cosh(r)
t4 = sinh(r)
t5 = t4 ^ 2
t6 = t1 ^ 2
t9 = 0.1e1 / (t5 + t6) * t2 * t1
Julia(evalc(Im(sech(r+i*sqrt(-1)))), optimize = true);
t1 = sin(i)
t2 = sinh(r)
t4 = t2 ^ 2
t5 = cos(i)
t6 = t5 ^ 2
t10 = -0.1e1 / (t4 + t6) * t2 * t1
=#
function sech(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    den = cos(2*ima) + cosh(2*rea)
    rea, ima = 2*(cos(ima)*cosh(rea)), 2*(sin(ima)*sinh(rea))
    rea = rea / den
    ima = -ima / den
    return Complex{DoubleFloat{T}}(rea, ima)
end

#=
Julia(evalc(Re(coth(r+i*sqrt(-1)))), optimize = true);
t1 = sinh(r)
t2 = cosh(r)
t4 = t1 ^ 2
t5 = sin(i)
t6 = t5 ^ 2
t9 = 0.1e1 / (t4 + t6) * t2 * t1
Julia(evalc(Im(coth(r+i*sqrt(-1)))), optimize = true);
t1 = sin(i)
t2 = cos(i)
t4 = sinh(r)
t5 = t4 ^ 2
t6 = t1 ^ 2
t10 = -0.1e1 / (t5 + t6) * t2 * t1
=#
function coth(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = 2*real(x), 2*imag(x)
    den = cos(ima) - cosh(rea)
    rea = -sinh(rea) / den
    ima = sin(ima) / den
    return Complex{DoubleFloat{T}}(rea, ima)
end

#=
Julia(evalc(Re(arcsin(r+i*sqrt(-1)))), optimize = true);
t2 = (r + 1) ^ 2
t3 = i ^ 2
t5 = sqrt(t2 + t3)
t7 = (r - 1) ^ 2
t9 = sqrt(t7 + t3)
t11 = asin(-t5 / 2 + t9 / 2)
t12 = -t11
Julia(evalc(Im(arcsin(r+i*sqrt(-1)))), optimize = true);
Warning, type signature [complex] for function sqrt is not recognized
Warning, the function names {csgn} are not recognized in the target language
t3 = csgn(-1im * r + i)
t5 = (r + 1) ^ 2
t6 = i ^ 2
t8 = sqrt(t5 + t6)
t11 = (r - 1) ^ 2
t13 = sqrt(t11 + t6)
t18 = sqrt((t8 + t13) ^ 2 / 4 - 1)
t20 = log(t8 / 2 + t13 / 2 + t18)
t21 = t20 * t3
=#
function asin(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    xp1sq = square(rea + 1)
    xm1sq = square(rea - 1)
    ysq   = square(ima)
    xp1py = sqrt(xp1sq + ysq)
    xm1py = sqrt(xm1sq + ysq)
    rea = asin((xp1py - xm1py) * 0.5)
    x = (xm1py + xp1py) * 0.5
    ima = sign(ima) * log(x + sqrt(square(x) - 1))
   return Complex{DoubleFloat{T}}(rea, ima)
end

#=
Julia(evalc(Re(arccos(r+i*sqrt(-1)))), optimize = true);
t2 = (r + 1) ^ 2
t3 = i ^ 2
t5 = sqrt(t2 + t3)
t7 = (r - 1) ^ 2
t9 = sqrt(t7 + t3)
t11 = acos(-t5 / 2 + t9 / 2)
t12 = pi - t11
Julia(evalc(Im(arccos(r+i*sqrt(-1)))), optimize = true);
Warning, type signature [complex] for function sqrt is not recognized
Warning, the function names {csgn} are not recognized in the target language
t3 = csgn(-1im * r + i)
t5 = (r + 1) ^ 2
t6 = i ^ 2
t8 = sqrt(t5 + t6)
t11 = (r - 1) ^ 2
t13 = sqrt(t11 + t6)
t18 = sqrt((t8 + t13) ^ 2 / 4 - 1)
t20 = log(t8 / 2 + t13 / 2 + t18)
t22 = -t20 * t3
=#
function acos(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    xp1sq = square(rea + 1)
    xm1sq = square(rea - 1)
    ysq   = square(ima)
    xp1py = sqrt(xp1sq + ysq)
    xm1py = sqrt(xm1sq + ysq)
    rea = acos((xp1py - xm1py) * 0.5)
    x = (xm1py + xp1py) * 0.5
    ima = -sign(ima) * log(x + sqrt(square(x) - 1))
   return Complex{DoubleFloat{T}}(rea, ima)
end

#=
Julia(evalc(Re(arctan(r+i*sqrt(-1)))), optimize = true);
t2 = atan2(r, 1 - i)
t4 = atan2(-r, i + 1)
t6 = t2 / 2 - t4 / 2
Julia(evalc(Im(arctan(r+i*sqrt(-1)))), optimize = true);
t1 = r ^ 2
t3 = (i + 1) ^ 2
t6 = (i - 1) ^ 2
t10 = log(1 / (t1 + t6) * (t1 + t3))
t11 = t10 / 4
=#
function atan(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    ima2 = square(ima)
    rea2 = square(rea)
    realpart = 0.5*(atan(2*rea/(1-rea2-ima2)) + (pi1o2(DoubleFloat{T})*(sign(rea2+ima2-1)+1)) * sign(rea))
    imagpart = (log(square(ima + 1) + rea2) - log(rea2 + square(1 - ima))) * 0.25
    return Complex{DoubleFloat{T}}(realpart, imagpart)
end

#=
Julia(evalc(Re(arccsc(r+i*sqrt(-1)))), optimize = true);
t1 = i ^ 2
t2 = r ^ 2
t3 = t1 + t2
t5 = 0.1e1 / t3 * r
t7 = (t5 + 1) ^ 2
t8 = t3 ^ 2
t10 = 0.1e1 / t8 * t1
t12 = sqrt(t7 + t10)
t14 = (t5 - 1) ^ 2
t16 = sqrt(t14 + t10)
t18 = asin(-t12 / 2 + t16 / 2)
t19 = -t18
Julia(evalc(Im(arccsc(r+i*sqrt(-1)))), optimize = true);
Warning, type signature [complex] for function sqrt is not recognized
Warning, the function names {csgn} are not recognized in the target language
t2 = csgn(1im * r + i)
t3 = i ^ 2
t4 = r ^ 2
t5 = t3 + t4
t7 = 0.1e1 / t5 * r
t9 = (t7 + 1) ^ 2
t10 = t5 ^ 2
t12 = 0.1e1 / t10 * t3
t14 = sqrt(t9 + t12)
t17 = (t7 - 1) ^ 2
t19 = sqrt(t17 + t12)
t24 = sqrt((t14 + t19) ^ 2 / 4 - 1)
t26 = log(t14 / 2 + t19 / 2 + t24)
t28 = -t26 * t2
=#
function acsc(z::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    return asin(inv(z))
end

#=
Julia(evalc(Re(arcsec(r+i*sqrt(-1)))), optimize = true);
t1 = i ^ 2
t2 = r ^ 2
t3 = t1 + t2
t5 = 0.1e1 / t3 * r
t7 = (t5 + 1) ^ 2
t8 = t3 ^ 2
t10 = 0.1e1 / t8 * t1
t12 = sqrt(t7 + t10)
t14 = (t5 - 1) ^ 2
t16 = sqrt(t14 + t10)
t18 = acos(-t12 / 2 + t16 / 2)
t19 = pi - t18
Julia(evalc(Im(arcsec(r+i*sqrt(-1)))), optimize = true);
Warning, type signature [complex] for function sqrt is not recognized
Warning, the function names {csgn} are not recognized in the target language
t2 = csgn(1im * r + i)
t3 = i ^ 2
t4 = r ^ 2
t5 = t3 + t4
t7 = 0.1e1 / t5 * r
t9 = (t7 + 1) ^ 2
t10 = t5 ^ 2
t12 = 0.1e1 / t10 * t3
t14 = sqrt(t9 + t12)
t17 = (t7 - 1) ^ 2
t19 = sqrt(t17 + t12)
t24 = sqrt((t14 + t19) ^ 2 / 4 - 1)
t26 = log(t14 / 2 + t19 / 2 + t24)
t27 = t26 * t2
=#
function asec(z::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    return acos(inv(z))
end

#=
Julia(evalc(Re(arccot(r+i*sqrt(-1)))), optimize = true);
t2 = atan2(r, 1 - i)
t4 = atan2(-r, i + 1)
t6 = pi / 2 - t2 / 2 + t4 / 2
Julia(evalc(Im(arccot(r+i*sqrt(-1)))), optimize = true);
t1 = r ^ 2
t3 = (i + 1) ^ 2
t6 = (i - 1) ^ 2
t10 = log(1 / (t1 + t6) * (t1 + t3))
t12 = -t10 / 4
=#
function acot(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    ima2 = square(ima)
    rea2 = square(rea)
    den = rea2 + ima2
    realpart = (atanxy((1 + (ima/den)), (rea/den)) - atanxy((1 - (ima/den)), -(rea/den))) * 0.5
    imagpart = (log((square(ima - 1) + rea2)/den) - log((square(ima +1) + rea2)/den)) * 0.25
    return Complex{DoubleFloat{T}}(realpart, imagpart)
end

#=
Julia(evalc(Re(arcsinh(r+i*sqrt(-1)))), optimize = true);
Warning, type signature [complex] for function sqrt is not recognized
Warning, the function names {csgn} are not recognized in the target language
t2 = csgn(r + 1im * i)
t3 = r ^ 2
t5 = (i + 1) ^ 2
t7 = sqrt(t3 + t5)
t10 = (i - 1) ^ 2
t12 = sqrt(t3 + t10)
t17 = sqrt((t7 + t12) ^ 2 / 4 - 1)
t19 = log(t7 / 2 + t12 / 2 + t17)
t20 = t19 * t2
Julia(evalc(Im(arcsinh(r+i*sqrt(-1)))), optimize = true);
t1 = r ^ 2
t3 = (i + 1) ^ 2
t5 = sqrt(t1 + t3)
t7 = (i - 1) ^ 2
t9 = sqrt(t1 + t7)
t11 = asin(-t5 / 2 + t9 / 2)
t12 = -t11
=#

function asinh(x::Complex{DoubleFloat{T}}) where {T}
    r, i = real(x), imag(x)
    signbit(r) && return - asinh(-x)
    r2 = r*r
    ip1sqr = (i + 1)^2
    sqrta  = sqrt(r2 + ip1sqr)
    im1sqr = (i - 1)^2
    sqrtb  = sqrt(r2 + im1sqr)
    t1 = sqrt((sqrta + sqrtb)^2 / 4 - 1)
    t2 = log(sqrta / 2 + sqrtb / 2 + t1)
    r = t2 * csgn(real(r + 1im * i), imag(r+1im * i))
    t3 = asin(-sqrta / 2 + sqrtb / 2)
    i  = - t3
    return Complex{DoubleFloat{T}}(r, i)
end

#=

Julia(evalc(Re(arccosh(r+i*sqrt(-1)))), optimize = true);
Warning, type signature [complex] for function sqrt is not recognized
Warning, the function names {csgn} are not recognized in the target language
t1 = -1im * r
t3 = csgn(t1 + 1im + i)
t5 = csgn(t1 + i)
t8 = (r + 1) ^ 2
t9 = i ^ 2
t11 = sqrt(t8 + t9)
t14 = (r - 1) ^ 2
t16 = sqrt(t14 + t9)
t21 = sqrt((t11 + t16) ^ 2 / 4 - 1)
t23 = log(t11 / 2 + t16 / 2 + t21)
t24 = t23 * t5 * t3
Julia(evalc(Im(arccosh(r+i*sqrt(-1)))), optimize = true);
Warning, type signature [complex] for function sqrt is not recognized
Warning, the function names {csgn} are not recognized in the target language
t3 = csgn(-1im * r + 1im + i)
t5 = (r + 1) ^ 2
t6 = i ^ 2
t8 = sqrt(t5 + t6)
t10 = (r - 1) ^ 2
t12 = sqrt(t10 + t6)
t14 = acos(-t8 / 2 + t12 / 2)
t16 = (pi - t14) * t3
=#
function acosh(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    t2 = (rea - 1) ^ 2
    t3 = ima ^ 2
    t5 = sqrt(t2 + t3)
    t8 = (rea + 1) ^ 2
    t10 = sqrt(t8 + t3)
    t15 = sqrt((t5 + t10) ^ 2 / 4 - 1)
    t17 = log(t5 / 2 + t10 / 2 + t15)
    t19 = acos(-t10 / 2 + t5 / 2)
    t19 = pi1o1(DoubleFloat{T}) - t19
    return Complex{DoubleFloat{T}}(t17, t19)
end

#=
Julia(evalc(Re(arctanh(r+i*sqrt(-1)))), optimize = true);
t2 = (r + 1) ^ 2
t3 = i ^ 2
t6 = (r - 1) ^ 2
t10 = log(1 / (t6 + t3) * (t2 + t3))
t11 = t10 / 4
Julia(evalc(Im(arctanh(r+i*sqrt(-1)))), optimize = true);
t2 = atan2(i, r + 1)
t4 = atan2(-i, 1 - r)
t6 = t2 / 2 - t4 / 2
=#

function atanh(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    ima2 = square(ima)
    rea2 = square(rea)
    realpart = (log(square(rea + 1) + ima2) - log(square(rea) - 2*rea + ima2 + 1)) * 0.25
    imagpart = 0.5*(atan(2*ima/(1-rea2-ima2)) + (pi1o2(DoubleFloat{T})*(sign(rea2+ima2-1)+1)) * sign(ima))
    return Complex{DoubleFloat{T}}(realpart, imagpart)
end

#=

Julia(evalc(Re(arccsch(r+i*sqrt(-1)))), optimize = true);
Warning, type signature [complex] for function sqrt is not recognized
Warning, the function names {csgn} are not recognized in the target language
t2 = csgn(1im * i - r)
t3 = r ^ 2
t4 = i ^ 2
t5 = t4 + t3
t6 = t5 ^ 2
t8 = 0.1e1 / t6 * t3
t10 = 0.1e1 / t5 * i
t12 = (-t10 + 1) ^ 2
t14 = sqrt(t8 + t12)
t17 = (-t10 - 1) ^ 2
t19 = sqrt(t8 + t17)
t24 = sqrt((t14 + t19) ^ 2 / 4 - 1)
t26 = log(t14 / 2 + t19 / 2 + t24)
t28 = -t26 * t2
Julia(evalc(Im(arccsch(r+i*sqrt(-1)))), optimize = true);
t1 = r ^ 2
t2 = i ^ 2
t3 = t2 + t1
t4 = t3 ^ 2
t6 = 0.1e1 / t4 * t1
t8 = 0.1e1 / t3 * i
t10 = (-t8 + 1) ^ 2
t12 = sqrt(t10 + t6)
t14 = (-t8 - 1) ^ 2
t16 = sqrt(t6 + t14)
t18 = asin(-t12 / 2 + t16 / 2)
t19 = -t18
=#
function acsch(z::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(z)
    t1 = square(ima)
    t2 = square(rea)
    t3 = t1 + t2
    t4 = inv(t3)
    t6 = square(t2)
    t8 = square(t1 - 1)
    t13 = square(t3)
    t14 = inv(t13)
    t16 = t14 * (t6 + t8 + 2 * (t1 + 1) * t2)
    t16 = sqrt(sqrt(t16))
    t23 = atanxy(1 + t14 * (-t1 + t2), -2 * t14 * ima * rea)
    t24 = t23 / 2
    t25 = cos(t24)
    t28 = square(t4 * rea + t25 * t16)
    t30 = sin(t24)
    t33 = square(t4 * ima - t30 * t16)
    t35 = log(t28 + t33)/2
    t37 = t16 * t3
    t44 = atanxy((t25 * t37 + rea) * t4, (t30 * t37 - ima) * t4)
    return Complex{DoubleFloat{T}}(t35, t44)
end

#=

Julia(evalc(Re(arcsech(r+i*sqrt(-1)))), optimize = true);
Warning, type signature [complex] for function sqrt is not recognized
Warning, the function names {csgn} are not recognized in the target language
t1 = i ^ 2
t2 = r ^ 2
t5 = csgn(1im * t1 + 1im * t2 + -1im * r - i)
t7 = csgn(1im * r + i)
t9 = t1 + t2
t11 = 0.1e1 / t9 * r
t13 = (t11 + 1) ^ 2
t14 = t9 ^ 2
t16 = 0.1e1 / t14 * t1
t18 = sqrt(t13 + t16)
t21 = (t11 - 1) ^ 2
t23 = sqrt(t21 + t16)
t28 = sqrt((t18 + t23) ^ 2 / 4 - 1)
t30 = log(t18 / 2 + t23 / 2 + t28)
t32 = -t30 * t7 * t5
Julia(evalc(Im(arcsech(r+i*sqrt(-1)))), optimize = true);
Warning, type signature [complex] for function sqrt is not recognized
Warning, the function names {csgn} are not recognized in the target language
t1 = i ^ 2
t2 = r ^ 2
t5 = csgn(1im * t1 + 1im * t2 + -1im * r - i)
t6 = t1 + t2
t8 = 0.1e1 / t6 * r
t10 = (t8 + 1) ^ 2
t11 = t6 ^ 2
t13 = 0.1e1 / t11 * t1
t15 = sqrt(t10 + t13)
t17 = (t8 - 1) ^ 2
t19 = sqrt(t17 + t13)
t21 = acos(-t15 / 2 + t19 / 2)
t23 = (pi - t21) * t5
=#
function asech(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    return log((sqrt(1-square(x))+1)/x)
end

#=
Julia(evalc(Re(arccoth(r+i*sqrt(-1)))), optimize = true);
t2 = (r + 1) ^ 2
t3 = i ^ 2
t6 = (r - 1) ^ 2
t10 = log(1 / (t6 + t3) * (t2 + t3))
t11 = t10 / 4
Julia(evalc(Im(arccoth(r+i*sqrt(-1)))), optimize = true);
t2 = atan2(i, r + 1)
t4 = atan2(i, r - 1)
t6 = t2 / 2 - t4 / 2
=#
function acoth(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    rea2 = square(rea)
    ima2 = square(ima)
    den = rea2 + ima2
    num = den + 1
    num1 = num + 2*rea
    num2 = num - 2*rea
    realpart = (log(num1/den) - log(num2/den)) * 0.25
    num1 = rea / den
    num2 = ima / den
    imagpart = (atanxy(num1 + 1, -num2) - atanxy(1 - num1, num2)) * 0.5
    return Complex{DoubleFloat{T}}(realpart, imagpart)
end


text:
##### Prognostic variable

export Mass, Momentum, Energy
export Moisture, TotalMoisture, LiquidMoisture, IceMoisture
export Precipitation, Rain, Snow
export Tracers

struct Mass <: PrognosticVariable end
struct Momentum <: PrognosticVariable end
struct Energy <: PrognosticVariable end

abstract type Moisture <: PrognosticVariable end
struct TotalMoisture <: Moisture end
struct LiquidMoisture <: Moisture end
struct IceMoisture <: Moisture end

abstract type Precipitation <: PrognosticVariable end
struct Rain <: Precipitation end
struct Snow <: Precipitation end

struct Tracers{N} <: PrognosticVariable end

ext: jl

max_forks_repo_name: jm-c/CLIMA


text:
module TimeOfDay_tests

export
  millisInOneSecond,
  microsInOnSecond,
  nanosInOneSecond,
  nanosInOneMilli,
  nanosInOneMicro,
  secondsInOneMinute,
  secondsInOneHour,
  secondsInOneDay,
  nanosInOneMinute,
  nanosInOneHour,
  nanosInOneDay,
  TimeOfDay,
  NegativeTimeOfDayError

const millisInOneSecond = Int64(1e3)
const microsInOnSecond = Int64(1e6)
const nanosInOneSecond = Int64(1e9)
const nanosInOneMilli = Int64(1e6)
const nanosInOneMicro = Int64(1e3)
const secondsInOneMinute = Int64(60)
const secondsInOneHour = Int64(3600)
const secondsInOneDay = secondsInOneHour * Int64(24)
const nanosInOneMinute = nanosInOneSecond * secondsInOneMinute
const nanosInOneHour = nanosInOneSecond * secondsInOneHour
const nanosInOneDay = nanosInOneSecond * secondsInOneDay

type NegativeTimeOfDayError{T} <: Exception
  val::T
end

immutable TimeOfDay <: Number
  nanosSinceMidnight::Int64

  function TimeOfDay{T<:Number}(val::T)
    rv = round(Int64, val)
    if (rv < 0)
      throw(NegativeTimeOfDayError(val))
    end
    new(rv)
  end
end

function TimeOfDay(str::AbstractString)
  tmp = split(str,":")
  minutes = 0
  seconds = 0
  nanos = 0

  hours = parse(Int64, tmp[1])
  if length(tmp) > 1
    minutes = parse(Int64, tmp[2])
  end

  if length(tmp) > 2
    tmp2 = split(tmp[3], ".")
    seconds = parse(Int64, tmp2[1])
    if length(tmp2) > 1
      nanos = parse(Int64, rpad(tmp2[2],9,"0")[1:9])
    end
  end

  TimeOfDay(hours * nanosInOneHour + minutes * nanosInOneMinute + seconds * nanosInOneSecond + nanos)
end

function Base.string(timeOfDay::TimeOfDay)
  secsSinceMidnight, nanos = divrem(timeOfDay.nanosSinceMidnight, nanosInOneSecond)
  hours, hourrem = divrem(secsSinceMidnight, secondsInOneHour)
  minutes, seconds = divrem(hourrem, secondsInOneMinute)
  microseconds = div(nanos, millisInOneSecond)
  string(lpad(hours,2,"0"), ":", lpad(minutes,2,"0"), ":", lpad(seconds,2,"0"), ".", lpad(microseconds,6,"0"))
end
Base.print(io::IO, time::TimeOfDay) = print(io, string(time))
Base.show(io::IO, time::TimeOfDay) = print(io, string(time))

# define how numbers get promoted to TimeOfDay
Base.promote_rule{T<:Real}(::Type{TimeOfDay}, ::Type{T}) = TimeOfDay
Base.convert(::Type{TimeOfDay}, str::AbstractString) = TimeOfDay(str)
Base.convert(::Type{TimeOfDay}, x::Real) = TimeOfDay(x)
Base.convert{T<:Number}(::Type{T}, x::TimeOfDay) = T(x.nanosSinceMidnight)

# now for each op, we define the op between TimeOfDay's, then the generic "promote" version of the op will handle conversion
Base.typemin(::Type{TimeOfDay}) = TimeOfDay(0)
Base.typemax(::Type{TimeOfDay}) = TimeOfDay(typemax(Int64))

for op in (:<, :>, :(==), :<=, :>=, :(Base.isless))
  @eval $op(t1::TimeOfDay, t2::TimeOfDay) = $op(t1.nanosSinceMidnight, t2.nanosSinceMidnight)
end
for op in (:+, :-, :(Base.rem))
  @eval $op(t1::TimeOfDay, t2::TimeOfDay) = TimeOfDay($op(t1.nanosSinceMidnight, t2.nanosSinceMidnight))
end

# some helpers for creating ranges of the form "timestart:timestep:timestop"
Base.colon(start::TimeOfDay, step::TimeOfDay, stop::TimeOfDay) = map(TimeOfDay, start.nanosSinceMidnight : step.nanosSinceMidnight : stop.nanosSinceMidnight)
Base.colon(start::TimeOfDay, step::Number, stop::TimeOfDay) = colon(start, TimeOfDay(step), stop)

iszero(t::TimeOfDay) = t.nanosSinceMidnight == 0

using FactCheck

facts("TimeOfDay") do
  @fact TimeOfDay("10") => 10 * nanosInOneHour
  @fact TimeOfDay("10:30") => 10 * nanosInOneHour + 30 * nanosInOneMinute
  @fact TimeOfDay("10:30:05") => 10 * nanosInOneHour + 30 * nanosInOneMinute + 5 * nanosInOneSecond
  @fact TimeOfDay("10:30:05.00001") => 10 * nanosInOneHour + 30 * nanosInOneMinute + 5 * nanosInOneSecond + 10 * nanosInOneMicro

  t = TimeOfDay("10:30:05")
  buf = IOBuffer()
  print(buf, t)
  @fact bytestring(buf) => "10:30:05.000000"

  tdiff = TimeOfDay(nanosInOneMinute)
  trange = t : tdiff : t + 5 * nanosInOneMinute
  @fact length(trange) => 6
  @fact trange[2] => TimeOfDay("10:31:05")

  @fact iszero(TimeOfDay(0)) => true
  @fact iszero(TimeOfDay(1)) => false
  @fact_throws TimeOfDay(-1)
end

end # module
tod = TimeOfDay_tests


text:
module GenieAuthorisation

using Reexport
using Genie
using SearchLight, SearchLight.Relationships
@reexport using GenieAuthentication

export Role, Permission
export assign_role, assign_permission
export has_role, has_permission, @authorised!

struct UnexpectedTypeException <: Exception
  expected_type::String
  received_type::String
end

struct UndefinedPermissionException <: Exception
  permission::String
end

Base.@kwdef mutable struct Role <: AbstractModel
  id::DbId = DbId()
  name::String = ""
end
Role(name::Union{String,Symbol}) = Role(name = string(name))

Base.@kwdef mutable struct Permission <: AbstractModel
  id::DbId = DbId()
  name::String = ""
end
Permission(name::Union{String,Symbol}) = Permission(name = string(name))

function assign_role(user::U, role::Role)::Bool where {U<:AbstractModel}
  Relationship!(user, role)

  true
end

function assign_permission(role::Role, permission::Permission) :: Bool
  Relationship!(role, permission)

  true
end

function has_role(user::U, role::Role)::Bool where {U<:AbstractModel}
  isrelated(user, role)
end

function has_role(user::U, role::Union{String,Symbol})::Bool where {U<:AbstractModel}
  has_role(user, findone(Role, name = string(role)))
end

function fetch_permission(permission::Union{Permission,String,Symbol}) :: Union{Permission,Nothing}
  isa(permission, Permission) || (permission = findone(Permission, name = string(permission)))

  permission
end

function has_permission(role::Role, permission::Union{Permission,String,Symbol})::Bool
  permission = fetch_permission(permission)
  permission === nothing && return false
  isrelated(role, permission)
end

function has_permission(user::U, permission::Union{Permission,String,Symbol})::Bool where {U<:AbstractModel}
  permission = fetch_permission(permission)
  permission === nothing && return false
  isrelated(user, permission, through = [Role])
end

function has_permission(u::Nothing, permission)::Bool
  false
end

macro authorised!(permission, exception = Genie.Exceptions.NotFoundException("Page"))
  :(has_permission($(esc( :( Main.UserApp.current_user() ) )), $(esc(permission))) || throw($exception))
end

"""
    install(dest::String; force = false, debug = false) :: Nothing

Copies the plugin's files into the host Genie application.
"""
function install(dest::String; force = false, debug = false) :: Nothing
  # automatically install the GenieAuthentication plugin -- however, do not force the install
  GenieAuthentication.install(dest; force = false, debug = debug)

  src = abspath(normpath(joinpath(pathof(@__MODULE__) |> dirname, "..", Genie.Plugins.FILES_FOLDER)))

  debug && @info "Preparing to install from $src into $dest"
  debug && @info "Found these to install $(readdir(src))"

  for f in readdir(src)
    debug && @info "Processing $(joinpath(src, f))"
    debug && @info "$(isdir(joinpath(src, f)))"

    isdir(joinpath(src, f)) || continue

    debug && "Installing from $(joinpath(src, f))"

    Genie.Plugins.install(joinpath(src, f), dest, force = force)
  end

  nothing
end

end


text:
# """
#     cse

# Common subexpression elimination.

# """
# function cse(x::Sym, subxs = Dict{Sym,Sym}())
#     if hashead(x, :call)
#         args = tailargs(x)
#     else
#         args = getargs(x)
#     end

#     for i in eachindex(args)
#         if args[i] in keys(subxs)
#             args[i] = subxs[args[i]]
#         else
#             args[i] = subxs[args[i]] = Sym{tagof(args[i])}(gensym())
#         end
#     end
#     return args, subxs
# end

"""
    normalorder

Sort by commuting adjacent elements.

"""
normalorder(x::Sym, op::Sym; kwargs...) = normalorder!(deepcopy(x), op; kwargs...)

normalorder!(x, op; kwargs...) = x
function normalorder!(x::Sym, op::Sym; context::Context = @__context__(), kwargs...)
    isatomic(x) && return x
    map(a -> normalorder!(a, op; context=context), getargs(x))
    (hashead(x, :call) && (firstarg(x) == op) === true) || return x
    sort!(getargs(x), alg=InsertionSort, lt=(a, b) -> _commute_lt(op, a, b; context=context))
    return query!(x, context=context)
end

function _commute_lt(op, a, b; context::Context)
    return @scope context begin
        (commuteswith(op, a, b) === true ||
         commuteswith(op, b, a) === true) &&
         isless(string(a), string(b))
    end
end

"""
    rle

Run length encode.

"""
rle(x::Sym, op::Sym, enc_op::Sym; kwargs...) = rle!(deepcopy(x), op, enc_op; kwargs...)

rle!(x, op, enc_op; kwargs...) = x
function rle!(x::Sym, op::Sym, enc_op::Sym; context::Context = @__context__(), kwargs...)
    isatomic(x) && return x
    map(a -> rle!(a, op, enc_op, context=context), getargs(x))
    (hashead(x, :call) && (firstarg(x) == op) === true) || return x

    N = length(getargs(x))
    args = []
    n = 1
    for i in 1:N
        if i < N && (getargs(x, i) == getargs(x, i + 1)) === true
            n += 1
        else
            if n > 1
                push!(args, @scope context firstarg(enc_op)(getargs(x, i), n))
            else
                push!(args, getargs(x, i))
            end
            n = 1
        end
    end
    return query!(setargs!(x, args), context=context)
end

"""
    remove_identities(x, op, isidentity)

Remove identities of the operator `op`.  `isidentity` is a function that
returns true, if the element is an identity element of `op`.

"""
remove_identities(x::Sym, op::Sym, isidentity::Sym; kwargs...) = remove_identities!(deepcopy(x), op, isidentity; kwargs...)

remove_identities!(x, op, isidentity; kwargs...) = x
function remove_identities!(x::Sym, op::Sym, isidentity::Sym; context::Context = @__context__)
    isatomic(x) && return x
    map(a -> remove_identities!(a, op, isidentity, context=context), getargs(x))
    (hashead(x, :call) && (firstarg(x) == op) === true) || return x
    args = getargs(x)
    deleteat!(
        args,
        findall(
            arg -> @scope(context, firstarg(isidentity)(arg)) === true,
            args
        )
    )
    return query!(x, context=context)
end

"""
    resolve_absorbing(x, op, isabsorbing)

Resolves absorbing elements of the operator `op`.  `isabsorbing` is a
function that returns true, if the element is an absorbing element of
`op`

"""
resolve_absorbing(x::Sym, op::Sym, isabsorbing::Sym; kwargs...) = resolve_absorbing!(deepcopy(x), op, isabsorbing; kwargs...)

resolve_absorbing!(x, op, isabsorbing; kwargs...) = x
function resolve_absorbing!(x::Sym, op::Sym, isabsorbing::Sym; context::Context = @__context__)
    isatomic(x) && return x
    map(a -> resolve_absorbing!(a, op, isabsorbing, context=context), getargs(x))
    (hashead(x, :call) && (firstarg(x) == op) === true) || return x
    args = getargs(x)
    absorber_index = findfirst(arg -> @scope(context, firstarg(isabsorbing)(arg)) === true, args)
    absorber_index === nothing && return x
    setargs!(x, Any[firstarg(x), getargs(x, absorber_index)])
    return query!(x, context=context)
end

split_reapply(x::Sym, op::Sym; kwargs...) = split_reapply!(deepcopy(x), op; kwargs...)

split_reapply!(x, op; kwargs...) = x
function split_reapply!(x::Sym{TAG}, op::Sym; context::Context = @__context__) where TAG
    isatomic(x) && return x
    map(a -> split_reapply!(a, op, context=context), getargs(x))
    (hashead(x, :call) && (firstarg(x) == op) === true) || return x
    args = []
    for arg in getargs(x)
        if hashead(arg, :call) && (firstarg(arg) == op) === true
            append!(args, tailargs(arg))
        else
            push!(args, arg)
        end
    end
    return setargs!(x, args)
end

##################################################
# Multi-algorithms
##################################################

const MAX_GATHER_ITERATIONS = 10

gather(x; kwargs...) = gather!(deepcopy(x); kwargs...)
function gather!(x; context::Context = @__context__)
    x_ = deepcopy(x)
    for _ in 1:MAX_GATHER_ITERATIONS
        _gather_one_iteration!(x_; context=context)
        if (x == x_) === true
            break
        else
            x, x_ = x_, deepcopy(x_)
        end
    end
    return query!(x_, context=context)
end

_gather_one_iteration!(x; kwargs...) = x
function _gather_one_iteration!(x::Sym; context::Context = @__context__)
    remove_identities!(x, Sym(+), Sym(iszero), context=context)
    remove_identities!(x, Sym(*), Sym(isone), context=context)
    resolve_absorbing!(x, Sym(*), Sym(iszero), context=context)
    split_reapply!(x, Sym(+), context=context)
    split_reapply!(x, Sym(*), context=context)
    normalorder!(x, Sym(+), context=context)
    normalorder!(x, Sym(*), context=context)
    rle!(x, Sym(+), Sym(*), context=context)
    rle!(x, Sym(*), Sym(^), context=context)

    remove_identities!(x, Sym(|), Sym(isfalse), context=context)
    remove_identities!(x, Sym(&), Sym(istrue), context=context)
    resolve_absorbing!(x, Sym(|), Sym(istrue), context=context)
    resolve_absorbing!(x, Sym(&), Sym(isfalse), context=context)
    normalorder!(x, Sym(&), context=context)
    normalorder!(x, Sym(|), context=context)

    return x
end


text:
int_rules_6_7_6 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 6.7.6*(c+d*x)^m*hyper(a+b*x)^n*hyper(a+b*x)^p =#
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sinh(~(a') + ~(b') * ~x) ^ ~(n') * cosh(~(a') + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ ~m * sinh(~a + ~b * ~x) ^ (~n + 1)) / (~b * (~n + 1)) - ((~d * ~m) / (~b * (~n + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * sinh(~a + ~b * ~x) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sinh(~(a') + ~(b') * ~x) * cosh(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => ((~c + ~d * ~x) ^ ~m * cosh(~a + ~b * ~x) ^ (~n + 1)) / (~b * (~n + 1)) - ((~d * ~m) / (~b * (~n + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * cosh(~a + ~b * ~x) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sinh(~(a') + ~(b') * ~x) ^ ~(n') * cosh(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => Antiderivative(ExpandTrigReduce((~c + ~d * ~x) ^ ~m, sinh(~a + ~b * ~x) ^ ~n * cosh(~a + ~b * ~x) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (IGtQ(~n, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sinh(~(a') + ~(b') * ~x) ^ ~(n') * tanh(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => Antiderivative((~c + ~d * ~x) ^ ~m * sinh(~a + ~b * ~x) ^ ~n * tanh(~a + ~b * ~x) ^ (~p - 2), ~x) - Antiderivative((~c + ~d * ~x) ^ ~m * sinh(~a + ~b * ~x) ^ (~n - 2) * tanh(~a + ~b * ~x) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (IGtQ(~n, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * cosh(~(a') + ~(b') * ~x) ^ ~(n') * coth(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => Antiderivative((~c + ~d * ~x) ^ ~m * cosh(~a + ~b * ~x) ^ ~n * coth(~a + ~b * ~x) ^ (~p - 2), ~x) + Antiderivative((~c + ~d * ~x) ^ ~m * cosh(~a + ~b * ~x) ^ (~n - 2) * coth(~a + ~b * ~x) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (IGtQ(~n, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sech(~(a') + ~(b') * ~x) ^ ~(n') * tanh(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => (-((~c + ~d * ~x) ^ ~m) * sech(~a + ~b * ~x) ^ ~n) / (~b * ~n) + ((~d * ~m) / (~b * ~n)) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * sech(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (EqQ(~p, 1) && GtQ(~m, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * csch(~(a') + ~(b') * ~x) ^ ~(n') * coth(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => (-((~c + ~d * ~x) ^ ~m) * csch(~a + ~b * ~x) ^ ~n) / (~b * ~n) + ((~d * ~m) / (~b * ~n)) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * csch(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (EqQ(~p, 1) && GtQ(~m, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sech(~(a') + ~(b') * ~x) ^ 2 * tanh(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => ((~c + ~d * ~x) ^ ~m * tanh(~a + ~b * ~x) ^ (~n + 1)) / (~b * (~n + 1)) - ((~d * ~m) / (~b * (~n + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * tanh(~a + ~b * ~x) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * csch(~(a') + ~(b') * ~x) ^ 2 * coth(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => (-((~c + ~d * ~x) ^ ~m) * coth(~a + ~b * ~x) ^ (~n + 1)) / (~b * (~n + 1)) + ((~d * ~m) / (~b * (~n + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * coth(~a + ~b * ~x) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sech(~(a') + ~(b') * ~x) * tanh(~(a') + ~(b') * ~x) ^ ~p, ~x) => Antiderivative((~c + ~d * ~x) ^ ~m * sech(~a + ~b * ~x) * tanh(~a + ~b * ~x) ^ (~p - 2), ~x) - Antiderivative((~c + ~d * ~x) ^ ~m * sech(~a + ~b * ~x) ^ 3 * tanh(~a + ~b * ~x) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && IGtQ(~p / 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sech(~(a') + ~(b') * ~x) ^ ~(n') * tanh(~(a') + ~(b') * ~x) ^ ~p, ~x) => Antiderivative((~c + ~d * ~x) ^ ~m * sech(~a + ~b * ~x) ^ ~n * tanh(~a + ~b * ~x) ^ (~p - 2), ~x) - Antiderivative((~c + ~d * ~x) ^ ~m * sech(~a + ~b * ~x) ^ (~n + 2) * tanh(~a + ~b * ~x) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && IGtQ(~p / 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * csch(~(a') + ~(b') * ~x) * coth(~(a') + ~(b') * ~x) ^ ~p, ~x) => Antiderivative((~c + ~d * ~x) ^ ~m * csch(~a + ~b * ~x) * coth(~a + ~b * ~x) ^ (~p - 2), ~x) + Antiderivative((~c + ~d * ~x) ^ ~m * csch(~a + ~b * ~x) ^ 3 * coth(~a + ~b * ~x) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && IGtQ(~p / 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * csch(~(a') + ~(b') * ~x) ^ ~(n') * coth(~(a') + ~(b') * ~x) ^ ~p, ~x) => Antiderivative((~c + ~d * ~x) ^ ~m * csch(~a + ~b * ~x) ^ ~n * coth(~a + ~b * ~x) ^ (~p - 2), ~x) + Antiderivative((~c + ~d * ~x) ^ ~m * csch(~a + ~b * ~x) ^ (~n + 2) * coth(~a + ~b * ~x) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && IGtQ(~p / 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sech(~(a') + ~(b') * ~x) ^ ~(n') * tanh(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => With([u = IntHide(sech(~a + ~b * ~x) ^ ~n * tanh(~a + ~b * ~x) ^ ~p, ~x)], Dist((~c + ~d * ~x) ^ ~m, u, ~x) - ~d * ~m * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * u, ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (IGtQ(~m, 0) && (IntegerQ(~n / 2) || IntegerQ((~p - 1) / 2)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * csch(~(a') + ~(b') * ~x) ^ ~(n') * coth(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => With([u = IntHide(csch(~a + ~b * ~x) ^ ~n * coth(~a + ~b * ~x) ^ ~p, ~x)], Dist((~c + ~d * ~x) ^ ~m, u, ~x) - ~d * ~m * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * u, ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (IGtQ(~m, 0) && (IntegerQ(~n / 2) || IntegerQ((~p - 1) / 2)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * csch(~(a') + ~(b') * ~x) ^ ~(n') * sech(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => 2 ^ ~n * Antiderivative((~c + ~d * ~x) ^ ~m * csch(2 * ~a + 2 * ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (RationalQ(~m) && IntegerQ(~n))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * csch(~(a') + ~(b') * ~x) ^ ~(n') * sech(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => With([u = IntHide(csch(~a + ~b * ~x) ^ ~n * sech(~a + ~b * ~x) ^ ~p, ~x)], Dist((~c + ~d * ~x) ^ ~m, u, ~x) - ~d * ~m * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * u, ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (IntegersQ(~n, ~p) && (GtQ(~m, 0) && NeQ(~n, ~p)))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~F)(~v) ^ ~(n') * (~G)(~w) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * (~F)(ExpandToSum(~v, ~x)) ^ ~n * (~G)(ExpandToSum(~v, ~x)) ^ ~p, ~x) <-- FreeQ([~m, ~n, ~p], ~x) && (HyperbolicQ(~F) && (HyperbolicQ(~G) && (EqQ(~v, ~w) && (LinearQ([~u, ~v, ~w], ~x) && Not(LinearMatchQ([~u, ~v, ~w], ~x))))))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * cosh(~(c') + ~(d') * ~x) * (~a + ~(b') * sinh(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => ((~e + ~f * ~x) ^ ~m * (~a + ~b * sinh(~c + ~d * ~x)) ^ (~n + 1)) / (~b * ~d * (~n + 1)) - ((~f * ~m) / (~b * ~d * (~n + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * (~a + ~b * sinh(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * sinh(~(c') + ~(d') * ~x) * (~a + ~(b') * cosh(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => ((~e + ~f * ~x) ^ ~m * (~a + ~b * cosh(~c + ~d * ~x)) ^ (~n + 1)) / (~b * ~d * (~n + 1)) - ((~f * ~m) / (~b * ~d * (~n + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * (~a + ~b * cosh(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * sech(~(c') + ~(d') * ~x) ^ 2 * (~a + ~(b') * tanh(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => ((~e + ~f * ~x) ^ ~m * (~a + ~b * tanh(~c + ~d * ~x)) ^ (~n + 1)) / (~b * ~d * (~n + 1)) - ((~f * ~m) / (~b * ~d * (~n + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * (~a + ~b * tanh(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * csch(~(c') + ~(d') * ~x) ^ 2 * (~a + ~(b') * coth(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => (-((~e + ~f * ~x) ^ ~m) * (~a + ~b * coth(~c + ~d * ~x)) ^ (~n + 1)) / (~b * ~d * (~n + 1)) + ((~f * ~m) / (~b * ~d * (~n + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * (~a + ~b * coth(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * sech(~(c') + ~(d') * ~x) * tanh(~(c') + ~(d') * ~x) * (~a + ~(b') * sech(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => (-((~e + ~f * ~x) ^ ~m) * (~a + ~b * sech(~c + ~d * ~x)) ^ (~n + 1)) / (~b * ~d * (~n + 1)) + ((~f * ~m) / (~b * ~d * (~n + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * (~a + ~b * sech(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * csch(~(c') + ~(d') * ~x) * coth(~(c') + ~(d') * ~x) * (~a + ~(b') * csch(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => (-((~e + ~f * ~x) ^ ~m) * (~a + ~b * csch(~c + ~d * ~x)) ^ (~n + 1)) / (~b * ~d * (~n + 1)) + ((~f * ~m) / (~b * ~d * (~n + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * (~a + ~b * csch(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * sinh(~(a') + ~(b') * ~x) ^ ~(p') * sinh(~(c') + ~(d') * ~x) ^ ~(q'), ~x) => Antiderivative(ExpandTrigReduce((~e + ~f * ~x) ^ ~m, sinh(~a + ~b * ~x) ^ ~p * sinh(~c + ~d * ~x) ^ ~q, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 0) && (IGtQ(~q, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * cosh(~(a') + ~(b') * ~x) ^ ~(p') * cosh(~(c') + ~(d') * ~x) ^ ~(q'), ~x) => Antiderivative(ExpandTrigReduce((~e + ~f * ~x) ^ ~m, cosh(~a + ~b * ~x) ^ ~p * cosh(~c + ~d * ~x) ^ ~q, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 0) && (IGtQ(~q, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * sinh(~(a') + ~(b') * ~x) ^ ~(p') * cosh(~(c') + ~(d') * ~x) ^ ~(q'), ~x) => Antiderivative(ExpandTrigReduce((~e + ~f * ~x) ^ ~m, sinh(~a + ~b * ~x) ^ ~p * cosh(~c + ~d * ~x) ^ ~q, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (IGtQ(~p, 0) && IGtQ(~q, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~F)(~(a') + ~(b') * ~x) ^ ~(p') * (~G)(~(c') + ~(d') * ~x) ^ ~(q'), ~x) => Antiderivative(ExpandTrigExpand((~e + ~f * ~x) ^ ~m * (~G)(~c + ~d * ~x) ^ ~q, ~F, ~c + ~d * ~x, ~p, ~b / ~d, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (MemberQ([sinh, cosh], ~F) && (MemberQ([sech, csch], ~G) && (IGtQ(~p, 0) && (IGtQ(~q, 0) && (EqQ(~b * ~c - ~a * ~d, 0) && IGtQ(~b / ~d, 1))))))
end

max_stars_repo_path: src/rules/6 Hyperbolic functions/6.7 Miscellaneous/6.7.6 (c+d x)^m hyper(a+b x)^n hyper(a+b x)^p.jl


text:
"""
    SynchronousExecution

Functor object for using synchronous execution in an L-shaped algorithm (assuming multiple Julia cores are available). Create by supplying a [`Synchronous`](@ref) object through `execution` in the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

"""
struct SynchronousExecution{H <: AbstractFeasibilityHandler,
                            T <: AbstractFloat,
                            A <: AbstractVector} <: AbstractLShapedExecution
    subworkers::Vector{SubWorker{H,T}}
    decisions::Vector{DecisionChannel}
    subobjectives::A
    model_objectives::A
    metadata::Vector{MetaData}
    cutqueue::CutQueue{T}

    function SynchronousExecution(structure::VerticalStructure{2, 1, <:Tuple{DistributedScenarioProblems}},
                                  ::Type{F}, ::Type{T}, ::Type{A}) where {F <: AbstractFeasibility,
                                                                          T <: AbstractFloat,
                                                                          A <: AbstractVector}
        H = HandlerType(F)
        return new{H,T,A}(Vector{SubWorker{H,T}}(undef, nworkers()),
                          scenarioproblems(structure).decisions,
                          A(),
                          A(),
                          Vector{MetaData}(undef, nworkers()),
                          RemoteChannel(() -> Channel{QCut{T}}(4 * nworkers() * num_scenarios(structure))))
    end
end

function initialize_subproblems!(execution::SynchronousExecution,
                                 scenarioproblems::DistributedScenarioProblems)
    load_subproblems!(execution.subworkers, scenarioproblems, execution.decisions)
    return nothing
end

function finish_initilization!(lshaped::AbstractLShaped, execution::SynchronousExecution)
    append!(execution.subobjectives, fill(1e10, num_thetas(lshaped)))
    append!(execution.model_objectives, fill(-1e10, num_thetas(lshaped)))
    for w in workers()
        execution.metadata[w-1] = RemoteChannel(() -> MetaChannel(), w)
    end
    return lshaped
end

function restore_subproblems!(::AbstractLShaped, execution::SynchronousExecution)
    restore_subproblems!(execution.subworkers)
    return nothing
end

function resolve_subproblems!(lshaped::AbstractLShaped, execution::SynchronousExecution{H,T}) where {H <: AbstractFeasibilityHandler, T <: AbstractFloat}
    # Update metadata
    for w in workers()
        put!(execution.metadata[w-1], timestamp(lshaped), :gap, gap(lshaped))
    end
    @sync begin
        for (i,w) in enumerate(workers())
            worker_aggregator = remote_aggregator(lshaped.aggregation, scenarioproblems(lshaped.structure), w)
            @async remotecall_fetch(resolve_subproblems!,
                                    w,
                                    execution.subworkers[w-1],
                                    execution.decisions[w-1],
                                    lshaped.x,
                                    execution.cutqueue,
                                    worker_aggregator,
                                    timestamp(lshaped),
                                    execution.metadata[w-1])
        end
    end
    # Assume no cuts are added
    added = false
    # Collect incoming cuts
    while isready(execution.cutqueue)
        _, cut::SparseHyperPlane{T} = take!(execution.cutqueue)
        added |= add_cut!(lshaped, cut)
    end
    # Return current objective value and cut_added flag
    return current_objective_value(lshaped), added
end

# API
# ------------------------------------------------------------
function (execution::Synchronous)(structure::VerticalStructure{2, 1, <:Tuple{DistributedScenarioProblems}},
                                  ::Type{F}, ::Type{T}, ::Type{A}) where {F <: AbstractFeasibility,
                                                                          T <: AbstractFloat,
                                                                          A <: AbstractVector}
    return SynchronousExecution(structure, F, T, A)
end

function str(::Synchronous)
    return "Synchronous "
end


text:
module Runoff

using LinearAlgebra
using DocStringExtensions

using ...VariableTemplates
using ...Land:
    SoilModel,
    pressure_head,
    hydraulic_conductivity,
    get_temperature,
    effective_saturation,
    impedance_factor,
    viscosity_factor,
    moisture_factor

export AbstractPrecipModel,
    DrivenConstantPrecip,
    AbstractSurfaceRunoffModel,
    NoRunoff,
    compute_surface_grad_bc,
    CoarseGridRunoff

"""
    AbstractPrecipModel{FT <: AbstractFloat}
"""
abstract type AbstractPrecipModel{FT <: AbstractFloat} end

"""
    DrivenConstantPrecip{FT, F} <: AbstractPrecipModel{FT}

Instance of a precipitation distribution where the precipication value
is constant across the domain. However, this value can change in time.

Precipitation is assumed to be aligned with the vertical, i.e. P⃗ = Pẑ,
with P<0.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct DrivenConstantPrecip{FT, F} <: AbstractPrecipModel{FT}
    "Mean precipitation in grid"
    mp::F
    function DrivenConstantPrecip{FT}(mp::F) where {FT, F}
        new{FT, F}(mp)
    end
end

function (dcp::DrivenConstantPrecip{FT})(t::Real) where {FT}
    return FT(dcp.mp(t))
end

"""
    AbstractSurfaceRunoffModel

Abstract type for different surface runoff models. Currently, only
`NoRunoff` is supported.
"""
abstract type AbstractSurfaceRunoffModel end

"""
    NoRunoff <: AbstractSurfaceRunoffModel

Chosen when no runoff is to be modeled.
"""
struct NoRunoff <: AbstractSurfaceRunoffModel end

"""
    CoarseGridRunoff{FT} <: AbstractSurfaceRunoffModel

Chosen when no subgrid effects are to be modeled.
"""
struct CoarseGridRunoff{FT} <: AbstractSurfaceRunoffModel
    "Mean vertical resolution at the surface"
    Δz::FT
end

"""
    function compute_surface_grad_bc(soil::SoilModel,
                                     runoff_model::CoarseGridRunoff,
                                     precip_model::AbstractPrecipModel,
                                     n̂,
                                     state⁻::Vars,
                                     diff⁻::Vars,
                                     aux⁻::Vars,
                                     t::Real
                                     )

Given a runoff model and a precipitation distribution function, compute
the surface water flux normal to the surface. The sign of the flux
(inwards or outwards) is determined by the magnitude of precipitation,
evaporation, and infiltration.
"""
function compute_surface_grad_bc(
    soil::SoilModel,
    runoff_model::CoarseGridRunoff,
    precip_model::AbstractPrecipModel,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    t::Real,
)
    FT = eltype(state⁻)
    precip_vector = (FT(0), FT(0), precip_model(t))
    incident_water_flux = dot(precip_vector, n̂)
    Δz = runoff_model.Δz
    water = soil.water
    param_functions = soil.param_functions
    hydraulics = water.hydraulics(aux⁻)
    ν = param_functions.porosity
    θ_r = param_functions.water.θ_r(aux⁻)
    S_s = param_functions.water.S_s(aux⁻)

    T = get_temperature(soil.heat, aux⁻, t)
    θ_i = state⁻.soil.water.θ_i
    # Ponding Dirichlet BC
    ϑ_bc = FT(ν - θ_i)
    # Value below surface
    ϑ_below = state⁻.soil.water.ϑ_l

    # Approximate derivative of hydraulic head with respect to z
    ∂h∂z =
        FT(1) +
        (
            pressure_head(hydraulics, ν, S_s, θ_r, ϑ_bc, θ_i) -
            pressure_head(hydraulics, ν, S_s, θ_r, ϑ_below, θ_i)
        ) / Δz

    S_l = effective_saturation(ν, ϑ_bc, θ_r)
    f_i = θ_i / (ϑ_bc + θ_i)
    impedance_f = impedance_factor(water.impedance_factor, f_i)
    viscosity_f = viscosity_factor(water.viscosity_factor, T)
    moisture_f = moisture_factor(water.moisture_factor, hydraulics, S_l)
    K = hydraulic_conductivity(
        param_functions.water.Ksat(aux⁻),
        impedance_f,
        viscosity_f,
        moisture_f,
    )

    i_c = n̂ * (K * ∂h∂z)
    if incident_water_flux < -norm(i_c) # More negative if both are negative,
        #ponding BC
        K∇h⁺ = i_c
    else

        K∇h⁺ = n̂ * (-FT(2) * incident_water_flux) - diff⁻.soil.water.K∇h
    end
    return K∇h⁺
end

"""
    function compute_surface_grad_bc(soil::SoilModel,
                                     runoff_model::NoRunoff,
                                     precip_model::AbstractPrecipModel,
                                     n̂,
                                     state⁻::Vars,
                                     aux⁻::Vars,
                                     diff⁻::Vars,
                                     t::Real
                                     )

Given a runoff model and a precipitation distribution function, compute
the surface water flux normal to the surface. In this case, no runoff is
assumed. The direction of the flux (inwards or outwards) depends on the
magnitude of evaporation, precipitation, and infiltration.
"""
function compute_surface_grad_bc(
    soil::SoilModel,
    runoff_model::NoRunoff,
    precip_model::AbstractPrecipModel,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    t::Real,
)
    FT = eltype(state⁻)
    precip_vector = (FT(0), FT(0), precip_model(t))
    incident_water_flux = dot(precip_vector, n̂)
    K∇h⁺ = n̂ * (-FT(2) * incident_water_flux) - diff⁻.soil.water.K∇h
    return K∇h⁺
end

end


text:
# include the include ...
include("Include.jl")

# ======================================================================================== #
# PHASE 1: SERVER CONFIGURATION
# load the configuration file -
configuration_dictionary = TOML.parsefile(joinpath(_PATH_TO_CONFIG, "Configuration.toml"))
port_number = configuration_dictionary["server"]["port"]
host_ip_address = Sockets.getipaddr()

# show message to user -
start_message = "Server started. address: $(host_ip_address) port: $(port_number) at $(now())";
@info start_message

# setup shared resources
const PALIQUANT_SERVER_ROUTER = HTTP.Router()       # router handles HTTP request/response traffic
const PALIQUANT_SERVER_SESSION = Dict{String,Any}() # session is a temp storage area

# setup logging -
path_to_log_file = joinpath(_PATH_TO_LOG, configuration_dictionary["log"]["log_file_name"])
io = open(path_to_log_file, "w+")
simpleLogger = SimpleLogger(io)
with_logger(simpleLogger) do
    @info(start_message)
    flush(io)
end
# ======================================================================================== #

# ======================================================================================== #
# PHASE 2: DATABASE CONFIGURATION
const dbname = configuration_dictionary["database"]["name"]

# connect to the database -
const path_to_database_file = joinpath(_PATH_TO_DATABASE, dbname)
const DB_CONNECTION = SQLite.DB(path_to_database_file)

with_logger(simpleLogger) do
    @info("Server startup ... database connection init")
    flush(io)
end
# ======================================================================================== #

# ======================================================================================== #
# PHASE 3: API ENDPOINTS AND SERVER START
try

    # endpoints -
    version = 1
    pq_base_url = "/paliquant/transaction/api/v$(version)"

    # test endpoint -
    HTTP.@register(PALIQUANT_SERVER_ROUTER, "POST", "$(pq_base_url)/test", test_api_endpoint)

    # start the server -
    HTTP.serve(PALIQUANT_SERVER_ROUTER, host_ip_address, port_number)

finally

    shutdown_msg = "Shutting down the pqserver at $(now())"
    with_logger(simpleLogger) do
        @info(shutdown_msg)
        flush(io)
        close(io)
    end
end
# ======================================================================================== #


text:
@inline function eval_results(
    df::DataFrame
)
    simple_data = df  |> @map({_.processor, _.algorithm, _.num_samples, _.num_ants, _.num_correlators, _.TrialObj}) |> @orderby_descending({_.num_samples}) |> DataFrame
    df_pretty = pretty_table(simple_data)
end

max_stars_repo_name: mfkiwl/GPUAcceleratedTracking

max_issues_repo_name: mfkiwl/GPUAcceleratedTracking

max_forks_repo_name: mfkiwl/GPUAcceleratedTracking


text:
using Gaston
g = Gaston;
set_terminal("x11");

c = g.CurveConf();
c.legend = "Random";
c.plotstyle = "errorbars";

y = exp(-(1:.1:4.9));
ylow = y - 0.05*rand(40);
yhigh = y + 0.05*rand(40);
g.addcoords(1:40,y,c);
g.adderror(0.1*rand(40));

a = g.AxesConf();
a.title = "Error bars (ydelta)";
g.addconf(a);
g.llplot();

max_stars_repo_path: Module 3/Chapter07/gaston-ex1.jl

max_issues_repo_path: Module 3/Chapter07/gaston-ex1.jl

max_forks_repo_path: Module 3/Chapter07/gaston-ex1.jl


text:
module TestSchedules

using  Test
using  NFLTables: Schedules

@testset "Valid Seasons" begin
    @test length(Schedules.SEASONS) == 12
end

@testset "Test schedule retrieval" begin
    df = Schedules.schedule(2019)
    @test size(df, 2) == 9
end

end  # module TestSchedules


text:
using DecomposedMDPSolver
using POMDPs, POMDPModels
using Test
using Random
using Flux

mdp = SimpleGridWorld(tprob = 1, discount = 1)
POMDPs.gen(mdp::SimpleGridWorld, s, a, rng = Random.GLOBAL_RNG) = (sp =rand(transition(mdp, s, a )), r=reward(mdp, s, a))

Na = 1
Np = 10
solutions = [(x) -> 0.5*ones(Na) for i=1:Np ]

## Weights network
base = Chain(Dense(2, 32, relu), Dense(32, Na))
attn = Chain(Dense(2, 32, relu), Dense(32, Np+1), softmax)
a2t_model = A2TNetwork(base, attn, solutions)

action_prob(mdp, s, a) = 0.25
p = ISPolicy(a2t_model, mdp, action_prob, Random.GLOBAL_RNG)
@test p.network == a2t_model
@test p.mdp isa SimpleGridWorld
@test p.action_probability(mdp, rand(initialstate_distribution(mdp)), :up) == 0.25
@test p.rng == Random.GLOBAL_RNG

@test length(value(p, GWPos(3,3))) == 1

@test value(p, GWPos(3,3), :up) == value(p, GWPos(3,4))

a, prob = action_and_probability(p, GWPos(2,2))
@test a in [:up, :left, :right, :down]
@test prob isa Float64
@test action(p, GWPos(2,2)) in [:up, :left, :right, :down]


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export ProductWeight

"""
`ProductWeight(x...)`

| Set                                               |           ``\\oplus``                                                     |   ``\\otimes``                                                              | ``\\bar{0}``                                                 | ``\\bar{1}``                                                 |
|:-------------------------------------------------:|:-------------------------------------------------------------------------:|:---------------------------------------------------------------------------:|:------------------------------------------------------------:|:------------------------------------------------------------:|
|``\\mathbb{W}_1\\times \\dots\\times\\mathbb{W}_N``|  ``\\oplus_{\\mathbb{W}_1} \\times \\dots\\times\\oplus_{\\mathbb{W}_N}`` |  ``\\otimes_{\\mathbb{W}_1} \\times \\dots\\times\\otimes_{\\mathbb{W}_N}`` |``(\\bar{0}_{\\mathbb{W}_1},\\dots,\\bar{0}_{\\mathbb{W}_N})``|``(\\bar{1}_{\\mathbb{W}_1},\\dots,\\bar{1}_{\\mathbb{W}_N})``|

"""
struct ProductWeight{N, T <:NTuple{N,Semiring} } <: Semiring
  x::T
end
ProductWeight(x...) = ProductWeight(x)
reversetype(::Type{S}) where {N,T,S <: ProductWeight{N,T}} =
# in future we can use fieldtypes instead of ntuple(i -> fieldtype(T, i), fieldcount(T))
ProductWeight{N,Tuple{reversetype.(ntuple(i -> fieldtype(T, i), fieldcount(T)))...}}
ProductWeight{N,T}(x::ProductWeight{N,T}) where {N,T <:NTuple{N,Semiring}} = x

zero(::Type{ProductWeight{N,T}}) where {N,T}  = ProductWeight{N,T}(zero.(ntuple(i -> fieldtype(T, i), fieldcount(T))))
one(::Type{ProductWeight{N,T}}) where {N,T}  = ProductWeight{N,T}(one.(ntuple(i -> fieldtype(T, i), fieldcount(T))))

*(a::ProductWeight{N,T}, b::ProductWeight{N,T}) where {N,T} = ProductWeight{N,T}(a.x .* b.x)
+(a::ProductWeight{N,T}, b::ProductWeight{N,T}) where {N,T} = ProductWeight{N,T}(a.x .+ b.x)
/(a::ProductWeight{N,T}, b::ProductWeight{N,T}) where {N,T} = ProductWeight{N,T}(a.x ./ b.x)
reverse(a::ProductWeight{N,T}) where {N,T} = ProductWeight{N,T}t(reverse.(a.x))

# properties
for p in [:isleft,
          :isright,
          :isweaklydivisible,
          :ispath,
          :isidempotent,
          :iscommulative,
          :iscomplete]
  @eval $p(::Type{ProductWeight{N,T}}) where {N,T} = all( $p.(ntuple(i -> fieldtype(T, i), fieldcount(T)))  )
end


text:
# Checks output from chernikova_general function an input matrix
# defines a specified cone. That is,
# { x: Ax ≥ 0 } == { B y + C z : y ≥ 0, z free }
#
# Arguments:
# A::Matrix - constraint matrix (rows correspond to individual constraints)
# bidrays::Matrix - matrix of unidirectional rays (columns correspond to individual rays)
# unirays::Matrix - matrix of (normalised) bidirectional rays (columns correspond to individual rays)
function test_chernikova_general(A::Matrix{T1}, bidrays::Matrix{T2}, rays::Matrix{T2}) where {T1<:Union{Integer, Rational}, T2<:Integer}
    norm_rays = Cones.normalise_rays(rays, bidrays)
    chern_bidrays, chern_rays = chernikova_general(A)
    norm_chern_rays = Cones.normalise_rays(chern_rays, chern_bidrays)
    @test Cones.check_columns_same(norm_rays, norm_chern_rays)
    @test Cones.check_column_span_same(bidrays, chern_bidrays)
end

# Pointed 2D cone
println("\tExample 1...")
A = [[-1//5 4//5];
     [4//5 -1//5]]
rays = [[1,4] [4, 1]]
bidrays = Array{Int}(undef, 2, 0)
test_chernikova_general(A, bidrays, rays)

# Hyperplane
println("\tExample 2...")
A = [[-1 -1 0]; [1 1 0]]
rays = Array{Int}(undef, 3, 0)
bidrays = [[1, -1, 0] [0, 0, 1]]
test_chernikova_general(A, bidrays, rays)

# 1D lineality space
println("\tExample 3...")
A = [[-3 3 0]; [2 2 0]]
rays = [[-1, 1, 1] [2, 2, -1]]
bidrays = reshape([0, 0, 1], 3, 1)
test_chernikova_general(A, bidrays, rays)

## function unirays{T<:Integer}(bidrays::Matrix{T}, unirays::Matrix{T})
##     d, n1 = size(bidrays)
##     n2 = size(unirays, 2)
##     rays = Array{T}(undef, d, 2*n1 + n2)
##     for i in 1:n1
##         rays[:, 2*i - 1] = bidrays[:,i]
##         rays[:, 2*i] = -bidrays[:,i]
##     end
##     rays[:, (2*n1 + 1):(2*n1 + n2)] = unirays
##     return rays
## end

## # Function which checks anything generated from Chernikova output is in Polyhedral cone
## function test_polyhedral_contains_finite_cone(A::Matrix{Int}, num_points::Int)
##     X = chernikova(A)
##     dim, gens = size(X)
##     dist = Uniform(0.0, 100.0)
##     B = [A; eye(Int,dim)]
##     for i in 1:num_points
##         λ = rand(dist, gens)
##         @test all(B*(X*λ) .>= 0.0)
##     end
## end


text:
lexi=ARGS[1]
lexi_out=ARGS[2]

function main(lexicon)
  # remove phones
  L = read(lexi, String);
  D = Dict()
  D["AA0"] = D["AA1"] = D["AA2"] = "AA"
  D["AE0"] = D["AE1"] = D["AE2"] = "AE"
  D["AH0"] = D["AH1"] = D["AH2"] = "AH"
  D["AO0"] = D["AO1"] = D["AO2"] = "AO"
  D["AW0"] = D["AW1"] = D["AW2"] = "AW"
  D["AY0"] = D["AY1"] = D["AY2"] = "AY"
  D["EH0"] = D["EH1"] = D["EH2"] = "EH"
  D["ER0"] = D["ER1"] = D["ER2"] = "ER"
  D["EY0"] = D["EY1"] = D["EY2"] = "EY"
  D["IH0"] = D["IH1"] = D["IH2"] = "IH"
  D["IY0"] = D["IY1"] = D["IY2"] = "IY"
  D["OW0"] = D["OW1"] = D["OW2"] = "OW"
  D["OY0"] = D["OY1"] = D["OY2"] = "OY"
  D["SH0"] = D["SH1"] = D["SH2"] = "SH"
  D["UH0"] = D["UH1"] = D["UH2"] = "UH"
  D["UW0"] = D["UW1"] = D["UW2"] = "UW"

  for d in D
    L = replace(L,d)
  end
  return L
end

L = main(lexi)
write(lexi_out, L)


text:
macro def(name, definition)
  esc(quote
    macro $name()
      esc($(Expr(:quote, definition)))
    end
  end)
end

@def add_linesearch_fields begin
    x_ls::Tx
    alpha::T
end

@def initial_linesearch begin
    (similar(initial_x), # Buffer of x for line search in state.x_ls
    real(one(T)))             # Keep track of step size in state.alpha
end


text:
"""
    initialize(csv::String) -> Nothing

If the file along the `csv` path does not exist, then create it and write the header;
if it exists but is empty, do the same; if it exists and is not empty, check structural
correctness.

# Arguments
- `csv::String`: $(CSV_ARGUMENT)

# Throws
- [`NotCSV`](@ref): $(NOT_CSV_EXCEPTION)
- Exceptions from: [`_check`](@ref)

# Example
```jldoctest; output = false
using PDFHighlights
HEADER = PDFHighlights.Internal.Constants.HEADER

_file, _ = mktemp()
file = _file * ".csv"
mv(_file, file)

initialize(file)

open(file, "r") do io
    readlines(io) == [HEADER]
end

# output

true
```
"""
function initialize(csv::String)::Nothing

    # Check that the path has the correct extension
    !endswith(csv, ".csv") && throw(NotCSV(csv))

    # Check if the csv exists
    if isfile(csv)

        if length(readlines(csv)) == 0

            # Write the header
            open(csv, "w") do io
                println(io, HEADER)
            end

        else

            # Check if the structure of the table is correct
            _check(csv)

        end

    else

        # Write the header
        open(csv, "w") do io
            println(io, HEADER)
        end

    end

    return nothing

end


text:
using PyCall
using FlightSims
# using FSimBase
using FSimZoo
using FSimROS
using FSimPlots
using Plots
ENV["GKSwstype"]="nul"
theme(:lime)
using UnPack
using OnlineStats
using DataFrames
using ComponentArrays

struct Multicopter_ZOH_Input <: AbstractEnv
    multicopter::Multicopter
end

function FSimBase.State(env::Multicopter_ZOH_Input)
    State(env.multicopter)
end

function FSimBase.Dynamics!(env::Multicopter_ZOH_Input)
    @Loggable function dynamics!(dx, x, input, t)
        @nested_log FSimZoo._Dynamics!(env.multicopter)(dx, x, nothing, t; u=input)
    end
end

struct ControllerEnv <: AbstractEnv
    controller::BacksteppingPositionController
    allocator::StaticAllocator
    multicopter::Multicopter  # TODO: not use the dummy multicopter; receive data from other nodes
end

function ControllerEnv(pos_cmd_func=nothing)
    multicopter = LeeHexacopter()
    @unpack m, B = multicopter
    controller = BacksteppingPositionController(m; pos_cmd_func=pos_cmd_func)
    allocator = PseudoInverseAllocator(B)
    env = ControllerEnv(controller, allocator, multicopter)
end

function FSimBase.State(env::ControllerEnv)
    @unpack controller, allocator, multicopter = env  # TODO: receive multicopter data from other nodes
    @unpack m, g = multicopter
    function state(; args_multicopter=())
        x_multicopter = State(multicopter)(args_multicopter...)
        pos0 = copy(x_multicopter.p)
        x_controller = State(controller)(pos0, m, g)
        ComponentArray(multicopter=x_multicopter, controller=x_controller)
    end
end

function FSimBase.Dynamics!(env::ControllerEnv)
    @unpack controller, allocator, multicopter = env
    @unpack m, J, g = multicopter
    @Loggable function dynamics!(dx, x, input, t; pos_cmd=nothing)
        @unpack p, v, R, ω = input
        # @unpack p, v, R, ω = x.multicopter
        @unpack ref_model, Td = x.controller
        xd, vd, ad, ȧd, äd = ref_model.x_0, ref_model.x_1, ref_model.x_2, ref_model.x_3, ref_model.x_4
        νd, Ṫd, _... = Command(controller)(
                                           p, v, R, ω,
                                           xd, vd, ad, ȧd, äd, Td,
                                           m, J, g,
                                          )
        @nested_log :controller Dynamics!(controller)(dx.controller, x.controller, (), t; pos_cmd=pos_cmd, Ṫd=Ṫd)
        @nested_log :controller νd
    end
end

function plot_stat(x)
	fig_stat = plot()
        stat_E = MovingWindow(Float64, 500)
        stat_N = MovingWindow(Float64, 500)
        stat_U = MovingWindow(Float64, 500)
	fit!(stat_E, x.p[2])
	fit!(stat_N, x.p[1])
	fit!(stat_U, -x.p[3])
	fig_stat_E = plot(1:minimum([stat_E.b, stat_E.n]), stat_E.value; title="E [m]", label=nothing)
	fig_stat_N = plot(1:minimum([stat_N.b, stat_N.n]), stat_N.value; title="N [m]", label=nothing)
	fig_stat_U = plot(1:minimum([stat_U.b, stat_U.n]), stat_U.value; title="U [m]", label=nothing)
	fig_stat = plot(fig_stat_E, fig_stat_N, fig_stat_U; layout=(3, 1))
	fig_stat
end

function plot_traj(multicopter, x)
	fig_multicopter = plot()
	plot!(fig_multicopter, multicopter, x; xlim=(-2, 2), ylim=(-2, 2), zlim=(-1, 10))
	fig_multicopter
end

function plot_figs(multicopter, x)
	fig_multicopter = plot_traj(multicopter, x)
	fig_stat = plot_stat(x)
	# display(plot(fig_multicopter, fig_stat, layout=(1, 2)))
end

function main()
	t = 0.01
	multicopter = LeeHexacopter()
	env = Multicopter_ZOH_Input(multicopter)
	state0 = State(env)()
        @unpack m, g = env.multicopter
        u0 = (m*g/6) * ones(6)  # initial control input (actually meaningless; for initialisation of simulator)
        simulator = Simulator(state0, Dynamics!(env), u0; tf=100)  # second
        pos_cmd_func = t -> [sin(t), cos(t), 0.0]
        env_ctrl = ControllerEnv(pos_cmd_func)
        state0_ctrl = State(env_ctrl.multicopter)()  # dummy initial state of multicopter
        simulator_ctrl = Simulator(State(env_ctrl)(), Dynamics!(env_ctrl), state0_ctrl; tf=100)  # second

	# state update
	println("state update test...")
	step_until!(simulator, t)
	# control update
	println("control update test...")
	df = DataFrame()
	push!(simulator_ctrl, df)
	νd = df.sol[end].controller.νd  # virtual control
	step_until!(simulator_ctrl, t)

	# plot figures
	println("plotting test...")
	plot_figs(multicopter, state0)
end

main()


text:
POMDPs.states(p::LaserTagPOMDP) = states(p.floor)

function states(f::Floor) # inefficient?
    vec = Array{LTState}(undef, n_pos(f)^2+1)
    for i in 1:f.n_cols, j in 1:f.n_rows, k in 1:f.n_cols, l in 1:f.n_rows
        s = LTState(Coord(i,j), Coord(k,l), false)
        vec[stateindex(f, s)] = s
    end
    vec[end] = LTState(Coord(1,1), Coord(1,1), true)
    return vec
end

n_pos(f::Floor) = f.n_rows*f.n_cols

# POMDPs.n_states(p::LaserTagPOMDP) = n_pos(p.floor)^2+1

POMDPs.stateindex(p::LaserTagPOMDP, s::LTState) = stateindex(p.floor, s)

function stateindex(f::Floor, s::LTState)
    nr = f.n_rows
    nc = f.n_cols
    n_pos = nr^2*nc^2
    if s.terminal
        return n_pos+1
    else
        rob = s.robot
        opp = s.opponent
        return LinearIndices((nc,nr,nc,nr))[rob[1], rob[2], opp[1], opp[2]]
    end
end


text:
#using Revise
using Test
using JDF, TimeZones, DataFrames

a = some_elm(ZonedDateTime)

ar = Vector{ZonedDateTime}(undef, 10)
ar .= a

meta = open("plsdel.io", "w") do io
  compress_then_write(ar, io)
end

ar_loaded = open("plsdel.io") do io
  column_loader(meta.type, io, meta)
end

@test ar == ar_loaded

df = DataFrame(ar = ar)

savejdf(df, "df.jdf")

df_loaded = loadjdf("df.jdf")

@test df == df_loaded

ar_timezone = [ar.timezone for ar in ar]
meta = open("plsdel.io", "w") do io
  compress_then_write(ar_timezone, io)
end

ar_timezone_copy = open("plsdel.io") do io
  column_loader(meta.type, io, meta)
end

@test ar_timezone == ar_timezone_copy

meta = open("plsdel.io", "w") do io
  compress_then_write(a.timezone.transitions, io)
end

a_timezone_transitions_loaded = open("plsdel.io") do io
  column_loader(meta.type, io, meta)
end

@test a.timezone.transitions == a_timezone_transitions_loaded

rm("plsdel.io", force=true)
rm("df.jdf", force=true, recursive=true)


text:
module Symbolics

using DocStringExtensions

using LinearAlgebra

using Reexport

@reexport using SymbolicUtils

import SymbolicUtils: Term, Add, Mul, Pow, Sym, symtype, to_symbolic,
                      FnType, @rule, Rewriters, substitute, similarterm,
                      promote_symtype, istree, operation, arguments

import SymbolicUtils.Rewriters: Chain, Prewalk, Postwalk, Fixpoint

# re-export

export simplify, substitute

export Num
include("num.jl")

include("utils.jl")

using MacroTools
include("register.jl")

using TreeViews
export @variables
include("variable.jl")

include("linearity.jl")

using DiffRules, SpecialFunctions, NaNMath

using SparseArrays

export Differential, expand_derivatives, gradient,
       jacobian, jacobian_sparsity, sparsejacobian,
       hessian, sparsehessian, hessian_sparsity

include("diff.jl")

include("linear_algebra.jl")

end # module


text:
function write_dot(path::AbstractString, pg::PropertyGraph)
    open(path, "w") do io
        println(io, "digraph G {")
        println(io, "rankdir=LR;")

        foreach(collect(pg.vmap)) do (label, code)
            println(io, "$code [label=$label];")
        end

        for e in edges(pg)
            u, v = src(e), dst(e)
            println(io, "$u -> $v;")
        end

        println(io, "}")
    end
end


text:
get_versions(res_name::AbstractString) =
    Core.find_versions(get_working_dir(), res_name);
function get_version(res_name::AbstractString, ver_num::Int)
    vers = get_versions(res_name);
    if 0 < ver_num <= length(vers)
        return vers[ver_num];
    else
        error("$(ver_num) is not a valid ver number, it must be "*
            "between 0 and $(length(vers))!!!")
    end
end

has_versions(res_name::AbstractString) = length(get_versions(res_name)) > 0;

get_last_version(res_name::AbstractString) = get_versions(res_name)[end];


text:
export @ham, FusedHamiltonian

import IterTools

struct FusedHamiltonian{D<:Tuple} <: AbstractHamiltonian
    data::D
end

FusedHamiltonian(h::LocalHamiltonian...) = FusedHamiltonian(h)
function (h::FusedHamiltonian)(lattice::AbstractLattice, rhs::AbstractSites)
    itrs = []
    for each in h.data
        push!(itrs, RegionIterator(each, lattice, rhs))
    end
    IterTools.chain(itrs...)
end

# TODO: refactor this
function make_hamiltonian(ex::Expr)
    expr = Expr(:call, :FusedHamiltonian)
    if Meta.isexpr(ex, :call)
        if ex.args[1] == :sum # only one region
            ex.args[1] = :LocalHamiltonian
            return ex
        elseif ex.args[1] == :+
            for each in ex.args[2:end]
                if each.args[1] == :sum
                    each.args[1] = :LocalHamiltonian
                    push!(expr.args, each)
                elseif each.args[1] == :*
                    factor = each.args[2]
                    each.args[3].args[1] = :LocalHamiltonian
                    each.args[3].args[end] = Expr(:call, :*, factor, each.args[3].args[end])
                    push!(expr.args, each.args[3])
                else
                    throw(ErrorException("Invalid Expression: $each"))
                end
            end
        elseif ex.args[1] == :-
            # preserve sign for the first one
            each = ex.args[2]
            if each.args[1] == :sum
                each.args[1] = :LocalHamiltonian
                push!(expr.args, each)
            elseif each.args[1] == :*
                factor = each.args[2]
                each.args[3].args[1] = :LocalHamiltonian
                each.args[3].args[end] = Expr(:call, :*, factor, each.args[3].args[end])
                push!(expr.args, each.args[3])
            else
                throw(ErrorException("Invalid Expression: $each"))
            end

            for each in ex.args[3:end]
                if each.args[1] == :sum
                    each.args[1] = :LocalHamiltonian
                    each.args[end] = Expr(:call, :*, -1, each.args[end]) # merge -1 to local matrix
                    push!(expr.args, each)
                elseif each.args[1] == :*
                    factor = each.args[2]
                    each.args[3].args[1] = :LocalHamiltonian
                    each.args[3].args[end] = Expr(:call, :*, -1, factor, each.args[3].args[end])
                    push!(expr.args, each.args[3])
                else
                    throw(ErrorException("Invalid Expression: $each"))
                end
            end
        end
    else
        throw(ErrorException("Invalid Expression: $each"))
    end
    return expr
end

macro ham(expr::Expr)
    esc(make_hamiltonian(expr))
end


text:
# MIT Mystery Hunt 2014 puzzle Venntersections
# http://www.mit.edu/~puzzle/2014/puzzle/venntersections/

@testset "ventersections" begin

    @testset "diagram 1" begin
        # Set 1
        f1 = best_feature(["lowered", "levitate", "inanimate", "paradise", "leveraged", "sizes", "tuxedo"])
        @test f1.description == "alternates consonant vowel"

        # Set 2
        f2 = best_feature(["leveraged", "sizes", "tuxedo", "lynx", "lightly", "crocodile", "triumph"])
        @test f2.description == "has scrabble score 14"

        # Set 3
        f3 = best_feature(["lowered", "levitate", "leveraged", "lynx", "lightly", "lengths", "legislator"])
        @test f3.description == "contains 'l' at index 1"

        # Set 4
        f4 = best_feature(["levitate", "inanimate", "sizes", "lightly", "crocodile", "legislator", "carousels"], 1)
        @test f4.description== "contains a repeated consonant"

        # Intersection
        checks = (f1, f2, f3, f4, word -> length(word) == 9)
        for word in words
            if all(c(word) for c in checks)
                @test word == "lakesides"
                break
            end
        end

    end

    @testset "diagram 2" begin
        # Set 2
        f2 = best_feature(["questionable", "businesswoman", "exhaustion", "discouraged", "communicated", "hallucinogen", "sequoia"])
        @test f2.description == "has 5 unique vowels"

        # Set 3
        f3 = best_feature(["grimaced", "formally", "questionable", "discouraged", "communicated", "chrysalis", "saccharin"])
        @test f3.description == "contains 'a' at index 4 from end"

        # Set 4
        f4 = best_feature(["formally", "thinnest", "businesswoman", "communicated", "hallucinogen", "saccharin", "cellophane"])
        @test f4.description == "contains 1 double letters"
    end

    @testset "diagram 3" begin
        # Set 1
        f1 = best_feature(["thumbtacks", "monologue", "frigidities", "statuesque", "testimony", "satirizing", "flawed"])
        @test f1.description == "contains a day of the week abbreviation"

        # Set 3
        f3 = best_feature(["thumbtacks", "monologue", "testimony", "camel", "meteorology", "trampoline", "achievement"])
        @test f3.description == "contains 1 'm's"

        # Set 4
        f4 = best_feature(["monologue", "frigidities", "satirizing", "meteorology", "avalance", "achievement", "constitute"])
        @test f4.description == "has 1 letters repeated at least 3 times each"
    end

    @testset "diagram 4" begin
        # Set 1
        f1 = best_feature(["philharmonic", "mischievous", "alphabet", "restaurant", "leeching", "mushroom", "pioneer"])
        @test f1.description == "contains a greek letter"

        # Set 2
        f2 = best_feature(["leeching", "mushroom", "pioneer", "loophole", "toothpaste", "seventeenth", "kneeling"])
        @test f2.description == "contains 1 double letters"

        # Set 3
        f3 = best_feature(["philharmonic", "mischievous", "leeching", "loophole", "toothpaste", "alcoholic", "narwhal"])
        @test f3.description == "contains 'h' at index 5"

        # Set 4
        # Note: it turns out that the "most unlikely" feature of all of the words
        # in this set of the puzzle is just that they contain an 'h'. However, the
        # puzzle actually expects you to choose "has 3 consonants in a row" instead.
        # So I've added a dummy word which does not contain an 'h' to the end of the
        # list.
        f4 = best_feature(["mischievous", "alphabet", "mushroom", "toothpaste", "seventeenth", "narwhal", "chromosome", "aardvark"])
        @test f4.description == "has 3 consonants in a row"

        # Intersection
        checks = (f1, f2, f3, f4, word -> length(word) == 12)
        for word in words
            if all(c(word) for c in checks)
                @test word == "neighborhood"
                break
            end
        end
    end
end


text:
using Test, DiscretePDEs
const DPE = DiscretePDEs
using PyCall: pyimport

@testset "remove_seams" begin
    points = [0,1,2,3,4,5,6,5,7,3,2,8,9]
    seam_end_inds = DPE.identify_seam_end_inds(points)
    @test seam_end_inds == [8,10,11]
    begin_inds, end_inds, recur_ind = DPE.find_latest_ending_seam(points, seam_end_inds)
    @test begin_inds == [3,4]
    @test end_inds == [11,10]
    @test recur_ind == 1
    outside = [points[1:begin_inds[1]]; points[(end_inds[1]+1):end]]
    inside = points[begin_inds[end]:(end_inds[end]-1)]
    seam_end_inds[1:recur_ind] .- (begin_inds[end]-1)
    @test DPE.remove_seams(points) == [[0,1,2,8,9], [3,4,5,7], [5,6]]
end

@testset "demarcate_self_intersections" begin
    points = [(1.65, -0.34), (0.79, -0.34), (0.79, -0.53), (1.38, -0.53),
             (1.38, -0.67), (1.22, -0.67), (1.22, -0.71), (1.48, -0.71), (1.48, -0.76), (1.22, -0.76), (1.22, -0.67),
             (1.1, -0.67), (1.1, -0.8), (1.38, -0.8), (1.38, -0.92), (0.79, -0.92), (0.79, -0.34), (0.55, -0.34), (0.55, -0.71),
             (0.67, -0.71), (0.67, -0.88), (0.55, -0.88), (0.55, -0.34), (0.48, -0.34), (0.48, -1.11), (1.65, -1.11)]
    demarcated_points = [points[1:10]; points[7]; points[11:16]; points[3]; points[17:22]; points[19]; points[23:end]]
    @test DPE.demarcate_self_intersections(points) == demarcated_points
end

@testset "GeoCode" begin
    gdspy = pyimport("gdspy")
    cell = gdspy.Cell("test", exclude_from_current=true)
    cell.add(gdspy.Rectangle((0,0), (4,8)))
    geo_code = GeoCode(cell, layer_extrusion=Dict(0 => 1),
        layer_surface_groups = Dict("surface" => [0]),
        layer_volume_groups = Dict("volume" => [0]))
    @test length(geo_code.points) == 4 # subtract 1 for leading comment in each section
    @test length(geo_code.lines) == 4
    @test length(geo_code.curves) == 1
    @test length(geo_code.planes) == 1
    @test length(geo_code.surface_differences) == 0
    @test length(geo_code.extrusions) == 1
    @test length(geo_code.volume_differences) == 0
    @test length(geo_code.surface_groups) == 1
    @test length(geo_code.volume_groups) == 1
end


text:
using Test
using LinearAlgebra
using WaveProp
using WaveProp.Geometry
using WaveProp.Integration

@testset "Trapezoidal quadrature" begin
    q = Trapezoidal{10}()
    D = domain(q)
    @test D == ReferenceLine()
    x,w = q()
    @test sum(w) ≈ 1
    # integrate a periodic function. Should be very accurate.
    @test isapprox(integrate(x->cos(2π*x),q),0,atol=1e-10)
    @test integrate(x->sin(2π*x[1])^2,q) ≈ 0.5
end

@testset "TrapezoidalP quadrature" begin
    q = TrapezoidalP{10}()
    D = domain(q)
    @test D == ReferenceLine()
    x,w = q()
    @test sum(w) ≈ 1
    # integrate a periodic function. Should be very accurate.
    @test isapprox(integrate(x->cos(2π*x[1]),q),0,atol=1e-10)
    @test integrate(x->sin(2π*x[1])^2,q) ≈ 0.5
end

@testset "Fejer quadrature" begin
    N = 5
    q = Fejer{N}()
    x,w = q()
    @test sum(w) ≈ 1
    # integrate all polynomial of degree N-1 exactly
    for n in 1:(N-1)
        @test integrate(x->x[1]^n,q) ≈ 1/(n+1)
    end
end

@testset "Gauss-Legendre quadrature" begin
    N = 5
    q = GaussLegendre{N}()
    x,w = q()
    @test sum(w) ≈ 1
    # integrate all polynomial of degree 2N-1 exactly
    for n in 1:(2*N-1)
        @test integrate(x->x[1]^n,q) ≈ 1/(n+1)
    end
end

@testset "Gauss quad on triangle" begin
    D= ReferenceTriangle
    for N in (1,3)
        q = Gauss{D,N}()
        x,w = q()
        @test sum(w) ≈ 1/2
    end
    # FIXME: check that we integrate all monomials up to `order` on the triangle
end

@testset "Gauss quad on tetrahedron" begin
    D= ReferenceTetrahedron
    for N in (1,4)
        q = Gauss{D,N}()
        x,w = q()
        @test sum(w) ≈ 1/6
    end
    # FIXME: check that we integrate all monomials up to `order`
end

@testset "Tensor product quad on square" begin
    D   = ReferenceSquare
    N,M = 10,12
    qx  = GaussLegendre{N}()
    qy  = Fejer{M}()
    q   = TensorProductQuadrature(qx,qy)
    a,b = 2*N-1,M-1 # maximum integration order of monomials
    @test integrate( x -> 1,q) ≈ 1
    f = x -> x[1]^a*x[2]^b
    @test integrate(f,q) ≈ 1/(a+1)*1/(b+1)
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
######################### Global Time Stepping ###########################
function forwardEuler(mesh::Mesh, fluxResidualFn, sln::SolutionState, boundaryConditions, fluid::Fluid, dt)
    sln.fluxResiduals = fluxResidualFn(mesh, sln, boundaryConditions, fluid)
    @fastmath sln.cellState .+= sln.fluxResiduals.*dt
    @fastmath decodeSolution_3D(sln, fluid)

    return sln
end

function RK2_Mid(mesh, fluxResidualFn, sln, boundaryConditions, fluid::Fluid, dt)
    fluxResiduals1 = fluxResidualFn(mesh, sln, boundaryConditions, fluid)
    halfwayEstimate = sln.cellState .+ fluxResiduals1.*dt/2
    solutionState2 = SolutionState(halfwayEstimate, sln.cellFluxes, sln.cellPrimitives, sln.fluxResiduals, sln.faceFluxes)
    decodeSolution_3D(solutionState2, fluid)

    sln.fluxResiduals = fluxResidualFn(mesh, solutionState2, boundaryConditions, fluid)
    sln.cellState .+= sln.fluxResiduals.*dt
    decodeSolution_3D(sln, fluid)

    return sln
end

function RK4(mesh, fluxResidualFn, sln, boundaryConditions, fluid::Fluid, dt)

    fluxResiduals1 = fluxResidualFn(mesh, sln, boundaryConditions, fluid)
    halfwayEstimate = sln.cellState .+ fluxResiduals1*dt/2
    lastSolutionState = SolutionState(halfwayEstimate, sln.cellFluxes, sln.cellPrimitives, sln.fluxResiduals, sln.faceFluxes)
    decodeSolution_3D(lastSolutionState, fluid)

    fluxResiduals2 = fluxResidualFn(mesh, lastSolutionState, boundaryConditions, fluid)
    halfwayEstimate2 = sln.cellState .+ fluxResiduals2*dt/2
    lastSolutionState.cellState = halfwayEstimate2
    decodeSolution_3D(lastSolutionState, fluid)

    fluxResiduals3 = fluxResidualFn(mesh, lastSolutionState, boundaryConditions, fluid)
    finalEstimate1 = sln.cellState .+ fluxResiduals3*dt
    lastSolutionState.cellState = finalEstimate1
    decodeSolution_3D(lastSolutionState, fluid)

    fluxResiduals4 = fluxResidualFn(mesh, lastSolutionState, boundaryConditions, fluid)
    sln.cellState .+= (fluxResiduals1 .+ 2*fluxResiduals2 .+ 2*fluxResiduals3 .+ fluxResiduals4 )*(dt/6)
    decodeSolution_3D(sln, fluid)

    return sln
end

function ShuOsher(mesh, fluxResidualFn, sln, boundaryConditions, fluid::Fluid, dt)

    fluxResiduals1 = fluxResidualFn(mesh, sln, boundaryConditions, fluid)
    endEstimate = sln.cellState .+ fluxResiduals1.*dt
    lastSolutionState = SolutionState(endEstimate, sln.cellFluxes, sln.cellPrimitives, sln.fluxResiduals, sln.faceFluxes)
    decodeSolution_3D(lastSolutionState, fluid)

    fluxResiduals2 = fluxResidualFn(mesh, lastSolutionState, boundaryConditions, fluid)
    estimate2 = (3/4).*sln.cellState .+ (1/4).*(endEstimate .+ fluxResiduals2.*dt)
    lastSolutionState.cellState = estimate2
    decodeSolution_3D(lastSolutionState, fluid)

    fluxResiduals3 = fluxResidualFn(mesh, lastSolutionState, boundaryConditions, fluid)
    sln.cellState .= (1/3).*sln.cellState .+ (2/3).*(estimate2 .+ dt.*fluxResiduals3)
    decodeSolution_3D(sln, fluid)

    return sln
end

######################### Local Time Stepping ###########################
# Incomplete, will be commented more fully once it produces nice solutions and the implementation is finalized
function LTSEuler(mesh, fluxResidualFn, sln, boundaryConditions, fluid::Fluid, dt)
    targetCFL = dt[1]

    fluxResiduals = fluxResidualFn(mesh, sln, boundaryConditions, fluid)

    CFL!(dt, mesh, sln, fluid, 1)
    dt .= targetCFL ./ dt
    smoothTimeStep!(dt, mesh, 0.1)
    smoothTimeStep!(dt, mesh, 0.1)
    sln.cellState .+= fluxResiduals .* dt
    decodeSolution_3D(sln, fluid)

    return sln
end

function smoothTimeStep!(dt, mesh::Mesh, diffusionCoefficient=0.2)
    nCells, nFaces, nBoundaries, nBdryFaces = unstructuredMeshInfo(mesh)

    timeFluxes = zeros(nCells)
    surfaceAreas = zeros(nCells)
    for f in 1:nFaces-nBdryFaces
        ownerCell = mesh.faces[f][1]
        neighbourCell = mesh.faces[f][2]
        timeFlux = (dt[ownerCell] - dt[neighbourCell]) * mag(mesh.fAVecs[f])
        surfaceAreas[ownerCell] += mag(mesh.fAVecs[f])
        surfaceAreas[neighbourCell] += mag(mesh.fAVecs[f])
        timeFluxes[ownerCell] -= timeFlux
        timeFluxes[neighbourCell] += timeFlux
    end

    timeFluxes .*= (diffusionCoefficient ./ surfaceAreas)

    for i in eachindex(timeFluxes)
        timeFluxes[i] = min(0, timeFluxes[i])
    end

    dt .+= timeFluxes
end

#TODO: For implicit methods, need to compute the flux Jacobians at each edge, instead of just the fluxes
    # Use Jacobians as coefficients in matrix representing timestepping equations
    # Then solve with GMRES or some other matrix solver

max_stars_repo_name: henrystoldt/fvCFD


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

JLLWrappers.@generate_main_file_header("Minuit2")
JLLWrappers.@generate_main_file("Minuit2", UUID("683d4edd-6ccf-5a74-985b-ebee5a2d6cb7"))
end  # module Minuit2_jll


text:
using QuNet
using DataFrames
using CSV
using Statistics

# Graph + sampling parameters
grid_size = 100
num_pairs = 50
num_trials = 100
max_paths = 1

# Get the "coordinates" associated with end users
#   that is, the associated efficiency and fidelity of the paths connecting them (if any)
edge_costs = Dict("loss"=>0.05, "Z"=>0.05)
net = GridNetwork(grid_size, grid_size, edge_costs=edge_costs)
#usage:
# (network::Union{QNetwork, QuNet.TemporalGraph},
#    num_trials::Int64, num_pairs::Int64; max_paths=3, src_layer::Int64=-1,
#    dst_layer::Int64=-1, edge_perc_rate=0.0)
coord_list = QuNet.heat_data(net, num_trials, num_pairs, max_paths=max_paths)

# List comprehension to extract data
e = [coord_list[i][1] for i in 1:length(coord_list)]
f = [coord_list[i][2] for i in 1:length(coord_list)]

# # Convert cood_list to data_frame and write to csv
df = DataFrame(Efficiency = e, Fidelity = f)
CSV.write("data/heatmap.csv", df)
# Continue on from heatmap.py to generate the plot

"""
This function gets the end-to-end failure rate (P₀) of the graph in the heatmap
plot. By my calculations, this is about 0.20
"""
function get_P0_for_graph()
    num_trials = 100
    max_paths = 4

    failure_counts = []
    for i in 1:num_trials
        println("num_trial $i")
        edge_costs = Dict("loss"=>0.05, "Z"=>0.05)
        net = GridNetwork(grid_size, grid_size, edge_costs=edge_costs)
        user_pairs = QuNet.make_user_pairs(net, 50)
        pathset, pur_paths, pathuse_count = QuNet.greedy_multi_path!(net, QuNet.purify, user_pairs, 4)
        push!(failure_counts, pathuse_count[1])
    end
    return(mean(failure_counts./num_pairs))
end

max_stars_repo_name: drpeterrohde/QuNet

max_issues_repo_name: SamuelMarks/QuNet

max_forks_repo_name: SamuelMarks/QuNet


text:
g = PathGraph(4)
add_vertices!(g,10)
add_edge!(g,5,6)
add_edge!(g,6,7)
add_edge!(g,8,9)
add_edge!(g,10,9)

@test !is_connected(g)
@test is_connected(g6)

cc = connected_components(g)
label = zeros(Int, nv(g))
LightGraphs.connected_components!(label, g)
@test label[1:10] == [1,1,1,1,5,5,5,8,8,8]
import LightGraphs: components, components_dict
cclab = components_dict(label)
@test cclab[1] == [1,2,3,4]
@test cclab[5] == [5,6,7]
@test cclab[8] == [8,9,10]
@test length(cc) >= 3 && sort(cc[3]) == [8,9,10]

g10 =DiGraph(4)
add_edge!(g10,1,3)
add_edge!(g10,2,4)
@test is_bipartite(g10) == true
add_edge!(g10,1,4)
@test is_bipartite(g10) == true

g10 = DiGraph(20)
for m=1:50
    i = rand(1:10)
    j = rand(11:20)
    if rand() < 0.5
        i, j = j, i
    end
    if !has_edge(g10, i, j)
        add_edge!(g10, i, j)
        @test is_bipartite(g10) == true
    end
end

# graph from https://en.wikipedia.org/wiki/Strongly_connected_component
h = DiGraph(8)
add_edge!(h,1,2); add_edge!(h,2,3); add_edge!(h,2,5);
add_edge!(h,2,6); add_edge!(h,3,4); add_edge!(h,3,7);
add_edge!(h,4,3); add_edge!(h,4,8); add_edge!(h,5,1);
add_edge!(h,5,6); add_edge!(h,6,7); add_edge!(h,7,6);
add_edge!(h,8,4); add_edge!(h,8,7)

@test is_connected(h)

scc = strongly_connected_components(h)
wcc = weakly_connected_components(h)

@test length(scc) == 3 && sort(scc[3]) == [1,2,5]
@test length(wcc) == 1 && length(wcc[1]) == nv(h)

function scc_ok(graph)
  """Check that all SCC really are strongly connected"""
  scc = strongly_connected_components(graph)
  scc_as_subgraphs = map(i -> induced_subgraph(graph, i), scc)
  return all(is_strongly_connected, scc_as_subgraphs)
end

# the two graphs below are isomorphic (exchange 2 <--> 4)
h = DiGraph(4);  add_edge!(h, 1, 4); add_edge!(h, 4, 2); add_edge!(h, 2, 3); add_edge!(h, 1, 3);
h2 = DiGraph(4); add_edge!(h2, 1, 2); add_edge!(h2, 2, 4); add_edge!(h2, 4, 3); add_edge!(h2, 1, 3);
@test scc_ok(h)
@test scc_ok(h2)

h = DiGraph(6)
add_edge!(h,1,3); add_edge!(h,3,4); add_edge!(h,4,2); add_edge!(h,2,1)
add_edge!(h,3,5); add_edge!(h,5,6); add_edge!(h,6,4)

scc = strongly_connected_components(h)

@test length(scc) == 1 && sort(scc[1]) == [1:6;]

# tests from Graphs.jl
h = DiGraph(4)
add_edge!(h,1,2); add_edge!(h,2,3); add_edge!(h,3,1); add_edge!(h,4,1)
scc = strongly_connected_components(h)
@test length(scc) == 2 && sort(scc[1]) == [1:3;] && sort(scc[2]) == [4]

h = DiGraph(12)
add_edge!(h,1,2); add_edge!(h,2,3); add_edge!(h,2,4); add_edge!(h,2,5);
add_edge!(h,3,6); add_edge!(h,4,5); add_edge!(h,4,7); add_edge!(h,5,2);
add_edge!(h,5,6); add_edge!(h,5,7); add_edge!(h,6,3); add_edge!(h,6,8);
add_edge!(h,7,8); add_edge!(h,7,10); add_edge!(h,8,7); add_edge!(h,9,7);
add_edge!(h,10,9); add_edge!(h,10,11); add_edge!(h,11,12); add_edge!(h,12,10)

scc = strongly_connected_components(h)
@test length(scc) == 4
@test sort(scc[1]) == [7,8,9,10,11,12]
@test sort(scc[2]) == [3, 6]
@test sort(scc[3]) == [2, 4, 5]
@test scc[4] == [1]

# Test examples with self-loops from
# Graph-Theoretic Analysis of Finite Markov Chains by J.P. Jarvis & D. R. Shier

# figure 1 example
fig1 = spzeros(5,5)
fig1[[3,4,9,10,11,13,18,19,22,24]] = [.5,.4,.1,1.,1.,.2,.3,.2,1.,.3]
fig1 = DiGraph(fig1)
scc_fig1 = Vector[[2,5],[1,3,4]]

# figure 2 example
fig2 = spzeros(5,5)
fig2[[3, 10, 11, 13, 14, 17, 18, 19, 22]] = 1
fig2 = DiGraph(fig2)

# figure 3 example
fig3 = spzeros(8,8)
fig3[[1,7,9,13,14,15,18,20,23,27,28,31,33,34,37,45,46,49,57,63,64]] = 1
fig3 = DiGraph(fig3)
scc_fig3 = Vector[[3,4],[2,5,6],[8],[1,7]]
fig3_cond = DiGraph(4);
add_edge!(fig3_cond,4,3); add_edge!(fig3_cond,2,1)
add_edge!(fig3_cond,4,1); add_edge!(fig3_cond,4,2)

# construct a n-number edge ring graph (period = n)
n = 10
n_ring_m = spdiagm(ones(n-1),1,n,n); n_ring_m[end,1] = 1
n_ring = DiGraph(n_ring_m)

n_ring_shortcut = copy(n_ring); add_edge!(n_ring_shortcut,1,4)

# figure 8 example
fig8 = spzeros(6,6)
fig8[[2,10,13,21,24,27,35]] = 1
fig8 = DiGraph(fig8)

@test Set(strongly_connected_components(fig1)) == Set(scc_fig1)
@test Set(strongly_connected_components(fig3)) == Set(scc_fig3)

@test period(n_ring) == n
@test period(n_ring_shortcut) == 2

@test condensation(fig3) == fig3_cond

@test attracting_components(fig1) == Vector[[2,5]]
@test attracting_components(fig3) == Vector[[3,4],[8]]

g10 = StarGraph(10)
@test neighborhood(g10, 1 , 0) == [1]
@test length(neighborhood(g10, 1, 1)) == 10
@test length(neighborhood(g10, 2, 1)) == 2
@test length(neighborhood(g10, 1, 2)) == 10
@test length(neighborhood(g10, 2, 2)) == 10

g10 = StarDiGraph(10)
@test neighborhood(g10, 1 , 0, dir=:out) == [1]
@test length(neighborhood(g10, 1, 1, dir=:out)) == 10
@test length(neighborhood(g10, 2, 1, dir=:out)) == 1
@test length(neighborhood(g10, 1, 2, dir=:out)) == 10
@test length(neighborhood(g10, 2, 2, dir=:out)) == 1
@test neighborhood(g10, 1 , 0, dir=:in) == [1]
@test length(neighborhood(g10, 1, 1, dir=:in)) == 1
@test length(neighborhood(g10, 2, 1, dir=:in)) == 2
@test length(neighborhood(g10, 1, 2, dir=:in)) == 1
@test length(neighborhood(g10, 2, 2, dir=:in)) == 2

g10 = StarGraph(10)
@test egonet(g10, 1, 0) == Graph(1,0)
@test egonet(g10, 1, 1) == g10

@test !isgraphical([1,1,1])
@test isgraphical([2,2,2])
@test isgraphical(fill(3,10))


text:
export viewSlice2D

"""
function jInv.Vis.getGridSlice

generates a slice of a 3D grid

Required Inputs:

    M::RegularMesh - mesh that describes geometry
    idx::Int       - slice index

Optional Inputs:

    view::Symbol   - possible choices :xy (default), :yz, :xz

"""
function getGridSlice(M::RegularMesh,idx,view=:xy)
    if view==:xy
        return getRegularMesh(M.domain[1:4], M.n[1:2])
    elseif view==:yz
        return getRegularMesh(M.domain[3:6], M.n[2:3])
    elseif view==:xz
        return getRegularMesh(M.domain[[1:2;5:6]], M.n[[1;3]])
    end
end

"""
function jInv.Vis.viewSlice2D

plots a slice of a 3D data set using viewImage2D and PyPlot.pcolormesh

Required Inputs:

    Img            - image data to be visualized (assumed to be cell-centered)
    M::RegularMesh - mesh that describes geometry
    idx::Int       - slice index

Optional Inputs:

    view::Symbol   - possible choices :xy (default), :yz, :xz
    addLabel::Bool - print axis labels (default=false)
    kwargs         - keyword arguments for PyPlot.pcolormesh

"""
function viewSlice2D(Img,M::RegularMesh,idx::Int;view=:xy,addLabel::Bool=false,kwargs...)

    if idx<1; error(BoundsError); end
    Img = reshape(Img,tuple(M.n...))
    if view==:xy
        if idx>M.n[3]; error(BoundsError); end;
        I2D = Img[:,:,idx]
        putLabels = () -> (xlabel("x"); ylabel("y"))
    elseif view==:yz
        if idx>M.n[1]; error(BoundsError); end;
        I2D = Img[idx,:,:]
        putLabels = () -> (xlabel("y"); ylabel("z"))
    elseif view==:xz
        if idx>M.n[2]; error(BoundsError); end;
        I2D = Img[:,idx,:]
        putLabels = () -> (xlabel("x"); ylabel("z"))
    end

    M2D = getGridSlice(M,idx,view)
    viewImage2D(I2D,M2D;kwargs...)

    addLabel && putLabels()
end


text:
using DelimitedFiles

function writemesh(p, t, filename, mode)

nd,np = size(p);
nve,ne = size(t);
tmp = [nd, np, nve, ne];

fileID = open(filename, "w");
if mode==0 # binary
    write(fileID,Float64.(tmp));
    write(fileID,Float64.(p));
    write(fileID,Float64.(t));
else # text
    writedlm(fileID, tmp');
    writedlm(fileID, p');
    writedlm(fileID, t');
end
close(fileID);

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
s = "I am a string"
if ismatch(r"string$", s)
    println("'$s' ends with 'string'")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode