text:
int_rules_7_5_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 7.5.1*u*(a+b*arcsech(c*x))^n =#
    @apply_utils Antiderivative(asech(~(c') * ~x), ~x) => ~x * asech(~c * ~x) + sqrt(1 + ~c * ~x) * sqrt(1 / (1 + ~c * ~x)) * Antiderivative(1 / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ(~c, ~x)
    @apply_utils Antiderivative(acsch(~(c') * ~x), ~x) => ~x * acsch(~c * ~x) + (1 / ~c) * Antiderivative(1 / (~x * sqrt(1 + 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x) <-- FreeQ(~c, ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * asech(~(c') * ~x)) ^ ~n, ~x) => (-1 / ~c) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * sech(~x) * tanh(~x), ~x), ~x, asech(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * acsch(~(c') * ~x)) ^ ~n, ~x) => (-1 / ~c) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * csch(~x) * coth(~x), ~x), ~x, acsch(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * asech(~(c') * ~x)) / ~x, ~x) => -(Subst(Antiderivative((~a + ~b * acosh(~x / ~c)) / ~x, ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * acsch(~(c') * ~x)) / ~x, ~x) => -(Subst(Antiderivative((~a + ~b * asinh(~x / ~c)) / ~x, ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * asech(~(c') * ~x)), ~x) => ((~d * ~x) ^ (~m + 1) * (~a + ~b * asech(~c * ~x))) / (~d * (~m + 1)) + ((~b * sqrt(1 + ~c * ~x)) / (~m + 1)) * sqrt(1 / (1 + ~c * ~x)) * Antiderivative((~d * ~x) ^ ~m / (sqrt(1 - ~c * ~x) * sqrt(1 + ~c * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * acsch(~(c') * ~x)), ~x) => ((~d * ~x) ^ (~m + 1) * (~a + ~b * acsch(~c * ~x))) / (~d * (~m + 1)) + ((~b * ~d) / (~c * (~m + 1))) * Antiderivative((~d * ~x) ^ (~m - 1) / sqrt(1 + 1 / ((~c) ^ 2 * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * asech(~(c') * ~x)) ^ ~n, ~x) => (-1 / (~c) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * sech(~x) ^ (~m + 1) * tanh(~x), ~x), ~x, asech(~c * ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ(~n) && (IntegerQ(~m) && (GtQ(~n, 0) || LtQ(~m, -1))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * acsch(~(c') * ~x)) ^ ~n, ~x) => (-1 / (~c) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * csch(~x) ^ (~m + 1) * coth(~x), ~x), ~x, acsch(~c * ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IntegerQ(~n) && (IntegerQ(~m) && (GtQ(~n, 0) || LtQ(~m, -1))))
    @apply_utils Antiderivative((~(a') + ~(b') * asech(~(c') * ~x)) / (~(d') + ~(e') * ~x), ~x) => (((((~a + ~b * asech(~c * ~x)) * log(1 + (~e - sqrt(-((~c) ^ 2) * (~d) ^ 2 + (~e) ^ 2)) / (~c * ~d * ℯ ^ asech(~c * ~x)))) / ~e + ((~a + ~b * asech(~c * ~x)) * log(1 + (~e + sqrt(-((~c) ^ 2) * (~d) ^ 2 + (~e) ^ 2)) / (~c * ~d * ℯ ^ asech(~c * ~x)))) / ~e) - ((~a + ~b * asech(~c * ~x)) * log(1 + 1 / ℯ ^ (2 * asech(~c * ~x)))) / ~e) + (~b / ~e) * Antiderivative((sqrt((1 - ~c * ~x) / (1 + ~c * ~x)) * log(1 + (~e - sqrt(-((~c) ^ 2) * (~d) ^ 2 + (~e) ^ 2)) / (~c * ~d * ℯ ^ asech(~c * ~x)))) / (~x * (1 - ~c * ~x)), ~x) + (~b / ~e) * Antiderivative((sqrt((1 - ~c * ~x) / (1 + ~c * ~x)) * log(1 + (~e + sqrt(-((~c) ^ 2) * (~d) ^ 2 + (~e) ^ 2)) / (~c * ~d * ℯ ^ asech(~c * ~x)))) / (~x * (1 - ~c * ~x)), ~x)) - (~b / ~e) * Antiderivative((sqrt((1 - ~c * ~x) / (1 + ~c * ~x)) * log(1 + 1 / ℯ ^ (2 * asech(~c * ~x)))) / (~x * (1 - ~c * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * asech(~(c') * ~x)), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * asech(~c * ~x))) / (~e * (~m + 1)) + ((~b * sqrt(1 + ~c * ~x)) / (~e * (~m + 1))) * sqrt(1 / (1 + ~c * ~x)) * Antiderivative((~d + ~e * ~x) ^ (~m + 1) / (~x * sqrt(1 - (~c) ^ 2 * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(a') + ~(b') * acsch(~(c') * ~x)) / (~(d') + ~(e') * ~x), ~x) => (((((~a + ~b * acsch(~c * ~x)) * log(1 - ((~e - sqrt((~c) ^ 2 * (~d) ^ 2 + (~e) ^ 2)) * ℯ ^ acsch(~c * ~x)) / (~c * ~d))) / ~e + ((~a + ~b * acsch(~c * ~x)) * log(1 - ((~e + sqrt((~c) ^ 2 * (~d) ^ 2 + (~e) ^ 2)) * ℯ ^ acsch(~c * ~x)) / (~c * ~d))) / ~e) - ((~a + ~b * acsch(~c * ~x)) * log(1 - ℯ ^ (2 * acsch(~c * ~x)))) / ~e) + (~b / (~c * ~e)) * Antiderivative(log(1 - ((~e - sqrt((~c) ^ 2 * (~d) ^ 2 + (~e) ^ 2)) * ℯ ^ acsch(~c * ~x)) / (~c * ~d)) / ((~x) ^ 2 * sqrt(1 + 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x) + (~b / (~c * ~e)) * Antiderivative(log(1 - ((~e + sqrt((~c) ^ 2 * (~d) ^ 2 + (~e) ^ 2)) * ℯ ^ acsch(~c * ~x)) / (~c * ~d)) / ((~x) ^ 2 * sqrt(1 + 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x)) - (~b / (~c * ~e)) * Antiderivative(log(1 - ℯ ^ (2 * acsch(~c * ~x))) / ((~x) ^ 2 * sqrt(1 + 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * acsch(~(c') * ~x)), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * acsch(~c * ~x))) / (~e * (~m + 1)) + (~b / (~c * ~e * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 1) / ((~x) ^ 2 * sqrt(1 + 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asech(~(c') * ~x)), ~x) => With([u = IntHide((~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * asech(~c * ~x), u, ~x) + ~b * sqrt(1 + ~c * ~x) * sqrt(1 / (1 + ~c * ~x)) * Antiderivative(SimplifyIntegrand(u / (~x * sqrt(1 - ~c * ~x) * sqrt(1 + ~c * ~x)), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~p, 0) || ILtQ(~p + 1 / 2, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acsch(~(c') * ~x)), ~x) => With([u = IntHide((~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * acsch(~c * ~x), u, ~x) - ((~b * ~c * ~x) / sqrt(-((~c) ^ 2) * (~x) ^ 2)) * Antiderivative(SimplifyIntegrand(u / (~x * sqrt(-1 - (~c) ^ 2 * (~x) ^ 2)), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~p, 0) || ILtQ(~p + 1 / 2, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asech(~(c') * ~x)) ^ ~(n'), ~x) => -(Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * acosh(~x / ~c)) ^ ~n) / (~x) ^ (2 * (~p + 1)), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acsch(~(c') * ~x)) ^ ~(n'), ~x) => -(Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * asinh(~x / ~c)) ^ ~n) / (~x) ^ (2 * (~p + 1)), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asech(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt((~x) ^ 2)) / ~x) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * acosh(~x / ~c)) ^ ~n) / (~x) ^ (2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~p + 1 / 2) && (GtQ(~e, 0) && LtQ(~d, 0)))))
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acsch(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt((~x) ^ 2)) / ~x) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * asinh(~x / ~c)) ^ ~n) / (~x) ^ (2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ(~e - (~c) ^ 2 * ~d, 0) && (IntegerQ(~p + 1 / 2) && (GtQ(~e, 0) && LtQ(~d, 0)))))
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asech(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt(~d + ~e * (~x) ^ 2)) / (~x * sqrt(~e + ~d / (~x) ^ 2))) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * acosh(~x / ~c)) ^ ~n) / (~x) ^ (2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~p + 1 / 2) && Not(GtQ(~e, 0) && LtQ(~d, 0)))))
    @apply_utils Antiderivative((~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acsch(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt(~d + ~e * (~x) ^ 2)) / (~x * sqrt(~e + ~d / (~x) ^ 2))) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * asinh(~x / ~c)) ^ ~n) / (~x) ^ (2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ(~e - (~c) ^ 2 * ~d, 0) && (IntegerQ(~p + 1 / 2) && Not(GtQ(~e, 0) && LtQ(~d, 0)))))
    @apply_utils Antiderivative(~x * (~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asech(~(c') * ~x)), ~x) => ((~d + ~e * (~x) ^ 2) ^ (~p + 1) * (~a + ~b * asech(~c * ~x))) / (2 * ~e * (~p + 1)) + ((~b * sqrt(1 + ~c * ~x)) / (2 * ~e * (~p + 1))) * sqrt(1 / (1 + ~c * ~x)) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~p + 1) / (~x * sqrt(1 - ~c * ~x) * sqrt(1 + ~c * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && NeQ(~p, -1)
    @apply_utils Antiderivative(~x * (~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acsch(~(c') * ~x)), ~x) => ((~d + ~e * (~x) ^ 2) ^ (~p + 1) * (~a + ~b * acsch(~c * ~x))) / (2 * ~e * (~p + 1)) - ((~b * ~c * ~x) / (2 * ~e * (~p + 1) * sqrt(-((~c) ^ 2) * (~x) ^ 2))) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~p + 1) / (~x * sqrt(-1 - (~c) ^ 2 * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && NeQ(~p, -1)
    @apply_utils Antiderivative((~(f') * ~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asech(~(c') * ~x)), ~x) => With([u = IntHide((~f * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * asech(~c * ~x), u, ~x) + ~b * sqrt(1 + ~c * ~x) * sqrt(1 / (1 + ~c * ~x)) * Antiderivative(SimplifyIntegrand(u / (~x * sqrt(1 - ~c * ~x) * sqrt(1 + ~c * ~x)), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~p], ~x) && (IGtQ(~p, 0) && Not(ILtQ((~m - 1) / 2, 0) && GtQ(~m + 2 * ~p + 3, 0)) || (IGtQ((~m + 1) / 2, 0) && Not(ILtQ(~p, 0) && GtQ(~m + 2 * ~p + 3, 0)) || ILtQ((~m + 2 * ~p + 1) / 2, 0) && Not(ILtQ((~m - 1) / 2, 0))))
    @apply_utils Antiderivative((~(f') * ~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acsch(~(c') * ~x)), ~x) => With([u = IntHide((~f * ~x) ^ ~m * (~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * acsch(~c * ~x), u, ~x) - ((~b * ~c * ~x) / sqrt(-((~c) ^ 2) * (~x) ^ 2)) * Antiderivative(SimplifyIntegrand(u / (~x * sqrt(-1 - (~c) ^ 2 * (~x) ^ 2)), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~p], ~x) && (IGtQ(~p, 0) && Not(ILtQ((~m - 1) / 2, 0) && GtQ(~m + 2 * ~p + 3, 0)) || (IGtQ((~m + 1) / 2, 0) && Not(ILtQ(~p, 0) && GtQ(~m + 2 * ~p + 3, 0)) || ILtQ((~m + 2 * ~p + 1) / 2, 0) && Not(ILtQ((~m - 1) / 2, 0))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asech(~(c') * ~x)) ^ ~(n'), ~x) => -(Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * acosh(~x / ~c)) ^ ~n) / (~x) ^ (~m + 2 * (~p + 1)), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && IntegersQ(~m, ~p))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acsch(~(c') * ~x)) ^ ~(n'), ~x) => -(Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * asinh(~x / ~c)) ^ ~n) / (~x) ^ (~m + 2 * (~p + 1)), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && IntegersQ(~m, ~p))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asech(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt((~x) ^ 2)) / ~x) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * acosh(~x / ~c)) ^ ~n) / (~x) ^ (~m + 2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~m) && (IntegerQ(~p + 1 / 2) && (GtQ(~e, 0) && LtQ(~d, 0))))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acsch(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt((~x) ^ 2)) / ~x) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * asinh(~x / ~c)) ^ ~n) / (~x) ^ (~m + 2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ(~e - (~c) ^ 2 * ~d, 0) && (IntegerQ(~m) && (IntegerQ(~p + 1 / 2) && (GtQ(~e, 0) && LtQ(~d, 0))))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * asech(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt(~d + ~e * (~x) ^ 2)) / (~x * sqrt(~e + ~d / (~x) ^ 2))) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * acosh(~x / ~c)) ^ ~n) / (~x) ^ (~m + 2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~m) && (IntegerQ(~p + 1 / 2) && Not(GtQ(~e, 0) && LtQ(~d, 0))))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(d') + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acsch(~(c') * ~x)) ^ ~(n'), ~x) => (-(sqrt(~d + ~e * (~x) ^ 2)) / (~x * sqrt(~e + ~d / (~x) ^ 2))) * Subst(Antiderivative(((~e + ~d * (~x) ^ 2) ^ ~p * (~a + ~b * asinh(~x / ~c)) ^ ~n) / (~x) ^ (~m + 2 * (~p + 1)), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (IGtQ(~n, 0) && (EqQ(~e - (~c) ^ 2 * ~d, 0) && (IntegerQ(~m) && (IntegerQ(~p + 1 / 2) && Not(GtQ(~e, 0) && LtQ(~d, 0))))))
    @apply_utils Antiderivative(~u * (~(a') + ~(b') * asech(~(c') * ~x)), ~x) => With([v = IntHide(~u, ~x)], Dist(~a + ~b * asech(~c * ~x), v, ~x) + ((~b * sqrt(1 - (~c) ^ 2 * (~x) ^ 2)) / (~c * ~x * sqrt(-1 + 1 / (~c * ~x)) * sqrt(1 + 1 / (~c * ~x)))) * Antiderivative(SimplifyIntegrand(v / (~x * sqrt(1 - (~c) ^ 2 * (~x) ^ 2)), ~x), ~x) <-- InverseFunctionFreeQ(v, ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~u * (~(a') + ~(b') * acsch(~(c') * ~x)), ~x) => With([v = IntHide(~u, ~x)], Dist(~a + ~b * acsch(~c * ~x), v, ~x) + (~b / ~c) * Antiderivative(SimplifyIntegrand(v / ((~x) ^ 2 * sqrt(1 + 1 / ((~c) ^ 2 * (~x) ^ 2))), ~x), ~x) <-- InverseFunctionFreeQ(v, ~x)) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * asech(~(c') * ~x)) ^ ~(n'), ~x) => Unintegrable(~u * (~a + ~b * asech(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
    @apply_utils Antiderivative(~(u') * (~(a') + ~(b') * acsch(~(c') * ~x)) ^ ~(n'), ~x) => Unintegrable(~u * (~a + ~b * acsch(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
end

max_stars_repo_path: src/rules/7 Inverse hyperbolic functions/7.5 Inverse hyperbolic secant/7.5.1 u (a+b arcsech(c x))^n.jl


text:
#################################################################
######################## BELIEF #################################
#################################################################
# For tools supporting belief updates see POMDPToolbox.jl

"""
Abstract type for an object that defines how the belief should be updated

A belief is a general construct that represents the knowledge an agent has
about the state of the system. This can be a probability distribution, an
action observation history or a more general representation.
"""
abstract type Updater end

"""
    update(updater::Updater, belief_old, action, observation)

Return a new instance of an updated belief given `belief_old` and the latest action and observation.
"""
function update end

"""
    initialize_belief(updater::Updater,
                         state_distribution::Any)
    initialize_belief(updater::Updater, belief::Any)

Returns a belief that can be updated using `updater` that has similar
distribution to `state_distribution` or `belief`.

The conversion may be lossy. This function is also idempotent, i.e. there is a
default implementation that passes the belief through when it is already the
correct type: `initialize_belief(updater::Updater, belief) = belief`
"""
function initialize_belief end

# default implementation if the input is the same type as the output
initialize_belief(updater::Updater, belief) = belief

ext: jl


text:
"""
    ascertain_pose(camera::CameraModel, 𝐞₁, 𝐞₂, 𝐞₃)

Determines the rotation and translation of the camera with respect to the
origin of a world coordinate system with axes {𝐞₁, 𝐞₂, 𝐞₃}.
"""
function ascertain_pose(camera::CameraModel, 𝐞₁, 𝐞₂, 𝐞₃)
    𝐭 = camera.𝐜
    𝐑 = inv(hcat(𝐞₁, 𝐞₂, 𝐞₃)) * hcat(camera.𝐞₁, camera.𝐞₂, camera.𝐞₃)
    𝐑, 𝐭
end

function obtain_intrinsics(camera::CameraModel, system::RasterSystem)
    # Convention I: 𝐭 = Point(-camera.image_width / 2, camera.image_height / 2)
    𝐭 = Point(-camera.image_width / 2, -camera.image_height / 2)
    𝐑 = inv(hcat(system.𝐞₁, system.𝐞₂)) * hcat(camera.𝐞₁′, camera.𝐞₂′)
    f = camera.focal_length
    𝐊 = SMatrix{3,3,Float64,9}(f, 0.0, 0.0, 0.0, f, 0.0, 0.0, 0.0 , 1)
    𝐊′ =vcat(hcat(𝐑', -𝐑'*𝐭), SMatrix{1,3,Float64}(0,0,1) )*𝐊
end

function obtain_intrinsics(camera::CameraModel, system::CartesianSystem)
    # Convention I: 𝐭 = Point(-camera.image_width / 2, camera.image_height / 2)
    𝐭 = Point(-camera.image_width / 2, camera.image_height / 2)
    𝐑 = inv(hcat(system.𝐞₁, system.𝐞₂)) * hcat(camera.𝐞₁′, camera.𝐞₂′)
    #𝐑 = vcat(hcat(𝐑₂₂, 0), SVector(0,0,0))
    #𝐑 =vcat(hcat(𝐑₂₂, SVector(0,0)), SMatrix{1,3,Float64}(0,0,1))
    f = camera.focal_length
    𝐊 = SMatrix{3,3,Float64,9}(f, 0.0, 0.0, 0.0, f, 0.0, 0.0, 0.0 , 1)
    𝐊′ =vcat(hcat(𝐑', -𝐑'*𝐭), SMatrix{1,3,Float64}(0,0,1))*𝐊
end

function obtain_intrinsics(camera::CameraModel, system::OpticalSystem)
    𝐭 = Point(0, 0)
    𝐑 = inv(hcat(system.𝐞₁, system.𝐞₂)) * hcat(camera.𝐞₁′, camera.𝐞₂′)
    @show 𝐑
    f = camera.focal_length
    𝐊 = SMatrix{3,3,Float64,9}(f, 0.0, 0.0, 0.0, f, 0.0, 0.0, 0.0 , 1)
    𝐊′ =vcat(hcat(𝐑', -𝐑'*𝐭), SMatrix{1,3,Float64}(0,0,1))*𝐊
end


text:
struct UnsupportedArgumentType <: Exception
    key::String
end

struct ArbArgTypes
    supported::Dict{String,DataType}
    unsupported::Set{String}
    supported_reversed::Dict{DataType,String}
end

function ArbArgTypes(supported, unsupported)
    supported_reversed = Dict(value => key for (key, value) in supported)
    return ArbArgTypes(supported, unsupported, supported_reversed)
end

function Base.getindex(arbargtypes::ArbArgTypes, key::AbstractString)
    haskey(arbargtypes.supported, key) && return arbargtypes.supported[key]
    key in arbargtypes.unsupported && throw(UnsupportedArgumentType(key))
    throw(KeyError(key))
end

const arbargtypes = ArbArgTypes(
    Dict{String,DataType}(
        "void" => Cvoid,
        "void *" => Ptr{Cvoid},
        "int" => Cint,
        "slong" => Int,
        "ulong" => UInt,
        "double" => Cdouble,
        "arf_t" => Arf,
        "arb_t" => Arb,
        "acb_t" => Acb,
        "mag_t" => Mag,
        "arb_srcptr" => ArbVector,
        "arb_ptr" => ArbVector,
        "acb_srcptr" => AcbVector,
        "acb_ptr" => AcbVector,
        "arb_poly_t" => ArbPoly,
        "acb_poly_t" => AcbPoly,
        "arb_mat_t" => ArbMatrix,
        "acb_mat_t" => AcbMatrix,
        "arf_rnd_t" => arb_rnd,
        "mpfr_t" => BigFloat,
        "mpfr_rnd_t" => Base.MPFR.MPFRRoundingMode,
        "mpz_t" => BigInt,
        "char *" => Cstring,
        "slong *" => Vector{Int},
        "ulong *" => Vector{UInt},
    ),
    Set(["FILE *", "fmpr_t", "fmpr_rnd_t", "flint_rand_t", "bool_mat_t"]),
    Dict{DataType,String}(
        Cvoid => "void",
        Ptr{Cvoid} => "void *",
        Cint => "int",
        Int => "slong",
        UInt => "ulong",
        Cdouble => "double",
        Arf => "arf_t",
        Arb => "arb_t",
        Acb => "acb_t",
        Mag => "mag_t",
        ArbVector => "arb_ptr",
        AcbVector => "acb_ptr",
        ArbPoly => "arb_poly_t",
        AcbPoly => "acb_poly_t",
        ArbMatrix => "arb_mat_t",
        AcbMatrix => "acb_mat_t",
        arb_rnd => "arf_rnd_t",
        BigFloat => "mpfr_t",
        Base.MPFR.MPFRRoundingMode => "mpfr_rnd_t",
        BigInt => "mpz_t",
        Cstring => "char *",
        Vector{Int} => "slong *",
        Vector{UInt} => "ulong *",
    ),
)

struct Carg{ArgT}
    name::Symbol
    isconst::Bool
end

Carg{T}(n::AbstractString, isconst::Bool) where {T} = Carg{T}(Symbol(n), isconst)

Base.:(==)(a::Carg{T}, b::Carg{S}) where {T,S} =
    T == S && name(a) == name(b) && isconst(a) == isconst(b)

function Carg(str)
    m = match(r"(?<const>const)?\s*(?<type>\w+(\s\*)?)\s+(?<name>\w+)", str)
    isnothing(m) && throw(ArgumentError("string doesn't match c-argument pattern"))
    cnst =
        !isnothing(m[:const]) ||
        (!isnothing(m[:type]) && (m[:type] == "arb_srcptr" || m[:type] == "acb_srcptr"))
    return Carg{arbargtypes[m[:type]]}(m[:name], cnst)
end

name(ca::Carg) = ca.name
isconst(ca::Carg) = ca.isconst

rawtype(::Carg{T}) where {T} = T

jltype(ca::Carg) = rawtype(ca)
jltype(ca::Carg{Cint}) = Integer
jltype(ca::Carg{Int}) = Integer
jltype(ca::Carg{UInt}) = Unsigned
jltype(ca::Carg{Cdouble}) = Base.GMP.CdoubleMax
jltype(ca::Carg{arb_rnd}) = Union{arb_rnd,RoundingMode}
jltype(ca::Carg{Base.MPFR.MPFRRoundingMode}) =
    Union{Base.MPFR.MPFRRoundingMode,RoundingMode}
jltype(ca::Carg{Cstring}) = AbstractString
jltype(ca::Carg{Vector{Int}}) = Vector{<:Integer}
jltype(ca::Carg{Vector{UInt}}) = Vector{<:Unsigned}
jltype(::Carg{Mag}) = MagLike
jltype(::Carg{Arf}) = ArfLike
jltype(::Carg{Arb}) = ArbLike
jltype(::Carg{Acb}) = AcbLike
jltype(::Carg{ArbVector}) = ArbVectorLike
jltype(::Carg{AcbVector}) = AcbVectorLike
jltype(::Carg{ArbMatrix}) = ArbMatrixLike
jltype(::Carg{AcbMatrix}) = AcbMatrixLike
jltype(::Carg{ArbPoly}) = ArbPolyLike
jltype(::Carg{AcbPoly}) = AcbPolyLike

ctype(ca::Carg) = rawtype(ca)
ctype(::Carg{T}) where {T<:Union{ArbVector,arb_vec_struct}} = Ptr{arb_struct}
ctype(::Carg{T}) where {T<:Union{AcbVector,acb_vec_struct}} = Ptr{acb_struct}
ctype(::Carg{T}) where {T<:Union{Mag,Arf,Arb,Acb,ArbPoly,AcbPoly,ArbMatrix,AcbMatrix}} =
    Ref{cstructtype(T)}
ctype(::Carg{T}) where {T<:Union{BigFloat,BigInt}} = Ref{T}
ctype(::Carg{Vector{T}}) where {T} = Ref{T}

struct Arbfunction{ReturnT}
    fname::String
    args::Vector{Carg}
end

function Arbfunction(str)
    m = match(r"(?<returntype>\w+(\s\*)?)\s+(?<arbfunction>[\w_]+)\((?<args>.*)\)", str)
    isnothing(m) &&
        throw(ArgumentError("string doesn't match arblib function signature pattern"))

    args = Carg.(strip.(split(m[:args], ",")))

    return Arbfunction{arbargtypes[m[:returntype]]}(m[:arbfunction], args)
end

function jlfname(
    arbfname,
    prefixes = ("arf", "arb", "acb", "mag", "mat", "vec", "poly", "scalar"),
    suffixes = ("si", "ui", "d", "mag", "arf", "arb", "acb", "mpz", "mpfr", "str");
    inplace = false,
)
    strs = filter(!isempty, split(arbfname, "_"))
    k = findfirst(s -> s ∉ prefixes, strs)
    l = findfirst(s -> s ∉ suffixes, reverse(strs))
    fname = join(strs[k:end-l+1], "_")
    return inplace ? Symbol(fname, "!") : Symbol(fname)
end

arbfname(af::Arbfunction) = af.fname
returntype(af::Arbfunction{ReturnT}) where {ReturnT} = ReturnT
arguments(af::Arbfunction) = af.args

function inplace(af::Arbfunction)
    firstarg = first(arguments(af))
    return !isconst(firstarg) &&
           (ctype(firstarg) <: Ref || ctype(firstarg) <: AbstractArray)
end

function jlfname(
    af::Arbfunction,
    prefixes = ("arf", "arb", "acb", "mag", "mat", "vec", "poly", "scalar"),
    suffixes = ("si", "ui", "d", "mag", "arf", "arb", "acb", "mpz", "mpfr", "str");
    inplace = inplace(af),
)
    return jlfname(arbfname(af), prefixes, suffixes, inplace = inplace)
end

function is_length_argument(carg, prev_carg, len_keywords)
    (startswith(string(name(carg)), "len") || carg == Carg{Int}(:n, false)) &&
        rawtype(prev_carg) ∈ (ArbVector, AcbVector) &&
        name(carg) ∉ len_keywords
end
function extract_length_argument!(kwargs, len_keywords, carg, prev_carg)
    vec_name = name(prev_carg)
    push!(kwargs, Expr(:kw, :($(name(carg))::Integer), :(length($vec_name))))
    push!(len_keywords, name(carg))
end
function ispredicate(af::Arbfunction)
    return isconst(first(arguments(af))) &&
           returntype(af) == Cint &&
           (
               any(s -> startswith(string(jlfname(af)), s), ("is_",)) ||
               any(
                   s -> occursin(s, string(jlfname(af))),
                   ("_is_", "contains", "can_", "check_", "validate_"),
               ) ||
               any(==(jlfname(af)), (:eq, :ne, :lt, :le, :gt, :ge, :overlaps, :equal))
           )
end

function jlargs(af::Arbfunction; argument_detection::Bool = true)
    cargs = arguments(af)

    jl_arg_names_types = Tuple{Symbol,Any}[]
    kwargs = Expr[]

    prec_kwarg = false
    rnd_kwarg = false
    len_keywords = Set{Symbol}()
    for (i, carg) in enumerate(cargs)
        if !argument_detection
            push!(jl_arg_names_types, (name(carg), jltype(carg)))
            continue
        end

        # Automatic detection of precision argument
        if carg == Carg{Int}(:prec, false)
            @assert !prec_kwarg
            prec_kwarg = true

            # If the first argument has a precision,
            # then use this otherwise make it a mandatory kwarg
            if rawtype(cargs[1]) <: ArbTypes && rawtype(cargs[1]) != Mag
                push!(kwargs, Expr(:kw, :(prec::Integer), :(_precision($(name(cargs[1]))))))
            else
                push!(kwargs, :(prec::Integer))
            end

            # Automatic detection of rounding mode argument
        elseif carg == Carg{arb_rnd}(:rnd, false)
            @assert !rnd_kwarg
            rnd_kwarg = true

            push!(
                kwargs,
                Expr(:kw, :(rnd::Union{$(arb_rnd),RoundingMode}), :(RoundNearest)),
            )
        elseif carg == Carg{Base.MPFR.MPFRRoundingMode}(:rnd, false)
            @assert !rnd_kwarg
            rnd_kwarg = true
            push!(
                kwargs,
                Expr(
                    :kw,
                    :(rnd::Union{Base.MPFR.MPFRRoundingMode,RoundingMode}),
                    :(RoundNearest),
                ),
            )
            # Automatic detection of length arguments for vectors
        elseif i > 1 && is_length_argument(carg, cargs[i-1], len_keywords)
            extract_length_argument!(kwargs, len_keywords, carg, cargs[i-1])
        else
            push!(jl_arg_names_types, (name(carg), jltype(carg)))
        end
    end

    args = [:($a::$T) for (a, T) in jl_arg_names_types]

    return args, kwargs
end

function arbsignature(af::Arbfunction)
    creturnT = arbargtypes.supported_reversed[returntype(af)]
    args = arguments(af)

    arg_consts = isconst.(args)
    arg_ctypes = [arbargtypes.supported_reversed[rawtype(arg)] for arg in args]
    arg_names = name.(args)

    c_args = join(
        [
            ifelse(
                isconst && (type == "arb_ptr" || type == "acb_ptr"),
                "$(split(type, "_")[1])_srcptr $name",
                ifelse(isconst, "const ", "") * "$type $name",
            ) for (isconst, type, name) in zip(arg_consts, arg_ctypes, arg_names)
        ],
        ", ",
    )

    "$creturnT $(arbfname(af))($c_args)"
end

function jlcode(af::Arbfunction, jl_fname = jlfname(af))
    jl_args, jl_kwargs = jlargs(af; argument_detection = true)
    jl_full_args, _ = jlargs(af; argument_detection = false)

    returnT = returntype(af)
    cargs = arguments(af)

    func_full_args = :(
        function $jl_fname($(jl_full_args...))
            __ret = ccall(
                Arblib.@libarb($(arbfname(af))),
                $returnT,
                $(Expr(:tuple, ctype.(cargs)...)),
                $(name.(cargs)...),
            )
            $(
                if returnT === Nothing && inplace(af)
                    name(first(arguments(af)))
                elseif ispredicate(af)
                    :(!iszero(__ret))
                else
                    :__ret
                end
            )
        end
    )

    if isempty(jl_kwargs)
        func_full_args
    else
        quote
            $func_full_args
            $jl_fname($(jl_args...); $(jl_kwargs...)) = $jl_fname($(name.(cargs)...))
        end
    end
end

macro arbcall_str(str)
    af = Arbfunction(str)
    return esc(jlcode(af))
end


text:
module CbiomesProcessing

using YAML, FortranFiles, JLD, Dates
using Distributed, DistributedArrays
using SparseArrays, MeshArrays
using MITgcmTools

include("utilities.jl")
include("fileloops.jl")
include("examples.jl")

export StartWorkers, TaskDriver
export cbioproc_example1, cbioproc_example2

#export AverageYear

end # module


text:
module TestKernels

using CUDA
using Test
using UnionArrays

function kernel_test_shmem!(ys)
    iglobal = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    T = Union{Nothing,Float32}
    S = UnionArrays.buffereltypefor(T)
    data = @cuDynamicSharedMem(S, (blockDim().x,))
    typeids = @cuDynamicSharedMem(UInt8, (blockDim().x,), sizeof(data))
    xs = UnionVector(T, data, typeids)
    if isodd(threadIdx().x)
        xs[threadIdx().x] = Float32(iglobal)
    else
        xs[threadIdx().x] = nothing
    end
    sync_threads()
    i = mod1(threadIdx().x + 1, blockDim().x)
    x = xs[i]
    sync_threads()
    xs[threadIdx().x] = x
    sync_threads()
    ys[iglobal] = xs[i]
    return
end

function host_test_shmem()
    threads = 8
    blocks = 3
    ys = UnionVector(undef, CuVector, Union{Float32, Nothing}, threads * blocks)
    fill!(ys, 0)
    shmem = (sizeof(UnionArrays.buffereltypeof(ys)) + sizeof(UInt8)) * threads
    CUDA.@sync @cuda threads=threads blocks=blocks shmem=shmem kernel_test_shmem!(ys)
    return ys
end

function test_shmem()
    xs = collect(Union{Nothing,Float32}, 1:8 * 3)
    xs[2:2:end] .= nothing
    for block in Iterators.partition(eachindex(xs), 8)
        xs[block] = circshift(xs[block], -2)
    end
    @test collect(host_test_shmem()) == xs
end

should_test_module() = VERSION >= v"1.6-"

end  # module


text:
@testset "GridIndex" begin
    Random.seed!(42)
    points = rand(SVector{3, Float64}, 1_000)

    ps = accelerate(points, GridIndex; spacing = 0.1)

    for p in points
    	radius = rand()
    	@test count(in(Sphere(p, radius)), ps) == count(in(Sphere(p, radius)), points)
    	@test issetequal(findall(in(Sphere(p, radius)), ps), findall(x -> x in Sphere(p, radius) #= Base bug... wants `length` of `Sphere`? =#, points))
    	@test issetequal(filter(in(Sphere(p, radius)), ps), filter(in(Sphere(p, radius)), points))
    end

    ps2 = accelerate!(points, GridIndex; spacing = 0.1)

    for p in points
    	radius = rand()
    	@test count(in(Sphere(p, radius)), ps2) == count(in(Sphere(p, radius)), points)
    	@test issetequal(findall(in(Sphere(p, radius)), ps2), findall(x -> x in Sphere(p, radius) #= Base bug... wants `length` of `Sphere`? =#, points))
    	@test issetequal(filter(in(Sphere(p, radius)), ps2), filter(in(Sphere(p, radius)), points))
    end
end


text:
# # Entries
#
# The key ingredient for data representations are `AxisEntries`.
#
# ## The `AxisEntries` type
#
# An `AxisEntries` object is
# made of three components:
# - axis,
# - entries,
# - scales.

using AlgebraOfGraphics, CairoMakie
using AlgebraOfGraphics: CategoricalScale, fitscale
using AlgebraOfGraphics: arguments, namedarguments
resolution = (600, 600)
fig = Figure(; resolution)
N = 11
rg = range(1, 2, length=N)
markerpalette = [:circle, :utriangle, :dtriangle, :rect]
ae = AxisEntries(
    Axis(fig[1, 1]),
    [
        Entry(
            plottype=Scatter,
            positional=arguments((rg, cosh.(rg))),
            named=namedarguments((color=1:N, marker=fill("b", N)));
            attributes=Dict(:markersize => 15)
        ),
        Entry(
            plottype=Scatter,
            positional=arguments((rg, sinh.(rg))),
            named=namedarguments((color=1:N, marker=fill("c", N)));
            attributes=Dict(:markersize => 15)
        ),
    ],
    Dict(
        :marker => fitscale(CategoricalScale(["a", "b", "c"], markerpalette, "class")),
    ), # scales
)
plot!(ae)
fig


text:
function readme(t::Template)
    name = t.pkg_name
    org = t.org
    bors_id = t.bors.id
    return "# $name

$(t.description)

| **Documentation**    | [![dev][docs-dev-img]][docs-dev-url]          |
|----------------------|-----------------------------------------------|
| **Docs Build**       | [![docs build][docs-bld-img]][docs-bld-url]   |
| **Azure Build**      | [![azure][azure-img]][azure-url]              |
| **Code Coverage**    | [![codecov][codecov-img]][codecov-url]        |
| **Bors**             | [![Bors enabled][bors-img]][bors-url]         |

[docs-bld-img]: https://github.com/$(org)/$name/workflows/Documentation/badge.svg
[docs-bld-url]: https://github.com/$(org)/$name/actions?query=workflow%3ADocumentation

[docs-dev-img]: https://img.shields.io/badge/docs-dev-blue.svg
[docs-dev-url]: https://$(org).github.io/$name/dev/

[azure-img]: https://dev.azure.com/$(org)/$name/_apis/build/status/$(org).$name?branchName=master
[azure-url]: https://dev.azure.com/$(org)/$name/_build/latest?definitionId=1&branchName=master

[codecov-img]: https://codecov.io/gh/$(org)/$name/branch/master/graph/badge.svg
[codecov-url]: https://codecov.io/gh/$(org)/$name

[bors-img]: https://bors.tech/images/badge_small.svg
[bors-url]: https://app.bors.tech/repositories/$bors_id

"

end

max_stars_repo_name: charleskawczynski/PkgTemplates

max_issues_repo_name: charleskawczynski/PkgTemplates

max_forks_repo_name: charleskawczynski/PkgTemplates


text:
# https://github.com/Evizero/UnicodePlots.jl

# don't warn on unsupported... there's just too many warnings!!
warnOnUnsupported_args(::UnicodePlotsBackend, plotattributes::KW) = nothing

# --------------------------------------------------------------------------------------

const _canvas_type = Ref(:auto)

function _canvas_map()
    KW(
        :braille => UnicodePlots.BrailleCanvas,
        :ascii => UnicodePlots.AsciiCanvas,
        :block => UnicodePlots.BlockCanvas,
        :dot => UnicodePlots.DotCanvas,
        :density => UnicodePlots.DensityCanvas,
    )
end

# do all the magic here... build it all at once, since we need to know about all the series at the very beginning
function rebuildUnicodePlot!(plt::Plot, width, height)
    plt.o = []

    for sp in plt.subplots
        xaxis = sp[:xaxis]
        yaxis = sp[:yaxis]
        xlim =  axis_limits(xaxis)
        ylim =  axis_limits(yaxis)

        # make vectors
        xlim = [xlim[1], xlim[2]]
        ylim = [ylim[1], ylim[2]]

        # we set x/y to have a single point, since we need to create the plot with some data.
        # since this point is at the bottom left corner of the plot, it shouldn't actually be shown
        x = Float64[xlim[1]]
        y = Float64[ylim[1]]

        # create a plot window with xlim/ylim set, but the X/Y vectors are outside the bounds
        ct = _canvas_type[]
        canvas_type = if ct == :auto
            isijulia() ? UnicodePlots.AsciiCanvas : UnicodePlots.BrailleCanvas
        else
            _canvas_map()[ct]
        end

        # special handling for spy
        if length(sp.series_list) == 1
            series = sp.series_list[1]
            if series[:seriestype] == :spy
                push!(plt.o, UnicodePlots.spy(
                    series[:z].surf,
                    width = width,
                    height = height,
                    title = sp[:title],
                    canvas = canvas_type
                ))
                continue
            end
        end

        # # make it a bar canvas if plotting bar
        # if any(series -> series[:seriestype] == :bar, series_list(sp))
        #     canvas_type = UnicodePlots.BarplotGraphics
        # end

        o = UnicodePlots.Plot(x, y, canvas_type;
            width = width,
            height = height,
            title = sp[:title],
            xlim = xlim,
            ylim = ylim,
            border = isijulia() ? :ascii : :solid
        )

        # set the axis labels
        UnicodePlots.xlabel!(o, xaxis[:guide])
        UnicodePlots.ylabel!(o, yaxis[:guide])

        # now use the ! functions to add to the plot
        for series in series_list(sp)
            addUnicodeSeries!(o, series.plotattributes, sp[:legend] != :none, xlim, ylim)
        end

        # save the object
        push!(plt.o, o)
    end
end

# add a single series
function addUnicodeSeries!(o, plotattributes::KW, addlegend::Bool, xlim, ylim)
    # get the function, or special handling for step/bar/hist
    st = plotattributes[:seriestype]
    if st == :histogram2d
        UnicodePlots.densityplot!(o, plotattributes[:x], plotattributes[:y])
        return
    end

    if st in (:path, :straightline)
        func = UnicodePlots.lineplot!
    elseif st == :scatter || plotattributes[:markershape] != :none
        func = UnicodePlots.scatterplot!
    # elseif st == :bar
    #     func = UnicodePlots.barplot!
    elseif st == :shape
        func = UnicodePlots.lineplot!
    else
        error("Linestyle $st not supported by UnicodePlots")
    end

    # get the series data and label
    x, y = if st == :straightline
        straightline_data(plotattributes)
    elseif st == :shape
        shape_data(series)
    else
        [collect(float(plotattributes[s])) for s in (:x, :y)]
    end
    label = addlegend ? plotattributes[:label] : ""

    # if we happen to pass in allowed color symbols, great... otherwise let UnicodePlots decide
    color = plotattributes[:linecolor] in UnicodePlots.color_cycle ? plotattributes[:linecolor] : :auto

    # add the series
    x, y = Plots.unzip(collect(Base.Iterators.filter(xy->isfinite(xy[1])&&isfinite(xy[2]), zip(x,y))))
    func(o, x, y; color = color, name = label)
end

# -------------------------------

# since this is such a hack, it's only callable using `png`... should error during normal `show`
function png(plt::AbstractPlot{UnicodePlotsBackend}, fn::AbstractString)
    fn = addExtension(fn, "png")

    # make some whitespace and show the plot
    println("\n\n\n\n\n\n")
    gui(plt)

    # @osx_only begin
    @static if Sys.isapple()
        # BEGIN HACK

        # wait while the plot gets drawn
        sleep(0.5)

        # use osx screen capture when my terminal is maximized and cursor starts at the bottom (I know, right?)
        # TODO: compute size of plot to adjust these numbers (or maybe implement something good??)
        run(`screencapture -R50,600,700,420 $fn`)

        # END HACK (phew)
        return
    end

    error("Can only savepng on osx with UnicodePlots (though even then I wouldn't do it)")
end

# -------------------------------

# we don't do very much for subplots... just stack them vertically

function unicodeplots_rebuild(plt::Plot{UnicodePlotsBackend})
    w, h = plt[:size]
    plt.attr[:color_palette] = [RGB(0,0,0)]
    rebuildUnicodePlot!(plt, div(w, 10), div(h, 20))
end

function _show(io::IO, ::MIME"text/plain", plt::Plot{UnicodePlotsBackend})
    unicodeplots_rebuild(plt)
    foreach(x -> show(io, x), plt.o)
    nothing
end

function _display(plt::Plot{UnicodePlotsBackend})
    unicodeplots_rebuild(plt)
    map(show, plt.o)
    nothing
end


text:
using Interpolations
using Base.Threads
using Printf
using FFTW
using NetCDF
using Base.Filesystem

const G = begin
    #              ∞     (-1)^n * (2n + 1)
    # G(r) ≈ 1/π * Σ ——————————————————————————
    #             n=0 [(2n + 1)^2 + r^2]^(3/2)
    #
    # r ∈ [0; √(800)]
    Δr = .01
    l_r = 0:Δr:(sqrt(800) + Δr)

    function Σ(r::Float64)
        Σ = 0.
        for n ∈ 10000:-1:0
            frac = (2n + 1) / ((2n + 1)^2 + r^2)^1.5
            if n & 0x1 != 0 # n is odd, subtract
                Σ -= frac
            else # n is even, add
                Σ += frac
            end
        end

        return Σ
    end

    l_G = [1/Float64(π) * Σ(r) for r in l_r]

    itp = interpolate(l_G, BSpline(Constant()))
    extrapolate(scale(itp, l_r), Flat())
end

function precalculate_σ(h_min :: Float64, h_max :: Float64, Δx :: Float64, Δy :: Float64; n_h :: Float64 = 20.)
    Δh = (h_max - h_min) / 10
    if Δh == 0.; Δh = 1.; end
    l_h = max(0.0001, h_min):Δh:(h_max + Δh)

    #                       h^2
    # σ = ———————————————————————————————————————
    #       Σ   Σ  G(√[(nΔx)^2 + (mΔy)^2]/h)ΔxΔy
    #      n∈N m∈M
    #
    # N = {n ∈ Z| |nΔx| ≤ n_h * h}
    # M = {m ∈ Z| |mΔy| ≤ n_h * h}
    #
    # As the domain of n and m are point symmetric around 0, only ~1/4th of the sum element have to be evaluated:
    # * One quadrant (n > 0, m > 0), multiplied by 4
    # * One half of the x-axis (n > 0), multiplied by 2
    # * One hald of the y-axis (m > 0), multiplied by 2
    # * The point at the origin (n, m) = (0, 0)

    function σ(h)
        σ_inv = 0.
        n_max = ceil(Int, (n_h * h / Δx))
        m_max = ceil(Int, (n_h * h / Δy))

        N = n_max:-1:1
        M = m_max:-1:1

        for n ∈ N, m ∈ M
            σ_inv += 4*G(sqrt((n*Δx)^2 + (m*Δy)^2) / h)
        end

        for n ∈ N
            σ_inv += 2*G(n*Δx / h)
        end

        for m ∈ M
            σ_inv += 2*G(m*Δy / h)
        end

        σ_inv += G(0)

        return h^2 / (σ_inv * Δx * Δy)
    end

    l_σ = [σ(h) for h ∈ l_h]

    itp = interpolate(l_σ, BSpline(Cubic(Flat(OnGrid()))))
    return extrapolate(scale(itp, l_h), Flat())
end

function apply_kajiura!(b::AbstractArray{Float64, 2}, d::AbstractArray{Float64, 2}, η::AbstractArray{Float64, 2},
                h_min :: Float64, h_max :: Float64, Δx :: Float64, Δy :: Float64; water_level :: Float64 = 0., n_h :: Float64 = 20.,
                σ = precalculate_σ(h_min, h_max, Δx, Δy; n_h=n_h))

    filter_nx_half = ceil(Int, n_h * h_max / Δx / 2)
    filter_ny_half = ceil(Int, n_h * h_max / Δy / 2)

    nx = size(η, 2)
    ny = size(η, 1)

    filter = zeros(Float64, (ny, nx))

    @assert filter_nx_half < nx
    @assert filter_ny_half < ny

    println("  Computing filter matrix...")
    Threads.@threads for x ∈ -filter_nx_half:filter_nx_half
        for y ∈ -filter_ny_half:filter_ny_half
            filter[((ny + y) % ny) + 1, ((nx + x) % nx) + 1] = G(sqrt((x * Δx)^2 + (y * Δy)^2) / h_max)
        end
    end

    println("  Computing displacement matrix...")
    Threads.@threads for x ∈ 1:size(η, 2)
        for y ∈ 1:size(η, 1)
            h_yx = max(0., water_level - b[y, x]) # height 0 on land

            η[y, x] =
                if h_yx != 0.
                    σ(h_yx) * Δx * Δy / h_yx^2 * d[y, x]
                else
                    0. # No displacement where no water is (Kajiura not applicable)
                end
        end
    end

    println("  Computing filter FFT...")
    Filter = fft(filter)

    println("  Computing displacement FFT...")
    Η = fft(η)

    println("  Computing product...")
    Η .*= Filter

    println("  Computing IFFT...")
    η_complex = reinterpret(Float64, ifft(Η))
    copyto!(η, @view η_complex[1:2:end-1])
end

function main()
    if length(ARGS) != 3
        println("Usage: julia ./kajiura.jl <in_file.nc> <out_file.nc> <timestep_end>")
        return
    end

    println("Using $(nthreads()) threads.")

    in_filename = ARGS[1]
    out_filename = ARGS[2]
    t_end = parse(Int, ARGS[3])
    nc = NetCDF.open(in_filename, mode=NC_NOWRITE)
    d = nc["d"]

    lx = nc["x"]
    ly = nc["y"]

    Δx = lx[2] - lx[1]
    Δy = ly[2] - ly[1]

    nx = length(lx)
    ny = length(ly)

    l_times = nc["time"]
    n_times = length(l_times)
    last_timestamp = l_times[end]

    if t_end == -1 || t_end > n_times - 1
	    t_end = n_times - 1
    end

    println("Processing $nx × $ny cells over $(t_end + 1) timesteps.")

    b = Array{Float64, 2}(undef, (ny, nx))
    ncread!(in_filename, "b", b)
    current_disp = zeros(ny, nx)
    current_η_diff = Array{Float64, 2}(undef, (ny, nx))
    current_d_diff = Array{Float64, 2}(undef, (ny, nx))

    println("Creating output file")

    isfile(out_filename) && rm(out_filename)
    timeatts = Dict("units" => "seconds")
    xatts = Dict("units" => "m")
    yatts = Dict("units" => "m")

    nccreate(out_filename, "eta_diff", "y", ly, yatts, "x", lx, xatts, "time", l_times[1:t_end+1], timeatts)
    nccreate(out_filename, "d_diff", "y", "x", "time")

    nccreate(out_filename, "b", "y", "x")
    ncwrite(b, out_filename, "b")

    for t in 1:t_end + 1
        println("Working on timestep $(t - 1) of $(t_end)")
        current_η_diff .= 0.
        current_d_diff .= d[:,:,t] .- current_disp
        apply_kajiura!(b .+ current_disp, current_d_diff, current_η_diff, -maximum(b), -minimum(b), Δx, Δy)
        current_disp = d[:,:,t]

        println("  Writing output for timestep")
        ncwrite(current_d_diff, out_filename, "d_diff", start=[1,1,t], count=[-1,-1,1])
        ncwrite(current_η_diff, out_filename, "eta_diff", start=[1,1,t], count=[-1,-1,1])
    end
end

println("Using $(nthreads()) threads.")

Pkg.precompile()
println("Done.")

main()

max_stars_repo_name: SeisSol/SAMPLER

max_issues_repo_name: SeisSol/SAMPLER

max_forks_repo_name: SeisSol/SAMPLER


text:
#This is a structure that hold all of the information needed about the Model

mutable struct ModelInformation
  Parameters::Array{Union{Symbol,String},1}
  States::Array{Union{Symbol,String},1}
end

mutable struct MCMCSetup
  MI::ModelInformation
  ODE::ODEProblem
  Alg::OrdinaryDiffEqAlgorithm
  Data::Array
  Samples::Int64
  Priors::Array{UnivariateDistribution,1}
end

function ConvertData(data::DataFrame)
  #Convert Data into an array of arrays for iterating
  timePoints = [data[data.Experiment .==i,2] for i ∈ unique(data.Experiment)]
  experiments = [data[data.Experiment .==i,3:end] for i ∈ unique(data.Experiment)]
  #Convert into matrices and vectors
  timeVec = convert.(Vector,timePoints)
  expArray = convert.(Matrix,experiments)

  return [i for i in zip(timeVec,expArray)]
end

function BestParSet(result::Chains,parNames)
  m = mean(result)
  bestParTup = NamedTuple{tuple(m[:parameters]...)}(m[:mean])

  return [bestParTup[p] for p in parNames]
end

max_stars_repo_name: RobertGregg/TuringODE

max_issues_repo_name: RobertGregg/TuringODE

max_forks_repo_name: RobertGregg/TuringODE


text:
#=
Copyright 2017 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

# Convergence is checked starting from the simulated third cardiac cycle.
# Waveforms are first loaded and subdivided depending on the number of cycles
# already passed. This task is accomplished by function
# [`extractWaveform`](check_convergence.html#extractWaveform). Waveforms are
# then resampled in order to contain a smaller number of points. This would
# speed up the following error calculation. These two steps are taken by
# [`resampleAndCalculateError`](check_convergence.html#resampleAndCalculateError)
# function. Loading, re-sampling, and error calculation operations are done
# for all the quantities calculated by `openBF`. All these functions are
# wrapped into
# [`checkAllQuantities`](check_convergence.html#checkAllQuantities) which is
# called inside [main.jl](main.html#check_convergence) loop.

# *function* __`checkAllQuantities`__ $\rightarrow$ `max_error::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# --------------- ------------------------------------------------------------
# `v`             `::Vessel` vessel of which the error is being checked.
#
# `passed_cycles` `::Int` number of cardiac cycles already passed.
#
# `n_pts`         `::Int` number of points in the resampled waveform.
# ----------------------------------------------------------------------------
#
function checkAllQuantities(v             :: Vessel,
                            passed_cycles :: Int64,
                            n_pts         :: Int64)
  # --------------------------------------------------------------------------
  # Functioning
  # --------------------------------------------------------------------------
  # `er` is an array which will contain the percentage error for each quantity
  # in the current vessel. `qs` collection contains all the suffixes needed to
  # build `.out` filenames for each quantity.
  # --------------------------------------------------------------------------
  er = zeros(Float64, 5)

  qs = ["_P.out", "_Q.out", "_A.out", "_c.out", "_u.out"]
  # For each quantity `q` listed in `qs`, the `.out` filename is created by
  # joining the `label` inside `v` and the suffix `q`. Then the desired
  # waveforms are loaded by
  # [`extractWaveform`](check_convergence.html#extractWaveform), re-sampled,
  # and the error is calculated by
  # [`resampleAndCalculateError`](check_convergence.html#resampleAndCalculateError).
  # Eventually, the maximum error scored among all the quantities is returned.
  #
  # --------------------------------------------------------------------------
  # Returns:
  # ------------- ------------------------------------------------------------
  # `max_error`   `::Float` maximum error within all the quantities in the
  #               current vessel.
  # --------------------------------------------------------------------------
  idx = 1
  for q in qs
    filename = join([v.label, q])

    W = extractWaveform(filename, passed_cycles)

    er[idx] = resampleAndCalculateError(v, W, passed_cycles, n_pts)

    idx += 1
  end

  return maximum(er)
end

# *function* __`checkAllQuantities`__ $\rightarrow$ `max_error::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# --------------- ------------------------------------------------------------
# `vs`             `::Array{Vessel, 1}` collection of vessels.
#
# `passed_cycles` `::Int` number of cardiac cycles already passed.
#
# `n_pts`         `::Int` number of points in the resampled waveform.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# To compute the error for all the vessels, `checkAllQuantities` is
# re-defined to handle also a vector of `Vessel`s. `checkAllQuantities` is
# called recursively.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ------------- --------------------------------------------------------------
# `max_error`   `::Float` maximum error among all the vessels.
# ----------------------------------------------------------------------------
#
function checkAllQuantities(vs            :: Array{Vessel, 1},
                            passed_cycles :: Int64,
                            n_pts         :: Int64)

  ers = zeros(Float64, length(vs))

  idx = 1
  for v in vs

    ers[idx] = checkAllQuantities(v, passed_cycles, n_pts)
    idx += 1

  end

  return maximum(ers)
end

# *function* __`extractWaveform`__ $\rightarrow$ `Fx::Array{Array{Float, 2}}`
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `filename`    `::String` current time in seconds.
#
# `cycles`      `::Int` $\Delta t$, current time step.
# ----------------------------------------------------------------------------
#
function extractWaveform(filename :: String, cycles :: Int64)
  # --------------------------------------------------------------------------
  # Functioning
  # --------------------------------------------------------------------------
  # `filename` is a string referring to the `.out` file containing the
  # waveform to be imported. This is read by `readdlm` function and stored
  # into
  # `F` array. `F` is a timesteps $\times$ 6 array; the first column contains
  # the time variable, the remaining columns are one for each monitored node
  # along the current vessel.
  # --------------------------------------------------------------------------
  F = readdlm(filename)
  # These three lines find the cardiac period `t` by knowing how many
  # cycles have already passed.
  sf1 = int(length(F[:,1])*(cycles-1)/cycles)
  T1 = F[sf1:end,1] .- F[sf1,1]
  t = T1[end]
  # `Fx` is an empty array which will contain one array for each cardiac
  # cycle.
  Fx = Array{Float64, 2}[]
  # This loop skims through `F`, separates waveforms, and add them to `Fx`.
  # `si` is the waveform source index, the index from which the waveform
  # starts. `ti` is the waveform terminal index, the index at which the
  # waveforms ends because a new cardiac cycle is due to start.
  si = 1
  ti = 1
  idx = 1
  for i in 1:length(F[:,1])-1
    if F[i,1] < t*idx && F[i+1, 1] >= t*idx
      ti = i
      push!(Fx, F[si:ti,:])
      si = ti
      idx += 1
    end
  end
  # --------------------------------------------------------------------------
  # Returns:
  # ------------- ------------------------------------------------------------
  # `Fx`           `::Array{Array{Float, 2}}`
  # --------------------------------------------------------------------------
  return Fx
end

# *function* __`resampleAndCalculateError`__ $\rightarrow$
# `max_error::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ----------------- ----------------------------------------------------------
# `W`               `::Array{Array{Float, 2}}` collection of all waveforms
#                   during the passed cardiac cycles for the current
#                   quantity.
#
# `passed_cycles`   `::Int` number of cardiac cycles already simulated.
#
# `n_pts`           `::Int` number of points to be contained in the
#                   re-sampled waveform.
# ----------------------------------------------------------------------------
#
function resampleAndCalculateError(v, W :: Array{Array{Float64, 2}, 1},
                                   passed_cycles :: Int64,
                                   n_pts         :: Int64)
  # --------------------------------------------------------------------------
  # Functioning
  # --------------------------------------------------------------------------
  # `F` will contain the re-sampled waveform, `er` is an empty variable
  # for the maximum error.
  # --------------------------------------------------------------------------
  F = zeros(Float64, n_pts, 2)
  er = 0
  # The error is computed for all the monitor nodes. Their values are stored
  # in `W`'s columns from 2 to 6. Each waveform is re-sampled by using
  # [`resample`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.resample.html#scipy.signal.resample)
  # function from `scipy.signal` library. In order to compute the error, only
  # the last and the previous waveforms are analysed.
  for j in 2:6

    F[:,1], t = signal.resample(W[end-1][:,j], n_pts, W[end-1][:,1])
    F[:,2], t = signal.resample(W[end  ][:,j], n_pts, W[end  ][:,1])

    er = abs(F[:,1] .- F[:,2])./F[:,1]*100
  end
  # --------------------------------------------------------------------------
  # Returns:
  # ----------- --------------------------------------------------------------
  # `max_error` `::Float` maximum error on the current quantity along the
  #             the current vessel.
  # --------------------------------------------------------------------------
  return maximum(er)
end

function checkConvergence(vessels :: Array{Vessel, 1},
  convergence_tollerance :: Float64)

  for v in vessels
    err = checkConvergence(v, convergence_tollerance)
    if err > convergence_tollerance
      return err
    end
  end
  return err
end

function checkConvergence(v :: Vessel,
  convergence_tollerance :: Float64)

  qs = ["_P", "_Q"]

  for q in qs
    filename_last = join([v.label, q, ".last"])
    filename_temp = join([v.label, q, ".temp"])

    err = checkQuantity(filename_last, filename_temp, convergence_tollerance)
    if err > convergence_tollerance
      return err
    end
  end
  return err
end

function checkQuantity(last, temp, convergence_tollerance)

  s1 = readdlm(last)
  s2 = readdlm(temp)
  n_pts = 1000

  r1 = n_pts/length(s1[:,2])
  r2 = n_pts/length(s2[:,2])
  for j in 2:5
    sr1 = resample(s1[:,j], r1)
    sr2 = resample(s2[:,j], r2)
    for i in 20:length(sr1)-20
      err = abs(100*(sr2[i]-sr1[i])/sr1[i])
      if err > convergence_tollerance
        # println("\n",err)
        return err
      end
    end
  end
  return 9.99
end

function checkConvergence(edge_list, vessels :: Array{Vessel, 1}, passed_cycles :: Int64)
    qs = "_A"
    err = zeros(size(edge_list)[1]).+100
    if passed_cycles >= 2

        for i in 1:size(edge_list)[1]
            v = vessels[i]
            lbl = v.label

            filename_last = join([lbl, qs, ".last"])
            filename_temp = join([lbl, qs, ".temp"])
            w_last = readdlm(filename_last)
            w_temp = readdlm(filename_temp)

            if length(w_last[:,1]) == length(w_temp[:,1])
                err[i] = maximum(abs.((w_last[2:end,:].-w_temp[2:end,:])./w_last[2:end,:])*100)
            end

        end
    end
    return maximum(err)
end


text:
@testset "Configurations" begin
    @testset "Construction" begin
        config = Configuration([o"1s", o"2s", o"2p", o"3s", o"3p"], [2,2,6,2,6], [:closed])
        rconfig = Configuration([ro"1s", ro"2s", ro"2p", ro"3s", ro"3p"], [2,2,6,2,6], [:closed])

        @test config.orbitals ==
            [o"1s", o"2s", o"2p", o"3s", o"3p"]
        @test rconfig.orbitals ==
            [ro"1s", ro"2s", ro"2p-", ro"2p", ro"3s", ro"3p-", ro"3p"]

        @test config.occupancy == [2, 2, 6, 2, 6]
        @test rconfig.occupancy == [2, 2, 2, 4, 2, 2, 4]

        @test config.states == [:closed, :open, :open, :open, :open]
        @test rconfig.states == [:closed, :open, :open, :open, :open, :open, :open]

        @test c"1s2c 2s2 2p6 3s2 3p6" == config
        @test c"1s2c.2s2.2p6.3s2.3p6" == config
        @test c"[He]c 2s2 2p6 3s2 3p6" == config

        @test rc"1s2c 2s2 2p6 3s2 3p6" == rconfig
        @test rc"1s2c.2s2.2p6.3s2.3p6" == rconfig
        @test rc"[He]c 2s2 2p6 3s2 3p6" == rconfig

        @test c"[Kr]c 5s2" == Configuration([o"1s", o"2s", o"2p", o"3s", o"3p", o"3d", o"4s", o"4p", o"5s"],
                                            [2,2,6,2,6,10,2,6,2],
                                            [:closed,:closed,:closed,:closed,:closed,:closed,:closed,:closed,:open])

        @test length(c"") == 0
        @test length(rc"") == 0

        @test rc"1s ld-2 kp6" == Configuration([ro"1s", ro"kp-", ro"kp", ro"ld-"], [1, 2, 4, 2])

        @test_throws ArgumentError parse(Configuration{Orbital}, "1sc")
        @test_throws ArgumentError parse(Configuration{Orbital}, "1s 1s")
        @test_throws ArgumentError parse(Configuration{Orbital}, "[He]c 1s")
        @test_throws ArgumentError parse(Configuration{Orbital}, "1s3")
        @test_throws ArgumentError parse(Configuration{RelativisticOrbital}, "1s3")
        @test_throws ArgumentError parse(Configuration{Orbital}, "1s2 2p-2")

        @test fill(c"1s 2s 2p") == c"1s2 2s2 2p6"
        @test fill(rc"1s 2s 2p- 2p") == rc"1s2 2s2 2p-2 2p4"
        @test close(c"1s2") == c"1s2c"
        @test_throws ArgumentError close(c"1s")

        # Tests for #19
        @test c"10s2" == Configuration([o"10s"], [2], [:open])
        @test c"9999l32" == Configuration([o"9999l"], [32], [:open])
    end

    @testset "Number of electrons" begin
        @test num_electrons(c"[He]") == 2
        @test num_electrons(rc"[He]") == 2
        @test num_electrons(c"[Xe]") == 54
        @test num_electrons(rc"[Xe]") == 54

        @test num_electrons(c"[Xe]", o"1s") == 2
        @test num_electrons(rc"[Xe]", ro"1s") == 2
        @test num_electrons(c"[Ne] 3s 3p", o"3p") == 1
        @test num_electrons(rc"[Kr] Af-2 Bf5", ro"Bf") == 5
        @test num_electrons(c"[Kr]", ro"5s") == 0
        @test num_electrons(rc"[Rn]", ro"As") == 0
    end

    @testset "Access subsets" begin
        @test core(c"[Kr]c 5s2") == c"[Kr]c"
        @test peel(c"[Kr]c 5s2") == c"5s2"
        @test core(c"[Kr]c 5s2c 5p6") == c"[Kr]c 5s2c"
        @test peel(c"[Kr]c 5s2c 5p6") == c"5p6"
        @test active(c"[Kr]c 5s2") == c"5s2"
        @test inactive(c"[Kr]c 5s2i") == c"5s2i"
        @test inactive(c"[Kr]c 5s2") == c""

        @test bound(c"[Kr] 5s2 5p5 ks") == c"[Kr] 5s2 5p5"
        @test continuum(c"[Kr] 5s2 5p5 ks") == c"ks"
        @test bound(c"[Kr] 5s2 5p4 ks ld") == c"[Kr] 5s2 5p4"
        @test continuum(c"[Kr] 5s2 5p4 ks ld") == c"ks ld"

        @test bound(rc"[Kr] 5s2 5p-2 5p3 ks") == rc"[Kr] 5s2 5p-2 5p3"
        @test continuum(rc"[Kr] 5s2 5p-2 5p3 ks") == rc"ks"
        @test bound(rc"[Kr] 5s2 5p-2 5p2 ks ld") == rc"[Kr] 5s2 5p-2 5p2"
        @test continuum(rc"[Kr] 5s2 5p-2 5p2 ks ld") == rc"ks ld"

        @test c"[Ne]"[1] == (o"1s",2,:closed)
        @test c"[Ne]"[1:2] == c"1s2c 2s2c"
        @test c"[Ne]"[end-1:end] == c"2s2c 2p6c"

        @test c"1s2 2p kp"[2:3] == c"2p kp"

        @test o"1s" ∈ c"[He]"
        @test ro"1s" ∈ rc"[He]"
    end

    @testset "Sorting" begin
        @test c"1s2" < c"1s 2s"
        @test c"1s2" < c"1s 2p"
        @test c"1s2" < c"ks2"
        # @test c"kp2" < c"kp kd" # Ideally this would be true, but too
        #                         # complicated to implement
    end

    @testset "Parity" begin
        @test iseven(parity(c"1s"))
        @test iseven(parity(c"1s2"))
        @test iseven(parity(c"1s2 2s"))
        @test iseven(parity(c"1s2 2s2"))
        @test iseven(parity(c"[He]c 2s"))
        @test iseven(parity(c"[He]c 2s2"))
        @test isodd(parity(c"[He]c 2s 2p"))
        @test isodd(parity(c"[He]c 2s2 2p"))
        @test iseven(parity(c"[He]c 2s 2p2"))
        @test iseven(parity(c"[He]c 2s2 2p2"))
        @test iseven(parity(c"[He]c 2s2 2p2 3d"))
        @test isodd(parity(c"[He]c 2s kp"))
    end

    @testset "Number of electrons" begin
        @test count(c"1s") == 1
        @test count(c"[He]") == 2
        @test count(c"[Xe]") == 54
        @test count(peel(c"[Kr]c 5s2 5p6")) == 8
    end

    @testset "Pretty printing" begin
        Xe⁺ = rc"[Kr]c 5s2 5p-2 5p3"
        map([c"1s" => "1s",
             c"1s2" => "1s²",
             c"1s2 2s2" => "1s² 2s²",
             c"1s2 2s2 2p" => "1s² 2s² 2p",
             c"1s2c 2s2 2p" => "[He]ᶜ 2s² 2p",
             c"[He]* 2s2" => "1s² 2s²",
             c"[He]c 2s2" => "[He]ᶜ 2s²",
             c"[He]i 2s2" => "1s²ⁱ 2s²",
             c"[Kr]c 5s2" => "[Kr]ᶜ 5s²",
             Xe⁺ => "[Kr]ᶜ 5s² 5p⁻² 5p³",
             core(Xe⁺) => "[Kr]ᶜ",
             peel(Xe⁺) => "5s² 5p⁻² 5p³",
             rc"[Kr] 5s2c 5p6" => "[Kr]ᶜ 5s²ᶜ 5p⁻² 5p⁴",
             c"[Ne]"[end:end] => "2p⁶ᶜ",
             rc"[Ne]"[end-1:end] => "2p⁻²ᶜ 2p⁴ᶜ",
             c"5s2" => "5s²",
             rc"[Kr]*" => "1s² 2s² 2p⁻² 2p⁴ 3s² 3p⁻² 3p⁴ 3d⁻⁴ 3d⁶ 4s² 4p⁻² 4p⁴",
             c"[Kr]c" =>"[Kr]ᶜ",
             c"1s2 kp" => "1s² kp",
             c"" => "∅"]) do (c,s)
                 @test "$(c)" == s
             end
    end

    @testset "Orbital substitutions" begin
        @test replace(c"1s2", o"1s"=>o"2p") == c"1s 2p"
        @test replace(c"1s2", o"1s"=>o"kp") == c"1s kp"
        @test replace(c"1s kp", o"kp"=>o"ld") == c"1s ld"
        @test_throws ArgumentError replace(c"1s2", o"2p"=>o"3p")
        @test_throws ArgumentError replace(c"1s2 2s", o"2s"=>o"1s")
    end

    @testset "Orbital removal" begin
        @test c"1s2" - o"1s" == c"1s"
        @test c"1s" - o"1s" == c""
        @test c"1s 2s" - o"1s" == c"2s"
        @test c"[Ne]" - o"2s" == c"[He] 2s 2p6c"
        @test -(c"1s2 2s", o"1s", 2) == c"2s"

        @test delete!(c"1s 2s", o"1s") == c"2s"
        @test delete!(c"1s2 2s", o"1s") == c"2s"
        @test delete!(c"[Ne]", o"2p") == c"1s2c 2s2c"
        @test delete!(rc"[Ne]", ro"2p-") == rc"1s2c 2s2c 2p4c"
    end

    @testset "Configuration additions" begin
        @test c"1s" + c"1s" == c"[He]*"
        @test c"1s" + c"2p" == c"1s 2p"
        @test c"1s" + c"ld" == c"1s ld"
        @test c"[Kr]*" + c"4d10" + c"5s2" + c"5p6" == c"[Xe]*"
        @test_throws ArgumentError c"[He]*" + c"1s"
        @test_throws ArgumentError c"1si" + c"1s"
    end

    @testset "Configuration juxtapositions" begin
        @test c"" ⊗ c"" == [c""]
        @test c"1s" ⊗ c"" == [c"1s"]
        @test c"" ⊗ c"1s" == [c"1s"]
        @test c"1s" ⊗ c"2s" == [c"1s 2s"]
        @test rc"1s" ⊗ [rc"2s", rc"2p-", rc"2p"] == [rc"1s 2s", rc"1s 2p-", rc"1s 2p"]
        @test [rc"1s", rc"2s"] ⊗ rc"2p-" == [rc"1s 2p-", rc"2s 2p-"]
        @test [rc"1s", rc"2s"] ⊗ [rc"2p-", rc"2p"] == [rc"1s 2p-", rc"1s 2p", rc"2s 2p-", rc"2s 2p"]
        @test [rc"1s 2s"] ⊗ [rc"2p-2", rc"2p4"] == [rc"1s 2s 2p-2", rc"1s 2s 2p4"]
        @test [rc"1s 2s"] ⊗ [rc"kp-2", rc"lp4"] == [rc"1s 2s kp-2", rc"1s 2s lp4"]
    end

    @testset "Non-relativistic orbitals" begin
        import AtomicLevels: rconfigurations_from_orbital
        @test rconfigurations_from_orbital(1, 0, 1) == [rc"1s"]
        @test rconfigurations_from_orbital(1, 0, 2) == [rc"1s2"]
        @test rconfigurations_from_orbital(2, 0, 2) == [rc"2s2"]

        @test rconfigurations_from_orbital(2, 1, 1) == [rc"2p-", rc"2p"]
        @test rconfigurations_from_orbital(2, 1, 2) == [rc"2p-2", rc"2p- 2p", rc"2p2"]
        @test rconfigurations_from_orbital(2, 1, 3) == [rc"2p-2 2p1", rc"2p- 2p2", rc"2p3"]
        @test rconfigurations_from_orbital(2, 1, 4) == [rc"2p-2 2p2", rc"2p- 2p3", rc"2p4"]
        @test rconfigurations_from_orbital(2, 1, 5) == [rc"2p-2 2p3", rc"2p- 2p4"]
        @test rconfigurations_from_orbital(2, 1, 6) == [rc"2p-2 2p4"]

        @test rconfigurations_from_orbital(4, 2, 10) == [rc"4d-4 4d6"]
        @test rconfigurations_from_orbital(4, 2, 5) == [rc"4d-4 4d1", rc"4d-3 4d2", rc"4d-2 4d3", rc"4d-1 4d4", rc"4d5"]

        @test rconfigurations_from_orbital(:k, 2, 5) == [rc"kd-4 kd1", rc"kd-3 kd2", rc"kd-2 kd3", rc"kd-1 kd4", rc"kd5"]

        @test_throws ArgumentError rconfigurations_from_orbital(1, 2, 1)
        @test_throws ArgumentError rconfigurations_from_orbital(1, 0, 3)

        @test rconfigurations_from_orbital(o"1s", 2) == [rc"1s2"]
        @test rconfigurations_from_orbital(o"2p", 2) == [rc"2p-2", rc"2p- 2p", rc"2p2"]
        @test_throws ArgumentError rconfigurations_from_orbital(o"3d", 20)

        @test rcs"1s" == [rc"1s"]
        @test rcs"2s2" == [rc"2s2"]
        @test rcs"2p4" == [rc"2p-2 2p2", rc"2p- 2p3", rc"2p4"]
        @test rcs"4d5" == [rc"4d-4 4d1", rc"4d-3 4d2", rc"4d-2 4d3", rc"4d-1 4d4", rc"4d5"]

        @test rc"1s2" ⊗ rcs"2p" == [rc"1s2 2p-", rc"1s2 2p"]
        @test rc"1s2" ⊗ rcs"kp" == [rc"1s2 kp-", rc"1s2 kp"]
    end

    @testset "Spin-orbitals" begin
        @test_throws ArgumentError Configuration(spin_orbitals(o"1s"), [2,1])
        @test_throws ArgumentError Configuration(spin_orbitals(o"1s"), [1,2])

        @test spin_configurations(c"1s2") ==
            [Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"1s",0,false)], [1, 1])]

        @test spin_configurations(c"1s 2p") ==
            [Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"2p",-1,true)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"2p",-1,false)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"2p",0,true)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"2p",0,false)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"2p",1,true)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"2p",1,false)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,false), SpinOrbital(o"2p",-1,true)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,false), SpinOrbital(o"2p",-1,false)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,false), SpinOrbital(o"2p",0,true)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,false), SpinOrbital(o"2p",0,false)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,false), SpinOrbital(o"2p",1,true)], [1, 1]),
             Configuration([SpinOrbital(o"1s",0,false), SpinOrbital(o"2p",1,false)], [1, 1])]

        @test spin_configurations(
            excited_configurations(c"1s2", os"k[s-p]"..., max_excitations=:singles, keep_parity=false)) ==
                [Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"1s",0,false)], [1, 1]),

                 Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"ks",0,true)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"ks",0,false)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"kp",-1,true)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"kp",-1,false)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"kp",0,true)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"kp",0,false)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"kp",1,true)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"kp",1,false)], [1, 1]),

                 Configuration([SpinOrbital(o"1s",0,false), SpinOrbital(o"ks",0,true)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,false), SpinOrbital(o"ks",0,false)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,false), SpinOrbital(o"kp",-1,true)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,false), SpinOrbital(o"kp",-1,false)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,false), SpinOrbital(o"kp",0,true)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,false), SpinOrbital(o"kp",0,false)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,false), SpinOrbital(o"kp",1,true)], [1, 1]),
                 Configuration([SpinOrbital(o"1s",0,false), SpinOrbital(o"kp",1,false)], [1, 1])]

        @test bound(Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"ks",0,true)], [1, 1])) ==
            Configuration([SpinOrbital(o"1s",0,true),], [1,])

        @test continuum(Configuration([SpinOrbital(o"1s",0,true), SpinOrbital(o"ks",0,true)], [1, 1])) ==
            Configuration([SpinOrbital(o"ks",0,true),], [1,])

        @test all([o ∈ spin_configurations(c"1s2")[1]
                   for o in spin_orbitals(o"1s")])

        @test map(spin_configurations(c"[Kr] 5s2 5p6")[1]) do (orb,occ,state)
            orb.orb.n < 5 && state == :closed || orb.orb.n == 5 && state == :open
        end |> all

        @test string.(spin_configurations(c"2s2 2p")) ==
            ["2s² 2p₋₁α",
             "2s² 2p₋₁β",
             "2s² 2p₀α",
             "2s² 2p₀β",
             "2s² 2p₁α",
             "2s² 2p₁β"]

        @test string.(spin_configurations(c"[Kr] 5s2 5p5 ks")) ==
            ["[Kr]ᶜ 5s² 5p₋₁² 5p₀² 5p₁α ks₀α",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀² 5p₁α ks₀β",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀² 5p₁β ks₀α",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀² 5p₁β ks₀β",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀α 5p₁² ks₀α",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀α 5p₁² ks₀β",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀β 5p₁² ks₀α",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀β 5p₁² ks₀β",
             "[Kr]ᶜ 5s² 5p₋₁α 5p₀² 5p₁² ks₀α",
             "[Kr]ᶜ 5s² 5p₋₁α 5p₀² 5p₁² ks₀β",
             "[Kr]ᶜ 5s² 5p₋₁β 5p₀² 5p₁² ks₀α",
             "[Kr]ᶜ 5s² 5p₋₁β 5p₀² 5p₁² ks₀β"]

        @test string.(spin_configurations(c"[Kr] 5s2 5p4")) ==
            ["[Kr]ᶜ 5s² 5p₋₁² 5p₀²",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀α 5p₁α",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀α 5p₁β",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀β 5p₁α",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₀β 5p₁β",
             "[Kr]ᶜ 5s² 5p₋₁² 5p₁²",
             "[Kr]ᶜ 5s² 5p₋₁α 5p₀² 5p₁α",
             "[Kr]ᶜ 5s² 5p₋₁α 5p₀² 5p₁β",
             "[Kr]ᶜ 5s² 5p₋₁α 5p₀α 5p₁²",
             "[Kr]ᶜ 5s² 5p₋₁α 5p₀β 5p₁²",
             "[Kr]ᶜ 5s² 5p₋₁β 5p₀² 5p₁α",
             "[Kr]ᶜ 5s² 5p₋₁β 5p₀² 5p₁β",
             "[Kr]ᶜ 5s² 5p₋₁β 5p₀α 5p₁²",
             "[Kr]ᶜ 5s² 5p₋₁β 5p₀β 5p₁²",
             "[Kr]ᶜ 5s² 5p₀² 5p₁²"]

        @test substitutions(spin_configurations(c"1s2")[1],
                            spin_configurations(c"1s2")[1]) == []

        @test substitutions(spin_configurations(c"1s2")[1],
                            spin_configurations(c"1s ks")[3]) ==
                                [SpinOrbital(o"1s",0,true)=>SpinOrbital(o"ks",0,true)]
    end

    @testset "Internal utilities" begin
        @test AtomicLevels.get_noble_core_name(c"1s2 2s2 2p6") === nothing
        @test AtomicLevels.get_noble_core_name(c"1s2c 2s2 2p6") === "He"
        @test AtomicLevels.get_noble_core_name(c"1s2c 2s2c 2p6c") === "Ne"
        @test AtomicLevels.get_noble_core_name(c"[He] 2s2 2p6c 3s2c 3p6c") === "He"
        for gas in ["Rn", "Xe", "Kr", "Ar", "Ne", "He"]
            c = parse(Configuration{Orbital}, "[$gas]")
            @test AtomicLevels.get_noble_core_name(c) === gas
            rc = parse(Configuration{RelativisticOrbital}, "[$gas]")
            @test AtomicLevels.get_noble_core_name(rc) === gas
        end
    end
end


text:
# Thanks @tlienart

using Franklin, JSON

using Markdown, Dates

include("youtube_videos.jl")

const DATEFMT = dateformat"yyyy-mm-dd HH:MMp"
const TZ = "America/New_York"

function hfun_doc(params)
    fname = join(params[1:max(1, length(params)-2)], " ")
    head = params[end-1]
    type = params[end]
    doc = eval(Meta.parse("@doc $fname"))
    txt = Markdown.plain(doc)
    # possibly further processing here
    body = Franklin.fd2html(txt, internal=true)
    return """

            $head

$type

$body

    """
end

function hfun_youtube(params)
    id = params[1]
    return """

    """
end

function hfun_showtime(params)
    id = params[1]
    str = locvar(id)
    if isnothing(str)
        @warn "Unknown datetime variable $str"
        return ""
    end
    try
        DateTime(str, DATEFMT)
    catch err
        @warn "There was an error parsing date $str, the format is yyyy-mm-dd HH:MMp (see ?DateFormat)"
        rethrow(err)
    end
end

function parse_duration(str)
    str = replace(str, r"^PT"=>"")
    hrex, mrex, srex = Regex.(string.("^([0-9]+)", ["H","M","S"]))

    t = 0
    hmatch = match(hrex, str)
    if !isnothing(hmatch)
        h = parse(Int, hmatch[1])
        t += 60*60*h
        str = replace(str, hrex=>"")
    end

    mmatch = match(mrex, str)
    if !isnothing(mmatch)
        m = parse(Int, mmatch[1])
        t += 60*m
        str = replace(str, mrex=>"")
    end

    smatch = match(srex, str)
    if !isnothing(smatch)
        s = parse(Int, smatch[1])
        t += s
        str = replace(str, srex=>"")
    end

    t
end

function hfun_go_live()
    seq = locvar("sequence")
    airtime = locvar("airtime")

    if isnothing(seq)
        @warn "airtime set, but no `sequence` variable not defined." *
        "sequence is an array of video IDs to play in order on this page"
    end

    vid_ids = [get(videos, s, s) for s in seq]

    f = tempname()
    # Get the duration of each video
    download("https://www.googleapis.com/youtube/v3/videos?id=$(join(vid_ids, ","))&part=contentDetails&key=AIzaSyDZhbWHc2PTEFTx173MaTgddnWCGPqdbB8", f)
    dict = JSON.parse(String(read(f)))

    durations = [parse_duration(video["contentDetails"]["duration"])
                 for video in dict["items"]]

    jrepr(x) = sprint(io->JSON.print(io, x))
    """

    """
end

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
## TypeSortedCollections addendum
# `foreach_with_extra_args` below is a hack to avoid allocations associated with creating closures over
# heap-allocated variables. Hopefully this will not be necessary in a future version of Julia.
for num_extra_args = 1 : 5
    extra_arg_syms = [Symbol("arg", i) for i = 1 : num_extra_args]
    @eval begin
        @generated function foreach_with_extra_args(f, $(extra_arg_syms...), A1::TypeSortedCollection{<:Any, N}, As::Union{<:TypeSortedCollection{<:Any, N}, AbstractVector}...) where {N}
            extra_args = $extra_arg_syms
            expr = Expr(:block)
            push!(expr.args, :(Base.@_inline_meta)) # required to achieve zero allocation
            push!(expr.args, :(leading_tsc = A1))
            push!(expr.args, :(@boundscheck TypeSortedCollections.lengths_match(A1, As...) || TypeSortedCollections.lengths_match_fail()))
            for i = 1 : N
                vali = Val(i)
                push!(expr.args, quote
                    let inds = leading_tsc.indices[$i]
                        @boundscheck TypeSortedCollections.indices_match($vali, inds, A1, As...) || TypeSortedCollections.indices_match_fail()
                        @inbounds for j in LinearIndices(inds)
                            vecindex = inds[j]
                            f($(extra_args...), TypeSortedCollections._getindex_all($vali, j, vecindex, A1, As...)...)
                        end
                    end
                end)
            end
            quote
                $expr
                nothing
            end
        end
    end
end


text:
"""
    unsafe_swap!(v::Vector, i::Int, j::Int) -> Vector

swap `i`th and `j`th elements of a vector.
"""
@inline function unsafe_swap!(v::Vector, i::Int, j::Int)
    @inbounds temp = v[i]
    @inbounds v[i] = v[j]
    @inbounds v[j] = temp
end

"""
    minx(f, vec)

Find the element in `vec` that gives minimum `f(x)`.
"""
function minx(f, vec)
    local xmin = vec[1]
    fmin = f(xmin)
    for j=2:length(vec)
        @inbounds x = vec[j]
        fmin_ = f(x)
        fmin_ < fmin && (fmin = fmin_; xmin=x)
    end
    return xmin
end

"""
    copyltu!(A::AbstractMatrix) -> AbstractMatrix

copy the lower triangular to upper triangular.
"""
function copyltu!(A::AbstractMatrix)
    m, n = size(A)
    for i=1:m
        A[i,i] = real(A[i,i])
        for j=i+1:n
            @inbounds A[i,j] = conj(A[j,i])
        end
    end
    A
end


text:
# details on status codes

const statusses = Dict(
    :ok     => "successful",
    :skip   => "skipped",
    :fail   => "unsuccessful",
    :kill   => "interrupted",
)
const reasons = Dict(
    missing                 => missing,
    # skip
    :explicit               => "package was blacklisted",
    :jll                    => "package is a untestable wrapper package",
    :unsupported            => "package is not supported by this Julia version",
    # fail
    :unsatisfiable          => "package could not be installed",
    :untestable             => "package does not have any tests",
    :binary_dependency      => "package requires a missing binary dependency",
    :missing_dependency     => "package is missing a package dependency",
    :missing_package        => "package is using an unknown package",
    :test_failures          => "package has test failures",
    :syntax                 => "package has syntax issues",
    :gc_corruption          => "GC corruption detected",
    :segfault               => "a segmentation fault happened",
    :abort                  => "the process was aborted",
    :unreachable            => "an unreachable instruction was executed",
    :unknown                => "there were unidentified errors",
    # kill
    :time_limit             => "test duration exceeded the time limit",
    :log_limit              => "test log exceeded the size limit",
    :inactivity             => "tests became inactive",
)

# simple comparison of two versions

print_status(status, val=status) = print_status(stdout, status, val)

function print_status(io::IO, status, val=status)
    if status == :ok
        printstyled(io, val; color = :green)
    elseif status == :fail || status == :kill
        printstyled(io, val; color = Base.error_color())
    elseif status == :skip
        printstyled(io, val; color = Base.warn_color())
    else
        error("Unknown status $status")
    end
end

function compare(result; rev::Bool=false)
    pkg_names = unique(result.name)

    primary, against = groupby(result, [:julia])
    if rev  # hack
        primary, against = against, primary
    end

    # overview
    for df in (primary, against)
        o = count(==(:ok),      df[!, :status])
        s = count(==(:skip),    df[!, :status])
        f = count(==(:fail),    df[!, :status])
        k = count(==(:kill),    df[!, :status])
        x = o + s + k + f
        @assert x == nrow(df)

        print("On v$(first(df.julia)), out of $x packages ")
        print_status(:ok, o)
        print(" passed, ")
        print_status(:fail, f)
        print(" failed, ")
        print_status(:kill, k)
        print(" got killed and ")
        print_status(:skip, s)
        println(" were skipped.")
    end

    println()

    # list of differences
    println("Comparing $(first(primary.julia)) against $(first(against.julia)):")
    new_failures = 0
    new_successes = 0
    for current in eachrow(primary)
        pkg_name = current[:name]

        previous = against[against[!, :name] .== pkg_name, :]
        nrow(previous) == 0 && continue
        previous = first(previous)

        if current[:status] != previous[:status]
            print("- $pkg_name status was $(previous[:status])")
            ismissing(previous[:reason]) || print(" (reason: $(previous[:reason]))")
            print(", now ")
            print_status(current[:status])
            ismissing(current[:reason]) || print(" (reason: $(current[:reason]))")
            println()
            if current.status == :fail || current.status == :kill
                new_failures += 1
            elseif current.status == :ok
                new_successes += 1
            end
        end
    end

    println()

    # summary of differences
    print("In summary, ")
    print_status(:ok, new_successes)
    print(" packages now succeed, while ")
    print_status(:fail, new_failures)
    println(" have started to fail.")

    return
end


text:
using JML
using PrettyPrint

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    inputfilename = ARGS[1]
    midfilename = ARGS[2]
    outputfilename = ARGS[3]
    python = open(inputfilename) do f
        codegen(read(f, String), inputfilename, :py)
    end
    open(midfilename, "w") do f
        println(f, "from rupy import Rupy")
        println(f, "from rupy.pythonize import compile_to")
        print(f, "dexp = ")
        println(f, python)
        println(f, "filename = $(repr(inputfilename))")
        println(f, "compile_to(dexp, filename, $(repr(outputfilename)))")
    end
    return 0
end

max_forks_repo_name: thautwarm/Rupy


text:
using DashBootstrapComponents, DashHtmlComponents

inputs = html_div([
    dbc_input(placeholder = "A large input...", size = "lg", className = "mb-3"),
    dbc_input(placeholder = "A medium input...", size = "md", className = "mb-3"),
    dbc_input(placeholder = "A small input...", size = "sm"),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
module Atmos

export AtmosModel,
       AtmosAcousticLinearModel, AtmosAcousticGravityLinearModel,
       RemainderModel,
       AtmosLESConfiguration,
       AtmosGCMConfiguration

using LinearAlgebra, StaticArrays
using ..VariableTemplates
using ..MoistThermodynamics
using ..PlanetParameters
import ..MoistThermodynamics: internal_energy
using ..SubgridScaleParameters
using GPUifyLoops
using ..MPIStateArrays: MPIStateArray
using ..Mesh.Grids: VerticalDirection, HorizontalDirection, min_node_distance

import CLIMA.DGmethods: BalanceLaw, vars_aux, vars_state, vars_gradient,
                        vars_diffusive, flux_nondiffusive!,
                        flux_diffusive!, source!, wavespeed, boundary_state!,
                        gradvariables!, diffusive!, init_aux!, init_state!,
                        update_aux!, LocalGeometry, lengthscale,
                        resolutionmetric, DGModel, calculate_dt,
                        nodal_update_aux!, num_state,
                        num_integrals,
                        vars_integrals, vars_reverse_integrals,
                        indefinite_stack_integral!,
                        reverse_indefinite_stack_integral!,
                        integral_load_aux!, integral_set_aux!,
                        reverse_integral_load_aux!,
                        reverse_integral_set_aux!
import ..DGmethods.NumericalFluxes: boundary_state!,
                                    boundary_flux_diffusive!,
                                    NumericalFluxNonDiffusive,
                                    NumericalFluxGradient,
                                    NumericalFluxDiffusive

"""
    AtmosModel <: BalanceLaw

A `BalanceLaw` for atmosphere modeling.

# Usage

    AtmosModel(orientation, ref_state, turbulence, moisture, radiation, source,
               boundarycondition, init_state)

"""
struct AtmosModel{FT,O,RS,T,M,P,R,S,BC,IS} <: BalanceLaw
  orientation::O
  ref_state::RS
  turbulence::T
  moisture::M
  precipitation::P
  radiation::R
  source::S
  # TODO: Probably want to have different bc for state and diffusion...
  boundarycondition::BC
  init_state::IS
end

function calculate_dt(grid, model::AtmosModel, Courant_number)
    T = 290.0
    return Courant_number * min_node_distance(grid, VerticalDirection()) / soundspeed_air(T)
end

abstract type AtmosConfiguration end
struct AtmosLESConfiguration <: AtmosConfiguration end
struct AtmosGCMConfiguration <: AtmosConfiguration end

function AtmosModel{FT}(::Type{AtmosLESConfiguration};
                         orientation::O=FlatOrientation(),
                         ref_state::RS=HydrostaticState(LinearTemperatureProfile(FT(200),
                                                                                FT(280),
                                                                                FT(grav) / FT(cp_d)),
                                                                                FT(0)),
                         turbulence::T=SmagorinskyLilly{FT}(0.21),
                         moisture::M=EquilMoist(),
                         precipitation::P=NoPrecipitation(),
                         radiation::R=NoRadiation(),
                         source::S=( Gravity(),
                                     Coriolis(),
                                     GeostrophicForcing{FT}(7.62e-5, 0, 0)),
                         # TODO: Probably want to have different bc for state and diffusion...
                         boundarycondition::BC=NoFluxBC(),
                         init_state::IS=nothing) where {FT<:AbstractFloat,O,RS,T,M,P,R,S,BC,IS}
  @assert init_state ≠ nothing

  atmos = (
        orientation,
        ref_state,
        turbulence,
        moisture,
        precipitation,
        radiation,
        source,
        boundarycondition,
        init_state,
       )

  return AtmosModel{FT,typeof.(atmos)...}(atmos...)
end
function AtmosModel{FT}(::Type{AtmosGCMConfiguration};
                         orientation::O        = SphericalOrientation(),
                         ref_state::RS         = HydrostaticState(
                                                   LinearTemperatureProfile(
                                                    FT(200),
                                                    FT(280),
                                                    FT(grav) / FT(cp_d)),
                                                  FT(0)),
                         turbulence::T         = SmagorinskyLilly{FT}(0.21),
                         moisture::M           = EquilMoist(),
                         precipitation::P      = NoPrecipitation(),
                         radiation::R          = NoRadiation(),
                         source::S             = (Gravity(), Coriolis()),
                         boundarycondition::BC = NoFluxBC(),
                         init_state::IS=nothing) where {FT<:AbstractFloat,O,RS,T,M,P,R,SU,S,BC,IS}
  @assert init_state ≠ nothing
  atmos = (
        orientation,
        ref_state,
        turbulence,
        moisture,
        precipitation,
        radiation,
        source,
        boundarycondition,
        init_state,
       )

  return AtmosModel{FT,typeof.(atmos)...}(atmos...)
end

function vars_state(m::AtmosModel, FT)
  @vars begin
    ρ::FT
    ρu::SVector{3,FT}
    ρe::FT
    turbulence::vars_state(m.turbulence, FT)
    moisture::vars_state(m.moisture, FT)
    radiation::vars_state(m.radiation, FT)
  end
end
function vars_gradient(m::AtmosModel, FT)
  @vars begin
    u::SVector{3,FT}
    h_tot::FT
    turbulence::vars_gradient(m.turbulence,FT)
    moisture::vars_gradient(m.moisture,FT)
  end
end
function vars_diffusive(m::AtmosModel, FT)
  @vars begin
    ∇h_tot::SVector{3,FT}
    turbulence::vars_diffusive(m.turbulence,FT)
    moisture::vars_diffusive(m.moisture,FT)
  end
end

function vars_aux(m::AtmosModel, FT)
  @vars begin
    ∫dz::vars_integrals(m, FT)
    ∫dnz::vars_reverse_integrals(m, FT)
    coord::SVector{3,FT}
    orientation::vars_aux(m.orientation, FT)
    ref_state::vars_aux(m.ref_state,FT)
    turbulence::vars_aux(m.turbulence,FT)
    moisture::vars_aux(m.moisture,FT)
    radiation::vars_aux(m.radiation,FT)
  end
end
function vars_integrals(m::AtmosModel,FT)
  @vars begin
    radiation::vars_integrals(m.radiation,FT)
  end
end
function vars_reverse_integrals(m::AtmosModel,FT)
  @vars begin
    radiation::vars_reverse_integrals(m.radiation,FT)
  end
end

include("orientation.jl")
include("ref_state.jl")
include("turbulence.jl")
include("moisture.jl")
include("precipitation.jl")
include("radiation.jl")
include("source.jl")
include("boundaryconditions.jl")
include("linear.jl")
include("remainder.jl")

"""
    flux_nondiffusive!(m::AtmosModel, flux::Grad, state::Vars, aux::Vars,
                       t::Real)

Computes flux non-diffusive flux portion of `F` in:

```
∂Y
-- = - ∇ • (F_{adv} + F_{press} + F_{nondiff} + F_{diff}) + S(Y)
∂t
```
Where

 - `F_{adv}`      Advective flux             ; see [`flux_advective!`]@ref()
 - `F_{press}`    Pressure flux              ; see [`flux_pressure!`]@ref()
 - `F_{diff}`     Fluxes that state gradients; see [`flux_diffusive!`]@ref()
"""
@inline function flux_nondiffusive!(m::AtmosModel, flux::Grad, state::Vars,
                                    aux::Vars, t::Real)
  ρ = state.ρ
  ρinv = 1/ρ
  ρu = state.ρu
  u = ρinv * ρu

  # advective terms
  flux.ρ   = ρ * u
  flux.ρu  = ρ * u .* u'
  flux.ρe  = u * state.ρe

  # pressure terms
  p = pressure(m.moisture, m.orientation, state, aux)

  if m.ref_state isa HydrostaticState
    flux.ρu += (p-aux.ref_state.p)*I
  else
    flux.ρu += p*I
  end
  flux.ρe += u*p
  flux_radiation!(m.radiation, m, flux, state, aux, t)
  flux_moisture!(m.moisture, m, flux, state, aux, t)
end

function gradvariables!(atmos::AtmosModel, transform::Vars, state::Vars, aux::Vars, t::Real)
  ρinv = 1/state.ρ
  transform.u = ρinv * state.ρu
  transform.h_tot = total_specific_enthalpy(atmos.moisture, atmos.orientation, state, aux)

  gradvariables!(atmos.moisture, transform, state, aux, t)
  gradvariables!(atmos.turbulence, transform, state, aux, t)
end

function diffusive!(atmos::AtmosModel, diffusive::Vars, ∇transform::Grad, state::Vars, aux::Vars, t::Real)
  diffusive.∇h_tot = ∇transform.h_tot

  # diffusion terms required for SGS turbulence computations
  diffusive!(atmos.turbulence, atmos.orientation, diffusive, ∇transform, state, aux, t)
  # diffusivity of moisture components
  diffusive!(atmos.moisture, diffusive, ∇transform, state, aux, t)
end

@inline function flux_diffusive!(atmos::AtmosModel, flux::Grad, state::Vars,
                                 diffusive::Vars, aux::Vars, t::Real)
  ν, τ = turbulence_tensors(atmos.turbulence, state, diffusive, aux, t)
  D_t = (ν isa Real ? ν : diag(ν)) * inv_Pr_turb
  d_h_tot = -D_t .* diffusive.∇h_tot
  flux_diffusive!(atmos, flux, state, τ, d_h_tot)
  flux_diffusive!(atmos.moisture, flux, state, diffusive, aux, t, D_t)
end

#TODO: Consider whether to not pass ρ and ρu (not state), foc BCs reasons
@inline function flux_diffusive!(atmos::AtmosModel, flux::Grad, state::Vars,
                                 τ, d_h_tot)
  flux.ρu += τ * state.ρ
  flux.ρe += τ * state.ρu
  flux.ρe += d_h_tot * state.ρ
end

@inline function wavespeed(m::AtmosModel, nM, state::Vars, aux::Vars, t::Real)
  ρinv = 1/state.ρ
  u = ρinv * state.ρu
  return abs(dot(nM, u)) + soundspeed(m.moisture, m.orientation, state, aux)
end

function update_aux!(dg::DGModel, m::AtmosModel, Q::MPIStateArray, t::Real)
  FT = eltype(Q)
  auxstate = dg.auxstate

  if num_integrals(m, FT) > 0
    indefinite_stack_integral!(dg, m, Q, auxstate, t)
    reverse_indefinite_stack_integral!(dg, m, Q, auxstate, t)
  end

  nodal_update_aux!(atmos_nodal_update_aux!, dg, m, Q, t)

  return true
end

function atmos_nodal_update_aux!(m::AtmosModel, state::Vars, aux::Vars,
                                 t::Real)
  atmos_nodal_update_aux!(m.moisture, m, state, aux, t)
  atmos_nodal_update_aux!(m.radiation, m, state, aux, t)
  atmos_nodal_update_aux!(m.turbulence, m, state, aux, t)
end

function integral_load_aux!(m::AtmosModel, integ::Vars, state::Vars, aux::Vars)
  integral_load_aux!(m.radiation, integ, state, aux)
end

function integral_set_aux!(m::AtmosModel, aux::Vars, integ::Vars)
  integral_set_aux!(m.radiation, aux, integ)
end

function reverse_integral_load_aux!(m::AtmosModel, integ::Vars, state::Vars, aux::Vars)
  reverse_integral_load_aux!(m.radiation, integ, state, aux)
end

function reverse_integral_set_aux!(m::AtmosModel, aux::Vars, integ::Vars)
  reverse_integral_set_aux!(m.radiation, aux, integ)
end

# TODO: figure out a nice way to handle this
function init_aux!(m::AtmosModel, aux::Vars, geom::LocalGeometry)
  aux.coord = geom.coord
  atmos_init_aux!(m.orientation, m, aux, geom)
  atmos_init_aux!(m.ref_state, m, aux, geom)
  atmos_init_aux!(m.turbulence, m, aux, geom)
end

"""
    source!(m::AtmosModel, source::Vars, state::Vars, diffusive::Vars, aux::Vars, t::Real)
Computes flux `S(Y)` in:
```
∂Y
-- = - ∇ • F + S(Y)
∂t
```
"""
function source!(m::AtmosModel, source::Vars, state::Vars, diffusive::Vars, aux::Vars, t::Real)
  atmos_source!(m.source, m, source, state, diffusive, aux, t)
end

boundary_state!(nf, m::AtmosModel, x...) =
  atmos_boundary_state!(nf, m.boundarycondition, m, x...)

function init_state!(m::AtmosModel, state::Vars, aux::Vars, coords, t, args...)
  m.init_state(m, state, aux, coords, t, args...)
end

boundary_flux_diffusive!(nf::NumericalFluxDiffusive,
                         atmos::AtmosModel,
                         F,
                         state⁺, diff⁺, aux⁺, n⁻,
                         state⁻, diff⁻, aux⁻,
                         bctype, t,
                         state1⁻, diff1⁻, aux1⁻) =
  atmos_boundary_flux_diffusive!(nf, atmos.boundarycondition, atmos,
                                 F,
                                 state⁺, diff⁺, aux⁺, n⁻,
                                 state⁻, diff⁻, aux⁻,
                                 bctype, t,
                                 state1⁻, diff1⁻, aux1⁻)

end # module

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
# comparing _searchindex from base.search.jl with variants
# Author Robert Rudolph, Rudolph Consulting, Ahrensburg, Germany
import Base.ByteArray, Base._search_bloom_mask, Base._nthbyte

"""
data types supported for pattern and serarchable data.

We are investigating with respect to speed, so we restrict
to data structures with efficient byte access via position
and intentionally do not support AbstractVector and AbstractString
"""
const SearchSequence = Union{String,Vector{UInt8},Vector{Int8}}

#=
"structure used by many bloom-search-variants"
struct BloomSearch{T<:SearchSequence} <: SearchStructure
    pattern::T
    bloom_mask::UInt64
    bloom_skip::Int
    skip::Int # skip if last byte matches but another byte not
end
=#

@enum StatsFields begin
    SFinit=1        # time (ns) for initializing the search structure
    SFkernel        # time (ns) of search kernel execution
    SFloops         # no. of main loop runs of the search kernel
    SFtests         # no. of performed higher order (eg bloom) tests.
    SFskips         # no. of performed higher order (eg bloom) skips
    SFbits          # no. of bits set in the bloom filter
    SFskipgood      # no. of bytes to skip if last byte matches, others not
end

"Benchmark measurements. nothing ==> compiler strips measurement code"
const MaybeVector = Union{Nothing,Vector{Int}}

function bitcount(bitset::UInt64)
    count = 0
    while  bitset>0
        count += bitset&(1%UInt64)
        bitset = bitset >>>1
    end
    count
end

"benchmark measurement recording"
mutable struct Benchmark
    dict:: IdDict{Function,Array{Int,2}} # measurements
    size:: Int # no. of measurements per function in benchmark
    source:: String ## either UInt8,Int (alphabet end and seed) or a filename
    result::String
    function Benchmark(size::Int, source::String, result::String)
        d = IdDict{Function,Array{Int,2}}()
        new(d,size,source,result)
    end
end

# write benchmark measurements as CSV file
function Base.write(io::IO,s::Benchmark)
    for p in s.dict
        for j in Int(typemin(StatsFields)):Int(typemax(StatsFields))
            print(io,string(p.first),",",StatsFields(j))
            for v in p.second[:,j]
                print(io,", ",v)
            end
            println(io)
        end
    end
end

const nullmat = zeros(Int,0)

"record a line of measurements in the benchmark"
function record(stats::Benchmark, f::Function, i::Int, sv::MaybeVector)
    mat = get(stats.dict,f,nullmat)
    if mat===nullmat
        mat = zeros(Int,Int(stats.size),Int(typemax(StatsFields)))
        push!(stats.dict,f=>mat)
    end
    mat[i,:] = sv
end

"put statistic values from a search run into vector"
function recordcase(sv::Vector{Int},loops::Integer,bloomtests::Integer,bloomskips::Integer, bloom_bits::Integer, skipgood::Integer)
    sv[Int(SFloops)] += loops
    sv[Int(SFtests)] += bloomtests
    sv[Int(SFskips)] += bloomskips
    sv[Int(SFbits)] = bloom_bits
    sv[Int(SFskipgood)] = skipgood
end

"return stats series for function and field."
function line(stats::Benchmark,f::Function, field::StatsFields)
    m = get(stats.dict,f,nullmat)
    m[:,Int(field)]
end

"return all search functions in this package (prefix 'bloom' or suffix 'search')"
function all_search_functions()
    searches = Function[]
    for nm in names(SearchBenchmarks)
        f = eval(nm)
        if f isa Function && (endswith(string(nm),"search") || startswith(string(nm),"bloom"))
            push!(searches,f)
        end
    end
    searches
end


text:
using Scanivalve
using Test

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