text:
int_rules_5_1_6 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 5.1.6*Miscellaneous*inverse*sine =#
    @apply_utils Antiderivative((~(a') + ~(b') * asin(~c + ~(d') * ~x)) ^ ~(n'), ~x) => (1 / ~d) * Subst(Antiderivative((~a + ~b * asin(~x)) ^ ~n, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * acos(~c + ~(d') * ~x)) ^ ~(n'), ~x) => (1 / ~d) * Subst(Antiderivative((~a + ~b * acos(~x)) ^ ~n, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * asin(~c + ~(d') * ~x)) ^ ~(n'), ~x) => (1 / ~d) * Subst(Antiderivative(((~d * ~e - ~c * ~f) / ~d + (~f * ~x) / ~d) ^ ~m * (~a + ~b * asin(~x)) ^ ~n, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(a') + ~(b') * acos(~c + ~(d') * ~x)) ^ ~(n'), ~x) => (1 / ~d) * Subst(Antiderivative(((~d * ~e - ~c * ~f) / ~d + (~f * ~x) / ~d) ^ ~m * (~a + ~b * acos(~x)) ^ ~n, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(A') + ~(B') * ~x + ~(C') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asin(~c + ~(d') * ~x)) ^ ~(n'), ~x) => (1 / ~d) * Subst(Antiderivative((-(~C) / (~d) ^ 2 + (~C / (~d) ^ 2) * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~x)) ^ ~n, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~A, ~B, ~C, ~n, ~p], ~x) && (EqQ(~B * (1 - (~c) ^ 2) + 2 * ~A * ~c * ~d, 0) && EqQ(2 * ~c * ~C - ~B * ~d, 0))
    @apply_utils Antiderivative((~(A') + ~(B') * ~x + ~(C') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acos(~c + ~(d') * ~x)) ^ ~(n'), ~x) => (1 / ~d) * Subst(Antiderivative((-(~C) / (~d) ^ 2 + (~C / (~d) ^ 2) * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~x)) ^ ~n, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~A, ~B, ~C, ~n, ~p], ~x) && (EqQ(~B * (1 - (~c) ^ 2) + 2 * ~A * ~c * ~d, 0) && EqQ(2 * ~c * ~C - ~B * ~d, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(A') + ~(B') * ~x + ~(C') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * asin(~c + ~(d') * ~x)) ^ ~(n'), ~x) => (1 / ~d) * Subst(Antiderivative(((~d * ~e - ~c * ~f) / ~d + (~f * ~x) / ~d) ^ ~m * (-(~C) / (~d) ^ 2 + (~C / (~d) ^ 2) * (~x) ^ 2) ^ ~p * (~a + ~b * asin(~x)) ^ ~n, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C, ~m, ~n, ~p], ~x) && (EqQ(~B * (1 - (~c) ^ 2) + 2 * ~A * ~c * ~d, 0) && EqQ(2 * ~c * ~C - ~B * ~d, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~(A') + ~(B') * ~x + ~(C') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acos(~c + ~(d') * ~x)) ^ ~(n'), ~x) => (1 / ~d) * Subst(Antiderivative(((~d * ~e - ~c * ~f) / ~d + (~f * ~x) / ~d) ^ ~m * (-(~C) / (~d) ^ 2 + (~C / (~d) ^ 2) * (~x) ^ 2) ^ ~p * (~a + ~b * acos(~x)) ^ ~n, ~x), ~x, ~c + ~d * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C, ~m, ~n, ~p], ~x) && (EqQ(~B * (1 - (~c) ^ 2) + 2 * ~A * ~c * ~d, 0) && EqQ(2 * ~c * ~C - ~B * ~d, 0))
    @apply_utils Antiderivative(sqrt(~(a') + ~(b') * asin(~c + ~(d') * (~x) ^ 2)), ~x) => (~x * sqrt(~a + ~b * asin(~c + ~d * (~x) ^ 2)) - (sqrt(pi) * ~x * (cos(~a / (2 * ~b)) + ~c * sin(~a / (2 * ~b))) * FresnelC(sqrt(~c / (pi * ~b)) * sqrt(~a + ~b * asin(~c + ~d * (~x) ^ 2)))) / (sqrt(~c / ~b) * (cos(asin(~c + ~d * (~x) ^ 2) / 2) - ~c * sin(asin(~c + ~d * (~x) ^ 2) / 2)))) + (sqrt(pi) * ~x * (cos(~a / (2 * ~b)) - ~c * sin(~a / (2 * ~b))) * FresnelS(sqrt(~c / (pi * ~b)) * sqrt(~a + ~b * asin(~c + ~d * (~x) ^ 2)))) / (sqrt(~c / ~b) * (cos(asin(~c + ~d * (~x) ^ 2) / 2) - ~c * sin(asin(~c + ~d * (~x) ^ 2) / 2))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c) ^ 2, 1)
    @apply_utils Antiderivative(sqrt(~(a') + ~(b') * acos(1 + ~(d') * (~x) ^ 2)), ~x) => ((-2 * sqrt(~a + ~b * acos(1 + ~d * (~x) ^ 2)) * sin(acos(1 + ~d * (~x) ^ 2) / 2) ^ 2) / (~d * ~x) - (2 * sqrt(pi) * sin(~a / (2 * ~b)) * sin(acos(1 + ~d * (~x) ^ 2) / 2) * FresnelC(sqrt(1 / (pi * ~b)) * sqrt(~a + ~b * acos(1 + ~d * (~x) ^ 2)))) / (sqrt(1 / ~b) * ~d * ~x)) + (2 * sqrt(pi) * cos(~a / (2 * ~b)) * sin(acos(1 + ~d * (~x) ^ 2) / 2) * FresnelS(sqrt(1 / (pi * ~b)) * sqrt(~a + ~b * acos(1 + ~d * (~x) ^ 2)))) / (sqrt(1 / ~b) * ~d * ~x) <-- FreeQ([~a, ~b, ~d], ~x)
    @apply_utils Antiderivative(sqrt(~(a') + ~(b') * acos(-1 + ~(d') * (~x) ^ 2)), ~x) => ((2 * sqrt(~a + ~b * acos(-1 + ~d * (~x) ^ 2)) * cos((1 / 2) * acos(-1 + ~d * (~x) ^ 2)) ^ 2) / (~d * ~x) - (2 * sqrt(pi) * cos(~a / (2 * ~b)) * cos(acos(-1 + ~d * (~x) ^ 2) / 2) * FresnelC(sqrt(1 / (pi * ~b)) * sqrt(~a + ~b * acos(-1 + ~d * (~x) ^ 2)))) / (sqrt(1 / ~b) * ~d * ~x)) - (2 * sqrt(pi) * sin(~a / (2 * ~b)) * cos(acos(-1 + ~d * (~x) ^ 2) / 2) * FresnelS(sqrt(1 / (pi * ~b)) * sqrt(~a + ~b * acos(-1 + ~d * (~x) ^ 2)))) / (sqrt(1 / ~b) * ~d * ~x) <-- FreeQ([~a, ~b, ~d], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * asin(~c + ~(d') * (~x) ^ 2)) ^ ~n, ~x) => (~x * (~a + ~b * asin(~c + ~d * (~x) ^ 2)) ^ ~n + (2 * ~b * ~n * sqrt(-2 * ~c * ~d * (~x) ^ 2 - (~d) ^ 2 * (~x) ^ 4) * (~a + ~b * asin(~c + ~d * (~x) ^ 2)) ^ (~n - 1)) / (~d * ~x)) - 4 * (~b) ^ 2 * ~n * (~n - 1) * Antiderivative((~a + ~b * asin(~c + ~d * (~x) ^ 2)) ^ (~n - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ((~c) ^ 2, 1) && GtQ(~n, 1))
    @apply_utils Antiderivative((~(a') + ~(b') * acos(~c + ~(d') * (~x) ^ 2)) ^ ~n, ~x) => (~x * (~a + ~b * acos(~c + ~d * (~x) ^ 2)) ^ ~n - (2 * ~b * ~n * sqrt(-2 * ~c * ~d * (~x) ^ 2 - (~d) ^ 2 * (~x) ^ 4) * (~a + ~b * acos(~c + ~d * (~x) ^ 2)) ^ (~n - 1)) / (~d * ~x)) - 4 * (~b) ^ 2 * ~n * (~n - 1) * Antiderivative((~a + ~b * acos(~c + ~d * (~x) ^ 2)) ^ (~n - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ((~c) ^ 2, 1) && GtQ(~n, 1))
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * asin(~c + ~(d') * (~x) ^ 2)), ~x) => (-(~x) * (~c * cos(~a / (2 * ~b)) - sin(~a / (2 * ~b))) * CosIntegral((~c / (2 * ~b)) * (~a + ~b * asin(~c + ~d * (~x) ^ 2)))) / (2 * ~b * (cos(asin(~c + ~d * (~x) ^ 2) / 2) - ~c * sin(asin(~c + ~d * (~x) ^ 2) / 2))) - (~x * (~c * cos(~a / (2 * ~b)) + sin(~a / (2 * ~b))) * SinIntegral((~c / (2 * ~b)) * (~a + ~b * asin(~c + ~d * (~x) ^ 2)))) / (2 * ~b * (cos(asin(~c + ~d * (~x) ^ 2) / 2) - ~c * sin(asin(~c + ~d * (~x) ^ 2) / 2))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c) ^ 2, 1)
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * acos(1 + ~(d') * (~x) ^ 2)), ~x) => (~x * cos(~a / (2 * ~b)) * CosIntegral((~a + ~b * acos(1 + ~d * (~x) ^ 2)) / (2 * ~b))) / (sqrt(2) * ~b * sqrt(-(~d) * (~x) ^ 2)) + (~x * sin(~a / (2 * ~b)) * SinIntegral((~a + ~b * acos(1 + ~d * (~x) ^ 2)) / (2 * ~b))) / (sqrt(2) * ~b * sqrt(-(~d) * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~d], ~x)
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * acos(-1 + ~(d') * (~x) ^ 2)), ~x) => (~x * sin(~a / (2 * ~b)) * CosIntegral((~a + ~b * acos(-1 + ~d * (~x) ^ 2)) / (2 * ~b))) / (sqrt(2) * ~b * sqrt(~d * (~x) ^ 2)) - (~x * cos(~a / (2 * ~b)) * SinIntegral((~a + ~b * acos(-1 + ~d * (~x) ^ 2)) / (2 * ~b))) / (sqrt(2) * ~b * sqrt(~d * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~d], ~x)
    @apply_utils Antiderivative(1 / sqrt(~(a') + ~(b') * asin(~c + ~(d') * (~x) ^ 2)), ~x) => (-(sqrt(pi)) * ~x * (cos(~a / (2 * ~b)) - ~c * sin(~a / (2 * ~b))) * FresnelC((1 / (sqrt(~b * ~c) * sqrt(pi))) * sqrt(~a + ~b * asin(~c + ~d * (~x) ^ 2)))) / (sqrt(~b * ~c) * (cos(asin(~c + ~d * (~x) ^ 2) / 2) - ~c * sin(asin(~c + ~d * (~x) ^ 2) / 2))) - (sqrt(pi) * ~x * (cos(~a / (2 * ~b)) + ~c * sin(~a / (2 * ~b))) * FresnelS((1 / (sqrt(~b * ~c) * sqrt(pi))) * sqrt(~a + ~b * asin(~c + ~d * (~x) ^ 2)))) / (sqrt(~b * ~c) * (cos(asin(~c + ~d * (~x) ^ 2) / 2) - ~c * sin(asin(~c + ~d * (~x) ^ 2) / 2))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c) ^ 2, 1)
    @apply_utils Antiderivative(1 / sqrt(~(a') + ~(b') * acos(1 + ~(d') * (~x) ^ 2)), ~x) => (-2 * sqrt(pi / ~b) * cos(~a / (2 * ~b)) * sin(acos(1 + ~d * (~x) ^ 2) / 2) * FresnelC(sqrt(1 / (pi * ~b)) * sqrt(~a + ~b * acos(1 + ~d * (~x) ^ 2)))) / (~d * ~x) - (2 * sqrt(pi / ~b) * sin(~a / (2 * ~b)) * sin(acos(1 + ~d * (~x) ^ 2) / 2) * FresnelS(sqrt(1 / (pi * ~b)) * sqrt(~a + ~b * acos(1 + ~d * (~x) ^ 2)))) / (~d * ~x) <-- FreeQ([~a, ~b, ~d], ~x)
    @apply_utils Antiderivative(1 / sqrt(~(a') + ~(b') * acos(-1 + ~(d') * (~x) ^ 2)), ~x) => (2 * sqrt(pi / ~b) * sin(~a / (2 * ~b)) * cos(acos(-1 + ~d * (~x) ^ 2) / 2) * FresnelC(sqrt(1 / (pi * ~b)) * sqrt(~a + ~b * acos(-1 + ~d * (~x) ^ 2)))) / (~d * ~x) - (2 * sqrt(pi / ~b) * cos(~a / (2 * ~b)) * cos(acos(-1 + ~d * (~x) ^ 2) / 2) * FresnelS(sqrt(1 / (pi * ~b)) * sqrt(~a + ~b * acos(-1 + ~d * (~x) ^ 2)))) / (~d * ~x) <-- FreeQ([~a, ~b, ~d], ~x)
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * asin(~c + ~(d') * (~x) ^ 2)) ^ (3 / 2), ~x) => (-(sqrt(-2 * ~c * ~d * (~x) ^ 2 - (~d) ^ 2 * (~x) ^ 4)) / (~b * ~d * ~x * sqrt(~a + ~b * asin(~c + ~d * (~x) ^ 2))) - ((~c / ~b) ^ (3 / 2) * sqrt(pi) * ~x * (cos(~a / (2 * ~b)) + ~c * sin(~a / (2 * ~b))) * FresnelC(sqrt(~c / (pi * ~b)) * sqrt(~a + ~b * asin(~c + ~d * (~x) ^ 2)))) / (cos((1 / 2) * asin(~c + ~d * (~x) ^ 2)) - ~c * sin(asin(~c + ~d * (~x) ^ 2) / 2))) + ((~c / ~b) ^ (3 / 2) * sqrt(pi) * ~x * (cos(~a / (2 * ~b)) - ~c * sin(~a / (2 * ~b))) * FresnelS(sqrt(~c / (pi * ~b)) * sqrt(~a + ~b * asin(~c + ~d * (~x) ^ 2)))) / (cos((1 / 2) * asin(~c + ~d * (~x) ^ 2)) - ~c * sin(asin(~c + ~d * (~x) ^ 2) / 2)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c) ^ 2, 1)
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * acos(1 + ~(d') * (~x) ^ 2)) ^ (3 / 2), ~x) => (sqrt(-2 * ~d * (~x) ^ 2 - (~d) ^ 2 * (~x) ^ 4) / (~b * ~d * ~x * sqrt(~a + ~b * acos(1 + ~d * (~x) ^ 2))) - (2 * (1 / ~b) ^ (3 / 2) * sqrt(pi) * sin(~a / (2 * ~b)) * sin(acos(1 + ~d * (~x) ^ 2) / 2) * FresnelC(sqrt(1 / (pi * ~b)) * sqrt(~a + ~b * acos(1 + ~d * (~x) ^ 2)))) / (~d * ~x)) + (2 * (1 / ~b) ^ (3 / 2) * sqrt(pi) * cos(~a / (2 * ~b)) * sin(acos(1 + ~d * (~x) ^ 2) / 2) * FresnelS(sqrt(1 / (pi * ~b)) * sqrt(~a + ~b * acos(1 + ~d * (~x) ^ 2)))) / (~d * ~x) <-- FreeQ([~a, ~b, ~d], ~x)
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * acos(-1 + ~(d') * (~x) ^ 2)) ^ (3 / 2), ~x) => (sqrt(2 * ~d * (~x) ^ 2 - (~d) ^ 2 * (~x) ^ 4) / (~b * ~d * ~x * sqrt(~a + ~b * acos(-1 + ~d * (~x) ^ 2))) - (2 * (1 / ~b) ^ (3 / 2) * sqrt(pi) * cos(~a / (2 * ~b)) * cos(acos(-1 + ~d * (~x) ^ 2) / 2) * FresnelC(sqrt(1 / (pi * ~b)) * sqrt(~a + ~b * acos(-1 + ~d * (~x) ^ 2)))) / (~d * ~x)) - (2 * (1 / ~b) ^ (3 / 2) * sqrt(pi) * sin(~a / (2 * ~b)) * cos(acos(-1 + ~d * (~x) ^ 2) / 2) * FresnelS(sqrt(1 / (pi * ~b)) * sqrt(~a + ~b * acos(-1 + ~d * (~x) ^ 2)))) / (~d * ~x) <-- FreeQ([~a, ~b, ~d], ~x)
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * asin(~c + ~(d') * (~x) ^ 2)) ^ 2, ~x) => (-(sqrt(-2 * ~c * ~d * (~x) ^ 2 - (~d) ^ 2 * (~x) ^ 4)) / (2 * ~b * ~d * ~x * (~a + ~b * asin(~c + ~d * (~x) ^ 2))) - (~x * (cos(~a / (2 * ~b)) + ~c * sin(~a / (2 * ~b))) * CosIntegral((~c / (2 * ~b)) * (~a + ~b * asin(~c + ~d * (~x) ^ 2)))) / (4 * (~b) ^ 2 * (cos(asin(~c + ~d * (~x) ^ 2) / 2) - ~c * sin(asin(~c + ~d * (~x) ^ 2) / 2)))) + (~x * (cos(~a / (2 * ~b)) - ~c * sin(~a / (2 * ~b))) * SinIntegral((~c / (2 * ~b)) * (~a + ~b * asin(~c + ~d * (~x) ^ 2)))) / (4 * (~b) ^ 2 * (cos(asin(~c + ~d * (~x) ^ 2) / 2) - ~c * sin(asin(~c + ~d * (~x) ^ 2) / 2))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c) ^ 2, 1)
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * acos(1 + ~(d') * (~x) ^ 2)) ^ 2, ~x) => (sqrt(-2 * ~d * (~x) ^ 2 - (~d) ^ 2 * (~x) ^ 4) / (2 * ~b * ~d * ~x * (~a + ~b * acos(1 + ~d * (~x) ^ 2))) + (~x * sin(~a / (2 * ~b)) * CosIntegral((~a + ~b * acos(1 + ~d * (~x) ^ 2)) / (2 * ~b))) / (2 * sqrt(2) * (~b) ^ 2 * sqrt(-(~d) * (~x) ^ 2))) - (~x * cos(~a / (2 * ~b)) * SinIntegral((~a + ~b * acos(1 + ~d * (~x) ^ 2)) / (2 * ~b))) / (2 * sqrt(2) * (~b) ^ 2 * sqrt(-(~d) * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~d], ~x)
    @apply_utils Antiderivative(1 / (~(a') + ~(b') * acos(-1 + ~(d') * (~x) ^ 2)) ^ 2, ~x) => (sqrt(2 * ~d * (~x) ^ 2 - (~d) ^ 2 * (~x) ^ 4) / (2 * ~b * ~d * ~x * (~a + ~b * acos(-1 + ~d * (~x) ^ 2))) - (~x * cos(~a / (2 * ~b)) * CosIntegral((~a + ~b * acos(-1 + ~d * (~x) ^ 2)) / (2 * ~b))) / (2 * sqrt(2) * (~b) ^ 2 * sqrt(~d * (~x) ^ 2))) - (~x * sin(~a / (2 * ~b)) * SinIntegral((~a + ~b * acos(-1 + ~d * (~x) ^ 2)) / (2 * ~b))) / (2 * sqrt(2) * (~b) ^ 2 * sqrt(~d * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~d], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * asin(~c + ~(d') * (~x) ^ 2)) ^ ~n, ~x) => ((~x * (~a + ~b * asin(~c + ~d * (~x) ^ 2)) ^ (~n + 2)) / (4 * (~b) ^ 2 * (~n + 1) * (~n + 2)) + (sqrt(-2 * ~c * ~d * (~x) ^ 2 - (~d) ^ 2 * (~x) ^ 4) * (~a + ~b * asin(~c + ~d * (~x) ^ 2)) ^ (~n + 1)) / (2 * ~b * ~d * (~n + 1) * ~x)) - (1 / (4 * (~b) ^ 2 * (~n + 1) * (~n + 2))) * Antiderivative((~a + ~b * asin(~c + ~d * (~x) ^ 2)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ((~c) ^ 2, 1) && (LtQ(~n, -1) && NeQ(~n, -2)))
    @apply_utils Antiderivative((~(a') + ~(b') * acos(~c + ~(d') * (~x) ^ 2)) ^ ~n, ~x) => ((~x * (~a + ~b * acos(~c + ~d * (~x) ^ 2)) ^ (~n + 2)) / (4 * (~b) ^ 2 * (~n + 1) * (~n + 2)) - (sqrt(-2 * ~c * ~d * (~x) ^ 2 - (~d) ^ 2 * (~x) ^ 4) * (~a + ~b * acos(~c + ~d * (~x) ^ 2)) ^ (~n + 1)) / (2 * ~b * ~d * (~n + 1) * ~x)) - (1 / (4 * (~b) ^ 2 * (~n + 1) * (~n + 2))) * Antiderivative((~a + ~b * acos(~c + ~d * (~x) ^ 2)) ^ (~n + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ((~c) ^ 2, 1) && (LtQ(~n, -1) && NeQ(~n, -2)))
    @apply_utils Antiderivative(asin(~(a') * (~x) ^ ~p) ^ ~(n') / ~x, ~x) => (1 / ~p) * Subst(Antiderivative((~x) ^ ~n * cot(~x), ~x), ~x, asin(~a * (~x) ^ ~p)) <-- FreeQ([~a, ~p], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative(acos(~(a') * (~x) ^ ~p) ^ ~(n') / ~x, ~x) => (-1 / ~p) * Subst(Antiderivative((~x) ^ ~n * tan(~x), ~x), ~x, acos(~a * (~x) ^ ~p)) <-- FreeQ([~a, ~p], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative(~(u') * asin(~(c') / (~(a') + ~(b') * (~x) ^ ~(n'))) ^ ~(m'), ~x) => Antiderivative(~u * acsc(~a / ~c + (~b * (~x) ^ ~n) / ~c) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~m], ~x)
    @apply_utils Antiderivative(~(u') * acos(~(c') / (~(a') + ~(b') * (~x) ^ ~(n'))) ^ ~(m'), ~x) => Antiderivative(~u * asec(~a / ~c + (~b * (~x) ^ ~n) / ~c) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~m], ~x)
    @apply_utils Antiderivative(asin(sqrt(1 + ~(b') * (~x) ^ 2)) ^ ~(n') / sqrt(1 + ~(b') * (~x) ^ 2), ~x) => (sqrt(-(~b) * (~x) ^ 2) / (~b * ~x)) * Subst(Antiderivative(asin(~x) ^ ~n / sqrt(1 - (~x) ^ 2), ~x), ~x, sqrt(1 + ~b * (~x) ^ 2)) <-- FreeQ([~b, ~n], ~x)
    @apply_utils Antiderivative(acos(sqrt(1 + ~(b') * (~x) ^ 2)) ^ ~(n') / sqrt(1 + ~(b') * (~x) ^ 2), ~x) => (sqrt(-(~b) * (~x) ^ 2) / (~b * ~x)) * Subst(Antiderivative(acos(~x) ^ ~n / sqrt(1 - (~x) ^ 2), ~x), ~x, sqrt(1 + ~b * (~x) ^ 2)) <-- FreeQ([~b, ~n], ~x)
    @apply_utils Antiderivative(~(u') * (~f) ^ (~(c') * asin(~(a') + ~(b') * ~x) ^ ~(n')), ~x) => (1 / ~b) * Subst(Antiderivative(ReplaceAll(~u, (~x->begin
                                        #= none:1 =#
                                        #= none:3 =#
                                        -(~a) / ~b + sin(~x) / ~b
                                    end)) * (~f) ^ (~c * (~x) ^ ~n) * cos(~x), ~x), ~x, asin(~a + ~b * ~x)) <-- FreeQ([~a, ~b, ~c, ~f], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative(~(u') * (~f) ^ (~(c') * acos(~(a') + ~(b') * ~x) ^ ~(n')), ~x) => (-1 / ~b) * Subst(Antiderivative(ReplaceAll(~u, (~x->begin
                                        #= none:1 =#
                                        #= none:3 =#
                                        -(~a) / ~b + cos(~x) / ~b
                                    end)) * (~f) ^ (~c * (~x) ^ ~n) * sin(~x), ~x), ~x, acos(~a + ~b * ~x)) <-- FreeQ([~a, ~b, ~c, ~f], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative(asin(~(a') * (~x) ^ 2 + ~(b') * sqrt(~c + ~(d') * (~x) ^ 2)), ~x) => ~x * asin(~a * (~x) ^ 2 + ~b * sqrt(~c + ~d * (~x) ^ 2)) - ((~x * sqrt((~b) ^ 2 * ~d + (~a) ^ 2 * (~x) ^ 2 + 2 * ~a * ~b * sqrt(~c + ~d * (~x) ^ 2))) / sqrt(-((~x) ^ 2) * ((~b) ^ 2 * ~d + (~a) ^ 2 * (~x) ^ 2 + 2 * ~a * ~b * sqrt(~c + ~d * (~x) ^ 2)))) * Antiderivative((~x * (~b * ~d + 2 * ~a * sqrt(~c + ~d * (~x) ^ 2))) / (sqrt(~c + ~d * (~x) ^ 2) * sqrt((~b) ^ 2 * ~d + (~a) ^ 2 * (~x) ^ 2 + 2 * ~a * ~b * sqrt(~c + ~d * (~x) ^ 2))), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~b) ^ 2 * ~c, 1)
    @apply_utils Antiderivative(acos(~(a') * (~x) ^ 2 + ~(b') * sqrt(~c + ~(d') * (~x) ^ 2)), ~x) => ~x * acos(~a * (~x) ^ 2 + ~b * sqrt(~c + ~d * (~x) ^ 2)) + ((~x * sqrt((~b) ^ 2 * ~d + (~a) ^ 2 * (~x) ^ 2 + 2 * ~a * ~b * sqrt(~c + ~d * (~x) ^ 2))) / sqrt(-((~x) ^ 2) * ((~b) ^ 2 * ~d + (~a) ^ 2 * (~x) ^ 2 + 2 * ~a * ~b * sqrt(~c + ~d * (~x) ^ 2)))) * Antiderivative((~x * (~b * ~d + 2 * ~a * sqrt(~c + ~d * (~x) ^ 2))) / (sqrt(~c + ~d * (~x) ^ 2) * sqrt((~b) ^ 2 * ~d + (~a) ^ 2 * (~x) ^ 2 + 2 * ~a * ~b * sqrt(~c + ~d * (~x) ^ 2))), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~b) ^ 2 * ~c, 1)
    @apply_utils Antiderivative(asin(~u), ~x) => ~x * asin(~u) - Antiderivative(SimplifyIntegrand((~x * D(~u, ~x)) / sqrt(1 - (~u) ^ 2), ~x), ~x) <-- InverseFunctionFreeQ(~u, ~x) && Not(FunctionOfExponentialQ(~u, ~x))
    @apply_utils Antiderivative(acos(~u), ~x) => ~x * acos(~u) + Antiderivative(SimplifyIntegrand((~x * D(~u, ~x)) / sqrt(1 - (~u) ^ 2), ~x), ~x) <-- InverseFunctionFreeQ(~u, ~x) && Not(FunctionOfExponentialQ(~u, ~x))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * asin(~u)), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * (~a + ~b * asin(~u))) / (~d * (~m + 1)) - (~b / (~d * (~m + 1))) * Antiderivative(SimplifyIntegrand(((~c + ~d * ~x) ^ (~m + 1) * D(~u, ~x)) / sqrt(1 - (~u) ^ 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (NeQ(~m, -1) && (InverseFunctionFreeQ(~u, ~x) && (Not(FunctionOfQ((~c + ~d * ~x) ^ (~m + 1), ~u, ~x)) && Not(FunctionOfExponentialQ(~u, ~x)))))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * acos(~u)), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * (~a + ~b * acos(~u))) / (~d * (~m + 1)) + (~b / (~d * (~m + 1))) * Antiderivative(SimplifyIntegrand(((~c + ~d * ~x) ^ (~m + 1) * D(~u, ~x)) / sqrt(1 - (~u) ^ 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (NeQ(~m, -1) && (InverseFunctionFreeQ(~u, ~x) && (Not(FunctionOfQ((~c + ~d * ~x) ^ (~m + 1), ~u, ~x)) && Not(FunctionOfExponentialQ(~u, ~x)))))
    @apply_utils Antiderivative(~v * (~(a') + ~(b') * asin(~u)), ~x) => With([w = IntHide(~v, ~x)], Dist(~a + ~b * asin(~u), w, ~x) - ~b * Antiderivative(SimplifyIntegrand((w * D(~u, ~x)) / sqrt(1 - (~u) ^ 2), ~x), ~x) <-- InverseFunctionFreeQ(w, ~x)) <-- FreeQ([~a, ~b], ~x) && (InverseFunctionFreeQ(~u, ~x) && Not(MatchQ(~v, (~(c') + ~(d') * ~x) ^ ~(m') <-- FreeQ([c, d, m], ~x))))
    @apply_utils Antiderivative(~v * (~(a') + ~(b') * acos(~u)), ~x) => With([w = IntHide(~v, ~x)], Dist(~a + ~b * acos(~u), w, ~x) + ~b * Antiderivative(SimplifyIntegrand((w * D(~u, ~x)) / sqrt(1 - (~u) ^ 2), ~x), ~x) <-- InverseFunctionFreeQ(w, ~x)) <-- FreeQ([~a, ~b], ~x) && (InverseFunctionFreeQ(~u, ~x) && Not(MatchQ(~v, (~(c') + ~(d') * ~x) ^ ~(m') <-- FreeQ([c, d, m], ~x))))
end

max_stars_repo_path: src/rules/5 Inverse trig functions/5.1 Inverse sine/5.1.6 Miscellaneous inverse sine.jl


text:
using LanguageServer
using SymbolServer

server = LanguageServerInstance(stdin, stdout, true, expanduser("~/.julia/environments/v1.0"), "", Dict())
run(server)

ext: jl

max_stars_repo_name: pfitzseb/atom-julia-lsp-client

max_issues_repo_name: pfitzseb/atom-julia-lsp-client

max_forks_repo_name: pfitzseb/atom-julia-lsp-client


text:
abstract type AbstractExpansion{T} end

struct GradientExpansion{T,N,M} <: AbstractExpansion{T}
	x::SizedVector{N,T,Vector{T}}
	u::SizedVector{M,T,Vector{T}}
	function GradientExpansion{T}(n::Int,m::Int) where T
		new{T,n,m}(SizedVector{n}(zeros(T,n)), SizedVector{m}(zeros(T,m)))
	end
end

struct QExpansionMC#{T,N,M,P} <: AbstractExpansion{T}
	x
	u
	λ
	ux
	uu
	uλ
	xx
	xu
	xλ
	λx
	λu
	λλ
end

# TODO: Move to ALTRO

abstract type GeneralDynamicsExpansion end

"""
	DynamicsExpansion{T,N,N̄,M}

Stores the dynamics expansion for a single time instance.
For a `LieGroupModel`, it will provide access to both the state and state
error Jacobians.

# Constructors
```julia
DynamicsExpansion{T}(n0, n, m)
DynamicsExpansion{T}(n, m)
```
where `n0` is the size of the full state, and `n` is the size of the error state.

# Methods
To evaluate the dynamics Jacobians, use

	dynamics_expansion!(::Type{Q}, D::DynamicsExpansion, model, Z)

To compute the Jacobians for the error state, use

	error_expansion!(D::DynamicsExpansion, model, G)

where `G` is a vector of error-state Jacobians. These can be computed using
`RobotDynamics.state_diff_jacobian(G, model, Z)`.

# Extracting Jacobians
The Jacobians should be extracted using

	fdx, fdu = error_expansion(D::DynamicsExpansion, model)

This method will provide the error state Jacobians for `LieGroupModel`s, and
	the normal Jacobian otherwise. Both `fdx` and `fdu` are a `SizedMatrix`.
"""
struct DynamicsExpansion{T,N,N̄,M} <: GeneralDynamicsExpansion
	∇f::Matrix{T} # n × (n+m)
	∇²f::Matrix{T}  # (n+m) × (n+m)
	A_::SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}
	B_::SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}
	A::SizedMatrix{N̄,N̄,T,2,Matrix{T}}
	B::SizedMatrix{N̄,M,T,2,Matrix{T}}
	fxx::SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}
	fuu::SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}
	fux::SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}
	tmpA::SizedMatrix{N,N,T,2,Matrix{T}}
	tmpB::SizedMatrix{N,M,T,2,Matrix{T}}
	tmp::SizedMatrix{N,N̄,T,2,Matrix{T}}
	function DynamicsExpansion{T}(n0::Int, n::Int, m::Int) where T
		∇f = zeros(n0,n0+m)
		∇²f = zeros(n0+m,n0+m)
		ix = 1:n0
		iu = n0 .+ (1:m)
		A_ = view(∇f, ix, ix)
		B_ = view(∇f, ix, iu)
		A = SizedMatrix{n,n}(zeros(n,n))
		B = SizedMatrix{n,m}(zeros(n,m))
		fxx = view(∇²f, ix, ix)
		fuu = view(∇²f, iu, iu)
		fux = view(∇²f, iu, ix)
		tmpA = SizedMatrix{n0,n0}(zeros(n0,n0))
		tmpB = SizedMatrix{n0,m}(zeros(n0,m))
		tmp = zeros(n0,n)
		new{T,n0,n,m}(∇f,∇²f,A_,B_,A,B,fxx,fuu,fux,tmpA,tmpB,tmp)
	end
	function DynamicsExpansion{T}(n::Int, m::Int) where T
		∇f = zeros(n,n+m)
		∇²f = zeros(n+m,n+m)
		ix = 1:n
		iu = n .+ (1:m)
		A_ = view(∇f, ix, ix)
		B_ = view(∇f, ix, iu)
		A = SizedMatrix{n,n}(zeros(n,n))
		B = SizedMatrix{n,m}(zeros(n,m))
		fxx = view(∇²f, ix, ix)
		fuu = view(∇²f, iu, iu)
		fux = view(∇²f, iu, ix)
		tmpA = A
		tmpB = B
		tmp = zeros(n,n)
		new{T,n,n,m}(∇f,∇²f,A_,B_,A,B,fxx,fuu,fux,tmpA,tmpB,tmp)
	end
end

function save_tmp!(D::DynamicsExpansion)
	D.tmpA .= D.A_
	D.tmpB .= D.B_
end

function dynamics_expansion!(Q, D::Vector{<:DynamicsExpansion}, model::AbstractModel,
		Z::Traj)
	for k in eachindex(D)
		RobotDynamics.discrete_jacobian!(Q, D[k].∇f, model, Z[k])
		# save_tmp!(D[k])
		# D[k].tmpA .= D[k].A_  # avoids allocations later
		# D[k].tmpB .= D[k].B_
	end
end

# function dynamics_expansion!(D::Vector{<:DynamicsExpansion}, model::AbstractModel,
# 		Z::Traj, Q=RobotDynamics.RK3)
# 	for k in eachindex(D)
# 		RobotDynamics.discrete_jacobian!(Q, D[k].∇f, model, Z[k])
# 		D[k].tmpA .= D[k].A_  # avoids allocations later
# 		D[k].tmpB .= D[k].B_
# 	end
# end

function error_expansion!(D::DynamicsExpansion,G1,G2)
    mul!(D.tmp, D.tmpA, G1)
    mul!(D.A, Transpose(G2), D.tmp)
    mul!(D.B, Transpose(G2), D.tmpB)
end

@inline error_expansion(D::DynamicsExpansion, model::LieGroupModel) = D.A, D.B
@inline error_expansion(D::DynamicsExpansion, model::AbstractModel) = D.tmpA, D.tmpB

@inline DynamicsExpansion(model::AbstractModel) = DynamicsExpansion{Float64}(model)
@inline function DynamicsExpansion{T}(model::AbstractModel) where T
	n,m = size(model)
	n̄ = state_diff_size(model)
	DynamicsExpansion{T}(n,n̄,m)
end

function error_expansion!(D::Vector{<:DynamicsExpansion}, model::AbstractModel, G)
	for d in D
		save_tmp!(d)
	end
end

function error_expansion!(D::Vector{<:DynamicsExpansion}, model::LieGroupModel, G)
	for k in eachindex(D)
		save_tmp!(D[k])
		error_expansion!(D[k], G[k], G[k+1])
	end
end

# function linearize(::Type{Q}, model::AbstractModel, z::AbstractKnotPoint) where Q
# 	D = DynamicsExpansion(model)
# 	linearize!(Q, D, model, z)
# end
#
# function linearize!(::Type{Q}, D::DynamicsExpansion{<:Any,<:Any,N,M}, model::AbstractModel,
# 		z::AbstractKnotPoint) where {N,M,Q}
# 	discrete_jacobian!(Q, D.∇f, model, z)
# 	D.tmpA .= D.A_  # avoids allocations later
# 	D.tmpB .= D.B_
# 	return D.tmpA, D.tmpB
# end
#
# function linearize!(::Type{Q}, D::DynamicsExpansion, model::LieGroupModel) where Q
# 	discrete_jacobian!(Q, D.∇f, model, z)
# 	D.tmpA .= D.A_  # avoids allocations later
# 	D.tmpB .= D.B_
# 	G1 = state_diff_jacobian(model, state(z))
# 	G2 = state_diff_jacobian(model, x1)
# 	error_expansion!(D, G1, G2)
# 	return D.A, D.B
# end

"""

DynamicsExpansionMC{T,N,N̄,M,P}

Note: this is a very early implementation without thorough design thinking
Stores the dynamics expansion for a single time instance for a system represented in maximal coordinate
Originally we thought we could just expand DynamicsExpansion by adding more
matrices, however the size of additional matrices are pretty different

{T,N,N̄,M}
 T  data type
 N  size of state
 N̄  size of the error state (in case the state has Lie group elements)
 M  size of the control
 P  size of the constraint force
"""

# mutable struct DynamicsExpansionMC{T,N,N̄,M,P} <: GeneralDynamicsExpansion
# 	A::SizedMatrix{N̄,N̄,T,2,Matrix{T}}  # nxn
# 	B::SizedMatrix{N̄,M,T,2,Matrix{T}}  # nxm
# 	C::SizedMatrix{N̄,P,T,2,Matrix{T}}  # nxp
# 	G::SizedMatrix{P,N̄,T,2,Matrix{T}}  # pxn

# 	function DynamicsExpansionMC{T}(n::Int, m::Int, p::Int) where T
# 		A = SizedMatrix{n,n}(zeros(n,n))
# 		B = SizedMatrix{n,m}(zeros(n,m))
# 		C = SizedMatrix{n,p}(zeros(n,p))
# 		G = SizedMatrix{p,n}(zeros(p,n))
# 		new{T,n,n,m,p}(A,B,C,G)

# 	end

# end

mutable struct DynamicsExpansionMC{T,N,N̄,M,P} <: GeneralDynamicsExpansion
∇f::Matrix{T} # n × (n+m)
A_::SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}
B_::SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}
C_::SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}
A::SizedMatrix{N̄,N̄,T,2,Matrix{T}}
B::SizedMatrix{N̄,M,T,2,Matrix{T}}
C::SizedMatrix{N̄,P,T,2,Matrix{T}}
tmpA::SizedMatrix{N,N,T,2,Matrix{T}}
tmpB::SizedMatrix{N,M,T,2,Matrix{T}}
tmpC::SizedMatrix{N,P,T,2,Matrix{T}}
tmp::SizedMatrix{N,N̄,T,2,Matrix{T}}
G::SizedMatrix{P,N̄,T,2,Matrix{T}}  # pxn
all_partials::Matrix{T}
function DynamicsExpansionMC{T}(n0::Int, n::Int, m::Int, p::Int) where T
    ∇f = zeros(n0,n0+m+p)
    ix = 1:n0
    iu = n0 .+ (1:m)
    iλ = (n0+m) .+ (1:p)
    A_ = view(∇f, ix, ix)
    B_ = view(∇f, ix, iu)
    C_ = view(∇f, ix, iλ)
    A = SizedMatrix{n,n}(zeros(n,n))
    B = SizedMatrix{n,m}(zeros(n,m))
    C = SizedMatrix{n,p}(zeros(n,p))
    tmpA = SizedMatrix{n0,n0}(zeros(n0,n0))
    tmpB = SizedMatrix{n0,m}(zeros(n0,m))
    tmpC = SizedMatrix{n0,p}(zeros(n0,p))
    tmp = zeros(n0,n)
    G = SizedMatrix{p,n}(zeros(p,n))
	all_partials = zeros(n0, 2n0+m+p)
    new{T,n0,n,m,p}(∇f,A_,B_,C_,A,B,C,tmpA,tmpB,tmpC,tmp,G,all_partials)
end
function DynamicsExpansionMC{T}(n::Int, m::Int, p::Int) where T
    DynamicsExpansionMC{T}(n, n, m, p)
end
end


text:
using AbstractPlotting, WGLMakie, JSServe
using Colors, ImageTransformations, Markdown
using WGLMakie, JSServe, AbstractPlotting
using JSServe.DOM
using JSServe: styled_slider, @js_str

folder = joinpath(@__DIR__, "simulation")

volume_frames = map(readdir(folder)) do file
    open(joinpath(folder, file)) do io
        data = Array{Float32, 3}(undef, 130, 258, 65)
        read!(io, data)
        return data
    end
end

struct Flip <: JSServe.WidgetsBase.AbstractWidget{Bool}
    label::String
    value::Observable{Bool}
end
Flip(label) = Flip(label, Observable(false))
# Implement interface for Flip to work with offline mode
JSServe.is_widget(::Flip) = true
JSServe.value_range(flip::Flip) = (true, false)
JSServe.update_value!(flip::Flip, value) = (flip.value[] = value)

function JSServe.jsrender(flip::Flip)
    return DOM.input(
        type = "button",
        value = flip.label,
        onclick = js"{
            const value = this.last_value;
            JSServe.update_obs($(flip.value), !value);
            this.last_value = !value;
        }";
        class="p-1 rounded m-2"
    )
end

function handler(s, r)
    sl = JSServe.Slider(1:64)
    absorption = JSServe.Slider(range(1f0, stop=10f0, step=0.1))
    flip_colormap = Flip("flip colormap")
    absorption.value[] = 5
    v = map(sl) do idx
        # make volumes a bit smaller!
        restrict(volume_frames[idx])
    end
    clims = (8, 12)
    cmapa = RGBAf.(to_colormap(:thermal), 1.0)
    cmap1 = vcat(fill(RGBAf(0,0,0,0), 20), cmapa)
    cmap2 = vcat(cmapa, fill(RGBAf(0,0,0,0), 20))
    cmap = map(flip_colormap) do x
        return x ? copy(cmap1) : copy(cmap2)
    end

    scene = volume(v;
        colorrange=clims, algorithm=:absorption, absorption=absorption, colormap=cmap, show_axis=false,
        resolution=(600, 600))

    markdown = md"""

    # Ocean Simulation

    Simulation of instability of a horizontal density gradient in a rotating channel using 256x512x128 grid points running on a GPU. A similar process called baroclinic instability acting on basin-scale temperature gradients fills the oceans with eddies, especially in regions with large temperature gradients. These eddies are the primary way the ocean transports heat, carbon dioxide, organic matter, and nutrients on a large scale.

    Simulation created with [Oceananigans.jl](https://github.com/CliMA/Oceananigans.jl/) which is part of the [CliMA](https://github.com/CliMA) project.

    [source code](https://github.com/SimonDanisch/JSServe.jl/blob/master/examples/oceananigans.jl)

    ---

    absorption: $(styled_slider(absorption, absorption.value))
    time: $(styled_slider(sl, sl.value))

    $(flip_colormap)

    $(scene)

    """
    dom = DOM.div(JSServe.MarkdownCSS, JSServe.Styling, JSServe.TailwindCSS, markdown)
    # return dom
    return JSServe.record_state_map(s, dom).dom
end

# Either serve
# app = JSServe.Server(handler, "0.0.0.0", 8083)

export_path = "dev/WGLDemos/oceananigans/"
# or export to e.g. github IO
JSServe.export_standalone(handler, export_path, clear_folder=true)


text:
module NC
    using NetCDF
    using ..Args

    export get_netcdf, create_netcdf

    t_atts      = Dict("units" => "seconds")
    vel_atts    = Dict("units" => "m/s")
    dist_atts   = Dict("units" => "m")
    no_atts     = Dict()

    """
    Opens the specified NetCDF file in write mode and returns its handle.
    """
    function get_netcdf(filename :: AbstractString) :: NcFile
        nc :: NcFile = NetCDF.open(filename, mode=NC_WRITE)
        NetCDF.nc_set_fill(nc.ncid, NetCDF.NC_NOFILL, Ptr{Int32}(0))
        return nc
    end

    """
    Create the NetCDF file `filename` with dimension values `x_vals`, `y_vals` and `t_vals`.
    Also create empty output variables according to `static_var_mappings` and `dynamic_var_mappings`.
    For each key in those var mappings, a variable corresponding to the mapped output name of that key is created,
    and assigned the attributes (e.g. units) of the input variable.
    """
    function create_netcdf(filename :: AbstractString, x_vals, y_vals, t_vals,
                           static_var_mappings  :: Array{Args.VarMapping, 1},
                           dynamic_var_mappings :: Array{Args.VarMapping, 1}) :: NcFile

        isfile(filename) && rm(filename)

        if isempty(x_vals) || isempty(y_vals) || isempty(t_vals)
            throw(ArgumentError("Dimension values have to be passed when creating a NetCDF file!"))
        end

        # `var_name` is a name of an INPUT variable. Thus, we can know which units should be assigned.
        function get_units(var_name:: String)
            if match(r"[uvw]", var_name) !== nothing
                return vel_atts
            elseif match(r"[UVW]", var_name) !== nothing
                return dist_atts
            else
                return no_atts
            end
        end

        x_dim = NcDim("x",    x_vals, dist_atts)
        y_dim = NcDim("y",    y_vals, dist_atts)
        t_dim = NcDim("time", t_vals, t_atts)

        static_vars  = []
        dynamic_vars = []

        for stat_mapping ∈ static_var_mappings
            append!(static_vars,  map(in_var -> NcVar(stat_mapping[in_var], [y_dim, x_dim],        atts=get_units(in_var)), collect(keys(stat_mapping))))
        end

        for dyn_mapping ∈ dynamic_var_mappings
            append!(dynamic_vars, map(in_var -> NcVar(dyn_mapping[in_var],  [y_dim, x_dim, t_dim], atts=get_units(in_var)), collect(keys(dyn_mapping))))
        end

        nc :: NcFile = NetCDF.create(filename, static_vars..., dynamic_vars...)

        NetCDF.nc_set_fill(nc.ncid, NetCDF.NC_NOFILL, Ptr{Int32}(0))
        return nc
    end
end

max_stars_repo_name: SeisSol/SAMPLER

max_issues_repo_name: SeisSol/SAMPLER

max_forks_repo_name: SeisSol/SAMPLER


text:
include("MultirateInfinitesimalStepMethod_kernels.jl")

export MIS2, MIS3C, MIS4, MIS4a, TVDMISA, TVDMISB

"""
    TimeScaledRHS(a, b, rhs!)

When evaluate at time `t`, evaluates `rhs!` at time `a + bt`.
"""
mutable struct TimeScaledRHS{RT}
  a::RT
  b::RT
  rhs!
end

function (o::TimeScaledRHS)(dQ, Q, params, tau; increment)
  o.rhs!(dQ, Q, params, o.a + o.b*tau; increment=increment)
end

"""
MultirateInfinitesimalStep(slowrhs!, fastrhs!, fastmethod,
                           α, β, γ,
                           Q::AT; dt=0, t0=0) where {AT<:AbstractArray}

This is a time stepping object for explicitly time stepping the partitioned differential
equation given by right-hand-side functions `f_fast` and `f_slow` with the state `Q`, i.e.,

```math
  \\dot{Q} = f_fast(Q, t) + f_slow(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

The constructor builds a multirate infinitesimal step Runge-Kutta scheme
based on the provided `α`, `β` and `γ` tableaux and `fastmethod` for solving
the fast modes.

### References
    @article{KnothWensch2014,
      title={Generalized split-explicit Runge--Kutta methods for the compressible Euler equations},
      author={Knoth, Oswald and Wensch, Joerg},
      journal={Monthly Weather Review},
      volume={142},
      number={5},
      pages={2067--2081},
      year={2014}
    }
"""
mutable struct MultirateInfinitesimalStep{T, RT, AT, FS, Nstages, Nstagesm1, Nstagesm2, Nstages_sq} <: AbstractODESolver
  "time step"
  dt::RT
  "time"
  t::RT
  "storage for y_n"
  yn::AT
  "Storage for ``Y_nj - y_n``"
  ΔYnj::NTuple{Nstagesm2, AT}
  "Storage for ``f(Y_nj)``"
  fYnj::NTuple{Nstagesm1, AT}
  "Storage for offset"
  offset::AT
  "slow rhs function"
  slowrhs!
  "RHS for fast solver"
  tsfastrhs!::TimeScaledRHS{RT}
  "fast rhs method"
  fastsolver::FS
  "number of substeps per stage"
  nsubsteps::Int
  α::SArray{NTuple{2, Nstages}, RT, 2, Nstages_sq}
  β::SArray{NTuple{2, Nstages}, RT, 2, Nstages_sq}
  γ::SArray{NTuple{2, Nstages}, RT, 2, Nstages_sq}
  d::SArray{NTuple{1, Nstages}, RT, 1, Nstages}
  c::SArray{NTuple{1, Nstages}, RT, 1, Nstages}
  c̃::SArray{NTuple{1, Nstages}, RT, 1, Nstages}
  function MultirateInfinitesimalStep(slowrhs!, fastrhs!, fastmethod, nsubsteps,
                                      α, β, γ,
                                      Q::AT; dt=0, t0=0) where {AT<:AbstractArray}

    T = eltype(Q)
    RT = real(T)

    Nstages = size(α, 1)

    yn = similar(Q)
    ΔYnj = ntuple(_ -> similar(Q), Nstages-2)
    fYnj = ntuple(_ -> similar(Q), Nstages-1)
    offset = similar(Q)
    tsfastrhs! = TimeScaledRHS(RT(0), RT(0), fastrhs!)
    fastsolver = fastmethod(tsfastrhs!, Q)

    d = sum(β, dims=2)

    c = similar(d)
    for i = eachindex(c)
      c[i] = d[i]
      if i > 1
        c[i] += sum(j-> (α[i,j] + γ[i,j])*c[j], 1:i-1)
      end
    end
    c̃ = α*c

    new{T, RT, AT, typeof(fastsolver), Nstages, Nstages-1, Nstages-2, Nstages ^ 2}(RT(dt), RT(t0), yn, ΔYnj, fYnj, offset,
                                         slowrhs!, tsfastrhs!, fastsolver, nsubsteps,
                                         α, β, γ, d, c, c̃)
  end
end

# TODO almost identical functions seem to be defined for every ode solver,
# define a common one in ODEsolvers ?
function dostep!(Q, mis::MultirateInfinitesimalStep, param,
                      timeend::AbstractFloat, adjustfinalstep::Bool)
  time, dt = mis.t, mis.dt
  @assert dt > 0
  if adjustfinalstep && time + dt > timeend
    dt = timeend - time
    @assert dt > 0
  end

  dostep!(Q, mis, param, time, dt)

  if dt == mis.dt
    mis.t += dt
  else
    mis.t = timeend
  end
  return mis.t
end

function dostep!(Q, mis::MultirateInfinitesimalStep, p,
                      time::AbstractFloat, dt::AbstractFloat)
  FT = eltype(dt)
  α = mis.α
  β = mis.β
  γ = mis.γ
  yn = mis.yn
  ΔYnj = mis.ΔYnj
  fYnj = mis.fYnj
  offset = mis.offset
  d = mis.d
  c = mis.c
  c̃ = mis.c̃
  slowrhs! = mis.slowrhs!
  fastsolver = mis.fastsolver
  fastrhs! = mis.tsfastrhs!

  nstages = size(α, 1)

  copyto!(yn, Q) # first stage
  for i = 2:nstages
    slowrhs!(fYnj[i-1], Q, p, time + c[i-1]*dt, increment=false)

    threads = 256
    blocks = div(length(realview(Q)) + threads - 1, threads)
    @launch(device(Q), threads=threads, blocks=blocks,
            update!(realview(Q), realview(offset), Val(i), realview(yn), map(realview, ΔYnj[1:i-2]), map(realview, fYnj[1:i-1]),
            α[i,:], β[i,:], γ[i,:], d[i], dt))

    fastrhs!.a = time + c̃[i]*dt
    fastrhs!.b = (c[i] - c̃[i]) / d[i]

    τ = zero(FT)
    dτ = d[i] * dt / mis.nsubsteps
    updatetime!(fastsolver, τ)
    updatedt!(fastsolver, dτ)
    # TODO: we want to be able to write
    #   solve!(Q, fastsolver, p; numberofsteps = mis.nsubsteps)  #(1c)
    # especially if we want to use StormerVerlet, but need some way to pass in `offset`
    for k = 1:mis.nsubsteps
      dostep!(Q, fastsolver, p, τ, dτ, FT(1), realview(offset))  #(1c)
      τ += dτ
    end
  end
end

function MIS2(slowrhs!, fastrhs!, fastmethod, nsubsteps, Q::AT; dt=0, t0=0) where {AT <: AbstractArray}
  α = [0 0              0              0;
       0 0              0              0;
       0 0.536946566710 0              0;
       0 0.480892968551 0.500561163566 0]

  β = [ 0                0                0              0;
        0.126848494553   0                0              0;
       -0.784838278826   1.37442675268    0              0;
       -0.0456727081749 -0.00875082271190 0.524775788629 0]

  γ = [0  0               0              0;
       0  0               0              0;
       0  0.652465126004  0              0;
       0 -0.0732769849457 0.144902430420 0]

  MultirateInfinitesimalStep(slowrhs!, fastrhs!, fastmethod, nsubsteps,  α, β, γ, Q; dt=dt, t0=t0)
end

function MIS3C(slowrhs!, fastrhs!, fastmethod, nsubsteps,  Q::AT; dt=0, t0=0) where {AT <: AbstractArray}
  α = [0 0              0              0;
       0 0              0              0;
       0 0.589557277145 0              0;
       0 0.544036601551 0.565511042564 0]

  β = [ 0                 0                0              0;
        0.397525189225    0                0              0;
       -0.227036463644    0.624528794618   0              0;
       -0.00295238076840 -0.270971764284   0.671323159437 0]

  γ = [0  0               0               0;
       0  0               0               0;
       0  0.142798786398  0               0;
       0 -0.0428918957402 0.0202720980282 0]

  MultirateInfinitesimalStep(slowrhs!, fastrhs!, fastmethod, nsubsteps,  α, β, γ, Q; dt=dt, t0=t0)
end

function MIS4(slowrhs!, fastrhs!, fastmethod, nsubsteps,  Q::AT; dt=0, t0=0) where {AT <: AbstractArray}
  α = [0 0              0              0              0;
       0 0              0              0              0;
       0 0.914092810304 0              0              0;
       0 1.14274417397 -0.295211246188 0              0;
       0 0.112965282231 0.337369411296 0.503747183119 0]

  β = [ 0                0               0              0              0;
        0.136296478423   0               0              0              0;
        0.280462398979  -0.0160351333596 0              0              0;
        0.904713355208  -1.04011183154   0.652337563489 0              0;
        0.0671969845546 -0.365621862610 -0.154861470835 0.970362444469 0]

  γ = [0  0              0               0              0;
       0  0              0               0              0;
       0  0.678951983291 0               0              0;
       0 -1.38974164070  0.503864576302  0              0;
       0 -0.375328608282 0.320925021109 -0.158259688945 0]

  MultirateInfinitesimalStep(slowrhs!, fastrhs!, fastmethod, nsubsteps,  α, β, γ, Q; dt=dt, t0=t0)
end

function MIS4a(slowrhs!, fastrhs!, fastmethod, nsubsteps,  Q::AT; dt=0, t0=0) where {AT <: AbstractArray}
  α = [0 0                     0                   0                   0;
       0 0                     0                   0                   0;
       0 0.52349249922385610   0                   0                   0;
       0 1.1683374366893629   -0.75762080241712637 0                   0;
       0 -0.036477233846797109 0.56936148730740477 0.47746263002599681 0]

  # β[5,1] in the paper is incorrect
  # the correct value is used below (from authors)
  β = [ 0                    0                   0                   0                   0;
        0.38758444641450318  0                   0                   0                   0;
       -0.025318448354142823 0.38668943087310403 0                   0                   0;
        0.20899983523553325 -0.45856648476371231 0.43423187573425748 0                   0;
       -0.10048822195663100 -0.46186171956333327 0.83045062122462809 0.27014914900250392 0]

  γ = [0  0                    0                    0                    0;
       0  0                    0                    0                    0;
       0  0.13145089796226542  0                    0                    0;
       0 -0.36855857648747881  0.33159232636600550  0                    0;
       0 -0.065767130537473045 0.040591093109036858 0.064902111640806712 0]

  MultirateInfinitesimalStep(slowrhs!, fastrhs!, fastmethod, nsubsteps, α, β, γ, Q; dt=dt, t0=t0)
end

function TVDMISA(slowrhs!, fastrhs!, fastmethod, nsubsteps,  Q::AT; dt=0, t0=0) where {AT <: AbstractArray}
  FT = eltype(Q)
  RT = real(FT)

  α = [0 0                    0                  0;
       0 0                    0                  0;
       0 0.1946360605647457   0                  0;
       0 0.3971200136786614   0.2609434606211801 0]

  β = [ 0                    0                   0                   0;
        RT(2//3)             0                   0                   0;
       -0.28247174703488398  RT(4//9)            0                   0;
       -0.31198081960042401  0.18082737579913699 RT(9//16)           0]

  γ = [0  0                    0                    0;
       0  0                    0                    0;
       0  0.5624048933209129   0                    0;
       0  0.4408467475713277  -0.2459300561692391   0]

  MultirateInfinitesimalStep(slowrhs!, fastrhs!, fastmethod, nsubsteps, α, β, γ, Q; dt=dt, t0=t0)
end

function TVDMISB(slowrhs!, fastrhs!, fastmethod, nsubsteps,  Q::AT; dt=0, t0=0) where {AT <: AbstractArray}
  FT = eltype(Q)
  RT = real(FT)

  α = [0 0                     0                   0;
       0 0                     0                   0;
       0 0.42668232863311001   0                   0;
       0 0.26570779016173801   0.41489966891866698 0]

  β = [ 0                    0                   0                   0;
        RT(2//3)             0                   0                   0;
       -0.25492859100078202  RT(4//9)            0                   0;
       -0.26452517179288798  0.11424084424766399 RT(9//16)           0]

  γ = [0  0                    0                    0;
       0  0                    0                    0;
       0  0.28904389120139701  0                    0;
       0  0.45113560071334202 -0.25006656847591002  0]

  MultirateInfinitesimalStep(slowrhs!, fastrhs!, fastmethod, nsubsteps, α, β, γ, Q; dt=dt, t0=t0)
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
module SialiaCLI

using ArgParse

include("./DockerUtils.jl")
include("./Generator.jl")
include("./Builder.jl")
include("./Runner.jl")
include("./Deployer.jl")

function julia_main()::Cint
  try
    real_main()
  catch
      Base.invokelatest(Base.display_error, Base.catch_stack())
      return 1
  end
  return 0
end

function real_main()
  parsed_args = parse_cli_args()

  if parsed_args[:_COMMAND_] == :build
    Builder.execute(parsed_args)
  elseif parsed_args[:_COMMAND_] == :generate
    Generator.execute(parsed_args)
  elseif parsed_args[:_COMMAND_] == :run
    Runner.execute(parsed_args)
  elseif parsed_args[:_COMMAND_] == :deploy
    Deployer.execute(parsed_args)
  end
end

include("./parseargs.jl")

if abspath(PROGRAM_FILE) == @__FILE__
  real_main()
end

end # module


text:
# Parameters

r_III = 16.0 # radius of region III. Other regions are derived from this
crack_surface = [0, 0, 1] # y-diretion, the surface which is opened
crack_front = [1, -1, 0] # z-direction, crack front line
relax_elasticity = false # if true, C_ij matrix computed with internal relaxation
relax_surface = false # if true, surface energy computed with internal relaxation

eam = EAM("w_fs_eam_core.fs")

unitcell = bulk(:W, cubic=true)
variablecell!(unitcell)
set_calculator!(unitcell, eam)

res = optimize(x -> energy(unitcell, x), x -> gradient(unitcell, x),
               dofs(unitcell), inplace=false)

##

# relaxed lattice constant
alat = unitcell.cell[1, 1]
@show alat
unitcell = bulk(:W, cubic=true)
set_cell!(unitcell, alat * I(3))
set_calculator!(unitcell, eam);

##

# 6x6 elastic constant matrix
C = voigt_moduli(unitcell)
C11, C12, C44 = cubic_moduli(C, tol=1e-2)
@show C11, C12, C44;
##

shift = 2.0

make_bulk() = Atoms(ASEAtoms(ase_lattice_cubic.BodyCenteredCubic(symbol="W",
                            latticeconstant=alat,
                            directions=[[1,-1,0],[1,1,0],[0,0,1]]) * (1, 1, 10)))
bulk_at = make_bulk()
set_calculator!(bulk_at, eam)
surface = make_bulk() #copy(bulk_at) # must be a better way to do this in JuLIP
X = positions(surface) |> mat
X[3, :] .+= shift
set_positions!(surface, X)
wrap_pbc!(surface)
c = Matrix(surface.cell)
c[3, :] += [0.0, 0.0, 10.0]
set_cell!(surface, c)
fixedcell!(surface)
set_calculator!(surface, eam)
area = norm(cross(bulk_at.cell[:, 1], bulk_at.cell[:, 2]))

if relax_surface
    rattle!(surface, 0.1)
    res = optimize(x -> energy(surface, x), x -> gradient(surface, x),
                dofs(surface), inplace=false)
    print(res)
end
γ = (energy(surface) - energy(bulk_at)) / (2 * area)

@show γ, γ / (ase_units.J / ase_units.m^2)

##

# build the final crystal in correct orientation
r_cut = cutoff(eam) # cutoff distance for the potential
r_I = r_III - 1.0 * r_cut

n = [2 * floor(Int, (r_III + r_cut)/ alat),
     2 * floor(Int, (r_III + r_cut)/ alat) - 1, 1]
@show n

cryst = clusters.bcc("W", alat, n,
                      crack_surface=crack_surface,
                      crack_front=crack_front)

# mark regions and convert to JuLIP
clust = clusters.set_regions(cryst, r_I, r_cut, r_III)
region = clust.arrays["region"]
clust = Atoms(ASEAtoms(clust))

X = positions(clust) |> mat

#X0 = copy(X)

set_calculator!(clust, eam);

#### shifting the domain, finding atoms in each region
Xm = diag(clust.cell) ./ 2

set_positions!(clust,[x .- Xm for x in clust.X]);

#set_positions!(cluster,[x .- cluster.X[1] for x in cluster.X]);
;

X = clust.X
II = sortperm(norm.(X))
X = X[II]
set_positions!(clust,X)
region = region[II]

N1 = findfirst(region .== 2) - 1

#NL = neighbourlist(cluster);
NL = pairs(clust,cutoff(eam)).nlist;

I_R1 = 1:N1

I_R2 = []
for i in I_R1
    for jj in NL.first[i]:(NL.first[i+1] - 1)
    push!(I_R2,NL.j[jj])
    end
end
I_R2 = setdiff(sort(unique(I_R2)),1:N1);

I_R3 = []
for i in I_R2
    for jj in NL.first[i]:(NL.first[i+1] - 1)
    push!(I_R3,NL.j[jj])
    end
end
I_R3 = setdiff(sort(unique(I_R3)),[1:N1;I_R2]);

I_R4 = setdiff(1:length(clust.X),[I_R1;I_R2;I_R3])
;

clust["N1"] = N1
clust["I_R1"] = I_R1
clust["I_R2"] = I_R2
clust["I_R3"] = I_R3
clust["I_R4"] = I_R4
;

# object for computing the CLE displacments and gradients
C = voigt_moduli(C11, C12, C44)
Crot = rotate_elastic_moduli(C, rotation_matrix(crack_surface, crack_front))

# old Python code - we fix the rotated elastic constants
# to avoid the disrecpancy in Voigt conversion
crk = crack.CubicCrystalCrack(crack_surface, crack_front,
                              Crot=Crot)

# new Julia code
rac = RectilinearAnisotropicCrack(PlaneStrain(), C11, C12, C44,
                                   crack_surface, crack_front)

k_G1 = crk.k1g(γ)
k_G2 = k1g(rac, γ)
;

# # displacement fields with both approaches
# u1, v1 = crk.displacements(X[1, :], X[2, :], 0.0, 0.0, 1.0)
# u2, v2 = displacements(rac, X[1, :] .- 0.0, X[2, :] .- 0.0)

# @assert maximum(abs.(u1 - u2)) < 1e-6
# @assert maximum(abs.(v1 - v2)) < 1e-6

u, ∇u = u_CLE(rac, clust, 0.0, 0.0);


text:
include("../gcc_common.jl")
build_and_upload_gcc(v"12.1.0")


text:
UI = UncertainValueDataset(example_uidxs)
UV = UncertainValueDataset(example_uvals)
UIV = UncertainIndexValueDataset(UI, UV)

n = 3
@test resample(UIV) isa Tuple{Vector}
@test resample(UIV, n) isa Vector{Tuple{Vector}}


text:
using Syslogs
using Test
using Sockets
import Libdl
using Distributed
using Base: BufferStream

@eval Syslogs begin
    UDP_PORT = 8080
    TCP_PORT = 8080
end

include("helpers.jl")

@testset "Syslog" begin
    @testset "Local" begin
        @info("Local Tests")
        io = Syslog()
        println(io, "info", "foobar")
    end

    @testset "Remote" begin
        logs = filter(!isempty, split(read("scripts/output.log", String), "\n\0"))
        test_logs = map(logs) do s
            input = replace(s, r"^<\d+>" => "")
            level = last(split(first(split(s, ':'))))
            (level, input, s * "\0")
        end

        @testset "Invalid Stream" begin
            io = Syslog((ip"127.0.0.1", 8080), :user, BufferStream())
            @test_throws ArgumentError println(io, :info, "foobar")
        end

        @testset "UDP" begin
            @info("UDP Tests")
            @testset "Simple" begin
                r = udp_srv(8080)
                io = Syslog(ip"127.0.0.1", 8080; tcp=false)
                println(io, "info", "foobar")
                s = fetch(r)
                @test occursin("foobar", s)
                close(io)
            end

            @testset "SysLogHandler" begin
                @testset "$input" for (level, input, output) in test_logs
                    r = udp_srv(8080)
                    io = Syslog(ip"127.0.0.1", 8080; tcp=false)
                    @info("Sending (UDP): $level, $input")
                    println(io, level, input)
                    s = fetch(r)
                    @test strip(s) == strip(output)
                    close(io)

                    r = udp_srv(8080)
                    io = Syslog(ip"127.0.0.1"; tcp=false)
                    @info("Sending (UDP): $level, $input")
                    println(io, level, input)
                    s = fetch(r)
                    @test strip(s) == strip(output)
                    close(io)
                end
            end
        end

        @testset "TCP" begin
            @info("TCP Tests")

            @testset "Simple" begin
                r = tcp_srv(8080)
                io = Syslog(ip"127.0.0.1", 8080; tcp=true)
                println(io, "info", "foobar")
                s = fetch(r)
                @test occursin("foobar", s)
                flush(io)
                close(io)
            end

            @testset "SysLogHandler" begin
                @testset "$input" for (level, input, output) in test_logs
                    r = tcp_srv(8080)
                    @info("Sending (TCP): $level, $input")
                    io = Syslog(ip"127.0.0.1", 8080; tcp=true)
                    println(io, level, input)
                    s = fetch(r)
                    @test s == output
                    close(io)
                end
            end

            @testset "Disonnect" begin
                serv = listen(ip"127.0.0.1", 8080)
                io = Syslog(ip"127.0.0.1", 8080; tcp=true)
                close(io)
                log(io, "info", "foobar\n")
                flush(io)
                close(io)
                close(serv)
            end
        end
    end

    @testset "Libc" begin
        @info("Libc Tests")
        Syslogs.openlog("syslog", 0, Syslogs.LOG_USER)
        Syslogs.closelog()
    end
end


text:
function next_normalized_double(obj::NormalizedRandomGenerator)
    return jcall(obj, "nextNormalizedDouble", jdouble, ())
end


text:
"""
    RecombinedBSplineBasis{k, T}

Functional basis defined from the recombination of a [`BSplineBasis`](@ref)
in order to satisfy certain homogeneous boundary conditions (BCs).

# Extended help

The basis recombination technique is a common way of applying BCs in Galerkin
methods. It is described for instance in Boyd 2000 (ch. 6), in the context of
a Chebyshev basis. In this approach, the original basis, ``\\{b_i(x), 1 ≤ i ≤ N\\}``,
is "recombined" into a new basis, ``\\{ϕ_j(x), 1 ≤ j ≤ M\\}``, so that each
basis function ``ϕ_j`` individually satisfies the chosen BCs.

The length ``M`` of the recombined basis is always smaller than the length ``N``
of the original basis.
The difference, ``δ = N - M``, is equal to the number of boundary conditions.
In the simplest (and most common) case, a single BC is applied on each boundary,
leading to ``δ = 2``.
More generally, as described further below, it is possible to simultaneously
impose different BCs, which further decreases the number of degrees of freedom
(increasing ``δ``).

Thanks to the local support of B-splines, basis recombination involves only a
little portion of the original B-spline basis. For instance, since there is only
one B-spline that is non-zero at each boundary, removing that function from the
basis is enough to apply homogeneous Dirichlet BCs. Imposing BCs for derivatives
is slightly more complex, but still possible.

Note that, when combining basis recombination with [collocation methods](@ref
collocation-api), there must be **no** collocation points at the boundaries,
or the resulting collocation matrices may not be invertible.

## Order of the boundary condition

In this section, we consider the simplest case where a single homogeneous
boundary condition, ``\\mathrm{d}^n u / \\mathrm{d}x^n = 0``, is to be satisfied
by the basis.

The recombined basis requires the specification of a `Derivative` object
determining the order of the homogeneous BCs to be applied at the two
boundaries.
Linear combinations of `Derivative`s are also supported.
The order of the B-spline needs to be ``k ≥ n + 1``, since a B-spline
of order ``k`` is a ``C^{k - 1}``-continuous function (except on the knots
where it is ``C^{k - 1 - p}``, with ``p`` the knot multiplicity).

Some usual choices are:

- `Derivative(0)` sets homogeneous [Dirichlet
  BCs](https://en.wikipedia.org/wiki/Dirichlet_boundary_condition) (``u = 0`` at
  the boundaries) by removing the first and last B-splines, i.e. ``ϕ_1 = b_2``;

- `Derivative(1)` sets homogeneous [Neumann
  BCs](https://en.wikipedia.org/wiki/Neumann_boundary_condition) (``u' = 0`` at
  the boundaries) by adding the two first (and two last) B-splines, i.e. ``ϕ_1 =
  b_1 + b_2``.

- more generally, `α Derivative(0) + β Derivative(1)` sets homogeneous [Robin
  BCs](https://en.wikipedia.org/wiki/Robin_boundary_condition) by defining
  ``ϕ_1`` as a linear combination of ``b_1`` and ``b_2``.
  Here it's important to note that `Derivative(1)` denotes the [normal
  derivative](https://en.wikipedia.org/wiki/Directional_derivative#Normal_derivative)
  at the boundary, ``\\frac{∂u}{∂n}``, which is equal to ``-\\frac{∂u}{∂x}`` on
  the left boundary.

Higher order BCs are also possible.
For instance, `Derivative(2)` recombines the first three B-splines into two
basis functions that satisfy ``ϕ_1'' = ϕ_2'' = 0`` at the left boundary, while
ensuring that lower and higher-order derivatives keep degrees of freedom at the
boundary.
Note that simply adding the first three B-splines, as in ``ϕ_1 = b_1 + b_2 +
b_3``, makes the first derivative vanish as well as the second one, which is
unwanted.

For `Derivative(2)`, the chosen solution is to set ``ϕ_i = α_i b_i + β_i b_{i +
1}`` for ``i ∈ \\{1, 2\\}``.
The ``α_i`` and ``β_i`` coefficients are chosen such that ``ϕ_i'' = 0`` at the
boundary.
Moreover, they satisfy the (somewhat arbitrary) constraint ``α_i + β_i = 2`` for
each ``i``, for consistency with the Neumann case described above.
This generalises to higher order BCs.
Note that, since each boundary function ``ϕ_i`` is defined from only two
neighbouring B-splines, its local support stays minimal, hence preserving the
small bandwidth of the resulting matrices.

Finally, note that in the current implementation, it is not possible to impose
different boundary conditions on both boundaries.

## Multiple boundary conditions

As an option, the recombined basis may simultaneously satisfy homogeneous BCs of
different orders. In this case, a tuple of `Derivative`s must be passed.

For more details on the supported combinations of BCs, see the different
`RecombinedBSplineBasis` constructors documented further below.

---

    RecombinedBSplineBasis(op::AbstractDifferentialOp, B::BSplineBasis)

Construct `RecombinedBSplineBasis` from B-spline basis `B`, satisfying
homogeneous boundary conditions associated to the given differential operator.

For instance, `op = Derivative(0)` and `op = Derivative(1)` correspond to
homogeneous Dirichlet and Neumann BCs, respectively.

Linear combinations of differential operators are also supported.
For instance, `op = Derivative(0) + λ Derivative(1)` corresponds to homogeneous
Robin BCs.

Higher-order derivatives are also allowed, being only limited by the order of
the B-spline basis.

## Examples

```jldoctest RecombinedBSplineBasis
julia> B = BSplineBasis(BSplineOrder(4), -1:0.2:1)
13-element BSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]

julia> R_neumann = RecombinedBSplineBasis(Derivative(1), B)
11-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 BCs left:  (D{1},)
 BCs right: (D{1},)

julia> R_robin = RecombinedBSplineBasis(Derivative(0) + 3Derivative(1), B)
11-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 BCs left:  (D{0} + 3 * D{1},)
 BCs right: (D{0} + 3 * D{1},)
```

---

    RecombinedBSplineBasis(ops::Tuple{Vararg{AbstractDifferentialOp}},
                           B::BSplineBasis)

Construct `RecombinedBSplineBasis` simultaneously satisfying homogeneous BCs
associated to multiple differential operators.

Currently, the following specific combinations of differential operators are
supported:

1. all derivatives up to order `m`:

        ops = (Derivative(0), ..., Derivative(m))

   This boundary condition is obtained by removing the first `m + 1` B-splines
   from the original basis.

   For instance, if `(Derivative(0), Derivative(1))` is passed, then the basis
   simultaneously satisfies homogeneous Dirichlet and Neumann BCs at the two
   boundaries.
   The resulting basis is ``ϕ_1 = b_3, ϕ_2 = b_4, …, ϕ_{N - 4} = b_{N - 2}``.

2. an extension of the above, with an additional differential operator of order
   `n` at the end:

        ops = (Derivative(0), ..., Derivative(m), D(n))

   The operator `D(n)` may be a `Derivative`, or a linear combination of
   derivatives.
   The only restriction is that its maximum degree must satisfy `n ≥ m + 1`.

   One example is the combination of homogeneous Dirichlet BCs, ``u = 0`` on the
   boundaries, with Robin BCs for the derivative, ``u' + λ u'' = 0``, which
   corresponds to `ops = (Derivative(0), Derivative(1) + λ Derivative(2))`.

In all cases, the degrees of the differential operators must be in increasing
order. For instance, `ops = (Derivative(1), Derivative(0))` fails with an error.

## Examples

```jldoctest RecombinedBSplineBasis
julia> ops = (Derivative(0), Derivative(1));

julia> R1 = RecombinedBSplineBasis(ops, B)
9-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 BCs left:  (D{0}, D{1})
 BCs right: (D{0}, D{1})

julia> ops = (Derivative(0), Derivative(1) - 4Derivative(2));

julia> R2 = RecombinedBSplineBasis(ops, B)
9-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 BCs left:  (D{0}, D{1} + -4 * D{2})
 BCs right: (D{0}, D{1} + -4 * D{2})
```
"""
struct RecombinedBSplineBasis{
            k, T, Parent <: BSplineBasis{k,T},
            DiffOps <: Tuple{Vararg{AbstractDifferentialOp}},
            RMatrix <: RecombineMatrix{Q,DiffOps} where Q,
        } <: AbstractBSplineBasis{k,T}
    B :: Parent   # original B-spline basis
    ops :: DiffOps  # list of differential operators for BCs
    M :: RMatrix  # basis recombination matrix

    function RecombinedBSplineBasis(ops::Tuple{Vararg{AbstractDifferentialOp}},
                                    B::BSplineBasis{k,T}) where {k,T}
        Parent = typeof(B)
        M = RecombineMatrix(ops, B)
        RMatrix = typeof(M)
        Ops = typeof(ops)
        new{k,T,Parent,Ops,RMatrix}(B, ops, M)
    end

    RecombinedBSplineBasis(op::AbstractDifferentialOp, args...) =
        RecombinedBSplineBasis((op, ), args...)
end

Base.:(==)(A::RecombinedBSplineBasis, B::RecombinedBSplineBasis) =
    A === B ||
    constraints(A) == constraints(B) &&
    parent(A) == parent(B)

function Base.summary(io::IO, R::RecombinedBSplineBasis)
    summary_basis(io, R)
    cl, cr = constraints(R)
    print(io, ", BCs {left => ", cl, ", right => ", cr, "}")
    nothing
end

function Base.show(io::IO, R::RecombinedBSplineBasis)
    summary_basis(io, R)
    cl, cr = constraints(R)
    let io = IOContext(io, :compact => true, :limit => true)
        print(io, "\n knots: ", knots(R))
        print(io, "\n BCs left:  ", cl)
        print(io, "\n BCs right: ", cr)
    end
end

"""
    RecombinedBSplineBasis(ops, args...; kwargs...)

Construct [`RecombinedBSplineBasis`](@ref) from B-spline basis, satisfying
homogeneous boundary conditions associated one or more differential operators.

This variant does not require a previously constructed [`BSplineBasis`](@ref).
Arguments are passed to the `BSplineBasis` constructor.
"""
RecombinedBSplineBasis(order, args...; kwargs...) =
    RecombinedBSplineBasis(order, BSplineBasis(args...; kwargs...))

"""
    parent(R::RecombinedBSplineBasis)

Get original B-spline basis.
"""
Base.parent(R::RecombinedBSplineBasis) = R.B

"""
    recombination_matrix(R::AbstractBSplineBasis)

Get [`RecombineMatrix`](@ref) associated to the recombined basis.

For non-recombined bases such as [`BSplineBasis`](@ref), this returns the
identity matrix (`LinearAlgebra.I`).
"""
recombination_matrix(R::RecombinedBSplineBasis) = R.M
recombination_matrix(B::AbstractBSplineBasis) = LinearAlgebra.I

"""
    length(R::RecombinedBSplineBasis)

Returns the number of functions in the recombined basis.
"""
@inline Base.length(R::RecombinedBSplineBasis) =
    length(parent(R)) - sum(num_constraints(R))

boundaries(R::RecombinedBSplineBasis) = boundaries(parent(R))

knots(R::RecombinedBSplineBasis) = knots(parent(R))
order(::Type{<:RecombinedBSplineBasis{k}}) where {k} = k
Base.eltype(::Type{RecombinedBSplineBasis{k,T}}) where {k,T} = T

"""
    constraints(R::AbstractBSplineBasis) -> (left, right)
    constraints(A::RecombineMatrix) -> (left, right)

Return the constraints (homogeneous boundary conditions) that the basis
satisfies on each boundary.

Constraints are returned as a tuple `(left, right)` indicating the BCs that are
satisfied on each boundary. Each element is a tuple of differential operators
specifying the BCs.

For example, if both Dirichlet and Neumann BCs are satisfied on the left
boundary, then `left = (Derivative(0), Derivative(1))`.

For non-recombined bases such as [`BSplineBasis`](@ref), this returns a tuple of
empty tuples: `((), ())`, since no BCs are satisfied on either boundary.
"""
constraints(R::AbstractBSplineBasis) = constraints(recombination_matrix(R))

"""
    num_constraints(R::AbstractBSplineBasis) -> (Int, Int)
    num_constraints(A::RecombineMatrix) -> (Int, Int)

Returns the number of constraints (number of BCs to satisfy) on each boundary.

For instance, if `R` simultaneously satisfies Dirichlet and Neumann boundary
conditions on each boundary, this returns `(2, 2)`.

Note that for non-recombined bases such as [`BSplineBasis`](@ref), the number of
constraints is zero, and this returns `(0, 0)`.
"""
num_constraints(B::AbstractBSplineBasis) = num_constraints(recombination_matrix(B))

"""
    num_recombined(R::AbstractBSplineBasis) -> (Int, Int)
    num_recombined(A::RecombineMatrix) -> (Int, Int)

Returns the number of recombined functions in the recombined basis for each
boundary.

For instance, if `R` satisfies Neumann boundary conditions on both boundaries,
then only the first and last basis functions are different from the original
B-spline basis, e.g. ``ϕ_1 = b_1 + b_2``, and this returns `(1, 1)`.

For non-recombined bases such as [`BSplineBasis`](@ref), this returns zero.
"""
num_recombined(R::AbstractBSplineBasis) = num_recombined(recombination_matrix(R))

function nonzero_in_segment(R::RecombinedBSplineBasis, n)
    # Same as the original B-spline basis on the interval n - δl.
    δl = num_constraints(R)[1]
    nonzero_in_segment(parent(R), n - δl; N = length(R))
end

function support(R::RecombinedBSplineBasis, j)
    δl = num_constraints(R)[1]
    support(parent(R), j + δl)
end

@propagate_inbounds function evaluate(R::RecombinedBSplineBasis, j, args...)
    B = parent(R)
    A = recombination_matrix(R)
    is = nzrows(A, j)
    @assert !isempty(is)

    i1 = is[1]
    ϕ = evaluate(B, i1, args...)
    T = typeof(ϕ)

    if length(is) == 1  # avoid accessing A[i1, j] (we know it's 1)
        # @assert A[i1, j] == 1
        return ϕ
    end

    ϕ::T *= A[i1, j]

    for i ∈ is[2:end]
        ϕ::T += A[i, j] * evaluate(B, i, args...)
    end

    ϕ::T
end


text:
"""
    above_diagonal_loop(N, itersym, imax, args...)

Generate `N` Cartesian loops for the iteration variables
`itersym_{1:N}`, where `itersym_N ∈ 1:imax`, `itersym_{N-1} ∈
itersym_N+1:imax`, etc, i.e. above the hyper-diagonal of the
`N`-dimensional hypercube with the side `imax`. `args...` is passed on
to `Base.Cartesian._nloops`. `above_diagonal_loop` is nestable.

# Examples

```jldoctest
julia> @above_diagonal_loop 2 i 3 begin
           println("==================================")
           println("i = ", Base.Cartesian.@ntuple 2 i)
           @above_diagonal_loop 2 j 3 begin
               println("j = ", Base.Cartesian.@ntuple 2 j)
           end
       end
==================================
i = (2, 1)
j = (2, 1)
j = (3, 1)
j = (3, 2)
==================================
i = (3, 1)
j = (2, 1)
j = (3, 1)
j = (3, 2)
==================================
i = (3, 2)
j = (2, 1)
j = (3, 1)
j = (3, 2)
```
"""
macro above_diagonal_loop(N, itersym, imax, args...)
    lim = Expr(:call, :(:), Expr(:call, :(+), Expr(:curly, Symbol("$(itersym)_"), Expr(:call, :(+), :d, 1)), 1), imax)
    rng = Expr(:(->), :d, Expr(:block, Expr(:if, :(d == $N), :(1:$imax), lim)))
    Base.Cartesian._nloops(N, itersym, rng, args...)
end

"""
    anti_diagonal_loop(N, itersym, imax, args...)

Generate `N` Cartesian loops for the iteration variables
`itersym_{1:N}`, where `itersym_N ∈ 1:imax`, `itersym_{N-1} ∈
1:imax\\itersym_N`, etc, i.e. no two iteration variables have the same
values simultaneously. `args...` is passed on to
`Base.Cartesian._nloops`; however, `preexpr` is already used to skip
the diagonal elements. `anti_diagonal_loop` is nestable.

# Examples

```jldoctest
julia> @anti_diagonal_loop 3 i 3 begin
           println("-----------------------------")
           t = (Base.Cartesian.@ntuple 3 i)
           println("\$t: ", allunique(t))
           @anti_diagonal_loop 2 j 2 begin
               u = (Base.Cartesian.@ntuple 2 j)
               println("\$u: ", allunique(u))
           end
       end
-----------------------------
(3, 2, 1): true
(2, 1): true
(1, 2): true
-----------------------------
(2, 3, 1): true
(2, 1): true
(1, 2): true
-----------------------------
(3, 1, 2): true
(2, 1): true
(1, 2): true
-----------------------------
(1, 3, 2): true
(2, 1): true
(1, 2): true
-----------------------------
(2, 1, 3): true
(2, 1): true
(1, 2): true
-----------------------------
(1, 2, 3): true
(2, 1): true
(1, 2): true
```
"""
macro anti_diagonal_loop(N, itersym, imax, args...)
    rng = :(d -> 1:$imax)
    # The preexpr body generates tests for inner loops that are true
    # if the inner loop variable equals any of the outer ones; then
    # that iteration is skipped.
    prebody = Expr(:(->), :d,
                   Expr(:call, :(==),
                        Symbol("$(itersym)_e"),
                        Expr(:curly, Symbol("$(itersym)_"), :($N - d + 1))))
    pre = :(e -> (Base.Cartesian.@nany $N-e $prebody) && continue)

    Base.Cartesian._nloops(N, itersym, rng, pre, args...)
end


text:
macro inv(expr, cond)
    cond isa Expr && cond.head == :if || throw(ArgumentError("$cond is not an if expression"))
    cond.args[2] = expr
    return cond
end

@inv println("Wow! Lucky Guess!") if true else println("Not!") end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using PolynomialZeros
const AMVW = PolynomialZeros.AMVW
using Polynomials
using Test
using LinearAlgebra

# transformations

# # givens rotation
# @testset "Givens rotations" begin
#     a,b = complex(1.0, 2.0), complex(2.0, 3.0)
#     c,s,r = AMVW.givensrot(a,b)
#     @test norm(([c -conj(s); s conj(c)] * [a,b])[2]) <= 4eps(Float64)
#     a,b = complex(rand(2)...), complex(rand(2)...)
#     c,s,r = AMVW.givensrot(a,b)
#     @test norm(([c -conj(s); s conj(c)] * [a,b])[2]) <= 4eps(Float64)
# end

# # dflip
# @testset "D flip" begin
#     # t1 = pi/3;
#     # alpha = complex(cos(t1), sin(t1))
#     # d = AMVW.ComplexComplexRotator(alpha, complex(0.0, 0.0), 1)
#     # u = one(AMVW.ComplexComplexRotator{Float64})
#     # AMVW.vals!(u, complex(1.0, 2.0), complex(2.0, 3.0)); AMVW.idx!(u, 2)
#     # M = AMVW.as_full(u, 3) * AMVW.as_full(d,3)
#     # AMVW.Dflip(u, d)
#     # M1 = AMVW.as_full(d, 3) * AMVW.as_full(u,3)
#     # u = M - M1
#     # @test  maximum(norm.(u)) <= 4eps()

#     # complexrealrotator is different
#     #  U --> U D
#     # D         D
#     # r1,r2 = ones(AMVW.ComplexRealRotator{Float64},2)
#     # AMVW.vals!(r1, complex(1.0, 2.0), 2.0); AMVW.idx!(r1, 1)
#     # di = one(AMVW.ComplexRealRotator{Float64})
#     # AMVW.vals!(di, complex(cos(pi/3), sin(pi/3)), 0.0); AMVW.idx!(di, 2)
#     # M = AMVW.as_full(di, 3) * AMVW.as_full(r1, 3)
#     # AMVW.Dflip(di, r1)
#     # dic = copy(di); AMVW.idx!(dic, AMVW.idx(r1))
#     # M1 = AMVW.as_full(r1, 3) * AMVW.as_full(dic, 3) * AMVW.as_full(di, 3)
#     # @test maximum(norm.(M - M1)) <= 4eps()

#     ## D   -->      D
#     ##   U      U D
#     ##
#     # r1,r2 = ones(AMVW.ComplexRealRotator{Float64},2)
#     # AMVW.vals!(r1, complex(1.0, 2.0), 2.0); AMVW.idx!(r1, 2)
#     # di = one(AMVW.ComplexRealRotator{Float64})
#     # AMVW.vals!(di, complex(cos(pi/3), sin(pi/3)), 0.0); AMVW.idx!(di, 1)
#     # M = AMVW.as_full(di, 3) * AMVW.as_full(r1, 3)
#     # AMVW.Dflip(di, r1)
#     # dic = copy(di); AMVW.idx!(dic, AMVW.idx(r1))
#     # M1 = AMVW.as_full(r1, 3) * AMVW.as_full(dic, 3) * AMVW.as_full(di, 3)
#     # @test maximum(norm.(M - M1)) <= 4eps()

#     # ##
#     # ## Q   --> D   Q
#     # ##   D       D
#     # r1,r2 = ones(AMVW.ComplexRealRotator{Float64},2)
#     # AMVW.vals!(r1, complex(1.0, 2.0), 2.0); AMVW.idx!(r1, 1)
#     # di = one(AMVW.ComplexRealRotator{Float64})
#     # AMVW.vals!(di, complex(cos(pi/3), sin(pi/3)), 0.0); AMVW.idx!(di, 2)
#     # M = AMVW.as_full(r1, 3) * AMVW.as_full(di, 3)
#     # AMVW.Dflip(di, r1)
#     # dic = copy(di); AMVW.idx!(dic, AMVW.idx(r1))
#     # M1 = AMVW.as_full(dic, 3) * AMVW.as_full(di, 3) * AMVW.as_full(r1, 3)
#     # @test maximum(norm.(M - M1)) <= 4eps()

# end

# @testset "Fuse" begin
#     ##
#     # fuse
#     r1,r2 = ones(AMVW.ComplexComplexRotator{Float64},2)
#     AMVW.vals!(r1, complex(1.0, 2.0), complex(2.0, 3.0)); AMVW.idx!(r1, 1)
#     AMVW.vals!(r2, complex(3.0, 2.0), complex(5.0, 3.0)); AMVW.idx!(r2, 1)
#     M = AMVW.as_full(r1,2) * AMVW.as_full(r2, 2)
#     AMVW.fuse(r1, r2, Val{:left})
#     M1 = AMVW.as_full(r1, 2)
#     u = M - M1
#     @test maximum(norm.(u)) <= 4eps()

#     r1,r2 = ones(AMVW.ComplexRealRotator{Float64},2)
#     AMVW.vals!(r1, complex(1.0, 2.0), 2.0); AMVW.idx!(r1, 1)
#     AMVW.vals!(r2, complex(3.0, 2.0), 5.0); AMVW.idx!(r2, 1)
#     M = AMVW.as_full(r1,2) * AMVW.as_full(r2, 2)
#     alpha = AMVW.fuse(r1, r2, Val{:left})
#     M1 = AMVW.as_full(r1, 2) * diagm(0 => [alpha, conj(alpha)])
#     u = M - M1
#     @test maximum(norm.(u)) <= 4eps()

#     r1,r2 = ones(AMVW.ComplexRealRotator{Float64},2)
#     AMVW.vals!(r1, complex(1.0, 2.0), 2.0); AMVW.idx!(r1, 1)
#     AMVW.vals!(r2, complex(3.0, 2.0), 5.0); AMVW.idx!(r2, 1)
#     M = AMVW.as_full(r1,2) * AMVW.as_full(r2, 2)
#     alpha = AMVW.fuse(r1, r2, Val{:right})
#     M1 =  AMVW.as_full(r2, 2)  * diagm(0 => [alpha, conj(alpha)])
#     u = M - M1
#     @test maximum(norm.(u)) <= 4eps()

# end

# ## Cascade
# @testset "Cascade" begin
#     D1, Q1, Q2,Q3,Q4 =  AMVW._ones(AMVW.ComplexRealRotator{Float64},5)
#     alpha = complex(1.0, -1.0)
#     alpha = alpha/norm(alpha)
#     AMVW.vals!(D1, alpha, 0.0); AMVW.idx!(D1, 1)
#     AMVW.vals!(Q2, complex(1.0, 2.0), 2.0); AMVW.idx!(Q2, 2)
#     AMVW.vals!(Q3, complex(3.0, 2.0), 5.0); AMVW.idx!(Q3, 3)
#     AMVW.vals!(Q4, complex(2.0, 2.0), 3.0); AMVW.idx!(Q4, 4)

#     M1 = AMVW.as_full(D1, 5) * AMVW.as_full(Q2, 5) * AMVW.as_full(Q3, 5) * AMVW.as_full(Q4, 5)
#     D = ones(Complex{Float64}, 5)
#     Qs = [Q1, Q2, Q3, Q4]
#     AMVW.cascade(Qs, D, alpha, 1, 4)
#     M2 = AMVW.as_full(Q2, 5) * AMVW.as_full(Q3, 5) * AMVW.as_full(Q4, 5) * diagm(0 => D)

#     u = M1 - M2
#     @test maximum(norm.(u)) <= 4eps()

# end

# ##
# # turnover
# @testset "Turnover" begin
#     r1,r2,r3 = AMVW._ones(AMVW.ComplexComplexRotator{Float64}, 3)
#     AMVW.vals!(r1, complex(1.0, 2.0), complex(2.0, 3.0)); AMVW.idx!(r1, 1)
#     AMVW.vals!(r2, complex(3.0, 2.0), complex(5.0, 3.0)); AMVW.idx!(r2, 2)
#     AMVW.vals!(r3, complex(4.0, 2.0), complex(6.0, 3.0)); AMVW.idx!(r3, 1)

#     M = AMVW.as_full(r1,3) * AMVW.as_full(r2,3) * AMVW.as_full(r3,3)
#     AMVW.turnover(r1, r2, r3, Val{:right})
#     M1 =  AMVW.as_full(r3,3) * AMVW.as_full(r1,3) * AMVW.as_full(r2,3)
#     u = M - M1
#     @test maximum(norm.(u)) <= 4eps()

#     r1,r2,r3 = AMVW._ones(AMVW.ComplexRealRotator{Float64}, 3)
#     AMVW.vals!(r1, complex(1.0, 2.0), 2.0); AMVW.idx!(r1, 1)
#     AMVW.vals!(r2, complex(3.0, 2.0), 5.0); AMVW.idx!(r2, 2)
#     AMVW.vals!(r3, complex(4.0, 2.0), 6.0); AMVW.idx!(r3, 1)

#     M = AMVW.as_full(r1,3) * AMVW.as_full(r2,3) * AMVW.as_full(r3,3)
#     AMVW.turnover(r1, r2, r3, Val{:right})
#     M1 =  AMVW.as_full(r3,3) * AMVW.as_full(r1,3) * AMVW.as_full(r2,3)
#     u = M - M1
#     @test maximum(norm.(u)) <= 4eps()
#     AMVW.turnover(r3, r1, r2, Val{:left})
#     M2 = AMVW.as_full(r1,3) * AMVW.as_full(r2,3) * AMVW.as_full(r3,3)
#     u = M - M1
#     @test maximum(norm.(u)) <= 4eps()

# end

# # passthrough
# @testset "Passthrough" begin
#     r1,r2,r3 = AMVW._ones(AMVW.ComplexRealRotator{Float64}, 3)
#     AMVW.vals!(r1, complex(1.0, 2.0), 2.0); AMVW.idx!(r1, 1)
#     t1,t2,t3=pi/3,pi/4, pi/5
#     cplx(t) = complex(cos(t), sin(t))
#     D = [cplx(t) for t in [t1,t2,t3]]
#     M = diagm(0 =>D) * AMVW.as_full(r1, 3)
#     AMVW.passthrough(D, r1)
#     M1 = AMVW.as_full(r1, 3) * diagm(0 => D)
#     u = M - M1
#     @test maximum(norm.(u)) <= 4eps()

# end

@testset "poly_roots" begin

    ## real coeffs
    rs = [1.0, 2.0, 3.0]
    p = poly(rs)
    rts =  AMVW.poly_roots(p.a)
    sort!(rts, by=norm)
    @test maximum(norm.(rts - rs)) <= 1e-6

    # complex coeffs
    rs = [1.0, 2.0, 3.0, 4.0 + 1.0im]
    p = poly(rs)
    rts = AMVW.poly_roots(p.a)
    sort!(rts, by=norm)
    @test maximum(norm.(rts - rs)) <= 1e-6

    # # ComplexComplex
    # rs = [1.0, 2.0, 3.0, 4.0+0im]
    # p = poly(rs)
    # qs, k = AMVW.reverse_poly(p.a)
    # state = AMVW.ComplexComplexSingleShift(qs)
    # AMVW.init_state(state)
    # AMVW.AMVW_algorithm(state)
    # rts = complex.(state.REIGS, state.IEIGS)
    # sort!(rts, by=norm)
    # @test maximum(norm.(rts - rs)) <= 1e-6

    ## simple cases
    n=1
    rs = [1.0]
    p = poly(rs)
    rts = AMVW.poly_roots(p.a)
    @test maximum(norm.(rts - rs)) <= 1e-6

    rs = [1.0 + im]
    p = poly(rs)
    rts = AMVW.poly_roots(p.a)
    @test maximum(norm.(rts - rs)) <= 1e-6

    # n = 2
    rs = [1.0, 2.0]
    p = poly(rs)
    rts = AMVW.poly_roots(p.a)
    sort!(rts, by=norm)
    @test maximum(norm.(rts - rs)) <= 1e-6

    rs = [1.0, 2.0+im]
    p = poly(rs)
    rts = AMVW.poly_roots(p.a)
    sort!(rts, by=norm)
    @test maximum(norm.(rts - rs)) <= 1e-6

    # zeros
    rs = [1.0, 2.0, 3.0, 0.0, 0.0]
    sort!(rs, by=norm)
    p = poly(rs)
    rts = AMVW.poly_roots(p.a)
    sort!(rts, by=norm)
    @test maximum(norm.(rts - rs)) <= 1e-6

    rs = [1.0, 2.0, 3.0+im, 0.0, 0.0]
    sort!(rs, by=norm)
    p = poly(rs)
    rts = AMVW.poly_roots(p.a)
    sort!(rts, by=norm)
    @test maximum(norm.(rts - rs)) <= 1e-6

end

@testset "pencil factorization" begin
    pr = poly([1.0, 2.0, 3.0])
    rts = AMVW.poly_roots(coeffs(pr), pencil=AMVW.basic_pencil)
    @test maximum(norm.(pr.(rts))) <= 1e-12

    pi = poly([im+1.0, 2.0, 3.0])
    rts = AMVW.poly_roots(coeffs(pi),  pencil=AMVW.basic_pencil)
    @test maximum(norm.(pi.(rts))) <= 1e-12

end


text:
# load code
include("SDE_u1_dm.jl")
include("BSE_u1_dm_s.jl")
include("BSE_u1_dm_t.jl")
include("BSE_u1_dm_u.jl")


text:
using NavAbilitySDK

function solveGraph!(dfg::NavAbilityDFG; maxWaitSeconds::Int = 300)
  context = Client(dfg.userId,dfg.robotId,dfg.sessionId)
  solveSession(dfg.navabilityClient, context)
end


text:
"""運動学

"""
module Kinematics

"""Phi0のラッパー"""
function Phi0!(q::Vector{Float64}, ξ::Float64, out::Vector{Float64})
    ccall(
        (:Phi0, "o/soft_robot/derived/ikko_dake/eqs/c_so/Phi0.so"),
        Cvoid,
        (Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
        q[1], q[2], q[3], ξ, out
    )
end

"""グローバル位置"""
function Phi0(q::Vector{Float64}, ξ)
    Z = Vector{Float64}(undef, 3)
    Phi0!(q, ξ, Z)
    return Z
end

const Ξ = Vector(0:0.01:1)

"""アームのディスク位置を計算"""
function Arm(q::Vector{T},) where T
    Phis = Vector{typeof(q)}(undef, length(Ξ))
    for (i, ξ) in enumerate(Ξ)
        Phis[i] = Phi0(q, ξ)
    end
    return Phis
end

end

# using .Kinematics: Phi0
# println(Phi0([0.0, 0.0, 0.0], 1.0))


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
g = (k^2 for k in -3:3)

# %%
collect(g)

# %%
v = [k^2 for k in -3:3]

# %%
Tuple(g)

# %%
t = Tuple(k^2 for k in -3:3)

# %%
Set(g)

# %%
s = Set(k^2 for k in -3:3)

# %%
d = Dict(k => k^2 for k in -3:3)

# %%
:([k^2 for k in -3:3]) |> Meta.show_sexpr

# %%
:(Tuple(k^2 for k in -3:3)) |> Meta.show_sexpr

# %%
:(Set(k^2 for k in -3:3)) |> Meta.show_sexpr

# %%
:(Dict(k => k^2 for k in -3:3)) |> Meta.show_sexpr

# %%
f(x, y) = (d = x^2+y^2; d == 0 ? d : x^2*y/d)
X = Y = range(-1, 1; length=201)

z = [f(x, y) for y in Y, x in X]
Z = f.(X', Y)
z == Z

# %%
using Plots
surface(X, Y, Z; camera=(40, 70))

# %%
surface(X, Y, f; camera=(40, 70))

# %%

max_stars_repo_path: 0007/comprehension examples.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0007/comprehension examples.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0007/comprehension examples.jl

max_forks_repo_name: genkuroki/public


text:
"""
    mbsvplist(cmd0::String=""; kwargs...)

List water sound velocity profiles in swath sonar data files.

Parameters
----------

- **C** | **uniquesvp** :: [Type => Bool]

    Output the number of unique SVPs in each file.
- **F** | **format** :: [Type => Int]

    Sets the format for the input swath sonar data.
- **M** | **mode** :: [Type => Int]		``Arg = 1 or 2 or 3``

    Sets the SVP output mode..
- **S** | **ssv** :: [Type => Bool]

    Sets the minimum speed in km/hr (5.5 kts ~ 10 km/hr) allowed in the input data.
- **Z** | **firstiszero** :: [Type => Bool]

    Sets the style of the plot.
- $(GMT.opt_R)
- $(GMT.opt_V)
"""
function mbsvplist(cmd0::String=""; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("mbsvplist", cmd0)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode

	cmd = parse_common_opts(d, "", [:R :yx :V_params :o])[1]
	cmd = parse_these_opts(cmd, d, [[:C :uniquesvp], [:F :format], [:H :help], [:M :mode], [:S :ssv], [:Z :firstiszero], [:z :z_down]])
	((val = find_in_dict(d, [:z :z_down])[1]) === nothing) && (cmd *= " -z")	# Means here default is Z-up
	(!occursin("-:", cmd)) && (cmd = " -:" * cmd)		# Means that the default here (contrary to C version) is speed-dept

	cmd = "mbsvplist -I" * cmd0 * cmd				# In any case we need this
	finish_PS_module(d, cmd, "", true, false, false)
end


text:
using Logging, Test, HSM

struct CommonParentState <: HSM.AbstractHsmState
    state_info::HSM.HsmStateInfo

    CommonParentState(parent_state) = new(HSM.HsmStateInfo(parent_state))
end

"""
    basic_state_machine()

Builds and returns the following state machine:

```plantuml
@startuml
state A {
    state B
    state C {
        state D
        state E
    }
}
@enduml
```
"""
function common_parent_state_machine()
    A = CommonParentState(nothing)
    B = CommonParentState(A)
    C = CommonParentState(A)
    D = CommonParentState(C)
    E = CommonParentState(C)

    return A, B, C, D, E
end

@testset "common_parent()" begin
    A, B, C, D, E = common_parent_state_machine()

    # Common parent of A and any other state is A.
    @test HSM.common_parent(A, B) == A
    @test HSM.common_parent(D, A) == A

    # Common parent of B and any other state is A.
    @test HSM.common_parent(B, C) == A
    @test HSM.common_parent(D, B) == A

    # Common parent of C and its substates is C.
    @test HSM.common_parent(C, D) == C
    @test HSM.common_parent(E, C) == C

    # Common parent of D and E is C.
    @test HSM.common_parent(D, E) == C
end


text:
struct Segment{T,U}
    value::T
    tstart::U
    tend::U
end

@inline value(s::Segment) = s.value
@inline timestamp(s::Segment) = s.tstart
@inline duration(s::Segment) = s.tend - s.tstart


text:
# Use BLAS.gemm for strided matrix-matrix multiplication sensitivites. Don't bother with
# BLAS for matrix-vector stuff yet. Definitely an optimisation that we might want to
# consider at some point in the future though.
const RS = StridedMatrix{<:∇Scalar}
const RST = Transpose{<:∇Scalar, RS}
const RSA = Adjoint{<:∇Scalar, RS}
strided_matmul = [
    (RS,  RS,  'N', 'C', :Ȳ, :B, 'C', 'N', :A, :Ȳ),
    (RST, RS,  'N', 'T', :B, :Ȳ, 'N', 'N', :A, :Ȳ),
    (RS,  RST, 'N', 'N', :Ȳ, :B, 'T', 'N', :Ȳ, :A),
    (RST, RST, 'T', 'T', :B, :Ȳ, 'T', 'T', :Ȳ, :A),
    (RSA, RS,  'N', 'C', :B, :Ȳ, 'N', 'N', :A, :Ȳ),
    (RS,  RSA, 'N', 'N', :Ȳ, :B, 'C', 'N', :Ȳ, :A),
    (RSA, RSA, 'C', 'C', :B, :Ȳ, 'C', 'C', :Ȳ, :A),
]
import Base: *
for (TA, TB, tCA, tDA, CA, DA, tCB, tDB, CB, DB) in strided_matmul

    # Add intercepts and export names.
    @eval @explicit_intercepts $(Symbol("*")) Tuple{$TA, $TB}

    # Define allocating and non-allocating sensitivities for each output.
    alloc_Ā = :(LinearAlgebra.BLAS.gemm($tCA, $tDA, $CA, $DA))
    alloc_B̄ = :(LinearAlgebra.BLAS.gemm($tCB, $tDB, $CB, $DB))
    no_alloc_Ā = :(LinearAlgebra.BLAS.gemm!($tCA, $tDA, 1., $CA, $DA, 1., Ā))
    no_alloc_B̄ = :(LinearAlgebra.BLAS.gemm!($tCB, $tDB, 1., $CB, $DB, 1., B̄))

    # Add sensitivity definitions.
    @eval ∇(::typeof(*), ::Type{Arg{1}}, p, Y::RS, Ȳ::RS, A::$TA, B::$TB) = $alloc_Ā
    @eval ∇(::typeof(*), ::Type{Arg{2}}, p, Y::RS, Ȳ::RS, A::$TA, B::$TB) = $alloc_B̄
    @eval ∇(Ā, ::typeof(*), ::Type{Arg{1}}, p, Y::RS, Ȳ::RS, A::$TA, B::$TB) = $no_alloc_Ā
    @eval ∇(B̄, ::typeof(*), ::Type{Arg{2}}, p, Y::RS, Ȳ::RS, A::$TA, B::$TB) = $no_alloc_B̄
end

# # Not every permutation of transpositions makes sense for matrix-vector multiplication. This
# # list just includes those which make sense.
# strided_matvecmul = [
#     (:*,         'C', :ȳ, :b, 'C'),
#     (:At_mul_B,  'T', :b, :ȳ, 'N'),
#     (:Ac_mul_B,  'C', :b, :ȳ, 'N'),
# ]
# for (f, tdA, CA, dA, tCb) in strided_matvecmul
#     n_Ā, u_Ā = tdA == 'C' ? :(Ā = $CA * $dA') : :(Ā = $CA * $dA'), :(ger!(1., $CA, $dA, Ā))
#     n_b̄, u_b̄ = :(b̄ = gemv($tCb, A, ȳ)), :(b̄ = gemv!($tCb, 1., A, ȳ, 1., b̄))
#     generate_primitive(f, [:(T <: StridedMatrix), :(V <: StridedVector)],
#         [:A, :b], [:Ā, :b̄], [:T, :V], [true, true], :y, :ȳ, [n_Ā, n_b̄], [u_Ā, u_b̄])
# end

# # Operations of the for Y = A \ B
# strided_ldiv = [
#     (:\, :(C = At_ldiv_B(A, Ȳ)), 'N', 'T', :C, :Y),
#     (:At_ldiv_B, :(C = A \ Ȳ), 'N', 'T', :Y, :C),
#     (:A_ldiv_Bt, :(C = At_rdiv_B(Ȳ, A)), 'T', 'T', :C, :Y),
#     (:At_ldiv_Bt, :(C = At_rdiv_Bt(Ȳ, A)), 'N', 'N', :Y, :C),
#     (:Ac_ldiv_B, :(C = A \ Ȳ), 'N', 'C', :Y, :C),
#     (:A_ldiv_Bc, :(C = Ac_rdiv_B(Ȳ, A)), 'C', 'C', :C, :Y),
#     (:Ac_ldiv_Bc, :(C = Ac_rdiv_Bc(Ȳ, A)), 'N', 'N', :Y, :C),
# ]

# # Iterate through primitive definitions and add methods for each.
# for (f, C, tA, tB, arg1, arg2) in strided_ldiv
#     new_Ā = :(Ā = gemm($tA, $tB, -1.0, $arg1, $arg2))
#     update_Ā = :(gemm!($tA, $tB, -1.0, $arg1, $arg2, 1.0, Ā))
#     new_B̄ = :(B̄ = C)
#     update_B̄ = :(broadcast!((b̄, c)->b̄ + c, B̄, B̄, C))
#     generate_primitive(f, [:(T <: Any), :(V <: Any)],
#         [:A, :B], [:Ā, :B̄], [:T, :V], [true, true], :Y, :Ȳ,
#         [new_Ā, new_B̄], [update_Ā, update_B̄], C)
# end