text:
int_rules_5_3_6 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 5.3.6*Exponentials*of*inverse*tangent =#
    @apply_utils Antiderivative(ℯ ^ (~n * atan(~(a') * ~x)), ~x) => Antiderivative((1 - im * ~a * ~x) ^ ((im * ~n + 1) / 2) / ((1 + im * ~a * ~x) ^ ((im * ~n - 1) / 2) * sqrt(1 + (~a) ^ 2 * (~x) ^ 2)), ~x) <-- FreeQ(~a, ~x) && IntegerQ((im * ~n - 1) / 2)
    @apply_utils Antiderivative((~x) ^ ~(m') * ℯ ^ (~n * atan(~(a') * ~x)), ~x) => Antiderivative((~x) ^ ~m * ((1 - im * ~a * ~x) ^ ((im * ~n + 1) / 2) / ((1 + im * ~a * ~x) ^ ((im * ~n - 1) / 2) * sqrt(1 + (~a) ^ 2 * (~x) ^ 2))), ~x) <-- FreeQ([~a, ~m], ~x) && IntegerQ((im * ~n - 1) / 2)
    @apply_utils Antiderivative(ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => Antiderivative((1 - im * ~a * ~x) ^ ((im * ~n) / 2) / (1 + im * ~a * ~x) ^ ((im * ~n) / 2), ~x) <-- FreeQ([~a, ~n], ~x) && Not(IntegerQ((im * ~n - 1) / 2))
    @apply_utils Antiderivative((~x) ^ ~(m') * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => Antiderivative(((~x) ^ ~m * (1 - im * ~a * ~x) ^ ((im * ~n) / 2)) / (1 + im * ~a * ~x) ^ ((im * ~n) / 2), ~x) <-- FreeQ([~a, ~m, ~n], ~x) && Not(IntegerQ((im * ~n - 1) / 2))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x) ^ ~(p') * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative((~u * (1 + (~d * ~x) / ~c) ^ ~p * (1 - im * ~a * ~x) ^ ((im * ~n) / 2)) / (1 + im * ~a * ~x) ^ ((im * ~n) / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~a) ^ 2 * (~c) ^ 2 + (~d) ^ 2, 0) && (IntegerQ(~p) || GtQ(~c, 0)))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x) ^ ~(p') * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => Antiderivative((~u * (~c + ~d * ~x) ^ ~p * (1 - im * ~a * ~x) ^ ((im * ~n) / 2)) / (1 + im * ~a * ~x) ^ ((im * ~n) / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~a) ^ 2 * (~c) ^ 2 + (~d) ^ 2, 0) && Not(IntegerQ(~p) || GtQ(~c, 0)))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / ~x) ^ ~(p') * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => (~d) ^ ~p * Antiderivative((~u / (~x) ^ ~p) * (1 + (~c * ~x) / ~d) ^ ~p * ℯ ^ (~n * atan(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~c) ^ 2 + (~a) ^ 2 * (~d) ^ 2, 0) && IntegerQ(~p))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / ~x) ^ ~p * ℯ ^ (~n * atanh(~(a') * ~x)), ~x) => (-1) ^ (~n / 2) * (~c) ^ ~p * Antiderivative((~u * (1 + ~d / (~c * ~x)) ^ ~p * (1 - 1 / (im * ~a * ~x)) ^ ((im * ~n) / 2)) / (1 + 1 / (im * ~a * ~x)) ^ ((im * ~n) / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ((~c) ^ 2 + (~a) ^ 2 * (~d) ^ 2, 0) && (Not(IntegerQ(~p)) && (IntegerQ((im * ~n) / 2) && GtQ(~c, 0))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / ~x) ^ ~p * ℯ ^ (~n * atan(~(a') * ~x)), ~x) => Antiderivative((~u * (~c + ~d / ~x) ^ ~p * (1 - im * ~a * ~x) ^ ((im * ~n) / 2)) / (1 + im * ~a * ~x) ^ ((im * ~n) / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ((~c) ^ 2 + (~a) ^ 2 * (~d) ^ 2, 0) && (Not(IntegerQ(~p)) && (IntegerQ((im * ~n) / 2) && Not(GtQ(~c, 0)))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / ~x) ^ ~p * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => (((~x) ^ ~p * (~c + ~d / ~x) ^ ~p) / (1 + (~c * ~x) / ~d) ^ ~p) * Antiderivative((~u / (~x) ^ ~p) * (1 + (~c * ~x) / ~d) ^ ~p * ℯ ^ (~n * atan(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~c) ^ 2 + (~a) ^ 2 * (~d) ^ 2, 0) && Not(IntegerQ(~p)))
    @apply_utils Antiderivative(ℯ ^ (~(n') * atan(~(a') * ~x)) / (~c + ~(d') * (~x) ^ 2) ^ (3 / 2), ~x) => ((~n + ~a * ~x) * ℯ ^ (~n * atan(~a * ~x))) / (~a * ~c * ((~n) ^ 2 + 1) * sqrt(~c + ~d * (~x) ^ 2)) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && Not(IntegerQ(im * ~n)))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => ((~n - 2 * ~a * (~p + 1) * ~x) * (~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * atan(~a * ~x))) / (~a * ~c * ((~n) ^ 2 + 4 * (~p + 1) ^ 2)) + ((2 * (~p + 1) * (2 * ~p + 3)) / (~c * ((~n) ^ 2 + 4 * (~p + 1) ^ 2))) * Antiderivative((~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * atan(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (LtQ(~p, -1) && (Not(IntegerQ(im * ~n)) && (NeQ((~n) ^ 2 + 4 * (~p + 1) ^ 2, 0) && IntegerQ(2 * ~p)))))
    @apply_utils Antiderivative(ℯ ^ (~(n') * atan(~(a') * ~x)) / (~c + ~(d') * (~x) ^ 2), ~x) => ℯ ^ (~n * atan(~a * ~x)) / (~a * ~c * ~n) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && EqQ(~d, (~a) ^ 2 * ~c)
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~n * atan(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative((1 + (~a) ^ 2 * (~x) ^ 2) ^ (~p - (im * ~n) / 2) * (1 - im * ~a * ~x) ^ (im * ~n), ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (IntegerQ(~p) && (IntegerQ((im * ~n + 1) / 2) && Not(IntegerQ(~p - (im * ~n) / 2)))))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative((1 - im * ~a * ~x) ^ (~p + (im * ~n) / 2) * (1 + im * ~a * ~x) ^ (~p - (im * ~n) / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (IntegerQ(~p) || GtQ(~c, 0)))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atan(~(a') * ~x)), ~x) => (~c) ^ ((im * ~n) / 2) * Antiderivative((~c + ~d * (~x) ^ 2) ^ (~p - (im * ~n) / 2) * (1 - im * ~a * ~x) ^ (im * ~n), ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (Not(IntegerQ(~p) || GtQ(~c, 0)) && IGtQ((im * ~n) / 2, 0)))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atan(~(a') * ~x)), ~x) => (1 / (~c) ^ ((im * ~n) / 2)) * Antiderivative((~c + ~d * (~x) ^ 2) ^ (~p + (im * ~n) / 2) / (1 + im * ~a * ~x) ^ (im * ~n), ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (Not(IntegerQ(~p) || GtQ(~c, 0)) && ILtQ((im * ~n) / 2, 0)))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => (((~c) ^ IntPart(~p) * (~c + ~d * (~x) ^ 2) ^ FracPart(~p)) / (1 + (~a) ^ 2 * (~x) ^ 2) ^ FracPart(~p)) * Antiderivative((1 + (~a) ^ 2 * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atan(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && Not(IntegerQ(~p) || GtQ(~c, 0)))
    @apply_utils Antiderivative((~x * ℯ ^ (~(n') * atan(~(a') * ~x))) / (~c + ~(d') * (~x) ^ 2) ^ (3 / 2), ~x) => (-((1 - ~a * ~n * ~x)) * ℯ ^ (~n * atan(~a * ~x))) / (~d * ((~n) ^ 2 + 1) * sqrt(~c + ~d * (~x) ^ 2)) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && Not(IntegerQ(im * ~n)))
    @apply_utils Antiderivative(~x * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => ((~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * atan(~a * ~x))) / (2 * ~d * (~p + 1)) - ((~a * ~c * ~n) / (2 * ~d * (~p + 1))) * Antiderivative((~c + ~d * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atan(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (LtQ(~p, -1) && (Not(IntegerQ(im * ~n)) && IntegerQ(2 * ~p))))
    #= Antiderivative((~x)*((~c)+(~d')*(~x)^2)^(~p)*ℯ^((~n')*atan((~a')*(~x))),~x) :=  (2*(p+1)+a*n*x)*(c+d*x^2)^(p+1)*ℯ^(n*atan(a*x))/(a^2*c*(n^2+4*(p+1)^ 2)) - n*(2*p+3)/(a*c*(n^2+4*(p+1)^2))*Antiderivative((c+d*x^2)^(p+1)*ℯ^(n*atan(a*x) ),x) <-- FreeQ([a,c,d,n],x) && EqQ(d,a^2*c) && LtQ(p,-1) &&  NeQ(n^2+4*(p+1)^2,0) && Not(IntegerQ(im*n)) =#
    @apply_utils Antiderivative((~x) ^ 2 * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => (-((1 - ~a * ~n * ~x)) * (~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * atan(~a * ~x))) / (~a * ~d * ~n * ((~n) ^ 2 + 1)) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (EqQ((~n) ^ 2 - 2 * (~p + 1), 0) && Not(IntegerQ(im * ~n))))
    @apply_utils Antiderivative((~x) ^ 2 * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => (-((~n - 2 * (~p + 1) * ~a * ~x)) * (~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * atan(~a * ~x))) / (~a * ~d * ((~n) ^ 2 + 4 * (~p + 1) ^ 2)) + (((~n) ^ 2 - 2 * (~p + 1)) / (~d * ((~n) ^ 2 + 4 * (~p + 1) ^ 2))) * Antiderivative((~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * atan(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (LtQ(~p, -1) && (Not(IntegerQ(im * ~n)) && (NeQ((~n) ^ 2 + 4 * (~p + 1) ^ 2, 0) && IntegerQ(2 * ~p)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~n * atan(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative((~x) ^ ~m * (1 + (~a) ^ 2 * (~x) ^ 2) ^ (~p - (im * ~n) / 2) * (1 - im * ~a * ~x) ^ (im * ~n), ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~p], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && ((IntegerQ(~p) || GtQ(~c, 0)) && (IntegerQ((im * ~n + 1) / 2) && Not(IntegerQ(~p - (im * ~n) / 2)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative((~x) ^ ~m * (1 - im * ~a * ~x) ^ (~p + (im * ~n) / 2) * (1 + im * ~a * ~x) ^ (~p - (im * ~n) / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (IntegerQ(~p) || GtQ(~c, 0)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atan(~(a') * ~x)), ~x) => (~c) ^ ((im * ~n) / 2) * Antiderivative((~x) ^ ~m * (~c + ~d * (~x) ^ 2) ^ (~p - (im * ~n) / 2) * (1 - im * ~a * ~x) ^ (im * ~n), ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~p], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (Not(IntegerQ(~p) || GtQ(~c, 0)) && IGtQ((im * ~n) / 2, 0)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atan(~(a') * ~x)), ~x) => (1 / (~c) ^ ((im * ~n) / 2)) * Antiderivative(((~x) ^ ~m * (~c + ~d * (~x) ^ 2) ^ (~p + (im * ~n) / 2)) / (1 + im * ~a * ~x) ^ (im * ~n), ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~p], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (Not(IntegerQ(~p) || GtQ(~c, 0)) && ILtQ((im * ~n) / 2, 0)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => (((~c) ^ IntPart(~p) * (~c + ~d * (~x) ^ 2) ^ FracPart(~p)) / (1 + (~a) ^ 2 * (~x) ^ 2) ^ FracPart(~p)) * Antiderivative((~x) ^ ~m * (1 + (~a) ^ 2 * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atan(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && Not(IntegerQ(~p) || GtQ(~c, 0)))
    @apply_utils Antiderivative(~u * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative(~u * (1 - im * ~a * ~x) ^ (~p + (im * ~n) / 2) * (1 + im * ~a * ~x) ^ (~p - (im * ~n) / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (IntegerQ(~p) || GtQ(~c, 0)))
    @apply_utils Antiderivative(~u * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~n * atan(~(a') * ~x)), ~x) => (((~c) ^ IntPart(~p) * (~c + ~d * (~x) ^ 2) ^ FracPart(~p)) / ((1 - im * ~a * ~x) ^ FracPart(~p) * (1 + im * ~a * ~x) ^ FracPart(~p))) * Antiderivative(~u * (1 - im * ~a * ~x) ^ (~p + (im * ~n) / 2) * (1 + im * ~a * ~x) ^ (~p - (im * ~n) / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && ((IntegerQ(~p) || GtQ(~c, 0)) && IntegerQ((im * ~n) / 2)))
    @apply_utils Antiderivative(~u * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => (((~c) ^ IntPart(~p) * (~c + ~d * (~x) ^ 2) ^ FracPart(~p)) / (1 + (~a) ^ 2 * (~x) ^ 2) ^ FracPart(~p)) * Antiderivative(~u * (1 + (~a) ^ 2 * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atan(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (Not(IntegerQ(~p) || GtQ(~c, 0)) && Not(IntegerQ((im * ~n) / 2))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => (~d) ^ ~p * Antiderivative((~u / (~x) ^ (2 * ~p)) * (1 + (~a) ^ 2 * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atan(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~c - (~a) ^ 2 * ~d, 0) && IntegerQ(~p))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / (~x) ^ 2) ^ ~p * ℯ ^ (~n * atan(~(a') * ~x)), ~x) => (~c) ^ ~p * Antiderivative(~u * (1 - im / (~a * ~x)) ^ ~p * (1 + im / (~a * ~x)) ^ ~p * ℯ ^ (~n * atan(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~p], ~x) && (EqQ(~c - (~a) ^ 2 * ~d, 0) && (Not(IntegerQ(~p)) && (IntegerQ((im * ~n) / 2) && GtQ(~c, 0))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / (~x) ^ 2) ^ ~p * ℯ ^ (~n * atan(~(a') * ~x)), ~x) => (((~x) ^ (2 * ~p) * (~c + ~d / (~x) ^ 2) ^ ~p) / ((1 - im * ~a * ~x) ^ ~p * (1 + im * ~a * ~x) ^ ~p)) * Antiderivative((~u / (~x) ^ (2 * ~p)) * (1 - im * ~a * ~x) ^ ~p * (1 + im * ~a * ~x) ^ ~p * ℯ ^ (~n * atan(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~c - (~a) ^ 2 * ~d, 0) && (Not(IntegerQ(~p)) && (IntegerQ((im * ~n) / 2) && Not(GtQ(~c, 0)))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * atan(~(a') * ~x)), ~x) => (((~x) ^ (2 * ~p) * (~c + ~d / (~x) ^ 2) ^ ~p) / (1 + (~a) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative((~u / (~x) ^ (2 * ~p)) * (1 + (~a) ^ 2 * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atan(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~c - (~a) ^ 2 * ~d, 0) && (Not(IntegerQ(~p)) && Not(IntegerQ((im * ~n) / 2))))
    @apply_utils Antiderivative(ℯ ^ (~(n') * atan(~(c') * (~a + ~(b') * ~x))), ~x) => Antiderivative(((1 - im * ~a * ~c) - im * ~b * ~c * ~x) ^ ((im * ~n) / 2) / (1 + im * ~a * ~c + im * ~b * ~c * ~x) ^ ((im * ~n) / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~n * atan(~(c') * (~a + ~(b') * ~x))), ~x) => (4 / (im ^ ~m * ~n * (~b) ^ (~m + 1) * (~c) ^ (~m + 1))) * Subst(Antiderivative(((~x) ^ (2 / (im * ~n)) * ((1 - im * ~a * ~c) - (1 + im * ~a * ~c) * (~x) ^ (2 / (im * ~n))) ^ ~m) / (1 + (~x) ^ (2 / (im * ~n))) ^ (~m + 2), ~x), ~x, (1 - im * ~c * (~a + ~b * ~x)) ^ ((im * ~n) / 2) / (1 + im * ~c * (~a + ~b * ~x)) ^ ((im * ~n) / 2)) <-- FreeQ([~a, ~b, ~c], ~x) && (ILtQ(~m, 0) && LtQ(-1, im * ~n, 1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ℯ ^ (~(n') * atan(~(c') * (~a + ~(b') * ~x))), ~x) => Antiderivative(((~d + ~e * ~x) ^ ~m * ((1 - im * ~a * ~c) - im * ~b * ~c * ~x) ^ ((im * ~n) / 2)) / (1 + im * ~a * ~c + im * ~b * ~c * ~x) ^ ((im * ~n) / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x + ~(e') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atan(~a + ~(b') * ~x)), ~x) => (~c / (1 + (~a) ^ 2)) ^ ~p * Antiderivative(~u * ((1 - im * ~a) - im * ~b * ~x) ^ (~p + (im * ~n) / 2) * (1 + im * ~a + im * ~b * ~x) ^ (~p - (im * ~n) / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (EqQ(~b * ~d, 2 * ~a * ~e) && (EqQ((~b) ^ 2 * ~c - ~e * (1 + (~a) ^ 2), 0) && (IntegerQ(~p) || GtQ(~c / (1 + (~a) ^ 2), 0))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x + ~(e') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * atan(~a + ~(b') * ~x)), ~x) => ((~c + ~d * ~x + ~e * (~x) ^ 2) ^ ~p / (1 + (~a) ^ 2 + 2 * ~a * ~b * ~x + (~b) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative(~u * (1 + (~a) ^ 2 + 2 * ~a * ~b * ~x + (~b) ^ 2 * (~x) ^ 2) ^ ~p * ℯ ^ (~n * atan(~a * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (EqQ(~b * ~d, 2 * ~a * ~e) && (EqQ((~b) ^ 2 * ~c - ~e * (1 + (~a) ^ 2), 0) && Not(IntegerQ(~p) || GtQ(~c / (1 + (~a) ^ 2), 0))))
    @apply_utils Antiderivative(~(u') * ℯ ^ (~(n') * atan(~(c') / (~(a') + ~(b') * ~x))), ~x) => Antiderivative(~u * ℯ ^ (~n * acot(~a / ~c + (~b * ~x) / ~c)), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
    @apply_utils Antiderivative(~(u') * ℯ ^ (~n * acot(~(a') * ~x)), ~x) => (-1) ^ ((im * ~n) / 2) * Antiderivative(~u * ℯ ^ (-(~n) * atan(~a * ~x)), ~x) <-- FreeQ(~a, ~x) && IntegerQ((im * ~n) / 2)
    @apply_utils Antiderivative(ℯ ^ (~n * acot(~(a') * ~x)), ~x) => -(Subst(Antiderivative((1 - (im * ~x) / ~a) ^ ((im * ~n + 1) / 2) / ((~x) ^ 2 * (1 + (im * ~x) / ~a) ^ ((im * ~n - 1) / 2) * sqrt(1 + (~x) ^ 2 / (~a) ^ 2)), ~x), ~x, 1 / ~x)) <-- FreeQ(~a, ~x) && IntegerQ((im * ~n - 1) / 2)
    @apply_utils Antiderivative((~x) ^ ~(m') * ℯ ^ (~n * acot(~(a') * ~x)), ~x) => -(Subst(Antiderivative((1 - (im * ~x) / ~a) ^ ((im * ~n + 1) / 2) / ((~x) ^ (~m + 2) * (1 + (im * ~x) / ~a) ^ ((im * ~n - 1) / 2) * sqrt(1 + (~x) ^ 2 / (~a) ^ 2)), ~x), ~x, 1 / ~x)) <-- FreeQ(~a, ~x) && (IntegerQ((im * ~n - 1) / 2) && IntegerQ(~m))
    @apply_utils Antiderivative(ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => -(Subst(Antiderivative((1 - (im * ~x) / ~a) ^ ((im * ~n) / 2) / ((~x) ^ 2 * (1 + (im * ~x) / ~a) ^ ((im * ~n) / 2)), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~n], ~x) && Not(IntegerQ(im * ~n))
    @apply_utils Antiderivative((~x) ^ ~(m') * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => -(Subst(Antiderivative((1 - (im * ~x) / ~a) ^ (~n / 2) / ((~x) ^ (~m + 2) * (1 + (im * ~x) / ~a) ^ (~n / 2)), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~n], ~x) && (Not(IntegerQ(im * ~n)) && IntegerQ(~m))
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~n * acot(~(a') * ~x)), ~x) => -((~x) ^ ~m) * (1 / ~x) ^ ~m * Subst(Antiderivative((1 - (im * ~x) / ~a) ^ ((im * ~n + 1) / 2) / ((~x) ^ (~m + 2) * (1 + (im * ~x) / ~a) ^ ((im * ~n - 1) / 2) * sqrt(1 + (~x) ^ 2 / (~a) ^ 2)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~m], ~x) && (IntegerQ((im * ~n - 1) / 2) && Not(IntegerQ(~m)))
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => -(Subst(Antiderivative((1 - (im * ~x) / ~a) ^ (~n / 2) / ((~x) ^ (~m + 2) * (1 + (im * ~x) / ~a) ^ (~n / 2)), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~m, ~n], ~x) && (Not(IntegerQ((im * ~n) / 2)) && Not(IntegerQ(~m)))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x) ^ ~(p') * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => (~d) ^ ~p * Antiderivative(~u * (~x) ^ ~p * (1 + ~c / (~d * ~x)) ^ ~p * ℯ ^ (~n * acot(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ((~a) ^ 2 * (~c) ^ 2 + (~d) ^ 2, 0) && (Not(IntegerQ((im * ~n) / 2)) && IntegerQ(~p)))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x) ^ ~p * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => ((~c + ~d * ~x) ^ ~p / ((~x) ^ ~p * (1 + ~c / (~d * ~x)) ^ ~p)) * Antiderivative(~u * (~x) ^ ~p * (1 + ~c / (~d * ~x)) ^ ~p * ℯ ^ (~n * acot(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~a) ^ 2 * (~c) ^ 2 + (~d) ^ 2, 0) && (Not(IntegerQ((im * ~n) / 2)) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~c + ~(d') / ~x) ^ ~(p') * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => -((~c) ^ ~p) * Subst(Antiderivative(((1 + (~d * ~x) / ~c) ^ ~p * (1 - (im * ~x) / ~a) ^ ((im * ~n) / 2)) / ((~x) ^ 2 * (1 + (im * ~x) / ~a) ^ ((im * ~n) / 2)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~c) ^ 2 + (~a) ^ 2 * (~d) ^ 2, 0) && (Not(IntegerQ((im * ~n) / 2)) && (IntegerQ(~p) || GtQ(~c, 0))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') / ~x) ^ ~(p') * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => -((~c) ^ ~p) * Subst(Antiderivative(((1 + (~d * ~x) / ~c) ^ ~p * (1 - (im * ~x) / ~a) ^ ((im * ~n) / 2)) / ((~x) ^ (~m + 2) * (1 + (im * ~x) / ~a) ^ ((im * ~n) / 2)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ((~c) ^ 2 + (~a) ^ 2 * (~d) ^ 2, 0) && (Not(IntegerQ((im * ~n) / 2)) && ((IntegerQ(~p) || GtQ(~c, 0)) && IntegerQ(~m))))
    @apply_utils Antiderivative((~c + ~(d') / ~x) ^ ~p * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => ((~c + ~d / ~x) ^ ~p / (1 + ~d / (~c * ~x)) ^ ~p) * Antiderivative((1 + ~d / (~c * ~x)) ^ ~p * ℯ ^ (~n * acot(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~c) ^ 2 + (~a) ^ 2 * (~d) ^ 2, 0) && (Not(IntegerQ((im * ~n) / 2)) && Not(IntegerQ(~p) || GtQ(~c, 0))))
    @apply_utils Antiderivative((~x) ^ ~m * (~c + ~(d') / ~x) ^ ~(p') * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => -((~c) ^ ~p) * (~x) ^ ~m * (1 / ~x) ^ ~m * Subst(Antiderivative(((1 + (~d * ~x) / ~c) ^ ~p * (1 - (im * ~x) / ~a) ^ ((im * ~n) / 2)) / ((~x) ^ (~m + 2) * (1 + (im * ~x) / ~a) ^ ((im * ~n) / 2)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ((~c) ^ 2 + (~a) ^ 2 * (~d) ^ 2, 0) && (Not(IntegerQ((im * ~n) / 2)) && ((IntegerQ(~p) || GtQ(~c, 0)) && Not(IntegerQ(~m)))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / ~x) ^ ~p * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => ((~c + ~d / ~x) ^ ~p / (1 + ~d / (~c * ~x)) ^ ~p) * Antiderivative(~u * (1 + ~d / (~c * ~x)) ^ ~p * ℯ ^ (~n * acot(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ((~c) ^ 2 + (~a) ^ 2 * (~d) ^ 2, 0) && (Not(IntegerQ((im * ~n) / 2)) && Not(IntegerQ(~p) || GtQ(~c, 0))))
    @apply_utils Antiderivative(ℯ ^ (~(n') * acot(~(a') * ~x)) / (~c + ~(d') * (~x) ^ 2), ~x) => -(ℯ ^ (~n * acot(~a * ~x))) / (~a * ~c * ~n) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && EqQ(~d, (~a) ^ 2 * ~c)
    @apply_utils Antiderivative(ℯ ^ (~(n') * acot(~(a') * ~x)) / (~c + ~(d') * (~x) ^ 2) ^ (3 / 2), ~x) => (-((~n - ~a * ~x)) * ℯ ^ (~n * acot(~a * ~x))) / (~a * ~c * ((~n) ^ 2 + 1) * sqrt(~c + ~d * (~x) ^ 2)) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && Not(IntegerQ((im * ~n - 1) / 2)))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => (-((~n + 2 * ~a * (~p + 1) * ~x)) * (~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * acot(~a * ~x))) / (~a * ~c * ((~n) ^ 2 + 4 * (~p + 1) ^ 2)) + ((2 * (~p + 1) * (2 * ~p + 3)) / (~c * ((~n) ^ 2 + 4 * (~p + 1) ^ 2))) * Antiderivative((~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * acot(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (LtQ(~p, -1) && (NeQ(~p, -3 / 2) && (NeQ((~n) ^ 2 + 4 * (~p + 1) ^ 2, 0) && (Not(IntegerQ(~p) && IntegerQ((im * ~n) / 2)) && Not(Not(IntegerQ(~p)) && IntegerQ((im * ~n - 1) / 2)))))))
    @apply_utils Antiderivative((~x * ℯ ^ (~(n') * acot(~(a') * ~x))) / (~c + ~(d') * (~x) ^ 2) ^ (3 / 2), ~x) => (-((1 + ~a * ~n * ~x)) * ℯ ^ (~n * acot(~a * ~x))) / ((~a) ^ 2 * ~c * ((~n) ^ 2 + 1) * sqrt(~c + ~d * (~x) ^ 2)) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && Not(IntegerQ((im * ~n - 1) / 2)))
    @apply_utils Antiderivative(~x * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => ((2 * (~p + 1) - ~a * ~n * ~x) * (~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * acot(~a * ~x))) / ((~a) ^ 2 * ~c * ((~n) ^ 2 + 4 * (~p + 1) ^ 2)) + ((~n * (2 * ~p + 3)) / (~a * ~c * ((~n) ^ 2 + 4 * (~p + 1) ^ 2))) * Antiderivative((~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * acot(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (LeQ(~p, -1) && (NeQ(~p, -3 / 2) && (NeQ((~n) ^ 2 + 4 * (~p + 1) ^ 2, 0) && (Not(IntegerQ(~p) && IntegerQ((im * ~n) / 2)) && Not(Not(IntegerQ(~p)) && IntegerQ((im * ~n - 1) / 2)))))))
    @apply_utils Antiderivative((~x) ^ 2 * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => ((~n + 2 * (~p + 1) * ~a * ~x) * (~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * acot(~a * ~x))) / ((~a) ^ 3 * ~c * (~n) ^ 2 * ((~n) ^ 2 + 1)) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (EqQ((~n) ^ 2 - 2 * (~p + 1), 0) && NeQ((~n) ^ 2 + 1, 0)))
    @apply_utils Antiderivative((~x) ^ 2 * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => ((~n + 2 * (~p + 1) * ~a * ~x) * (~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * acot(~a * ~x))) / ((~a) ^ 3 * ~c * ((~n) ^ 2 + 4 * (~p + 1) ^ 2)) + (((~n) ^ 2 - 2 * (~p + 1)) / ((~a) ^ 2 * ~c * ((~n) ^ 2 + 4 * (~p + 1) ^ 2))) * Antiderivative((~c + ~d * (~x) ^ 2) ^ (~p + 1) * ℯ ^ (~n * acot(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (LeQ(~p, -1) && (NeQ((~n) ^ 2 - 2 * (~p + 1), 0) && (NeQ((~n) ^ 2 + 4 * (~p + 1) ^ 2, 0) && (Not(IntegerQ(~p) && IntegerQ((im * ~n) / 2)) && Not(Not(IntegerQ(~p)) && IntegerQ((im * ~n - 1) / 2)))))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => (-((~c) ^ ~p) / (~a) ^ (~m + 1)) * Subst(Antiderivative((ℯ ^ (~n * ~x) * cot(~x) ^ (~m + 2 * (~p + 1))) / cos(~x) ^ (2 * (~p + 1)), ~x), ~x, acot(~a * ~x)) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (IntegerQ(~m) && (LeQ(3, ~m, -2 * (~p + 1)) && IntegerQ(~p))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => (~d) ^ ~p * Antiderivative(~u * (~x) ^ (2 * ~p) * (1 + 1 / ((~a) ^ 2 * (~x) ^ 2)) ^ ~p * ℯ ^ (~n * acot(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (Not(IntegerQ((im * ~n) / 2)) && IntegerQ(~p)))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => ((~c + ~d * (~x) ^ 2) ^ ~p / ((~x) ^ (2 * ~p) * (1 + 1 / ((~a) ^ 2 * (~x) ^ 2)) ^ ~p)) * Antiderivative(~u * (~x) ^ (2 * ~p) * (1 + 1 / ((~a) ^ 2 * (~x) ^ 2)) ^ ~p * ℯ ^ (~n * acot(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~d, (~a) ^ 2 * ~c) && (Not(IntegerQ((im * ~n) / 2)) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => ((~c) ^ ~p / (im * ~a) ^ (2 * ~p)) * Antiderivative((~u / (~x) ^ (2 * ~p)) * (-1 + im * ~a * ~x) ^ (~p - (im * ~n) / 2) * (1 + im * ~a * ~x) ^ (~p + (im * ~n) / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~c, (~a) ^ 2 * ~d) && (Not(IntegerQ((im * ~n) / 2)) && ((IntegerQ(~p) || GtQ(~c, 0)) && IntegersQ(2 * ~p, ~p + (im * ~n) / 2))))
    @apply_utils Antiderivative((~c + ~(d') / (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => -((~c) ^ ~p) * Subst(Antiderivative(((1 - (im * ~x) / ~a) ^ (~p + (im * ~n) / 2) * (1 + (im * ~x) / ~a) ^ (~p - (im * ~n) / 2)) / (~x) ^ 2, ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~c, (~a) ^ 2 * ~d) && (Not(IntegerQ((im * ~n) / 2)) && ((IntegerQ(~p) || GtQ(~c, 0)) && Not(IntegerQ(2 * ~p) && IntegerQ(~p + (im * ~n) / 2)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~c + ~(d') / (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => -((~c) ^ ~p) * Subst(Antiderivative(((1 - (im * ~x) / ~a) ^ (~p + (im * ~n) / 2) * (1 + (im * ~x) / ~a) ^ (~p - (im * ~n) / 2)) / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~c, (~a) ^ 2 * ~d) && (Not(IntegerQ((im * ~n) / 2)) && ((IntegerQ(~p) || GtQ(~c, 0)) && (Not(IntegerQ(2 * ~p) && IntegerQ(~p + (im * ~n) / 2)) && IntegerQ(~m)))))
    @apply_utils Antiderivative((~x) ^ ~m * (~c + ~(d') / (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => -((~c) ^ ~p) * (~x) ^ ~m * (1 / ~x) ^ ~m * Subst(Antiderivative(((1 - (im * ~x) / ~a) ^ (~p + (im * ~n) / 2) * (1 + (im * ~x) / ~a) ^ (~p - (im * ~n) / 2)) / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ(~c, (~a) ^ 2 * ~d) && (Not(IntegerQ((im * ~n) / 2)) && ((IntegerQ(~p) || GtQ(~c, 0)) && (Not(IntegerQ(2 * ~p) && IntegerQ(~p + (im * ~n) / 2)) && Not(IntegerQ(~m))))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') / (~x) ^ 2) ^ ~p * ℯ ^ (~(n') * acot(~(a') * ~x)), ~x) => ((~c + ~d / (~x) ^ 2) ^ ~p / (1 + 1 / ((~a) ^ 2 * (~x) ^ 2)) ^ ~p) * Antiderivative(~u * (1 + 1 / ((~a) ^ 2 * (~x) ^ 2)) ^ ~p * ℯ ^ (~n * acot(~a * ~x)), ~x) <-- FreeQ([~a, ~c, ~d, ~n, ~p], ~x) && (EqQ(~c, (~a) ^ 2 * ~d) && (Not(IntegerQ((im * ~n) / 2)) && Not(IntegerQ(~p) || GtQ(~c, 0))))
    @apply_utils Antiderivative(~(u') * ℯ ^ (~n * acot(~(c') * (~a + ~(b') * ~x))), ~x) => (-1) ^ ((im * ~n) / 2) * Antiderivative(~u * ℯ ^ (-(~n) * atan(~c * (~a + ~b * ~x))), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && IntegerQ((im * ~n) / 2)
    @apply_utils Antiderivative(ℯ ^ (~(n') * acot(~(c') * (~a + ~(b') * ~x))), ~x) => (((im * ~c * (~a + ~b * ~x)) ^ ((im * ~n) / 2) * (1 + 1 / (im * ~c * (~a + ~b * ~x))) ^ ((im * ~n) / 2)) / (1 + im * ~a * ~c + im * ~b * ~c * ~x) ^ ((im * ~n) / 2)) * Antiderivative((1 + im * ~a * ~c + im * ~b * ~c * ~x) ^ ((im * ~n) / 2) / (-1 + im * ~a * ~c + im * ~b * ~c * ~x) ^ ((im * ~n) / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && Not(IntegerQ((im * ~n) / 2))
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~n * acoth(~(c') * (~a + ~(b') * ~x))), ~x) => (4 / (im ^ ~m * ~n * (~b) ^ (~m + 1) * (~c) ^ (~m + 1))) * Subst(Antiderivative(((~x) ^ (2 / (im * ~n)) * (1 + im * ~a * ~c + (1 - im * ~a * ~c) * (~x) ^ (2 / (im * ~n))) ^ ~m) / (-1 + (~x) ^ (2 / (im * ~n))) ^ (~m + 2), ~x), ~x, (1 + 1 / (im * ~c * (~a + ~b * ~x))) ^ ((im * ~n) / 2) / (1 - 1 / (im * ~c * (~a + ~b * ~x))) ^ ((im * ~n) / 2)) <-- FreeQ([~a, ~b, ~c], ~x) && (ILtQ(~m, 0) && LtQ(-1, im * ~n, 1))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * ℯ ^ (~(n') * acoth(~(c') * (~a + ~(b') * ~x))), ~x) => (((im * ~c * (~a + ~b * ~x)) ^ ((im * ~n) / 2) * (1 + 1 / (im * ~c * (~a + ~b * ~x))) ^ ((im * ~n) / 2)) / (1 + im * ~a * ~c + im * ~b * ~c * ~x) ^ ((im * ~n) / 2)) * Antiderivative(((~d + ~e * ~x) ^ ~m * (1 + im * ~a * ~c + im * ~b * ~c * ~x) ^ ((im * ~n) / 2)) / (-1 + im * ~a * ~c + im * ~b * ~c * ~x) ^ ((im * ~n) / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && Not(IntegerQ((im * ~n) / 2))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x + ~(e') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acot(~a + ~(b') * ~x)), ~x) => (~c / (1 + (~a) ^ 2)) ^ ~p * ((im * ~a + im * ~b * ~x) / (1 + im * ~a + im * ~b * ~x)) ^ ((im * ~n) / 2) * ((1 + im * ~a + im * ~b * ~x) / (im * ~a + im * ~b * ~x)) ^ ((im * ~n) / 2) * (((1 - im * ~a) - im * ~b * ~x) ^ ((im * ~n) / 2) / (-1 + im * ~a + im * ~b * ~x) ^ ((im * ~n) / 2)) * Antiderivative(~u * ((1 - im * ~a) - im * ~b * ~x) ^ (~p - (im * ~n) / 2) * (1 + im * ~a + im * ~b * ~x) ^ (~p + (im * ~n) / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (Not(IntegerQ((im * ~n) / 2)) && (EqQ(~b * ~d - 2 * ~a * ~e, 0) && (EqQ((~b) ^ 2 * ~c - ~e * (1 + (~a) ^ 2), 0) && (IntegerQ(~p) || GtQ(~c / (1 + (~a) ^ 2), 0)))))
    @apply_utils Antiderivative(~(u') * (~c + ~(d') * ~x + ~(e') * (~x) ^ 2) ^ ~(p') * ℯ ^ (~(n') * acot(~a + ~(b') * ~x)), ~x) => ((~c + ~d * ~x + ~e * (~x) ^ 2) ^ ~p / (1 + (~a) ^ 2 + 2 * ~a * ~b * ~x + (~b) ^ 2 * (~x) ^ 2) ^ ~p) * Antiderivative(~u * (1 + (~a) ^ 2 + 2 * ~a * ~b * ~x + (~b) ^ 2 * (~x) ^ 2) ^ ~p * ℯ ^ (~n * acot(~a * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (Not(IntegerQ((im * ~n) / 2)) && (EqQ(~b * ~d - 2 * ~a * ~e, 0) && (EqQ((~b) ^ 2 * ~c - ~e * (1 + (~a) ^ 2), 0) && Not(IntegerQ(~p) || GtQ(~c / (1 + (~a) ^ 2), 0)))))
    @apply_utils Antiderivative(~(u') * ℯ ^ (~(n') * acot(~(c') / (~(a') + ~(b') * ~x))), ~x) => Antiderivative(~u * ℯ ^ (~n * atan(~a / ~c + (~b * ~x) / ~c)), ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
end

max_stars_repo_path: src/rules/5 Inverse trig functions/5.3 Inverse tangent/5.3.6 Exponentials of inverse tangent.jl


text:
@doc raw"""
    Centered()

Used in conjunction with `SparseVariationalApproximation`.
States that the `q` field of [`SparseVariationalApproximation`](@ref) is to be interpreted
directly as the approximate posterior over the pseudo-points.

This is also known as the "unwhitened" parametrization [1].

See also [`NonCentered`](@ref).

[1] - https://en.wikipedia.org/wiki/Whitening_transformation
"""
struct Centered end

@doc raw"""
    NonCentered()

Used in conjunction with `SparseVariationalApproximation`.
States that the `q` field of [`SparseVariationalApproximation`](@ref) is to be interpreted
as the approximate posterior over `cholesky(cov(u)).L \ (u - mean(u))`, where `u` are the
pseudo-points.

This is also known as the "whitened" parametrization [1].

See also [`Centered`](@ref).

[1] - https://en.wikipedia.org/wiki/Whitening_transformation
"""
struct NonCentered end

struct SparseVariationalApproximation{Parametrization,Tfz<:FiniteGP,Tq<:AbstractMvNormal}
    fz::Tfz
    q::Tq
end

@doc raw"""
    SparseVariationalApproximation(
        ::Parametrization, fz::FiniteGP, q::AbstractMvNormal
    ) where {Parametrization}

Produce a `SparseVariationalApproximation{Parametrization}`, which packages the prior over
the pseudo-points, `fz`, and the approximate posterior at the pseudo-points, `q`, together
into a single object.

The `Parametrization` determines the precise manner in which `q` and `fz` are interpreted.
Existing parametrizations include [`Centered`](@ref) and [`NonCentered`](@ref).
"""
function SparseVariationalApproximation(
    ::Parametrization, fz::Tfz, q::Tq
) where {Parametrization,Tfz<:FiniteGP,Tq<:AbstractMvNormal}
    return SparseVariationalApproximation{Parametrization,Tfz,Tq}(fz, q)
end

"""
    SparseVariationalApproximation(fz::FiniteGP, q::AbstractMvNormal)

Packages the prior over the pseudo-points `fz`, and the approximate posterior at the
pseudo-points, which is `mean(fz) + cholesky(cov(fz)).U' * ε`, `ε ∼ q`.

Shorthand for
```julia
SparseVariationalApproximation(NonCentered(), fz, q)
```
"""
function SparseVariationalApproximation(fz::FiniteGP, q::AbstractMvNormal)
    return SparseVariationalApproximation(NonCentered(), fz, q)
end

@doc raw"""
    posterior(sva::SparseVariationalApproximation{Centered})

Compute the approximate posterior [1] over the process `f =
sva.fz.f`, given inducing inputs `z = sva.fz.x` and a variational
distribution over inducing points `sva.q` (which represents ``q(u)``
where `u = f(z)`). The approximate posterior at test points ``x^*``
where ``f^* = f(x^*)`` is then given by:

```math
q(f^*) = \int p(f | u) q(u) du
```
which can be found in closed form.

[1] - Hensman, James, Alexander Matthews, and Zoubin Ghahramani. "Scalable
variational Gaussian process classification." Artificial Intelligence and
Statistics. PMLR, 2015.
"""
function AbstractGPs.posterior(sva::SparseVariationalApproximation{Centered})
    q, fz = sva.q, sva.fz
    m, S = mean(q), _chol_cov(q)
    Kuu = _chol_cov(fz)
    B = Kuu.L \ S.L
    α = Kuu \ (m - mean(fz))
    data = (S=S, m=m, Kuu=Kuu, B=B, α=α)
    return ApproxPosteriorGP(sva, fz.f, data)
end

function AbstractGPs.posterior(
    sva::SparseVariationalApproximation, fx::FiniteGP, ::AbstractVector{<:Real}
)
    @assert sva.fz.f === fx.f
    return posterior(sva)
end

#
# Various methods implementing the Internal AbstractGPs API.
# See AbstractGPs.jl API docs for more info.
#

function Statistics.mean(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation{Centered}}, x::AbstractVector
)
    return mean(f.prior, x) + cov(f.prior, x, inducing_points(f)) * f.data.α
end

function Statistics.cov(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation{Centered}}, x::AbstractVector
)
    Cux = cov(f.prior, inducing_points(f), x)
    D = f.data.Kuu.L \ Cux
    return cov(f.prior, x) - At_A(D) + At_A(f.data.B' * D)
end

function Statistics.var(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation{Centered}}, x::AbstractVector
)
    Cux = cov(f.prior, inducing_points(f), x)
    D = f.data.Kuu.L \ Cux
    return var(f.prior, x) - diag_At_A(D) + diag_At_A(f.data.B' * D)
end

function Statistics.cov(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation{Centered}},
    x::AbstractVector,
    y::AbstractVector,
)
    B = f.data.B
    Cxu = cov(f.prior, x, inducing_points(f))
    Cuy = cov(f.prior, inducing_points(f), y)
    D = f.data.Kuu.L \ Cuy
    E = Cxu / f.data.Kuu.L'
    return cov(f.prior, x, y) - (E * D) + (E * B * B' * D)
end

function StatsBase.mean_and_cov(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation{Centered}}, x::AbstractVector
)
    Cux = cov(f.prior, inducing_points(f), x)
    D = f.data.Kuu.L \ Cux
    μ = Cux' * f.data.α
    Σ = cov(f.prior, x) - At_A(D) + At_A(f.data.B' * D)
    return μ, Σ
end

function StatsBase.mean_and_var(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation{Centered}}, x::AbstractVector
)
    Cux = cov(f.prior, inducing_points(f), x)
    D = f.data.Kuu.L \ Cux
    μ = Cux' * f.data.α
    Σ_diag = var(f.prior, x) - diag_At_A(D) + diag_At_A(f.data.B' * D)
    return μ, Σ_diag
end

#
# NonCentered Parametrization.
#

@doc raw"""
    posterior(sva::SparseVariationalApproximation{NonCentered})

Compute the approximate posterior [1] over the process `f =
sva.fz.f`, given inducing inputs `z = sva.fz.x` and a variational
distribution over inducing points `sva.q` (which represents ``q(ε)``
where `ε = cholesky(cov(fz)).U' \ (f(z) - mean(f(z)))`). The approximate posterior at test
points ``x^*`` where ``f^* = f(x^*)`` is then given by:

```math
q(f^*) = \int p(f | ε) q(ε) du
```
which can be found in closed form.

[1] - Hensman, James, Alexander Matthews, and Zoubin Ghahramani. "Scalable
variational Gaussian process classification." Artificial Intelligence and
Statistics. PMLR, 2015.
"""
function AbstractGPs.posterior(approx::SparseVariationalApproximation{NonCentered})
    fz = approx.fz
    data = (Cuu=_chol_cov(fz), C_ε=_chol_cov(approx.q))
    return ApproxPosteriorGP(approx, fz.f, data)
end

#
# Various methods implementing the Internal AbstractGPs API.
# See AbstractGPs.jl API docs for more info.
#

# Produces a matrix that is consistently referred to as A in this file. A more descriptive
# name is, unfortunately, not obvious. It's just an intermediate quantity that happens to
# get used a lot.
_A(f, x) = f.data.Cuu.U' \ cov(f.prior, inducing_points(f), x)

function Statistics.mean(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation{NonCentered}}, x::AbstractVector
)
    return mean(f.prior, x) + _A(f, x)' * mean(f.approx.q)
end

function Statistics.cov(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation{NonCentered}}, x::AbstractVector
)
    A = _A(f, x)
    return cov(f.prior, x) - At_A(A) + Xt_A_X(f.data.C_ε, A)
end

function Statistics.var(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation{NonCentered}}, x::AbstractVector
)
    A = _A(f, x)
    return var(f.prior, x) - diag_At_A(A) + diag_Xt_A_X(f.data.C_ε, A)
end

function Statistics.cov(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation{NonCentered}},
    x::AbstractVector,
    y::AbstractVector,
)
    Ax = _A(f, x)
    Ay = _A(f, y)
    return cov(f.prior, x, y) - Ax'Ay + Xt_A_Y(Ax, f.data.C_ε, Ay)
end

function StatsBase.mean_and_cov(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation{NonCentered}}, x::AbstractVector
)
    A = _A(f, x)
    μ = mean(f.prior, x) + A' * mean(f.approx.q)
    Σ = cov(f.prior, x) - At_A(A) + Xt_A_X(f.data.C_ε, A)
    return μ, Σ
end

function StatsBase.mean_and_var(
    f::ApproxPosteriorGP{<:SparseVariationalApproximation{NonCentered}}, x::AbstractVector
)
    A = _A(f, x)
    μ = mean(f.prior, x) + A' * mean(f.approx.q)
    Σ = var(f.prior, x) - diag_At_A(A) + diag_Xt_A_X(f.data.C_ε, A)
    return μ, Σ
end

#
# Misc utility.
#

inducing_points(f::ApproxPosteriorGP{<:SparseVariationalApproximation}) = f.approx.fz.x

#
# elbo
#

"""
    elbo(
        sva::SparseVariationalApproximation,
        fx::FiniteGP,
        y::AbstractVector{<:Real};
        num_data=length(y),
        quadrature=DefaultQuadrature(),
    )

Compute the Evidence Lower BOund from [1] for the process `f = fx.f ==
svgp.fz.f` where `y` are observations of `fx`, pseudo-inputs are given by `z =
svgp.fz.x` and `q(u)` is a variational distribution over inducing points `u =
f(z)`.

`quadrature` selects which method is used to calculate the expected loglikelihood in
the ELBO. The options are: `DefaultQuadrature()`, `Analytic()`, `GaussHermite()` and
`MonteCarlo()`. For likelihoods with an analytic solution, `DefaultQuadrature()` uses this
exact solution. If there is no such solution, `DefaultQuadrature()` either uses
`GaussHermite()` or `MonteCarlo()`, depending on the likelihood.

N.B. the likelihood is assumed to be Gaussian with observation noise `fx.Σy`.
Further, `fx.Σy` must be isotropic - i.e. `fx.Σy = α * I`.

[1] - Hensman, James, Alexander Matthews, and Zoubin Ghahramani. "Scalable
variational Gaussian process classification." Artificial Intelligence and
Statistics. PMLR, 2015.
"""
function AbstractGPs.elbo(
    sva::SparseVariationalApproximation,
    fx::FiniteGP{<:AbstractGP,<:AbstractVector,<:Diagonal{<:Real,<:Fill}},
    y::AbstractVector{<:Real};
    num_data=length(y),
    quadrature=DefaultQuadrature(),
)
    @assert sva.fz.f === fx.f
    return _elbo(quadrature, sva, fx, y, GaussianLikelihood(fx.Σy[1]), num_data)
end

function AbstractGPs.elbo(
    ::SparseVariationalApproximation, ::FiniteGP, ::AbstractVector; kwargs...
)
    return error(
        "The observation noise fx.Σy must be homoscedastic.\n",
        "To avoid this error, construct fx using: f = GP(kernel); fx = f(x, σ²)",
        ", where σ² is a positive Real.",
    )
end

"""
    elbo(
        sva::SparseVariationalApproximation,
        lfx::LatentFiniteGP,
        y::AbstractVector;
        num_data=length(y),
        quadrature=DefaultQuadrature(),
    )

Compute the ELBO for a LatentGP with a possibly non-conjugate likelihood.
"""
function AbstractGPs.elbo(
    sva::SparseVariationalApproximation,
    lfx::LatentFiniteGP,
    y::AbstractVector;
    num_data=length(y),
    quadrature=DefaultQuadrature(),
)
    @assert sva.fz.f === lfx.fx.f
    return _elbo(quadrature, sva, lfx.fx, y, lfx.lik, num_data)
end

# Compute the common elements of the ELBO
function _elbo(
    quadrature::QuadratureMethod,
    sva::SparseVariationalApproximation,
    fx::FiniteGP,
    y::AbstractVector,
    lik,
    num_data::Integer,
)
    @assert sva.fz.f === fx.f
    post = posterior(sva)
    q_f = marginals(post(fx.x))
    variational_exp = expected_loglik(quadrature, y, q_f, lik)

    n_batch = length(y)
    scale = num_data / n_batch
    return sum(variational_exp) * scale - kl_term(sva, post)
end

kl_term(sva::SparseVariationalApproximation{Centered}, post) = KL(sva.q, sva.fz)

function kl_term(sva::SparseVariationalApproximation{NonCentered}, post)
    m_ε = mean(sva.q)
    return (tr(cov(sva.q)) + m_ε'm_ε - length(m_ε) - logdet(post.data.C_ε)) / 2
end

ext: jl


text:
module TestPromises

using Base.Experimental: @sync
using Reagents
using Test

include("../../../examples/promises.jl")

function Base.show(io::IO, mime::MIME"text/plain", @nospecialize(p::Promise{T})) where {T}
    show(io, Promise)
    print(io, '{')
    show(io, T)
    print(io, '}')
    println(io, ":")
    print(io, "  value: ")
    show(IOContext(io, :typeinfo => typeof(p.value)), mime, p.value)
    println(io)
    print(io, "  send: ")
    show(io, mime, p.send)
    println(io)
    print(io, "  receive: ")
    show(io, mime, p.receive)
end

end  # module


text:
using BioSimulatorPetriNets
using Test

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


text:
using Pkg
Pkg.add(Pkg.PackageSpec(;name="LocalRegistry", version="0.3.2"))
using LocalRegistry
using Git
using TOML

# Read input parameters
registry_name = ARGS[1]
registry_url = ARGS[2]
registry_branch = ARGS[3]
dry_run = ARGS[4]

if lowercase(dry_run) != "true" && lowercase(dry_run) != "yes"
    # Add the private Julia registry
    Pkg.Registry.add(RegistrySpec(url=registry_url))
    pkg"develop ."

    # Read the Project.toml file in the package
    GITHUB_REPOSITORY = ENV["GITHUB_REPOSITORY"]
    TOKEN = ""
    URL = ""
    if haskey(ENV, "GITHUB_TOKEN")
        TOKEN = ENV["GITHUB_TOKEN"]
        URL = "https://x-access-token:$(TOKEN)@github.com/$(GITHUB_REPOSITORY).git"
    end

    fname = "Project.toml";
    dict_project = TOML.parsefile(fname)
    VERSION = dict_project["version"]
    run(`$(git()) tag -a -f "v$(VERSION)" -m "Update version to v$(VERSION)"`)
    if isempty(URL)
        run(`$(git()) push --tags`)
    else
        run(`$(git()) push --tags $(URL)`)
    end

    # Get the package location
    proj_path = pwd()

    # Checkout to the specified branch
    cd("/home/runner/.julia/registries/$(registry_name)/")
    run(`$(git()) fetch --all`)
    run(`$(git()) checkout $(registry_branch)`)

    # Register the new version in the Julia registry
    register(proj_path, registry_name)
    run(`$(git()) push`)
end


text:
using MieSeries;
using BoundaryElements;
using CompScienceMeshes;
using LinearForms;
using FixedSizeArrays;

	#################################################
	#  Xlm and Ulm testing

type VectorSphericalHarmonicX <: BoundaryElements.Functional
	l::Int
	m::Int
end

type VectorSphericalHarmonicU <: BoundaryElements.Functional
	l::Int
	m::Int
end

function (vshX::VectorSphericalHarmonicX)(p::MeshPointNM)
	r = cartesian(p);
	l = vshX.l;
	m = vshX.m;
	lm = l*(l+1)+m+1;
	Xx, Xy, Xz, Yx, Yy, Yz, Ux, Uy, Uz = MieSeries.Xlm(l+1, r);
	return Vec(Xx[lm], Xy[lm], Xz[lm]);
end

function (vshX::VectorSphericalHarmonicU)(p::MeshPointNM)
	r = cartesian(p);
	l = vshX.l;
	m = vshX.m;
	lm = l*(l+1)+m+1;
	Xx, Xy, Xz, Yx, Yy, Yz, Ux, Uy, Uz = MieSeries.Xlm(l+1, r);
	return Vec(Ux[lm], Uy[lm], Uz[lm]);
end

function BoundaryElements.integrand(vshX::VectorSphericalHarmonicX, tv, fv)
	return dot(tv[1], fv);
end

function BoundaryElements.integrand(vshU::VectorSphericalHarmonicU, tv, fv)
	return dot(tv[1], fv);
end

	#########################################
	# Simulation parameters

const epsilon0 = 8.8541878176e-12;
const mu0 = pi*4e-7;
const eta0 = sqrt(mu0/epsilon0);
const c0 = 1.0/sqrt(epsilon0*mu0);

wavelength = 5.0;
frequency = c0/wavelength;
omega = 2.0 * pi * frequency;
k = omega * sqrt(epsilon0*mu0);

Z = (0.7+0.3*im)*eta0;
delta = 1.2;

#l = 2;
l = 1
m = -1;

	####################################
	# Operators

R = 4.0;
h = 1.2;
gamma = meshsphere(R, h);
f = raviartthomas(gamma);
g = buffachristiansen(gamma);

#Ts = MWWeaklySingular(k);
#Th = MWHyperSingular(k);
Ts = MWWeaklySingular(im*k);
Th = MWHyperSingular(im*k);
K = MWDoubleLayer3D(im*k);
N = NCross();
I = Identity();
Tdelta = MWWeaklySingular(1/delta);

println("Ts_nxff"); Ts_nxff = 1.0*assemble(Ts,f,f);
println("Th_nxff"); Th_nxff = -1.0*assemble(Th,f,f); # T = -ikTs + Th/ik is preferred instead of -ikTs - Th/ik
println("K_nxff"); K_nxff = 1.0*assemble(K,f,f);
println("K_nxfg"); K_nxfg = assemble(K,f,g);
println("Gmix = I_nxfg"); Gmix = -assemble(N,f,g); # <nxf,g> = -<f,Ng>
println("IGmix"); IGmix = inv(Gmix);
println("I_ff = Gram"); I_ff = assemble(I,f,f);
println("IGram"); IGram = inv(I_ff);
println("N_ff = I_fnxf"); N_ff = assemble(N,f,f);
println("Nd_nxff"); Nd_nxff = (2.0/delta)*assemble(Tdelta,f,f);
# =#

	###########################################
	# Vector Spherical harmonics testing

VSH_Xlm = VectorSphericalHarmonicX(l, m);
VSH_Ulm = VectorSphericalHarmonicU(l, m);

println("Xlm_f"); Xlm_f = assemble(VSH_Xlm, f);
println("Ulm_f"); Ulm_f = assemble(VSH_Ulm, f);
println("Xlm_nxf"); Xlm_nxf = -Ulm_f;  # <nxf,X> = -<f,U>
println("Ulm_nxf"); Ulm_nxf = Xlm_f;   # <nxf,U> = <f,X>

	##########################################
	# Special functions

#   input :
# order  integer order of the functions
# arg    complex argument of the functions
#   output :
# j      Spherical Bessel
# h1     Spherical Hankel of 1st kind
# h2     Spherical Hankel of 2nd kind
# dj     Derivative of Spherical Bessel
# dh1    Derivative of Spherical Hankel of 1st kind
# dh2    Derivative of Spherical Hankel of 2nd kind
# J      Ricatti Bessel
# H1     Ricatti Hankel of 1st kind
# H2     Ricatti Hankel of 2nd kind
# dJ     Derivative of Ricatti Bessel
# dH1    Derivative of Ricatti Hankel of 1st kind
# dH2    Derivative of Ricatti Hankel of 2nd kind
function BesselFamilyFunctions(order, arg)
	j = MieSeries.sbessel(order, arg);
	h1 = MieSeries.shankel(order, arg); # j+im*y
	h2 = 2.0*j-h1;  # j-im*y equals conj(h1) only when arg is real
	dj = (order*MieSeries.sbessel(order-1, arg) - (order+1)*MieSeries.sbessel(order+1, arg)) / (2*order+1);
	dh1 = (order*MieSeries.shankel(order-1, arg) - (order+1)*MieSeries.shankel(order+1, arg)) / (2*order+1);
	dh2 = 2.0*dj-dh1;

	J = arg*j;
	H1 = arg*h1;
	H2 = 2.0*J-H1;
	dJ = j + arg*dj; # d(z*j(z))/dz
	dH1 = h1 + arg*dh1;
	dH2 = 2.0*dJ-dH1;

	return (j, h1, h2, dj, dh1, dh2, J, H1, H2, dJ, dH1, dH2);
end

(j, h1, h2, dj, dh1, dh2,
 J, H1, H2, dJ, dH1, dH2) = BesselFamilyFunctions(l, k*R);

Rid = R/(im*delta);
(j_Rid, h1_Rid, h2_Rid, dj_Rid, dh1_Rid, dh2_Rid,
 J_Rid, H1_Rid, H2_Rid, dJ_Rid, dH1_Rid, dH2_Rid) = BesselFamilyFunctions(l, Rid);

println("im*J*dH2-im*dJ*H2 = ($(im*J*dH2-im*dJ*H2)) should be 1");
println("im*J_Rid*dH2_Rid-im*dJ_Rid*H2_Rid = ($(im*J_Rid*dH2_Rid-im*dJ_Rid*H2_Rid)) should be 1");

	#########################################
	# Numerical/Analytical comparisons

# N

num_NX = N_ff*IGram*Xlm_f;
anal_NX = Ulm_f;
re_NX = norm(num_NX-anal_NX)./norm(anal_NX);
println("re_NX : $re_NX");

num_NU = N_ff*IGram*Ulm_f;
anal_NU = -Xlm_f;
re_NU = norm(num_NU-anal_NU)./norm(anal_NU);
println("re_NU : $re_NU");

# T

T_nxff = -im*k*Ts_nxff + 1.0/(im*k)*Th_nxff;

num_TX = T_nxff*IGram*Xlm_f;
anal_TX = -J*H2*Ulm_nxf;
re_TX = norm(num_TX-anal_TX)./norm(anal_TX);
println("re_TX : $re_TX");

num_TU = T_nxff*IGram*Ulm_f;
anal_TU = dJ*dH2*Xlm_nxf;
re_TU = norm(num_TU-anal_TU)./norm(anal_TU);
println("re_TU : $re_TU");

# Ts

num_TsX = Ts_nxff*IGram*Xlm_f;
anal_TsX = 1.0/(im*k)*J*H2*Ulm_nxf;
re_TsX = norm(num_TsX-anal_TsX)./norm(anal_TsX);
println("re_TsX : $re_TsX");

num_TsU = Ts_nxff*IGram*Ulm_f;
anal_TsU = 1.0/(-im*k)*(dJ*dH2+l*(l+1)*j*h2)*Xlm_nxf;
re_TsU = norm(num_TsU-anal_TsU)./norm(anal_TsU);
println("re_TsU : $re_TsU");

# Th

# Print the absolute error for the numerical Th(X) because it equals 0 analytically
# Compare it relatively to T(X)
num_ThX = Th_nxff*IGram*Xlm_f;
anal_ThX = 0*Ulm_nxf;
abserror_ThX = norm(num_ThX);
reToTX_ThX = abserror_ThX / norm(anal_TX);
println("abserror_ThX = $(abserror_ThX) (norm(anal_TX) = $(norm(anal_TX)), reToTX_ThX = $reToTX_ThX)");

num_ThU = Th_nxff*IGram*Ulm_f;
anal_ThU = -im*k*l*(l+1)*j*h2*Xlm_nxf;
re_ThU = norm(num_ThU-anal_ThU)./norm(anal_ThU);
println("re_ThU : $re_ThU");

# K

num_KX = K_nxff*IGram*Xlm_f;
anal_KX = (im*J*dH2-0.5)*Xlm_nxf;
re_KX = norm(num_KX-anal_KX)./norm(anal_KX);
println("re_KX : $re_KX");

num_KU = K_nxff*IGram*Ulm_f;
anal_KU = (-im*dJ*H2-0.5)*Ulm_nxf;
re_KU = norm(num_KU-anal_KU)./norm(anal_KU);
println("re_KU : $re_KU");

# (I/2-K)N

Nover2minusKN_nxff = (0.5*Gmix-K_nxfg)*IGmix*I_ff;

num_Nover2minusKNX = Nover2minusKN_nxff*IGram*Xlm_f;
anal_Nover2minusKNX = im*J*dH2*Ulm_nxf;
re_Nover2minusKNX = norm(num_Nover2minusKNX-anal_Nover2minusKNX)./norm(anal_Nover2minusKNX);
println("re_Nover2minusKNX : $re_Nover2minusKNX");

num_Nover2minusKNU = Nover2minusKN_nxff*IGram*Ulm_f;
anal_Nover2minusKNU = im*dJ*H2*Xlm_nxf;
re_Nover2minusKNU = norm(num_Nover2minusKNU-anal_Nover2minusKNU)./norm(anal_Nover2minusKNU);
println("re_Nover2minusKNU : $re_Nover2minusKNU");

# S = eta * T - Z * (I/2-K)N

S_nxff = eta0 * T_nxff - Z * Nover2minusKN_nxff;

num_SX = S_nxff*IGram*Xlm_f;
anal_SX = (-eta0*H2-im*Z*dH2)*J*Ulm_nxf;
re_SX = norm(num_SX-anal_SX)./norm(anal_SX);
println("re_SX : $re_SX");

num_SU = S_nxff*IGram*Ulm_f;
anal_SU = (eta0*dH2-im*Z*H2)*dJ*Xlm_nxf;
re_SU = norm(num_SU-anal_SU)./norm(anal_SU);
println("re_SU : $re_SU");

# N delta

num_NdX = Nd_nxff*IGram*Xlm_f;
anal_NdX = 2.0*J_Rid*H2_Rid*Ulm_nxf;
re_NdX = norm(num_NdX-anal_NdX)./norm(anal_NdX);
println("re_NdX : $re_NdX");

num_NdU = Nd_nxff*IGram*Ulm_f;
anal_NdU = -2.0*(dJ_Rid*dH2_Rid+l*(l+1)*j_Rid*h2_Rid)*Xlm_nxf;
re_NdU = norm(num_NdU-anal_NdU)./norm(anal_NdU);
println("re_NdU : $re_NdU");


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm

function extendlines(V::Lar.Points, EV::Lar.Cells, s=1.001)
	@assert size(V,1)==2 "extendlines() is for 2D models!"
	segments = [hcat(V[:,v1], V[:,v2]) for (v1,v2) in EV]
	n = length(EV)

	for segi in segments
		# renaming of segment vertices
		x1 = segi[1,1]; x2 = segi[1,2];
		y1 = segi[2,1]; y2 = segi[2,2];
		# middle segment point
		tx = sum(x1+x2)/2
		ty = sum(y1+y2)/2
		# composite translation and scaling transformation
		segi[1,1] = s*(x1-tx)+tx; segi[1,2] = s*(x2-tx)+tx;
		segi[2,1] = s*(y1-ty)+ty; segi[2,2] = s*(y2-ty)+ty;
	end
	V = zeros(Float64, 2, 2*length(EV))
	# new vertices for expanded line segments
	for (i,segi) in enumerate(segments)
		V[:,i] = segi[:,1]
		V[:,n+i] = segi[:,2]
	end
	EV = [[k,n+k] for k=1:length(segments)]
	return V, EV
end

n=100; t=0.6
V = zeros(Float64,2,2*n)
EV = [zeros(Int64,2) for k=1:n]

for k=1:n
	v1 = rand(Float64,2)
	v2 = rand(Float64,2)
	vm = (v1+v2)/2
	transl = rand(Float64,2)
	V[:,k] = (v1-vm)*t + transl
	V[:,n+k] = (v2-vm)*t + transl
	EV[k] = [k,n+k]
end

V = Plasm.normalize(V)
model = (V,EV)
Sigma = Lar.spaceindex(model)

Plasm.view(Plasm.numbering(.1)((V,[[[k] for k=1:size(V,2)], EV])))
Z,EZ = extendlines(V,EV,1.1)
Plasm.view(Plasm.numbering(.1)((Z,[[[k] for k=1:size(Z,2)], EZ])))

Plasm.view(V,EV)
model = V,EV;
model2 = Z,EZ;
W,EW = Lar.fragmentlines(model);
Z,EZ = Lar.fragmentlines(model2);

Plasm.viewexploded(W,EW)(1.2,1.2,1.2)
Plasm.viewexploded(Z,EZ)(1.2,1.2,1.2)
W,EW = Lar.fragmentlines(model);
Z,EZ = Lar.fragmentlines(model2);
U,EVs = Lar.biconnectedComponent((W,EW::Lar.Cells));
T,ETs = Lar.biconnectedComponent((Z,EZ::Lar.Cells));
# 2-connected components (H & T)

W = convert(Lar.Points, U')
T = convert(Lar.Points, T')
cop_EV = Lar.coboundary_0(EW::Lar.Cells)
cop_EW = convert(Lar.ChainOp, cop_EV)
cop_EZ = Lar.coboundary_0(EZ::Lar.Cells)
cop_EZ = convert(Lar.ChainOp, cop_EZ)
V, copEV, copFE = Lar.Arrangement.planar_arrangement(W::Lar.Points, cop_EW::Lar.ChainOp)
Z, copEZ, copFE2 = Lar.Arrangement.planar_arrangement(T::Lar.Points, cop_EZ::Lar.ChainOp)

triangulated_faces = Lar.triangulate2D(V, [copEV, copFE])
triangulated_faces2 = Lar.triangulate2D(Z, [copEZ, copFE2])
FVs = convert(Array{Lar.Cells}, triangulated_faces)
FZs = convert(Array{Lar.Cells}, triangulated_faces2)
W = convert(Lar.Points, V')
Z = convert(Lar.Points, Z')
Plasm.viewcolor(W::Lar.Points, FVs::Array{Lar.Cells})
Plasm.viewcolor(Z::Lar.Points, FZs::Array{Lar.Cells})

model = U,EVs;
model2 = T',ETs;
Plasm.view(Plasm.lar_exploded(model)(1.2,1.2,1.2))
Plasm.view(Plasm.lar_exploded(model2)(1.2,1.2,1.2))


text:
using Serialization

"""
Supertype for types encoding the family of perturbations allowed.
"""
abstract type PerturbationFamily end

struct UnrestrictedPerturbationFamily <: PerturbationFamily end
Base.show(io::IO, pp::UnrestrictedPerturbationFamily) = print(io, "unrestricted")

abstract type RestrictedPerturbationFamily <: PerturbationFamily end

"""
For blurring perturbations, we currently allow colors to "bleed" across color channels -
that is, the value of the output of channel 1 can depend on the input to all channels.
(This is something that is worth reconsidering if we are working on color input).
"""
struct BlurringPerturbationFamily <: RestrictedPerturbationFamily
    blur_kernel_size::NTuple{2}
end
Base.show(io::IO, pp::BlurringPerturbationFamily) =
    print(io, filter(x -> !isspace(x), "blur-$(pp.blur_kernel_size)"))

struct LInfNormBoundedPerturbationFamily <: RestrictedPerturbationFamily
    norm_bound::Real

    function LInfNormBoundedPerturbationFamily(norm_bound::Real)
        @assert(norm_bound > 0, "Norm bound $(norm_bound) should be positive")
        return new(norm_bound)
    end
end
Base.show(io::IO, pp::LInfNormBoundedPerturbationFamily) =
    print(io, "linf-norm-bounded-$(pp.norm_bound)")

struct BrightnessPerturbationFamily <: RestrictedPerturbationFamily
    bright_bound::Real

    function BrightnessPerturbationFamily(bright_bound::Real)
        @assert(bright_bound > 0, "Norm bound $(bright_bound) should be positive")
        return new(bright_bound)
    end
end
Base.show(io::IO, pp::BrightnessPerturbationFamily) =
    print(io, "brightness-bounded-$(pp.bright_bound)")

function get_model(
    nn::NeuralNet,
    input::Array{<:Real},
    pp::PerturbationFamily,
    optimizer,
    tightening_options::Dict,
    tightening_algorithm::TighteningAlgorithm,
)::Dict{Symbol,Any}
    notice(
        MIPVerify.LOGGER,
        "Determining upper and lower bounds for the input to each non-linear unit.",
    )
    m = Model(optimizer_with_attributes(optimizer, tightening_options...))
    m.ext[:MIPVerify] = MIPVerifyExt(tightening_algorithm)

    d_common = Dict(
        :Model => m,
        :PerturbationFamily => pp,
        :TighteningApproach => string(tightening_algorithm),
    )

    return merge(d_common, get_perturbation_specific_keys(nn, input, pp, m))
end

function get_model(
    nn::NeuralNet,
    input_size,
    global_predicted,
    global_confidence,
    target_indexes,
    pp::PerturbationFamily,
    optimizer,
    tightening_options::Dict,
    tightening_algorithm::TighteningAlgorithm,
)::Dict{Symbol,Any}
    notice(
        MIPVerify.LOGGER,
        "Determining upper and lower bounds for the REGULAR input to each non-linear unit.",
    )
    m = Model(optimizer_with_attributes(optimizer, tightening_options...))
    m.ext[:MIPVerify] = MIPVerifyExt(tightening_algorithm)

    d_common = Dict(
        :Model => m,
        :PerturbationFamily => pp,
        :TighteningApproach => string(tightening_algorithm),
    )

    input_range = CartesianIndices(input_size)
    # input is the input without perturbation added
    input = map(_ -> @variable(m, lower_bound = 0, upper_bound = 1), input_range)
    d_common[:Input] = input

    v_output = input |> nn

    target_vars = v_output[Bool[i ∈ target_indexes for i in 1:length(v_output)]]
    @constraint(m, target_vars .<= v_output[global_predicted] - global_confidence)

    notice(
        MIPVerify.LOGGER,
        "Determining upper and lower bounds for the PERTUBATED input to each non-linear unit.",
    )

    m.ext[:MIPVerify] = MIPVerifyExt(tightening_algorithm)

    return merge(d_common, get_perturbation_specific_keys(nn, input, pp, m))
end

function get_perturbation_specific_keys(
    nn::NeuralNet,
    input::Array{<:Real},
    pp::UnrestrictedPerturbationFamily,
    m::Model,
)::Dict{Symbol,Any}
    input_range = CartesianIndices(size(input))

    # v_x0 is the input with the perturbation added
    v_x0 = map(_ -> @variable(m, lower_bound = 0, upper_bound = 1), input_range)

    v_output = v_x0 |> nn

    return Dict(:PerturbedInput => v_x0, :Perturbation => v_x0 - input, :Output => v_output)
end

function get_perturbation_specific_keys(
    nn::NeuralNet,
    input::Array{<:Real},
    pp::BlurringPerturbationFamily,
    m::Model,
)::Dict{Symbol,Any}

    input_size = size(input)
    num_channels = size(input)[4]
    filter_size = (pp.blur_kernel_size..., num_channels, num_channels)

    v_f = map(_ -> @variable(m, lower_bound = 0, upper_bound = 1), CartesianIndices(filter_size))
    @constraint(m, sum(v_f) == num_channels)
    v_x0 = map(_ -> @variable(m, lower_bound = 0, upper_bound = 1), CartesianIndices(input_size))
    @constraint(m, v_x0 .== input |> Conv2d(v_f))

    v_output = v_x0 |> nn

    return Dict(
        :PerturbedInput => v_x0,
        :Perturbation => v_x0 - input,
        :Output => v_output,
        :BlurKernel => v_f,
    )
end

function get_perturbation_specific_keys(
    nn::NeuralNet,
    input::Array{<:Real},
    pp::LInfNormBoundedPerturbationFamily,
    m::Model,
)::Dict{Symbol,Any}

    input_range = CartesianIndices(size(input))
    # v_e is the perturbation added
    v_e = map(
        _ -> @variable(m, lower_bound = -pp.norm_bound, upper_bound = pp.norm_bound),
        input_range,
    )
    # v_x0 is the input with the perturbation added
    v_x0 = map(
        i -> @variable(
            m,
            lower_bound = max(0, input[i] - pp.norm_bound),
            upper_bound = min(1, input[i] + pp.norm_bound)
        ),
        input_range,
    )
    @constraint(m, v_x0 .== input + v_e)

    v_output = v_x0 |> nn

    return Dict(:PerturbedInput => v_x0, :Perturbation => v_e, :Output => v_output)
end

function get_perturbation_specific_keys(
    nn::NeuralNet,
    input::Array{<:Real},
    pp::BrightnessPerturbationFamily,
    m::Model,
)::Dict{Symbol,Any}

    input_range = CartesianIndices(size(input))
    # v_e is the perturbation added
    v_tmp = @variable(m, lower_bound = -pp.bright_bound, upper_bound = pp.bright_bound)
    v_e = map(
        _ -> v_tmp,
        input_range,
    )
    # v_x0 is the input with the perturbation added
    v_x0 = map(
        i -> @variable(
            m,
            lower_bound = max(0, input[i] - pp.bright_bound),
            upper_bound = min(1, input[i] + pp.bright_bound)
        ),
        input_range,
    )
    @constraint(m, v_x0 .== input + v_tmp)

    v_output = v_x0 |> nn

    return Dict(:PerturbedInput => v_x0, :Perturbation => v_e, :Output => v_output)
end

function get_perturbation_specific_keys(
    nn::NeuralNet,
    input::AbstractArray{<:JuMPLinearType},
    pp::LInfNormBoundedPerturbationFamily,
    m::Model,
)::Dict{Symbol,Any}

    input_range = CartesianIndices(size(input))
    # v_e is the perturbation added
    v_e = map(
        _ -> @variable(m, lower_bound = -pp.norm_bound, upper_bound = pp.norm_bound),
        input_range,
    )
    # v_x0 is the input with the perturbation added
    v_x0 = map(
        i -> @variable(
            m,
            lower_bound = 0,
            upper_bound = 1
        ),
        input_range,
    )
    @constraint(m, v_x0 .== input + v_e)

    v_output = v_x0 |> nn

    return Dict(:PerturbedInput => v_x0, :Perturbation => v_e, :Output => v_output)
end

function get_perturbation_specific_keys(
    nn::NeuralNet,
    input::AbstractArray{<:JuMPLinearType},
    pp::BrightnessPerturbationFamily,
    m::Model,
)::Dict{Symbol,Any}

    input_range = CartesianIndices(size(input))
    # v_e is the perturbation added
    v_tmp = @variable(m, lower_bound = -pp.bright_bound, upper_bound = pp.bright_bound)
    v_e = [v_tmp for i in input_range]
    # v_x0 is the input with the perturbation added
    v_x0 = map(
        i -> @variable(
            m,
            lower_bound = 0,
            upper_bound = 1,
        ),
        input_range,
    )
    @constraint(m, v_x0 .== input + v_e)

    v_output = v_x0 |> nn

    return Dict(:PerturbedInput => v_x0, :Perturbation => v_e, :Output => v_output)
end

struct MIPVerifyExt
    tightening_algorithm::TighteningAlgorithm
end


text:
using Test
using Syncer


text:
using DrWatson, GPUAcceleratedTracking, CUDA, Tracking, GNSSSignals, StructArrays, ProgressMeter;
import Tracking: Hz, ms;
@quickactivate "GPUAcceleratedTracking"

N = 2048:32:262_144
err_rel  = zeros(length(N))
@showprogress 0.5 for (idx, num_samples) in enumerate(N)
    # num_samples = 2_048
    num_ants = 1
    num_correlators = 3
    enable_gpu = Val(true)
    system = GPSL1(use_gpu = Val(true));
    # system_h = GPSL1(use_gpu = Val(false));
    codes = system.codes
    # codes_text_mem_simple = CuTexture(
        # CuTextureArray(codes)
    # )
    codes_text_mem = CuTexture(
        CuTextureArray(codes),
        address_mode = CUDA.ADDRESS_MODE_WRAP,
        interpolation = CUDA.NearestNeighbour(),
        normalized_coordinates = true
    )
    code_frequency = get_code_frequency(system)
    code_length = get_code_length(system)
    start_code_phase = 0.0f0
    carrier_phase = 0.0f0
    carrier_frequency = 1500Hz
    prn = 1

    # Generate the signal;
    signal, sampling_frequency = gen_signal(system, prn, carrier_frequency, num_samples, num_ants = NumAnts(num_ants), start_code_phase = start_code_phase, start_carrier_phase = carrier_phase)

    # Generate correlator;
    correlator = EarlyPromptLateCorrelator(NumAnts(num_ants), NumAccumulators(num_correlators))
    correlator_sample_shifts = get_correlator_sample_shifts(system, correlator, sampling_frequency, 0.5)
    num_of_shifts = correlator_sample_shifts[end] - correlator_sample_shifts[1]

    # Generate blank code and carrier replica, and downconverted signal;
    # code_replica_cpu = zeros(Float32, num_samples + num_of_shifts)
    code_replica = CUDA.zeros(Float32, num_samples + num_of_shifts)
    code_replica_text_mem = CUDA.zeros(Float32, num_samples + num_of_shifts)

    # Generate CUDA kernel tuning parameters;
    threads_per_block = 768
    blocks_per_grid = cld.(num_samples, threads_per_block)

    @cuda threads=threads_per_block blocks=blocks_per_grid gen_code_replica_texture_mem_kernel!(
        code_replica_text_mem,
        codes_text_mem, # texture memory codes
        code_frequency,
        sampling_frequency,
        start_code_phase,
        prn,
        num_samples + num_of_shifts,
        correlator_sample_shifts[1],
        code_length
    )
    @cuda threads=threads_per_block blocks=blocks_per_grid gen_code_replica_kernel!(
        code_replica,
        codes, # texture memory codes
        code_frequency,
        sampling_frequency,
        start_code_phase,
        prn,
        num_samples + num_of_shifts,
        correlator_sample_shifts[1],
        code_length
    )
    # Tracking.gen_code_replica!(code_replica_cpu, system, code_frequency, sampling_frequency, start_code_phase, 1, num_samples, correlator_sample_shifts, prn)

    # code_replica_h = Array(code_replica)
    # code_replica_text_mem_h = Array(code_replica_text_mem)
    # signal = StructArray{ComplexF32}((ones(Float32, num_samples), zeros(Float32, num_samples) ))
    # code_phases = get_code_frequency(system) / sampling_frequency .* (0:num_samples-1) .+ start_code_phase
    # spread_signal = StructArray(signal .* system_h.codes[1 .+ mod.(floor.(Int, code_phases), get_code_length(system)), prn])
    # accums_true = Tracking.correlate(correlator, spread_signal, code_replica_cpu, correlator_sample_shifts, 1, num_samples)
    # accums = Tracking.correlate(correlator, spread_signal, code_replica_h, correlator_sample_shifts, 1, num_samples)
    # accums_text_mem = Tracking.correlate(correlator, spread_signal, code_replica_text_mem_h, correlator_sample_shifts, 1, num_samples)

    # err_rel = sum(abs.(code_replica - code_replica_text_mem)) / num_samples
    err_rel[idx] = sum(abs.(code_replica - code_replica_text_mem)) / num_samples
end

x = vec(collect(N / 0.001)) # convert to Hz
data = vec(100 .* err_rel)
data_bar = mean(data)
data_med = median(data)
data_max = maximum(data)

using CairoMakie
fig = Figure(font = "Times New Roman")
ax = Axis(
    fig,
    xlabel = "Sampling Frequency [Hz]",
    ylabel = "Relative Code Phase Error [%]",
    xscale = log10,
    title = "Relative code phase error of the texture memory code replica generation for 1 ms GPS L1 C/A signal",
    # xlim = [0 400_000],
    xminorgridvisible = true,
    xminorticksvisible = true,
    xminorticks = IntervalsBetween(9),
    # yticks = (10.0 .^(-5:1:-3)),
        # yticks = (10.0 .^(-5:1:-3)),

    xticklabelsize = 18,
    yticklabelsize = 18
)
xlims!(ax, 10^6, 5*10^8)
# string_data_bar = "$(round(data_bar, sigdigits=3))%"
# string_data_max = "$(round(data_max, sigdigits=3))%"
# string_data_med = "$(round(data_med, sigdigits=3))%"

# # textmu = "μ = " * string_data_bar
# # textmax = "max = " * string_data_max
# textmed = "median = " * string_data_med
lines!(ax, x,data)
# hlines!(ax, data_bar, color = :dimgrey, linestyle = :dash)
# hlines!(ax, data_med, color = :dimgrey, linestyle = :dash)
# hlines!(ax, data_max, color = :dimgrey, linestyle = :dash)
# text!(textmu, position = (9*10^8, 0.1+data_bar), align = (:right, :baseline))
# text!(textmax, position = (9*10^8, data_max - 0.2), align = (:right, :baseline))
# text!(textmed, position = (5*10^8, 0.5+data_med), align = (:center, :baseline))

fig[1,1] = ax
fig

@quickactivate "GPUAcceleratedTracking"

raw_data_df = collect_results(datadir("benchmarks/codereplica"))

sort!(raw_data_df, :num_samples)
samples = unique(Vector{Int64}(raw_data_df[!, :num_samples]))
algorithm_names = unique(Vector{String}(raw_data_df[!, :algorithm]))
samples = unique(Vector{Int64}(raw_data_df[!, :num_samples]))
x = samples ./ 0.001 # convert to Hz
algorithm_names = unique(Vector{String}(raw_data_df[!, :algorithm]))

# fig = Figure(
#     # resolution = (1000, 700),
#     font = "Times New Roman"
# )
ax2 = Axis(
    fig,
    xlabel = "Sampling Frequency [Hz]",
    ylabel = "Generation Time [s]",
    xscale = log10,
    yscale = log10,
    title = "Comparison between global memory and texture memory code replica generation for 1 ms GPS L1 C/A signal",
    xminorgridvisible = true,
    xminorticksvisible = true,
    xminorticks = IntervalsBetween(9),
    # yticks = (10.0 .^(-5:1:-3)),
    xticklabelsize = 18,
    yticklabelsize = 18

)
xlims!(ax2, 10^6, 5*10^8)
ylims!(ax2, 1.0e-5, 3.0e-3)

lin = Array{Lines}(undef, length(algorithm_names));
sca = Array{Scatter}(undef, length(algorithm_names));
markers = [:circle, :rect]
for (idx, name) = enumerate(algorithm_names)
    time = 10 ^ (-9) * vec((
        raw_data_df |>
        @filter(
            _.algorithm == name
        ) |> DataFrame
    ).Minimum)

    lin[idx] = lines!(
        ax2,
        x,
        time
    )
    sca[idx] = scatter!(
        ax2,
        x,
        time,
        marker = markers[idx],
        markersize = 15
    )
end
realtime = lines!(ax2, [10^6, 5 * 10^8], [10 ^ (-3), 10 ^ (-3)], color=:grey, linestyle=:dashdot)

fig[2,1] = ax2
fig
elements = [[lin[1] sca[1]], [lin[2] sca[2]]]
labels = ["Global Memory", "Texture Memory"]
axislegend(ax2, elements, labels, "Code Replication Algorithms",  position = :lt)
fig
# save( plotsdir("benchmark_textmem.pdf"), fig)

save(plotsdir("code_phase.pdf"), fig)

max_stars_repo_name: mfkiwl/GPUAcceleratedTracking

max_issues_repo_name: mfkiwl/GPUAcceleratedTracking

max_forks_repo_name: mfkiwl/GPUAcceleratedTracking


text:
using MultipleViewGeometry, Base.Test
using BenchmarkTools
using StaticArrays

ℳ = [Point2DH(x,y,1)  for x=-1000:0.5:1000 for y=-1000:0.5:1000]
𝐓 = hartley_transformation(ℳ)
@time hartley_transformation(ℳ)
@time hartley_transformation(ℳ)

@time hartley_normalization(ℳ)
@time hartley_normalization(ℳ)

𝒪, 𝐓 = hartley_normalization(ℳ)

ℳ[1]
𝒪[1]

@time transform(HomogeneousCoordinates(),CanonicalToHartley(),ℳ)


text:
mutable struct Person
    age::Int64

    function Person(age::Int64)
        this = new()
        if age < 0
            println("Age is not valid")
            this.age = 0
            return this
        else
            this.age = age
            return this
        end
    end
end

function yearPasses(self::Person)
    self.age += 1
end

function amIOld(self::Person)
    if self.age < 13
        println("You are young.")
    elseif self.age < 18
        println("You are a teenager.")
    else
        println("You are old.")
    end
end

# a = Person(111)
# yearPasses(a)
# println(a.age)
# amIOld(a)

max_stars_repo_path: Hackerrank/30 Days of Code/Julia/day 04.jl

max_stars_repo_name: Next-Gen-UI/Code-Dynamics

max_issues_repo_name: Next-Gen-UI/Code-Dynamics

max_forks_repo_name: Next-Gen-UI/Code-Dynamics


text:
# include("MyGraph.jl")
include("DagFlowInstance.jl")

include("../ZeroHalfCuts/ZeroHalfCuts.jl")

module DAGIntFlow

using JuMP
using CPLEX
using Printf
using MathProgBase
using ..MyGraphModule
using ..DagFlowInstanceModule
using ..ZeroHalfChvatGom

#=
function createVarDict(graph::MyGraphOut)
    varDict = Dict{Tuple{Int32,Int32},Int64}()
    n = getNNodes(graph)
    id = 1
    for i=1:n
        for node in graph.outNodes[i]
            varDict[(i,node)] = id
            id += 1
        end
    end
    return varDict
end
=#

function setupIP(instance::DagFlowInstance)
    #varDict = createVarDict(g)
    g = instance.g
    n,m = MyGraphModule.size(g)
    # the nodes where flow conservation is active includes all ndoes except the source and the sink:
    flowConsNodes = setdiff(1:n, [instance.sourceNode instance.sinkNode])

    IP = Model(solver=CplexSolver())

    #gInOut = createMyGraph(graph)
    # @@@ low arcs should be general integers!
    @variable(IP, x[allArcs(g)], Bin)

    @objective(IP, Min, sum(a.cost*x[a] for a in allArcs(g)))

    @constraint(IP, [subset in instance.nodePartition], sum(x[a] for node in subset for a in outArcs(g,node)) == 1)
    #=for subset in instance.nodePartition
        println("subset = $subset")
        for node in subset
            println("Node = $node")
            for a in outArcs(g,node)
                println(a)
            end
            throw("stopping!!!")
        end
    end=#
    @constraint(IP, [node in flowConsNodes], sum(x[a] for a in outArcs(g,Int32(node))) == sum(x[a] for a in inArcs(g,Int32(node))))
    @constraint(IP, [outFlowLimit in instance.outFlowLimits], sum(x[a] for a in outArcs(g,outFlowLimit.node)) <= outFlowLimit.limit)
    writeLP(IP, "DagIntFlow.lp", genericnames=true);
    solve(IP)
    println("objective value: ", getobjectivevalue(IP))
end

function writeDOTFile(instance::DagFlowInstanceArcs, filename)
    open(filename, "w") do f
        write(f, "digraph G {\n")
        write(f, "\toverlap=scale;\n")
      	write(f, "\tsplines=true;\n")

        n,m = MyGraphModule.size(instance.g)
        for i=1:n
            println(f,"\t$i [label=\"$(instance.nodeInfoDict[i])\"]")
        end

        for i=1:length(instance.g.arcs)
            arc = instance.g.arcs[i]
            #println("(",src(e),",",dst(e),")" )
            style = ""
            if i in instance.lowArcs
                style = "style=dashed"
            end
            println(f, "\t", arc.from," -> ", arc.to, "[label = \"$i\" $style];")
        end
        write(f, "}\n")
     end
end

function writeDOTFileSol(instance::DagFlowInstanceArcs, xVal, filename)
    open(filename, "w") do f
        write(f, "digraph G {\n")
        write(f, "\toverlap=scale;\n")
      	write(f, "\tsplines=true;\n")

        n,m = MyGraphModule.size(instance.g)
        for i=1:n
            println(f,"\t$i [label=\"$(instance.nodeInfoDict[i])\"]")
        end

        for i=1:length(instance.g.arcs)
            arc = instance.g.arcs[i]
            #println("(",src(e),",",dst(e),")" )
            style = ""
            color = ""
            if i in instance.lowArcs
                style = "style=dashed"
            end
            if xVal[i] > 0.005
                color = "color=\"red\""
            end
            val = round(xVal[i],digits=2)
            println(f, "\t", arc.from," -> ", arc.to, "[label = \"$val\" $style $color];")
        end
        write(f, "}\n")
     end
end

function setupIP(instance::DagFlowInstanceArcs, file_name_instance::String)
    fileNameInstance = splitext(basename(file_name_instance))[1]
    writeDOTFile(instance, "$(fileNameInstance)_graph.dot")
    #varDict = createVarDict(g)
    g = instance.g
    n,m = MyGraphModule.size(g)
    # the nodes where flow conservation is active includes all nodes except the source and the sink:
    flowConsNodes = setdiff(1:n, [instance.sourceNode instance.sinkNode])

    IP = Model(solver=CplexSolver())

    #gInOut = createMyGraph(graph)
    @variable(IP, x[allArcs(g)] >= 0, Int)

    for i in 1:length(allArcs(g))
        arc = g.arcs[i]
        if i in instance.lowArcs
            setupperbound(x[arc], 2.0)
        else
            setupperbound(x[arc],1.0)
        end
    end

    @objective(IP, Min, sum(a.cost*x[a] for a in allArcs(g)))

    @constraint(IP, [subset in instance.arcPartition], sum(x[getArc(g,a)] for a in subset ) == 1)
    @constraint(IP, [node in flowConsNodes], sum(x[a] for a in outArcs(g,Int32(node))) == sum(x[a] for a in inArcs(g,Int32(node))))
    @constraint(IP, [outFlowLimit in instance.outFlowLimits], sum(x[a] for a in outArcs(g,outFlowLimit.node)) <= outFlowLimit.limit)
    writeLP(IP, "$(fileNameInstance)_formulation.lp", genericnames=false)

    JuMP.build(IP)
    MPB=MathProgBase
    iModel = internalmodel(IP)
    zhOrigSubSys = ZeroHalfChvatGom.findZeroHalfOrigSubSystem(MPB.getconstrmatrix(iModel), MPB.getconstrLB(iModel), MPB.getconstrUB(iModel), MPB.getvartype(iModel), false, 1)
    solve(IP, relaxation=true)
    println("objective value: ", getobjectivevalue(IP))

    xValOrig = getvalue(x)
    xVal = Vector{Float64}(undef,length(xValOrig))
    for i=1:length(xValOrig)
        xVal[i] =  xValOrig[getArc(g,i)]
        if xVal[i] > 0.001
            println("x[$i] = $(xVal[i])")
        end
    end

    strength=3
    maxIter = 40
    printCuts = true
    zhsPrintLevel = 0

    done = false
    iter = 1
    while !done && iter <= maxIter
        file_name_sol = "$(fileNameInstance)_sol_$(lpad(iter,2,"0")).dot"
        writeDOTFileSol(instance, xVal, file_name_sol)
        file_name_cuts = "$(fileNameInstance)_cuts_$(lpad(iter,2,"0")).txt"
        zhsParams = ZeroHalfChvatGom.ZeroHalfSepParams(strength,zhsPrintLevel)
        cuts = ZeroHalfChvatGom.ZeroHalfSepTranslateVars(zhOrigSubSys.ASub, zhOrigSubSys.bSub, xVal, 0.0001, zhsParams, zhOrigSubSys.varsSorted)
        f = open(file_name_cuts, "w+")
        for i = 1:length(cuts)
            #println("*** Cut $i***")
            cut = cuts[i]
            vars = cut[1]
            coeffs = cut[2]
            rhs = cut[3]
            if printCuts
                for j=1:length(vars)
                    preString = ""
                    if j > 1
                        preString = " +"
                    end
                    if coeffs[j] == 1.0
                        print(f,"$preString x$(lpad(vars[j],3,"0"))")
                    elseif coeffs[j] == -1.0
                        print(f, " - x$(lpad(vars[j],3,"0"))")
                    else
                        print(f, "$preString $(coeffs[j])*x$(lpad(vars[j],3,"0"))")
                    end
                end
                println(f, " <= $rhs")
            end
            @constraint(IP, sum(coeffs[j] * x[getArc(g,vars[j])] for j=1:length(vars)) <= rhs )
        end

        close(f)
        if length(cuts) > 0
            solve(IP, relaxation=true)

            println("Objective value (iter $iter): ", getobjectivevalue(IP))
            #println("x = ", getvalue(x))
            xValOrig = getvalue(x)
            for i=1:length(xValOrig)
                xVal[i] = xValOrig[getArc(g,i)]
            end
            iter += 1
        else
            done = true
        end
    end
end

end


text:
using MPIMeasurements
using Unitful

# Create Bruker BaseScanner object
bR = BrukerRobot("RobotServer")

################# Use case 1 Basic Move #######################################

# Move absolue LowLevel no safetyCheck!
movePark(bR)
moveCenter(bR)
getPos(bR)
moveAbs(bR, 1.0Unitful.mm, 0.0Unitful.mm, 0.0Unitful.mm)

# Move absolute MidLevel
posXYZ = [1.0Unitful.mm,0.0Unitful.mm,0.0Unitful.mm]
moveAbs(bS, posXYZ)

################# Use case 2 Move and Measure #################################


text:
"""
    Checkpointer{I, T, P, A} <: AbstractOutputWriter

An output writer for checkpointing models to a JLD2 file from which models can be restored.
"""
mutable struct Checkpointer{I, T, P, A} <: AbstractOutputWriter
         frequency :: I
          interval :: T
          previous :: Float64
               dir :: String
            prefix :: String
        properties :: P
    has_array_refs :: A
             force :: Bool
           verbose :: Bool
end

"""
    Checkpointer(model; frequency=nothing, interval=nothing, dir=".", prefix="checkpoint",
                        force=false, verbose=false,
                        properties = [:architecture, :boundary_conditions, :grid, :clock,
                                      :eos, :constants, :closure, :velocities, :tracers,
                                      :timestepper])

Construct a `Checkpointer` that checkpoints the model to a JLD2 file every so often as
specified by `frequency` or `interval`. The `model.clock.iteration` is included in the
filename to distinguish between multiple checkpoint files.

Note that extra model `properties` can be safely specified, but removing crucial properties
such as `:velocities` will make restoring from the checkpoint impossible.

The checkpoint file is generated by serializing model properties to JLD2. However,
functions cannot be serialized to disk (at least not with JLD2). So if a model property
contains a reference somewhere in its hierarchy it will not be included in the checkpoint
file (and you will have to manually restore them).

Keyword arguments
=================
- `frequency::Int`   : Save output every `n` model iterations.
- `interval::Int`    : Save output every `t` units of model clock time.
- `dir::String`      : Directory to save output to. Default: "." (current working directory).
- `prefix::String`   : Descriptive filename prefixed to all output files. Default: "checkpoint".
- `force::Bool`      : Remove existing files if their filenames conflict. Default: `false`.
- `verbose::Bool`    : Log what the output writer is doing with statistics on compute/write times and file sizes.
                       Default: `false`.
- `properties::Array`: List of model properties to checkpoint.
"""
function Checkpointer(model; frequency=nothing, interval=nothing, dir=".", prefix="checkpoint", force=false,
                      verbose=false, properties = [:architecture, :boundary_conditions, :grid, :clock, :coriolis,
                                                   :buoyancy, :closure, :velocities, :tracers, :timestepper])

    validate_interval(frequency, interval)

    # Grid needs to be checkpointed for restoring to work.
    :grid ∉ properties && push!(properties, :grid)

    has_array_refs = Symbol[]

    for p in properties
        isa(p, Symbol) || @error "Property $p to be checkpointed must be a Symbol."
        p ∉ propertynames(model) && @error "Cannot checkpoint $p, it is not a model property!"

        if has_reference(Function, getproperty(model, p))
            @warn "model.$p contains a function somewhere in its hierarchy and will not be checkpointed."
            filter!(e -> e != p, properties)
        end

        has_reference(AbstractField, getproperty(model, p)) && push!(has_array_refs, p)
    end

    mkpath(dir)
    return Checkpointer(frequency, interval, 0.0, dir, prefix, properties, has_array_refs, force, verbose)
end

function write_output(model, c::Checkpointer)
    filepath = joinpath(c.dir, c.prefix * string(model.clock.iteration) * ".jld2")
    c.verbose && @info "Checkpointing to file $filepath..."

    t0 = time_ns()
    jldopen(filepath, "w") do file
        file["checkpointed_properties"] = c.properties
        file["has_array_refs"] = c.has_array_refs

        serializeproperties!(file, model, filter(e -> !(e ∈ c.has_array_refs), c.properties))
        saveproperties!(file, model, filter(e -> e ∈ c.has_array_refs, c.properties))
    end

    t1, sz = time_ns(), filesize(filepath)
    c.verbose && @info "Checkpointing done: time=$(prettytime((t1-t0)/1e9)), size=$(pretty_filesize(sz))"
end

defaultname(::Checkpointer, nelems) = :checkpointer
_arr(::CPU, a) = a
_arr(::GPU, a) = CuArray(a)

function restore_fields!(model, file, arch, fieldset; location="$fieldset")
    if fieldset == :timestepper
        restore_fields!(model.timestepper, file, arch, :Gⁿ; location="timestepper/Gⁿ")
        restore_fields!(model.timestepper, file, arch, :G⁻; location="timestepper/G⁻")
    else
        for p in propertynames(getproperty(model, fieldset))
            getproperty(getproperty(model, fieldset), p).data.parent .= _arr(arch, file[location * "/$p"])
        end
    end
end

"""
    restore_from_checkpoint(filepath; kwargs=Dict())

Restore a model from the checkpoint file stored at `filepath`. `kwargs` can be passed
to the `Model` constructor, which can be especially useful if you need to manually
restore forcing functions or boundary conditions that rely on functions.
"""
function restore_from_checkpoint(filepath; kwargs=Dict())
    file = jldopen(filepath, "r")

    cps = file["checkpointed_properties"]
    has_array_refs = file["has_array_refs"]

    # Restore model properties that were just serialized.
    # We skip fields that contain structs and restore them after model creation.
    for p in cps
        if p ∉ has_array_refs
            kwargs[p] = file["$p"]
        end
    end

    model = Model(; kwargs...)

    # Now restore fields.
    for p in cps
        if p in has_array_refs
            restore_fields!(model, file, model.architecture, p)
        end
    end

    close(file)

    return model
end


text:
struct Boundary{T}
    points::Vector{NTuple{2,T}}

    function Boundary(points::Vector{NTuple{2,T}}) where T
        if first(points) != last(points)
            push!(points,first(points))
        end
        return new{T}(points)
    end
end

const PlasmaBoundary = Boundary
const FluxSurface = Boundary
const Wall = Boundary

function Base.show(io::IO, b::Boundary)
    print(io,"Boundary: npoints = $(length(b.points))")
end

Base.broadcastable(b::Boundary) = (b,)

function Boundary(x::Vector,y::Vector)
    Boundary(collect(zip(x,y)))
end

function Base.getproperty(b::Boundary,s::Symbol)
    if s == :points
        return getfield(b,s)
    end
    if s in (:r,:x)
        v = getfield(b,:points)
        return getindex.(v,1)
    end
    if s in (:z,:y)
        v = getfield(b,:points)
        return getindex.(v,2)
    end
    error("type $(typeof(b)) has no field $s")
end

in_boundary(b::Boundary, p) = PolygonOps.inpolygon(p,b.points,in=true,on=true,out=false)
in_boundary(b::Boundary,x,y) = in_boundary(b,(x,y))

const in_vessel = in_boundary
const in_plasma = in_boundary

function boundary(M::AbstractEquilibrium, psi, dx=0.01,dy=0.01)
    xlims, ylims = limits(M)
    x = range(xlims...,step=dx)
    y = range(ylims...,step=dy)
    Psi = [M(xx,yy) for xx in x, yy in y]
    cl = Contour.contour(x,y,Psi,psi)
    l = Contour.lines(cl)[1]
    xc, yc = Contour.coordinates(l)
    return Boundary(xc,yc)
end

const flux_surface = boundary

function circumference(b::Boundary)
    p = b.points
    return sum(@inbounds norm(p[i+1] .- p[i]) for i=1:(length(p)-1))
end

function average(b::Boundary, F)
    p = b.points
    return sum(@inbounds norm(p[i+1] .- p[i])*(F(p[i][1],p[i][2])) for i=1:(length(p)-1))/circumference(b)
end

area(s::Boundary) = PolygonOps.area(s.points)

function surface_area(b::Boundary)
    p = b.points
    return 2pi*sum(@inbounds norm(p[i+1] .- p[i])*(p[i+1][1] + p[i][1])/2 for i=1:(length(p)-1))
end

function surface_area_average(b::Boundary, F)
    A = surface_area(b)
    p = b.points
    return (2pi/A)*sum(@inbounds norm(p[i+1] .- p[i])*((p[i+1][1] + p[i][1])/2)*
                   (F(p[i+1][1],p[i+1][2]) + F(p[i][1],p[i][2]))/2 for i=1:(length(p)-1))
end

function surface_area_average(b::Boundary, F::Vector)
    A = surface_area(b)
    p = b.points
    return (2pi/A)*sum(@inbounds norm(p[i+1] .- p[i])*((p[i+1][1] + p[i][1])/2)*
                       (F[i+1] + F[i])/2 for i=1:(length(p)-1))
end

function area_average(b::Boundary, F; dx=0.01, dy=0.01)
    x = range(extrema(b.r)...,step=dx) .+ dx/2
    y = range(extrema(b.z)...,step=dy) .+ dy/2

    A = zero(dx)
    for xx in x, yy in y
        if in_boundary(b,(xx,yy))
            A = A + F(xx,yy)*dx*dy
        end
    end
    return A/area(b)
end

function area_average(b::Boundary, F::AbstractMatrix, x::AbstractRange, y::AbstractRange)
    dx = step(x)
    dy = step(y)
    A = zero(eltype(x))
    @inbounds for ix in eachindex(x), iy in eachindex(y)
        if in_boundary(b, (x[ix],y[iy]))
            A = A + F[ix,iy]*dx*dy
        end
    end
    return A/area(b)
end

function volume(b::Boundary; dx=0.01,dy=0.01)
    x = range(extrema(b.r)...,step=dx) .+ dx/2
    y = range(extrema(b.z)...,step=dy) .+ dy/2

    V = zero(eltype(x))
    for xx in x, yy in y
        if in_boundary(b, (xx,yy))
            V = V + xx*dx*dy
        end
    end
    return 2pi*V
end

function volume_average(b::Boundary, F; dx=0.01, dy=0.01)
    x = range(extrema(b.r)...,step=dx) .+ dx/2
    y = range(extrema(b.z)...,step=dy) .+ dy/2

    A = zero(dx)
    for xx in x, yy in y
        if in_boundary(b, (xx,yy))
            A = A + F(xx,yy)*xx*dx*dy
        end
    end
    return 2pi*A/volume(b)
end

function volume_average(b::Boundary, F::AbstractMatrix, x::AbstractRange, y::AbstractRange)
    dx = step(x)
    dy = step(y)
    A = zero(eltype(x))
    @inbounds for ix in eachindex(x), iy in eachindex(y)
        if in_boundary(b, (x[ix],y[iy]))
            A = A + F[ix,iy]*x[ix]*dx*dy
        end
    end
    return 2pi*A/volume(b)
end


text:
# run this before graphs.jl and rest.jl variants.

#include("../src/GISSMOReader.jl")
import GISSMOReader
import JLD

# get the all GISSMO entries.
tmp = GISSMOReader.getGISSMOentriesall()
GISSMO_entries = GISSMOReader.extractfields(tmp, "entry")
molecule_names = GISSMOReader.extractfields(tmp, "molecule_name")

# save foldder.
base_dir = "/home/roy/Documents/repo/NMRData/input/GISSMO_data"

molar_masses = collect( GISSMOReader.fetchentrymetadata(GISSMO_entries[i])[1] for i = 1:length( GISSMO_entries))

save_path = joinpath(base_dir, "molar_masses.jld")
JLD.save(save_path, "molecule_names", molecule_names,
                    "molar_masses", molar_masses)

# load.
load_path = save_path
dic = JLD.load(load_path)

@assert all(isfinite.(dic["molar_masses"]))


text:
export crossval

crossval(a::Tuple, args...; kargs...) = loocv(a..., args...; kargs...)

function crossval(data, labels, train; predict = predict, n = 4, predictdata = [], shuffle = true, jointmodel = false)
    if !isempty(predictdata) && len(predictdata) != len(data)
        error("len(predictdata)==$(len(predictdata)) is not equal len(data)==$(len(data))")
    end

    ind = shuffle ? randperm(1:len(data)) : collect(1:len(data))
    invind = invperm(ind)

    D = @p part data ind | partition n
    L = @p part labels ind | partition n

    predictdata = isempty(predictdata) ? D : @p part predictdata ind | partition n

    r = cell(n)
    for i in 1:n
        d1, _ = cut(D,i)
        _, d2 = cut(predictdata,i)
        l1, l2 = cut(L,i)
        m = train(flatten(d1),flatten(l1))
        r[i] = predict(m,flatten(d2))
    end
    r = @p flatten r | part invind
    if jointmodel
        return r, @p train data labels
    else
        return r
    end
end


text:
U = Unum{3,5}
x = Ubound(Unum{3,5}(0x0000000000000000, 0x0000000000000000, 0x0001, 0x0000, 0x0000), Unum{3,5}(0x0000000000000001, 0x8000000000000000, 0x0001, 0x0001, 0x0000))
Unums.udivide(x)

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
# This file is part of the TaylorIntegration.jl package; MIT licensed

using TaylorIntegration
using Test
using LinearAlgebra: norm

const _order = 28
const _abstol = 1.0E-20

function pendulum!(dx, x, p, t)
    dx[1] = x[2]
    dx[2] = -sin(x[1])
    nothing
end

g(dx, x, p, t) = x[2]

@testset "Test intersection of surfaces and root-finding: simple pendulum" begin
    t0 = 0.0
    x0 = [1.3, 0.0]
    T = 7.019250311844546

    p = set_variables("ξ", numvars=length(x0), order=2)
    x0N = x0 + p
    ξ = Taylor1(2)
    x01 = x0 + [ξ, ξ]

    #warm-up lap and preliminary tests
    tv, xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, t0, T, _order, _abstol, maxsteps=1)
    @test size(tv) == (2,)
    @test tv[1] == t0
    @test size(xv) == (2,2)
    @test xv[1,:] == x0
    @test size(tvS) == (0,)
    @test size(xvS) == (0,2)
    @test size(gvS) == (0,)

    tvN, xvN, tvSN, xvSN, gvSN = taylorinteg(pendulum!, g, x0N, t0, 3T, _order, _abstol, maxsteps=1)
    @test eltype(tvN) == Float64
    @test eltype(xvN) == TaylorN{Float64}
    @test eltype(tvSN) == TaylorN{Float64}
    @test eltype(xvSN) == TaylorN{Float64}
    @test eltype(gvSN) == TaylorN{Float64}
    @test size(tvN) == (2,)
    @test tvN[1] == t0
    @test size(xvN) == (2,2)
    @test xvN[1,:] == x0N
    @test size(tvSN) == (0,)
    @test size(xvSN) == (0,2)
    @test size(gvSN) == (0,)

    tv1, xv1, tvS1, xvS1, gvS1 = taylorinteg(pendulum!, g, x01, t0, 3T, _order, _abstol, maxsteps=1)
    @test eltype(tv1) == Float64
    @test eltype(xv1) == Taylor1{Float64}
    @test eltype(tvS1) == Taylor1{Float64}
    @test eltype(xvS1) == Taylor1{Float64}
    @test eltype(gvS1) == Taylor1{Float64}
    @test size(tv1) == (2,)
    @test tv1[1] == t0
    @test size(xv1) == (2,2)
    @test xv1[1,:] == x01
    @test size(tvS1) == (0,)
    @test size(xvS1) == (0,2)
    @test size(gvS1) == (0,)

    #testing 0-th order root-finding
    tv, xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, t0, 3T, _order, _abstol, maxsteps=1000)
    @test tv[1] == t0
    @test xv[1,:] == x0
    @test size(tvS) == (5,)
    @test norm(tvS-[T/2,T,3T/2,2T,5T/2],Inf) < 1E-13
    @test norm(gvS,Inf) < eps()

    #testing 0-th order root-finding with time ranges/vectors
    tvr = [t0, T/2, T, 3T/2, 2T, 5T/2, 3T]
    @test_throws AssertionError taylorinteg(pendulum!, g, x0, view(tvr, :), _order, _abstol, maxsteps=1000, eventorder=_order+1)
    xvr, tvSr, xvSr, gvSr = taylorinteg(pendulum!, g, x0, view(tvr, :), _order, _abstol, maxsteps=1000)
    @test xvr[1,:] == x0
    @test size(tvSr) == (5,)
    @test size(tvSr) == size(tvr[2:end-1])
    @test norm(tvSr-[T/2,T,3T/2,2T,5T/2],Inf) < 1E-13
    @test norm(tvr[2:end-1]-tvSr, Inf) < 1E-14
    @test norm(xvr[2:end-1,:]-xvSr, Inf) < 1E-14
    @test norm(gvSr[:]) < eps()
    @test norm(tvS-tvSr, Inf) < 5E-15

    #testing 0-th order root-finding + TaylorN jet transport
    tvN, xvN, tvSN, xvSN, gvSN = taylorinteg(pendulum!, g, x0N, t0, 3T, _order, _abstol, maxsteps=1000)
    @test size(tvSN) == size(tvS)
    @test size(xvSN) == size(xvS)
    @test size(gvSN) == size(gvS)
    @test norm(gvSN[:]) < 1E-15
    @test norm( tvSN()-[T/2,T,3T/2,2T,5T/2], Inf ) < 1E-13
    @test norm( xvSN()-xvS, Inf ) < 1E-14

    #testing 0-th root-finding + Taylor1 jet transport
    tv1, xv1, tvS1, xvS1, gvS1 = taylorinteg(pendulum!, g, x01, t0, 3T, _order, _abstol, maxsteps=1000)
    @test size(tvS1) == size(tvS)
    @test size(xvS1) == size(xvS)
    @test size(gvS1) == size(gvS)
    @test norm(gvS1[:]) < 1E-14
    @test norm( tvS1()-[T/2,T,3T/2,2T,5T/2], Inf ) < 1E-13
    @test norm( xvS1()-xvS, Inf ) < 1E-14

    #testing surface higher order crossing detections and root-finding
    @test_throws AssertionError taylorinteg(pendulum!, g, x0, t0, 3T, _order, _abstol, maxsteps=1000, eventorder=_order+1, newtoniter=2)
    tv, xv, tvS, xvS, gvS = taylorinteg(pendulum!, g, x0, t0, 3T, _order, _abstol, maxsteps=1000, eventorder=2, newtoniter=2)
    @test tv[1] == t0
    @test xv[1,:] == x0
    @test size(tvS) == (5,)
    @test norm(tvS-[T/2,T,3T/2,2T,5T/2],Inf) < 1E-13
    @test norm(gvS[:]) < 1E-15

    #testing higher order root-finding + TaylorN jet transport
    tvN, xvN, tvSN, xvSN, gvSN = taylorinteg(pendulum!, g, x0N, t0, 3T, _order, _abstol, maxsteps=1000, eventorder=2)
    @test size(tvSN) == size(tvS)
    @test size(xvSN) == size(xvS)
    @test size(gvSN) == size(gvS)
    @test norm(gvSN[:]) < 1E-13
    @test norm(tvSN()-[T/2,T,3T/2,2T,5T/2],Inf) < 1E-13
    @test norm(xvSN()-xvS, Inf) < 1E-14

    #testing higher root-finding + Taylor1 jet transport
    tv1, xv1, tvS1, xvS1, gvS1 = taylorinteg(pendulum!, g, x01, t0, 3T, _order, _abstol, maxsteps=1000, eventorder=2)
    @test size(tvS1) == size(tvS)
    @test size(xvS1) == size(xvS)
    @test size(gvS1) == size(gvS)
    @test norm(gvS1[:]) < 1E-14
    @test norm( tvS1()-[T/2,T,3T/2,2T,5T/2], Inf ) < 1E-13
    @test norm( xvS1()-xvS, Inf ) < 1E-14
end


text:
include("euler/euler.jl")
using .Calculus: fibonacci_index, get_fibonacci
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

compute(n::Int)::Int = trunc(Int, get_fibonacci(fibonacci_index(n) + 2) / 2)

compute(4_000_000)

@benchmark compute(4_000_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler