text:
int_rules_4_1_2_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.2.3*(g*sin)^p*(a+b*sin)^m*(c+d*sin)^n =#
    @apply_utils Antiderivative((sqrt(~(g') * sin(~(e') + ~(f') * ~x)) * sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x))) / (~c + ~(d') * sin(~(e') + ~(f') * ~x)), ~x) => (~g / ~d) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / sqrt(~g * sin(~e + ~f * ~x)), ~x) - ((~c * ~g) / ~d) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / (sqrt(~g * sin(~e + ~f * ~x)) * (~c + ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) || EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative((sqrt(~(g') * sin(~(e') + ~(f') * ~x)) * sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x))) / (~c + ~(d') * sin(~(e') + ~(f') * ~x)), ~x) => (~b / ~d) * Antiderivative(sqrt(~g * sin(~e + ~f * ~x)) / sqrt(~a + ~b * sin(~e + ~f * ~x)), ~x) - ((~b * ~c - ~a * ~d) / ~d) * Antiderivative(sqrt(~g * sin(~e + ~f * ~x)) / (sqrt(~a + ~b * sin(~e + ~f * ~x)) * (~c + ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) / (sqrt(~(g') * sin(~(e') + ~(f') * ~x)) * (~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((-2 * ~b) / ~f) * Subst(Antiderivative(1 / (~b * ~c + ~a * ~d + ~c * ~g * (~x) ^ 2), ~x), ~x, (~b * cos(~e + ~f * ~x)) / (sqrt(~g * sin(~e + ~f * ~x)) * sqrt(~a + ~b * sin(~e + ~f * ~x)))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) / (sqrt(sin(~(e') + ~(f') * ~x)) * (~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => (-(sqrt(~a + ~b)) / (~c * ~f)) * Elliptic.E(asin(cos(~e + ~f * ~x) / (1 + sin(~e + ~f * ~x))), -((~a - ~b)) / (~a + ~b)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ(~d, ~c) && (GtQ((~b) ^ 2 - (~a) ^ 2, 0) && GtQ(~b, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) / (sqrt(~(g') * sin(~(e') + ~(f') * ~x)) * (~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((-(sqrt(~a + ~b * sin(~e + ~f * ~x))) * sqrt((~d * sin(~e + ~f * ~x)) / (~c + ~d * sin(~e + ~f * ~x)))) / (~d * ~f * sqrt(~g * sin(~e + ~f * ~x)) * sqrt(((~c) ^ 2 * (~a + ~b * sin(~e + ~f * ~x))) / ((~a * ~c + ~b * ~d) * (~c + ~d * sin(~e + ~f * ~x)))))) * Elliptic.E(asin((~c * cos(~e + ~f * ~x)) / (~c + ~d * sin(~e + ~f * ~x))), (~b * ~c - ~a * ~d) / (~b * ~c + ~a * ~d)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) / (sqrt(~(g') * sin(~(e') + ~(f') * ~x)) * (~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => (~a / ~c) * Antiderivative(1 / (sqrt(~g * sin(~e + ~f * ~x)) * sqrt(~a + ~b * sin(~e + ~f * ~x))), ~x) + ((~b * ~c - ~a * ~d) / (~c * ~g)) * Antiderivative(sqrt(~g * sin(~e + ~f * ~x)) / (sqrt(~a + ~b * sin(~e + ~f * ~x)) * (~c + ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) / (sin(~(e') + ~(f') * ~x) * (~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => (1 / ~c) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / sin(~e + ~f * ~x), ~x) - (~d / ~c) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / (~c + ~d * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) / (sin(~(e') + ~(f') * ~x) * (~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => (~a / ~c) * Antiderivative(1 / (sin(~e + ~f * ~x) * sqrt(~a + ~b * sin(~e + ~f * ~x))), ~x) + ((~b * ~c - ~a * ~d) / ~c) * Antiderivative(1 / (sqrt(~a + ~b * sin(~e + ~f * ~x)) * (~c + ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && NeQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative(sqrt(~(g') * sin(~(e') + ~(f') * ~x)) / (sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * (~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((-(~a) * ~g) / (~b * ~c - ~a * ~d)) * Antiderivative(1 / (sqrt(~g * sin(~e + ~f * ~x)) * sqrt(~a + ~b * sin(~e + ~f * ~x))), ~x) + ((~c * ~g) / (~b * ~c - ~a * ~d)) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / (sqrt(~g * sin(~e + ~f * ~x)) * (~c + ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) || EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(sqrt(~(g') * sin(~(e') + ~(f') * ~x)) / (sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * (~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((2 * sqrt(-(cot(~e + ~f * ~x) ^ 2)) * sqrt(~g * sin(~e + ~f * ~x))) / (~f * (~c + ~d) * cot(~e + ~f * ~x) * sqrt(~a + ~b * sin(~e + ~f * ~x)))) * sqrt((~b + ~a * csc(~e + ~f * ~x)) / (~a + ~b)) * Elliptic.Pi((2 * ~c) / (~c + ~d), asin(sqrt(1 - csc(~e + ~f * ~x)) / sqrt(2)), (2 * ~a) / (~a + ~b)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(1 / (sqrt(~(g') * sin(~(e') + ~(f') * ~x)) * sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * (~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => (~b / (~b * ~c - ~a * ~d)) * Antiderivative(1 / (sqrt(~g * sin(~e + ~f * ~x)) * sqrt(~a + ~b * sin(~e + ~f * ~x))), ~x) - (~d / (~b * ~c - ~a * ~d)) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / (sqrt(~g * sin(~e + ~f * ~x)) * (~c + ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) || EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(1 / (sqrt(~(g') * sin(~(e') + ~(f') * ~x)) * sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * (~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => (1 / ~c) * Antiderivative(1 / (sqrt(~g * sin(~e + ~f * ~x)) * sqrt(~a + ~b * sin(~e + ~f * ~x))), ~x) - (~d / (~c * ~g)) * Antiderivative(sqrt(~g * sin(~e + ~f * ~x)) / (sqrt(~a + ~b * sin(~e + ~f * ~x)) * (~c + ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(1 / (sin(~(e') + ~(f') * ~x) * sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * (~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((~d) ^ 2 / (~c * (~b * ~c - ~a * ~d))) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / (~c + ~d * sin(~e + ~f * ~x)), ~x) + (1 / (~c * (~b * ~c - ~a * ~d))) * Antiderivative(((~b * ~c - ~a * ~d) - ~b * ~d * sin(~e + ~f * ~x)) / (sin(~e + ~f * ~x) * sqrt(~a + ~b * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative(1 / (sin(~(e') + ~(f') * ~x) * sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * (~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => (1 / ~c) * Antiderivative(1 / (sin(~e + ~f * ~x) * sqrt(~a + ~b * sin(~e + ~f * ~x))), ~x) - (~d / ~c) * Antiderivative(1 / (sqrt(~a + ~b * sin(~e + ~f * ~x)) * (~c + ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && NeQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) / (sin(~(e') + ~(f') * ~x) * sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => (-(~d) / ~c) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / sqrt(~c + ~d * sin(~e + ~f * ~x)), ~x) + (1 / ~c) * Antiderivative((sqrt(~a + ~b * sin(~e + ~f * ~x)) * sqrt(~c + ~d * sin(~e + ~f * ~x))) / sin(~e + ~f * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && EqQ(~b * ~c + ~a * ~d, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) / (sin(~(e') + ~(f') * ~x) * sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((-2 * ~a) / ~f) * Subst(Antiderivative(1 / (1 - ~a * ~c * (~x) ^ 2), ~x), ~x, cos(~e + ~f * ~x) / (sqrt(~a + ~b * sin(~e + ~f * ~x)) * sqrt(~c + ~d * sin(~e + ~f * ~x)))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ(~b * ~c + ~a * ~d, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) / (sin(~(e') + ~(f') * ~x) * sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((~b * ~c - ~a * ~d) / ~c) * Antiderivative(1 / (sqrt(~a + ~b * sin(~e + ~f * ~x)) * sqrt(~c + ~d * sin(~e + ~f * ~x))), ~x) + (~a / ~c) * Antiderivative(sqrt(~c + ~d * sin(~e + ~f * ~x)) / (sin(~e + ~f * ~x) * sqrt(~a + ~b * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) / (sin(~(e') + ~(f') * ~x) * sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => ((-2 * (~a + ~b * sin(~e + ~f * ~x))) / (~c * ~f * Rt((~a + ~b) / (~c + ~d), 2) * cos(~e + ~f * ~x))) * sqrt((-((~b * ~c - ~a * ~d)) * (1 - sin(~e + ~f * ~x))) / ((~c + ~d) * (~a + ~b * sin(~e + ~f * ~x)))) * sqrt(((~b * ~c - ~a * ~d) * (1 + sin(~e + ~f * ~x))) / ((~c - ~d) * (~a + ~b * sin(~e + ~f * ~x)))) * Elliptic.Pi((~a * (~c + ~d)) / (~c * (~a + ~b)), asin((Rt((~a + ~b) / (~c + ~d), 2) * sqrt(~c + ~d * sin(~e + ~f * ~x))) / sqrt(~a + ~b * sin(~e + ~f * ~x))), ((~a - ~b) * (~c + ~d)) / ((~a + ~b) * (~c - ~d))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(1 / (sin(~(e') + ~(f') * ~x) * sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => (cos(~e + ~f * ~x) / (sqrt(~a + ~b * sin(~e + ~f * ~x)) * sqrt(~c + ~d * sin(~e + ~f * ~x)))) * Antiderivative(1 / (cos(~e + ~f * ~x) * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(1 / (sin(~(e') + ~(f') * ~x) * sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x))), ~x) => (-(~b) / ~a) * Antiderivative(1 / (sqrt(~a + ~b * sin(~e + ~f * ~x)) * sqrt(~c + ~d * sin(~e + ~f * ~x))), ~x) + (1 / ~a) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / (sin(~e + ~f * ~x) * sqrt(~c + ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) || NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative((sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x))) / sin(~(e') + ~(f') * ~x), ~x) => ((sqrt(~a + ~b * sin(~e + ~f * ~x)) * sqrt(~c + ~d * sin(~e + ~f * ~x))) / cos(~e + ~f * ~x)) * Antiderivative(cot(~e + ~f * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative((sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) * sqrt(~c + ~(d') * sin(~(e') + ~(f') * ~x))) / sin(~(e') + ~(f') * ~x), ~x) => ~d * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / sqrt(~c + ~d * sin(~e + ~f * ~x)), ~x) + ~c * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / (sin(~e + ~f * ~x) * sqrt(~c + ~d * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) || NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (~a) ^ ~n * (~c) ^ ~n * Antiderivative(tan(~e + ~f * ~x) ^ ~p * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (EqQ(~p + 2 * ~n, 0) && IntegerQ(~n))))
    @apply_utils Antiderivative((~(g') * sin(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((sqrt(~a - ~b * sin(~e + ~f * ~x)) * sqrt(~a + ~b * sin(~e + ~f * ~x))) / (~f * cos(~e + ~f * ~x))) * Subst(Antiderivative(((~g * ~x) ^ ~p * (~a + ~b * ~x) ^ (~m - 1 / 2) * (~c + ~d * ~x) ^ ~n) / sqrt(~a - ~b * ~x), ~x), ~x, sin(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && IntegerQ(~m - 1 / 2))))
    @apply_utils Antiderivative((~(g') * sin(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandTrig((~g * sin(~e + ~f * ~x)) ^ ~p * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n, ~p], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && ((IntegersQ(~m, ~n) || (IntegersQ(~m, ~p) || IntegersQ(~n, ~p))) && NeQ(~p, 2)))
    @apply_utils Antiderivative((~(g') * sin(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~n, ~x) => Unintegrable((~g * sin(~e + ~f * ~x)) ^ ~p * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && NeQ(~p, 2)
    @apply_utils Antiderivative((~(g') * sin(~(e') + ~(f') * ~x)) ^ ~(p') * (~(a') + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (~g) ^ (~m + ~n) * Antiderivative((~g * sin(~e + ~f * ~x)) ^ ((~p - ~m) - ~n) * (~b + ~a * sin(~e + ~f * ~x)) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~p], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (Not(IntegerQ(~p)) && (IntegerQ(~m) && IntegerQ(~n))))
    @apply_utils Antiderivative((~(g') * sin(~(e') + ~(f') * ~x)) ^ ~(p') * (~(a') + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (~g * csc(~e + ~f * ~x)) ^ ~p * (~g * sin(~e + ~f * ~x)) ^ ~p * Antiderivative(((~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n) / (~g * csc(~e + ~f * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (Not(IntegerQ(~p)) && Not(IntegerQ(~m) && IntegerQ(~n))))
    @apply_utils Antiderivative((~(g') * sin(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (~g) ^ ~n * Antiderivative((~g * sin(~e + ~f * ~x)) ^ (~p - ~n) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~p], ~x) && IntegerQ(~n)
    @apply_utils Antiderivative(sin(~(e') + ~(f') * ~x) ^ ~(p') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => Antiderivative(((~b + ~a * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n) / csc(~e + ~f * ~x) ^ (~m + ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (Not(IntegerQ(~n)) && (IntegerQ(~m) && IntegerQ(~p)))
    @apply_utils Antiderivative((~(g') * sin(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => csc(~e + ~f * ~x) ^ ~p * (~g * sin(~e + ~f * ~x)) ^ ~p * Antiderivative(((~b + ~a * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n) / csc(~e + ~f * ~x) ^ (~m + ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n, ~p], ~x) && (Not(IntegerQ(~n)) && (IntegerQ(~m) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~(g') * sin(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (((~g * sin(~e + ~f * ~x)) ^ ~n * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n) / (~d + ~c * sin(~e + ~f * ~x)) ^ ~n) * Antiderivative((~g * sin(~e + ~f * ~x)) ^ (~p - ~n) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~n)) && Not(IntegerQ(~m)))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ ~(p') * (~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (~g) ^ (~m + ~n) * Antiderivative((~g * csc(~e + ~f * ~x)) ^ ((~p - ~m) - ~n) * (~b + ~a * csc(~e + ~f * ~x)) ^ ~m * (~d + ~c * csc(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~p], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (Not(IntegerQ(~p)) && (IntegerQ(~m) && IntegerQ(~n))))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ ~(p') * (~(a') + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * sin(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (~g * csc(~e + ~f * ~x)) ^ ~p * (~g * sin(~e + ~f * ~x)) ^ ~p * Antiderivative(((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~c + ~d * sin(~e + ~f * ~x)) ^ ~n) / (~g * sin(~e + ~f * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (Not(IntegerQ(~p)) && Not(IntegerQ(~m) && IntegerQ(~n))))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (~g) ^ ~m * Antiderivative((~g * csc(~e + ~f * ~x)) ^ (~p - ~m) * (~b + ~a * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n, ~p], ~x) && IntegerQ(~m)
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) ^ ~(p') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => Antiderivative(((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n) / sin(~e + ~f * ~x) ^ (~n + ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (Not(IntegerQ(~m)) && (IntegerQ(~n) && IntegerQ(~p)))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => sin(~e + ~f * ~x) ^ ~p * (~g * csc(~e + ~f * ~x)) ^ ~p * Antiderivative(((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n) / sin(~e + ~f * ~x) ^ (~n + ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~p], ~x) && (Not(IntegerQ(~m)) && (IntegerQ(~n) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~g * csc(~e + ~f * ~x)) ^ ~m) / (~b + ~a * csc(~e + ~f * ~x)) ^ ~m) * Antiderivative((~g * csc(~e + ~f * ~x)) ^ (~p - ~m) * (~b + ~a * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~m)) && Not(IntegerQ(~n)))
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.2.3 (g sin)^p (a+b sin)^m (c+d sin)^n.jl


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module CatPerf
import Perftests: @perf, meta

function hvcat_perf(a,b)
    return [a b; b a]
end

function hvcat_setind_perf(a, b)
    n = size(a,1)
    c = Array(Float64,2n,2n)
    c[  1:n,    1:n  ] = a
    c[  1:n,  n+1:end] = b
    c[n+1:end,  1:n  ] = b
    c[n+1:end,n+1:end] = a
    return c
end

function hcat_perf(a,b)
    return [a b b a]
end

function hcat_setind_perf(a,b)
    n = size(a,1)
    c = Array(Float64, n, 4n)
    c[:,    1:  n] = a
    c[:,  n+1: 2n] = b
    c[:, 2n+1: 3n] = b
    c[:, 3n+1:end] = a
    return c
end

function vcat_perf(a,b)
    return [a; b; b; a]
end

function vcat_setind_perf(a,b)
    n = size(a,1)
    c = Array(Float64, 4n, n)
    c[   1: n, :] = a
    c[ n+1:2n, :] = b
    c[2n+1:3n, :] = b
    c[3n+1:4n, :] = a
    return c
end

function catnd_perf(a,b)
    return cat(3, a, b, b, a)
end

function catnd_setind_perf(a, b)
    n = size(a,2)
    c = Array(Float64, 1, n, 4n, 1)
    c[1,:,   1: n,1] = a
    c[1,:, n+1:2n,1] = b
    c[1,:,2n+1:3n,1] = b
    c[1,:,3n+1:4n,1] = a
    return c
end

sizes = [("small", 5), ("large", 500)]
for (size_str, n) in sizes
    a = rand(n,n)
    b = rand(n,n)
    @perf hvcat_perf(a,b) meta("hvcat", size_str, "Horizontal/vertical matrix concatenation")
    @perf hvcat_setind_perf(a,b) meta("hvcat_setind", size_str, "Horizontal/vertical matrix concatenation using setindex")
    @perf hcat_perf(a,b) meta("hcat", size_str, "Horizontal matrix concatenation")
    @perf hcat_setind_perf(a,b) meta("hcat_setind", size_str, "Horizontal matrix concatenation using setindex")
    @perf vcat_perf(a,b) meta("vcat", size_str, "Vertical matrix concatenation")
    @perf vcat_setind_perf(a,b) meta("vcat_setind", size_str, "Vertical matrix concatenation using setindex")

    a = rand(1,n,n,1)
    b = rand(1,n,n)
    @perf catnd_perf(a,b) meta("catnd", size_str, "N-dimensional matrix concatenation")
    @perf catnd_setind_perf(a,b) meta("catnd_setind", size_str, "N-dimensional matrix concatenation using setindex")
end

end # module

ext: jl


text:
function BL_get_f_g(featuresTP::Matrix{Float64}, labels::MinusPlusOneVector, W::Matrix{Float64}, index::Int64)
	X = featuresTP[:,index:index]
	Y = labels[index]
	ym=Y*((X'*W)[1])
	(log(1 + exp(-ym)), X*(-Y / (1 + exp(ym))))
end

function BL_get_f(features::Matrix{Float64}, labels::MinusPlusOneVector, W::Matrix{Float64})
	ym=labels.field.*(features*W)
	mean(log(1 + exp(-ym)))
end

# Returns f, pcc, mcc, tp, tn, fp, fn
function BL_for_output(features::Matrix{Float64},labels::MinusPlusOneVector,W::Matrix{Float64})
	ym=labels.field.*(features*W)
	fp, fn = 0,0
	for i in 1:length(labels)
		if ym[i]<=0
			if labels.field[i]<0
				fp += 1
			elseif labels.field[i]>0
				fn += 1
			end
		end
	end
	tp = labels.numPlus - fn
	tn = length(labels) - labels.numPlus - fp
	pcc = (tp+tn)/(tp+tn+fp+fn)
	mcc = (tp*tn - fp*fn)/(sqrt(tp+fp)*sqrt(tp+fn)*sqrt(tn+fp)*sqrt(tn+fn))
	(mean(log(1 + exp(-ym))), pcc, mcc, tp, tn, fp, fn)
end


text:
using Clang
using Test

include("generators.jl")
include("test_bitfield.jl")


text:
using DataStructures, GraphPlot

function create_graph_database(graph::Any, graph_type::String)

    if isfile("graph_database.js")
        run(`rm graph_database.js`)
    end
    run(`touch graph_database.js`)

    graph_file = open("graph_database.js", "w")

    print(graph_file, string("var graph_edges = ["))

    added_pairs = Stack{Any}()

    for (i, vertex) in enumerate(collect(keys(graph.vprops)))
        neighbours = neighbors(graph, vertex)
        for (j, neighbour) in enumerate(neighbours)
            if Pair(vertex, neighbour) ∉ added_pairs && Pair(neighbour, vertex) ∉ added_pairs
                push!(added_pairs, Pair(vertex, neighbour))
            end
        end
    end
    while !isempty(added_pairs)
        edge = pop!(added_pairs)
        if isempty(added_pairs)
            println(graph_file, string("{from: ", edge[1], ", to: ", edge[2], "}];"))
        else
            print(graph_file, string("{from: ", edge[1], ", to: ", edge[2], "},"))
        end
    end

    ##############################

    print(graph_file, string("var graph_nodes = ["))

    center_vertex = 1
    if graph_type == "circles"
        nlist = [collect(i*20-19:i*20) for i=1:(Int(floor(length(vertices(graph))/20))+20)]
        x, y = shell_layout(graph, nlist)
        center_vertex = collect(filter(x->length(neighbors(g, x))==Δ(g), keys(g.vprops)))[1]
    elseif graph_type == "spring"
        x, y = spring_layout(graph)
    end
    x*=500
    y*=500

    for (i, vertex) in enumerate(collect(keys(graph.vprops)))
        v_info = props(graph, vertex)
        if graph_type == "circles" && vertex == center_vertex
            x[i] = y[i] = 100
        end
        msg = string("{id: ", vertex, ", label: '", vertex ,"', title: '<b>name:</b> ", string(v_info[:name]), "', x:", x[i], ", y:", y[i], ", color: {border: '", v_info[:colour], "', background: ", "'#999DA0", "'}}")
        if i!=length(collect(keys(graph.vprops)))
            print(graph_file, string(msg, ","))
        else
            print(graph_file, string(msg, "];"))
        end
    end

    close(graph_file)
end

max_stars_repo_name: kozvojtex/kozvojtexMGvisjs

max_issues_repo_name: kozvojtex/kozvojtexMGvisjs

max_forks_repo_name: kozvojtex/kozvojtexMGvisjs


text:
# Cost functions
"""
    mae(ŷ, y)

Return the mean of absolute error `sum(abs.(ŷ .- y)) / length(y)`
"""
mae(ŷ, y) = sum(abs.(ŷ .- y)) * 1 // length(y)

"""
    mse(ŷ, y)

Return the mean squared error `sum((ŷ .- y).^2) / length(y)`.
"""
mse(ŷ, y) = sum((ŷ .- y).^2) * 1 // length(y)

"""
    msle(ŷ, y; ϵ=eps(eltype(ŷ)))

Returns the mean of the squared logarithmic errors `sum((log.(ŷ .+ ϵ) .- log.(y .+ ϵ)).^2) / length(y)`.
The `ϵ` term provides numerical stability.

This error penalizes an under-predicted estimate greater than an over-predicted estimate.
"""
msle(ŷ, y; ϵ=eps(eltype(ŷ))) = sum((log.(ŷ .+ ϵ) .- log.(y .+ ϵ)).^2) * 1 // length(y)

"""
    huber_loss(ŷ, y; δ=1.0)

Computes the mean of the Huber loss given the prediction `ŷ` and true values `y`. By default, δ is set to 1.0.

                    | 0.5*|ŷ - y|,   for |ŷ - y| <= δ
      Hubber loss = |
                    |  δ*(|ŷ - y| - 0.5*δ),  otherwise

[`Huber Loss`](https://en.wikipedia.org/wiki/Huber_loss).
"""
function huber_loss(ŷ, y;  δ=eltype(ŷ)(1))
   abs_error = abs.(ŷ .- y)
   temp = abs_error .<  δ
   x = eltype(ŷ)(0.5)
   hub_loss = sum(((abs_error.^2) .* temp) .* x .+ δ*(abs_error .- x*δ) .* (1 .- temp)) * 1 // length(y)
end

function _crossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat, weight::Nothing)
  return -sum(y .* log.(ŷ)) * 1 // size(y, 2)
end

function _crossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat, weight::Number)
  return -sum(y .* log.(ŷ)) .* weight * 1 // size(y, 2)
end

function _crossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat, weight::AbstractVector)
  return -sum(y .* log.(ŷ) .* weight) * 1 // size(y, 2)
end

"""
    crossentropy(ŷ, y; weight=1)

Return the crossentropy computed as `-sum(y .* log.(ŷ) .* weight) / size(y, 2)`.

See also [`logitcrossentropy`](@ref), [`binarycrossentropy`](@ref).
"""
crossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat; weight=nothing) = _crossentropy(ŷ, y, weight)

"""
    logitcrossentropy(ŷ, y; weight=1)

Return the crossentropy computed after a [softmax](@ref) operation:

  -sum(y .* logsoftmax(ŷ) .* weight) / size(y, 2)

See also [`crossentropy`](@ref), [`binarycrossentropy`](@ref).
"""
function logitcrossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat; weight = 1)
  return -sum(y .* logsoftmax(ŷ) .* weight) * 1 // size(y, 2)
end

"""
    binarycrossentropy(ŷ, y; ϵ=eps(ŷ))

Return `-y*log(ŷ + ϵ) - (1-y)*log(1-ŷ + ϵ)`. The ϵ term provides numerical stability.

Typically, the prediction `ŷ` is given by the output of a [`sigmoid`](@ref) activation.
"""
binarycrossentropy(ŷ, y; ϵ=eps(ŷ)) = -y*log(ŷ + ϵ) - (1 - y)*log(1 - ŷ + ϵ)

# Re-definition to fix interaction with CuArrays.
CuArrays.@cufunc binarycrossentropy(ŷ, y; ϵ=eps(ŷ)) = -y*log(ŷ + ϵ) - (1 - y)*log(1 - ŷ + ϵ)

"""
    logitbinarycrossentropy(ŷ, y)

`logitbinarycrossentropy(ŷ, y)` is mathematically equivalent to `binarycrossentropy(σ(ŷ), y)`
but it is more numerically stable.

See also [`binarycrossentropy`](@ref), [`sigmoid`](@ref), [`logsigmoid`](@ref).
"""
logitbinarycrossentropy(ŷ, y) = (1 - y)*ŷ - logσ(ŷ)

# Re-definition to fix interaction with CuArrays.
CuArrays.@cufunc logitbinarycrossentropy(ŷ, y) = (1 - y)*ŷ - logσ(ŷ)

"""
    normalise(x; dims=1)

Normalises `x` to mean 0 and standard deviation 1, across the dimensions given by `dims`. Defaults to normalising over columns.

```julia-repl
julia> a = reshape(collect(1:9), 3, 3)
3×3 Array{Int64,2}:
  1  4  7
  2  5  8
  3  6  9

julia> normalise(a)
3×3 Array{Float64,2}:
  -1.22474  -1.22474  -1.22474
  0.0       0.0       0.0
  1.22474   1.22474   1.22474

julia> normalise(a, dims=2)
3×3 Array{Float64,2}:
  -1.22474  0.0  1.22474
  -1.22474  0.0  1.22474
  -1.22474  0.0  1.22474
```
"""
function normalise(x::AbstractArray; dims=1)
  μ′ = mean(x, dims = dims)
  σ′ = std(x, dims = dims, mean = μ′, corrected=false)
  return (x .- μ′) ./ σ′
end

"""
    kldivergence(ŷ, y)

KLDivergence is a measure of how much one probability distribution is different from the other.
It is always non-negative and zero only when both the distributions are equal everywhere.

[KL Divergence](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence).
"""
function kldivergence(ŷ, y)
  entropy = sum(y .* log.(y)) * 1 //size(y,2)
  cross_entropy = crossentropy(ŷ, y)
  return entropy + cross_entropy
end

"""
    poisson(ŷ, y)

Poisson loss function is a measure of how the predicted distribution diverges from the expected distribution.
Returns `sum(ŷ .- y .* log.(ŷ)) / size(y, 2)`

[Poisson Loss](https://peltarion.com/knowledge-center/documentation/modeling-view/build-an-ai-model/loss-functions/poisson).
"""
poisson(ŷ, y) = sum(ŷ .- y .* log.(ŷ)) * 1 // size(y,2)

"""
    hinge(ŷ, y)

Measures the loss given the prediction `ŷ` and true labels `y` (containing 1 or -1).
Returns `sum((max.(0, 1 .- ŷ .* y))) / size(y, 2)`

[Hinge Loss](https://en.wikipedia.org/wiki/Hinge_loss)
See also [`squared_hinge`](@ref).
"""
hinge(ŷ, y) = sum(max.(0, 1 .-  ŷ .* y)) * 1 // size(y, 2)

"""
    squared_hinge(ŷ, y)

Computes squared hinge loss given the prediction `ŷ` and true labels `y` (conatining 1 or -1).
Returns `sum((max.(0, 1 .- ŷ .* y)).^2) / size(y, 2)`

See also [`hinge`](@ref).
"""
squared_hinge(ŷ, y) = sum((max.(0, 1 .- ŷ .* y)).^2) * 1 // size(y, 2)

"""
    dice_coeff_loss(ŷ, y; smooth=1)

Loss function used in Image Segmentation. Calculates loss based on dice coefficient. Similar to F1_score.
Returns `1 - 2*sum(|ŷ .* y| + smooth) / (sum(ŷ.^2) + sum(y.^2) + smooth)`

[V-Net: Fully Convolutional Neural Networks forVolumetric Medical Image Segmentation](https://arxiv.org/pdf/1606.04797v1.pdf)
"""
dice_coeff_loss(ŷ, y; smooth=eltype(ŷ)(1.0)) = 1 - (2*sum(y .* ŷ) + smooth) / (sum(y.^2) + sum(ŷ.^2) + smooth)

"""
    tversky_loss(ŷ, y; β=0.7)

Used with imbalanced data to give more weightage to False negatives.
Larger β weigh recall higher than precision (by placing more emphasis on false negatives)
Returns `1 - sum(|y .* ŷ| + 1) / (sum(y .* ŷ + β*(1 .- y) .* ŷ + (1 - β)*y .* (1 .- ŷ)) + 1)`

[Tversky loss function for image segmentation using 3D fully convolutional deep networks](https://arxiv.org/pdf/1706.05721.pdf)
"""
tversky_loss(ŷ, y; β=eltype(ŷ)(0.7)) = 1 - (sum(y .* ŷ) + 1) / (sum(y .* ŷ + β*(1 .- y) .* ŷ + (1 - β)*y .* (1 .- ŷ)) + 1)

"""
    flatten(x::AbstractArray)

Transforms (w,h,c,b)-shaped input into (w x h x c,b)-shaped output,
by linearizing all values for each element in the batch.
"""
function flatten(x::AbstractArray)
  return reshape(x, :, size(x)[end])
end


text:
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(The BTG Program)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#################################################################################
#copy this box into a file to use the BTG program
using LinearAlgebra
using Test
using Distributions
using DataFrames
using CSV
using Polynomials
include("../range_selector.jl")
include("../validation/det_loocv.jl")
include("../validation/loocv.jl")
include("../quadrature/quadrature.jl")
include("../transforms/transforms.jl")
include("../priors/priors.jl")
include("../various_mle_things/incremental.jl")
include("../kernels/kernel.jl")
include("../datastructs.jl")
include("../computation/buffers0.jl") #datastruct, kernel, incremental, quadrature
include("../model0.jl") #buffers, datastructs, several auxiliary
include("../computation/finitedifference.jl")
include("../computation/tdist.jl") #model0 and buffer0
#################################################################################
using Plots

df = DataFrame(CSV.File("../datasets/abalone.csv"))
data = convert(Matrix, df[:,2:8]) #length, diameter, height, whole weight, shucked weight, viscera weight, shell weight
target = convert(Array, df[:, 9]) #age
normalizing_constant = maximum(target)
target = target/normalizing_constant #normalization

#it appears that only btg2 test is working
#~~~~~~~~~~~~~~~~~~~~~~~~~~(to test or not to test)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
test_btg1 = false
test_btg2 = true #used to test derivatives of cdf w.r.t augmented label-location vector [u, s]
test_btg3 = false
test_btg4 = false #weird finite difference behavior when n >1000

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~(btg1)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Attributes:
# - single length scale
# - 2-dimensional covariates
# - 3-multidimensional location vectors

ind = 120:125
posx = 1:3 #
posc = 1:3
x = data[ind, posx]
#choose a subset of variables to be regressors for the mean

Fx = data[ind, posc]
y = float(target[ind])
pind = 10:10 #prediction index
trainingData1 = trainingData(x, Fx, y) #training data used for testing various functions
#prediction data

d = getDimension(trainingData1); n = getNumPts(trainingData1); p = getCovDimension(trainingData1)
Fx0 = reshape(data[pind, posc], 1, length(posc))
x0 = reshape(data[pind, posx], 1, length(posx))
rangeθ = [100.0 200]
rangeλ = [0.5 5] #we will always used 1 range scale for lambda
btg1 = btg(trainingData1, rangeθ, rangeλ)
θ1 = btg1.nodesWeightsθ.nodes[1, 6] #pick some theta value which doubles as quadrature node
λ1 = btg1.nodesWeightsλ.nodes[3]
##################################################################################################

if (test_btg1)  #test comp_tdist
    (dpdf, pdf, cdf) = comp_tdist(btg1, θ1, λ1)
    dpdf_fixed = y0 -> dpdf(x0, Fx0, y0)
    pdf_fixed = y0 -> pdf(x0, Fx0, y0)
    cdf_fixed = y0 -> cdf(x0, Fx0, y0)
    a = pdf_fixed
    b = cdf_fixed
    c = dpdf_fixed
    #rr = 8:16
    #h = abs.(b.(.5 .+ [2.718 ^(-i) for i=rr]) .- b.(.5) - [2.718 ^(-i) for i = rr] .* a(.5))
    #xx = [2.718 ^(-i) for i=rr]
    #lxx = log.(xx)
    #lh = log.(h)
    #polyfit(lxx, lh, 1)

    #buf = btg1.test_buffer_dict[θ1]
    #println(buf.Eθ)
    #unpack(buf)

    #a = unpack(btg1.test_buffer_dict[θ1])
    #println("unpacked val in test_tdist: ", a)
    (_, _, _, pol1) = checkDerivative(a, c, 0.5, nothing, 5, 13, 10) #function first, then derivative
    (_, _, _, pol2) = checkDerivative(b, a, 0.5, nothing, 5, 13, 10) #function first, then derivative
    @testset "comp_tdist" begin
        @test coeffs(pol1)[end] ≈ 2 atol = 1e-1
        @test coeffs(pol2)[end] ≈ 2 atol = 1e-1
    end

if true #test bayesian predictive distribution (pdf, cdf, pdf_deriv)
    #rangeθ = [2.0 5; 4 7; 5 10]  #number of length scales is height of rangeθ

    (pdf, cdf, dpdf) = solve(btg1)
    a = y0 -> dpdf(x0, Fx0, y0)
    b = y0 -> pdf(x0, Fx0, y0)
    c = y0 -> cdf(x0, Fx0, y0)

    #println("Checking derivative of btg_pdf...")
    (_, _, plt1, pol1) = checkDerivative(b, a, 0.5, nothing, 4, 12, 10) #first arg is function, second arg is derivative
    (_, _, plt2, pol2) = checkDerivative(c, b, 0.5, nothing, 4, 12, 10) #first arg is function, second arg is derivative

    @testset "comp_btg_dist" begin
        @test coeffs(pol1)[end] ≈ 2 atol = 3e-1
        @test coeffs(pol2)[end] ≈ 2 atol = 3e-1
    end
    #display(plt1)
    #println(pol)
    #println("Plotting dpdf, pdf, and cdf...")
    if false
        plt(a, .1, 1, 100)
        plt(b, .1, 1, 100)
        plt!(c, .1, 1, 100, title = "PDF and CDF of Bayesian Predictive Distribution")
    end
end

end

if test_btg2 #test btg2
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~(btg2)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Attributes:
# - 2 length scales
# - 3-dimensional covariate vectors
# - 2-dimensional location vectors
# polynomial mean basis

ind = 120:140
posx = 1:7 #
posc = 1:7
x = data[ind, posx]
@info "recommended theta range: ", select_single_theta_range(x)
#choose a subset of variables to be regressors for the mean
POLYNOMIAL_BASIS_COVARIATES = true
if POLYNOMIAL_BASIS_COVARIATES
    Fx = hcat(ones(length(ind), 1), x)
else
    Fx = data[ind, posc]
end

y = float(target[ind])
pind = 10:10 #prediction index
trainingData1 = trainingData(x, Fx, y) #training data used for testing various functions

d = getDimension(trainingData1); n = getNumPts(trainingData1); p = getCovDimension(trainingData1)
x0 = reshape(data[pind, posx], 1, length(posx))
if POLYNOMIAL_BASIS_COVARIATES
    Fx0 = hcat(1, x0)
else
    Fx0 = reshape(data[pind, posc], 1, length(posc))
end
rangeθ = [10 1000.0]
rangeλ = [0.2 2] #we will always used 1 range scale for lambda
btg2 = btg(trainingData1, rangeθ, rangeλ);
θ2 = btg2.nodesWeightsθ.nodes[6] #pick some theta value which doubles as quadrature node
λ2 = btg2.nodesWeightsλ.nodes[1]
##################################################################################################

if false  #test comp_tdist for btg2
    (dpdf, pdf, cdf) = comp_tdist(btg2, θ2, λ2)
    dpdf_fixed = y0 -> dpdf(x0, Fx0, y0)
    pdf_fixed = y0 -> pdf(x0, Fx0, y0)
    cdf_fixed = y0 -> cdf(x0, Fx0, y0)
    a = pdf_fixed
    b = cdf_fixed
    c = dpdf_fixed
    #plot_multiple(a, b, c) #plot pdf, cdf, dpdf
    (_, _, plt1, pol1) = checkDerivative(a, c, 0.5, nothing, 3, 9, 10) #function first, then derivative
    (_, _, plt2, pol2) = checkDerivative(b, a, 0.5, nothing, 3, 9, 10) #function first, then derivative
    @testset "comp_tdist2" begin
        @test coeffs(pol1)[end] ≈ 2 atol = 3e-1
        @test coeffs(pol2)[end] ≈ 2 atol = 3e-1
    end
end

if false #test bayesian predictive distribution (pdf, cdf, pdf_deriv) for btg 2
    #rangeθ = [2.0 5; 4 7; 5 10]  #number of length scales is height of rangeθ
    (pdf, cdf, dpdf) = solve(btg2)
    a = y0 -> dpdf(x0, Fx0, y0)
    b = y0 -> pdf(x0, Fx0, y0)
    c = y0 -> cdf(x0, Fx0, y0)
    #println("Checking derivative of btg_pdf...")
    (_, _, plt1, pol1) = checkDerivative(b, a, 0.5, nothing, 8, 16, 10) #first arg is function, second arg is derivative
    (_, _, plt2, pol2) = checkDerivative(c, b, 0.5, nothing, 4, 10, 10) #first arg is function, second arg is derivative
    @testset "comp_btg_dist2" begin
        @test coeffs(pol1)[end] ≈ 2 atol = 3e-1
        @test coeffs(pol2)[end] ≈ 2 atol = 3e-1
    end
    #display(plt1)
    #println(pol)
    #println("Plotting dpdf, pdf, and cdf...")
    if false
        plt(a, .1, 1, 100)
        plt(b, .1, 1, 100)
        plt!(c, .1, 1, 100, title = "PDF and CDF of Bayesian Predictive Distribution")
    end
end

if true
    #if true #test derivatives of new function
        println("test derivative of cdf/cdf-related quantities w.r.t augmented vector [u, s]...\n")
        #v is augmented vector
            function linear_polynomial_basis(x0)#linear polynomial basis
                return hcat([1], reshape(x0, 1, length(x0)))
            end
            #val = [0.45, 0.375, 0.115, 0.4105]
            (dpdf, pdf, cdf, cdf_grad_us, cdf_hess_us) = comp_tdist(btg2, θ2, λ2)
           # A = v -> cdf((tmp = v[2:end]; reshape(tmp, 1, length(tmp))), (tmp = Fx0(v[2:end]); reshape(tmp, 1, length(tmp))), v[1])
            #B = v -> cdf_grad_us(v[2:end], Fx0(v[2:end]), v[1])
            #y=.26
            if false
            #@testset "augmented comp_tdist test1" begin #close to actual data point
                y0 = .26
                x0 = [0.35 0.275 0.1 0.2355 0.0895 0.0585 0.08]
                A = au -> cdf_grad_us(au[2:end], linear_polynomial_basis(au[2:end]), au[1]);
                B = au -> cdf(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                C = au -> cdf_hess_us(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                init = hcat(y0, x0)
                (_, _, plt1, pol1) = checkDerivative(B, A, init, nothing, 5, 11, 10); #first arg is function, second arg is derivative
                #@test coeffs(pol1)[end] > 2 - 3e-1
                println("pol1:",pol1)
                (_, _, plt2, pol2) = checkDerivative(B, A, init, C, 5, 11, 10) #first arg is function, second arg is derivative
                println("pol1:",pol2)
                #@test coeffs(pol2)[end]> 3 - 3e-1
            #end

            @testset "augmented comp_tdist test2 " begin #perturbation to y
                y0 = .8
                x0 = [0.35 0.275 0.1 0.2355 0.0895 0.0585 0.08]
                A = au -> cdf_grad_us(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                B = au -> cdf(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                C = au -> cdf_hess_us(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                init = hcat(y0, x0)
                (_, _, plt1, pol1) = checkDerivative(B, A, init, nothing, 5, 11, 10) #first arg is function, second arg is derivative
                @test coeffs(pol1)[end] > 2 - 3e-1
                println("pol1:",pol1)
                (_, _, plt2, pol2) = checkDerivative(B, A, init, C, 5, 11, 10) #first arg is function, second arg is derivative
                println("pol1:",pol2)
                @test coeffs(pol2)[end]> 3 - 3e-1
            end

            @testset "augmented comp_tdist test3 " begin #perturbation to y and x
                y0 = .8
                x0 = [0.3 0.4 0.05 0.3 0.03 0.05 0.15]
                A = au -> cdf_grad_us(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                B = au -> cdf(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                C = au -> cdf_hess_us(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                init = hcat(y0, x0)
                (_, _, plt1, pol1) = checkDerivative(B, A, init, nothing, 5, 11, 10) #first arg is function, second arg is derivative
                @test coeffs(pol1)[end] > 2 - 3e-1
                println("pol1:",pol1)
                (_, _, plt2, pol2) = checkDerivative(B, A, init, C, 5, 11, 10) #first arg is function, second arg is derivative
                println("pol1:",pol2)
                @test coeffs(pol2)[end]> 3 - 3e-1
            end

            @testset "augmented comp_tdist test4" begin #major perturbation (.8 in all x coords) to y and x; y outside [0, 1]
                y0 = 1.2
                x0 = [0.8 0.8 0.8 0.8 0.8 0.8 0.8]
                A = au -> cdf_grad_us(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                B = au -> cdf(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                C = au -> cdf_hess_us(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                init = hcat(y0, x0)
                (_, _, plt1, pol1) = checkDerivative(B, A, init, nothing, 5, 11, 10) #first arg is function, second arg is derivative
                @test coeffs(pol1)[end] > 2 - 3e-1
                println("pol1:",pol1)
                (_, _, plt2, pol2) = checkDerivative(B, A, init, C, 5, 11, 10) #first arg is function, second arg is derivative
                println("pol1:",pol2)
                @test coeffs(pol2)[end]> 3 - 3e-1
            end
        end
        if true
            println("test derivatives of average augmented cdf_deriv...\n ")
            @testset "btg augmented cdf_function test 1" begin
                y0 = 0.4
                x0 = [0.3 0.4 0.05 0.3 0.03 0.05 0.15]
                @timeit to "get pdf, cdf, dpdf" (pdf, cdf, dpdf, cdf_grad_us, cdf_hess_us, quantInfo) = solve(btg2; derivatives = true)
                A = au -> cdf_grad_us(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                B = au -> cdf(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                C = au -> cdf_hess_us(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                init = hcat(y0, x0)
                @timeit to "check derivative 1" (_, _, plt1, pol1) = checkDerivative(B, A, init, nothing, 5, 11, 10) #first arg is function, second arg is derivative
                @test coeffs(pol1)[end] > 2 - 3e-1
                println("pol1:",pol1)
                #@timeit to "check derivative 2" (_, _, plt2, pol2) = checkDerivative(B, A, init, C, 5, 11, 10) #first arg is function, second arg is derivative
                #println("pol1:",pol2)
                #@test coeffs(pol2)[end]> 3 - 3e-1
            end
            if true
            @testset "btg augmented cdf_function test 2" begin
                y0 =1.1
                x0 = [0.8 0.02 0.4 0.3 0.5 0.2 0.4]
                @timeit to "get pdf, cdf, dpdf" (pdf, cdf, dpdf, cdf_grad_us, cdf_hess_us, quantInfo) = solve(btg2; derivatives = true)
                A = au -> cdf_grad_us(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                B = au -> cdf(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                C = au -> cdf_hess_us(au[2:end], linear_polynomial_basis(au[2:end]), au[1])
                init = hcat(y0, x0)
                @timeit to "check derivative 3" (_, _, plt1, pol1) = checkDerivative(B, A, init, nothing, 5, 11, 10) #first arg is function, second arg is derivative
                @test coeffs(pol1)[end] > 2 - 3e-1
                println("pol1:",pol1)
                #@timeit to "check derivative 4" (_, _, plt2, pol2) = checkDerivative(B, A, init, C, 5, 11, 10) #first arg is function, second arg is derivative
                #println("pol1:",pol2)
                #@test coeffs(pol2)[end]> 3 - 3e-1
            end
        end
        end
        #end
    #end
    end
if false
#if true #test derivatives of new function
    println("test derivative of cdf/cdf-related quantities computed using comp_tdist w.r.t s")
    #v is augmented vector
        Fx0 = x0 -> hcat([1], reshape(x0, 1, length(x0))) #linear polynomial basis
        #val = [0.45, 0.375, 0.115, 0.4105]
        (dpdf, pdf, cdf, cdf_grad_us, cdf_hess_us) = comp_tdist(btg2, θ2, λ2)
       # A = v -> cdf((tmp = v[2:end]; reshape(tmp, 1, length(tmp))), (tmp = Fx0(v[2:end]); reshape(tmp, 1, length(tmp))), v[1])
        #B = v -> cdf_grad_us(v[2:end], Fx0(v[2:end]), v[1])
        y=.26
        A = x0 -> cdf_grad_us(x0, Fx0(x0), y)
        B = x0 -> cdf(x0, Fx0(x0), y)
        C = x0 -> cdf_hess_us(x0, Fx0(x0), y)
        #init_v = (tmp = hcat(0.5, x[1:1, :] + rand(1, length(posx)) .* 1e-1); reshape(tmp, 1, length(tmp)))
        #init_v = [0.32, 0.45, 0.375, 0.115, 0.4105]
        init_v =  [0.35 0.275 0.1 0.2355 0.0895 0.0585 0.08]
        if false
            #cdf_fixed = y -> cdf([0.45, 0.375, 0.115, 0.4105], Fx0([0.45, 0.375, 0.115, 0.4105]), y)
            #cdf_fixed
            #plt(cdf_fixed, .01, 1, 100)
        end
        (_, _, plt1, pol1) = checkDerivative(B, A, init_v, nothing, 4, 10, 10); #first arg is function, second arg is derivative
        #@test coeffs(pol1)[end] > 2 - 3e-1
        println(pol1)
        #(_, _, plt2, pol2) = checkDerivative(B, A, init_v, C, 3, 13, 10); #first arg is function, second arg is derivative
        #@test coeffs(pol1)[end] > 2 - 3e-1
        #println(pol2)

        if false #temporary test for second derivative of cdf w.r.t x
        A = y -> cdf_hess_us(init_v, Fx0(init_v), y)[1]
        B = y -> cdf_hess_us(init_v, Fx0(init_v), y)[2]
        C = y -> cdf_hess_us(init_v, Fx0(init_v), y)[3]
        (_, _, plt1, pol1) = checkDerivative(A, B, 0.26, nothing, 4, 10, 10) #first arg is function, second arg is derivative
        #@test coeffs(pol1)[end] > 2 - 3e-1
        println(pol1)
        (_, _, plt2, pol2) = checkDerivative(B, C, 0.26, nothing, 3, 13, 10); #first arg is function, second arg is derivative
        #@test coeffs(pol1)[end] > 2 - 3e-1
        println(pol2)
        end
    #end
#end
end
if false
#if true #test derivatives of new function
println("test derivative of cdf w.r.t value/label u")
#v is augmented vector
    Fx0 = x0 -> hcat([1], reshape(x0, 1, length(x0))) #linear polynomial basis
    val = [.425 .3 .1 .35 .15 .09 .15]
    (pdf, dpdf, cdf, cdf_grad_us) = comp_tdist(btg2, θ2, λ2)
   # A = v -> cdf((tmp = v[2:end]; reshape(tmp, 1, length(tmp))), (tmp = Fx0(v[2:end]); reshape(tmp, 1, length(tmp))), v[1])
    #B = v -> cdf_grad_us(v[2:end], Fx0(v[2:end]), v[1])
    A = y -> cdf_grad_us(val, Fx0(val), y)
    B = y -> cdf(val, Fx0(val), y)

    #init_v = (tmp = hcat(0.5, x[1:1, :] + rand(1, length(posx)) .* 1e-1); reshape(tmp, 1, length(tmp)))
    #init_v = [0.32, 0.45, 0.375, 0.115, 0.4105]
    init_v = 0.3
    if false
        #cdf_fixed = y -> cdf([0.45, 0.375, 0.115, 0.4105], Fx0([0.45, 0.375, 0.115, 0.4105]), y)
        plt(B, .01, 1, 100)
    end
    (_, _, plt1, pol1) = checkDerivative(B, A, init_v, nothing, 4, 15, 10) #first arg is function, second arg is derivative
    #@test coeffs(pol1)[end] > 2 - 3e-1
    println(pol1)
end

end
if test_btg3 #test btg3
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~(btg3)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Attributes:
# - 3 length scales
# - 3-D covariate vectors
# - 2-D location vectors

ind = 120:135
posx = [1;3;4] #
posc = 1:3
x = data[ind, posx]
#choose a subset of variables to be regressors for the mean
Fx = data[ind, posc]
y = float(target[ind])
pind = 10:10 #prediction index
trainingData1 = trainingData(x, Fx, y) #training data used for testing various functions

d = getDimension(trainingData1); n = getNumPts(trainingData1); p = getCovDimension(trainingData1)
Fx0 = reshape(data[pind, posc], 1, length(posc))
x0 = reshape(data[pind, posx], 1, length(posx))
rangeθ = [100.0 200; 200.0 400.0; 100.0 1000.0]
rangeλ = [0.5 5] #we will always used 1 range scale for lambda
btg3 = btg(trainingData1, rangeθ, rangeλ)
θ3 = btg3.nodesWeightsθ.nodes[1:3, 6] #pick some theta value which doubles as quadrature node
λ3 = btg3.nodesWeightsλ.nodes[3]
##################################################################################################

if true  #test comp_tdist for btg3
    (dpdf, pdf, cdf) = comp_tdist(btg3, θ3, λ3)
    dpdf_fixed = y0 -> dpdf(x0, Fx0, y0)
    pdf_fixed = y0 -> pdf(x0, Fx0, y0)
    cdf_fixed = y0 -> cdf(x0, Fx0, y0)
    a = pdf_fixed
    b = cdf_fixed
    c = dpdf_fixed
    (_, _, _, pol1) = checkDerivative(a, c, 0.5, nothing, 4, 12, 10) #function first, then derivative
    (_, _, _, pol2) = checkDerivative(b, a, 0.5, nothing, 4, 12, 10) #function first, then derivative
    @testset "comp_tdist3" begin
        @test coeffs(pol1)[end] ≈ 2 atol = 3e-1
        @test coeffs(pol2)[end] ≈ 2 atol = 3e-1
    end
end

if true #test bayesian predictive distribution (pdf, cdf, pdf_deriv) for btg 3
    #rangeθ = [2.0 5; 4 7; 5 10]  #number of length scales is height of rangeθ
    (pdf, cdf, dpdf) = solve(btg3)
    a = y0 -> dpdf(x0, Fx0, y0)
    b = y0 -> pdf(x0, Fx0, y0)
    c = y0 -> cdf(x0, Fx0, y0)
    #println("Checking derivative of btg_pdf...")
    (_, _, plt1, pol1) = checkDerivative(b, a, 0.5, nothing, 4, 12, 10) #first arg is function, second arg is derivative
    (_, _, plt2, pol2) = checkDerivative(c, b, 0.5, nothing, 4, 12, 10) #first arg is function, second arg is derivative
    @testset "comp_btg_dist3" begin
        @test coeffs(pol1)[end] ≈ 2 atol = 1e-1
        @test coeffs(pol2)[end] ≈ 2 atol = 1e-1
    end
    #display(plt1)
    #println(pol)
    #println("Plotting dpdf, pdf, and cdf...")
    if false
        plt(a, .1, 1, 100)
        plt(b, .1, 1, 100)
        plt!(c, .1, 1, 100, title = "PDF and CDF of Bayesian Predictive Distribution")
    end
end
end

if test_btg4 #test btg4
    println("Runnings tests for BTG4....")
    ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~(btg4)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #Attributes:
    # - 3 length scales
    # - 3-D covariate vectors
    # - 2-D location vectors
    #- use Monte-Carlo integration

    ind = 120:250
    posx = [1;3;4] #
    posc = 1:3
    x = data[ind, posx]
    #choose a subset of variables to be regressors for the mean
    Fx = data[ind, posc]
    y = float(target[ind])
    pind = 10:10 #prediction index
    trainingData1 = trainingData(x, Fx, y) #training data used for testing various functions

    d = getDimension(trainingData1); n = getNumPts(trainingData1); p = getCovDimension(trainingData1)
    Fx0 = reshape(data[pind, posc], 1, length(posc))
    x0 = reshape(data[pind, posx], 1, length(posx))
    rangeθ = [100.0 200; 200.0 400.0; 100.0 1000.0]
    rangeλ = [0.5 5] #we will always used 1 range scale for lambda
    btg4 = btg(trainingData1, rangeθ, rangeλ; quadtype = ["MonteCarlo", "Gaussian"])
    θ4 = btg4.nodesWeightsθ.nodes[1:3, 6] #pick some theta value which doubles as quadrature node
    λ4 = btg4.nodesWeightsλ.nodes[3]
    ##################################################################################################

    if true  #test comp_tdist for btg4
        (dpdf, pdf, cdf) = comp_tdist(btg4, θ4, λ4)
        dpdf_fixed = y0 -> dpdf(x0, Fx0, y0)
        pdf_fixed = y0 -> pdf(x0, Fx0, y0)
        cdf_fixed = y0 -> cdf(x0, Fx0, y0)
        a = pdf_fixed
        b = cdf_fixed
        c = dpdf_fixed
        (_, _, plt1, pol1) = checkDerivative(a, c, 0.5, nothing, 4, 12, 20) #function first, then derivative
        (_, _, plt2, pol2) = checkDerivative(b, a, 0.5, nothing, 4, 12, 20) #function first, then derivative
        @testset "comp_tdist4" begin
            @test coeffs(pol1)[end] ≈ 2 atol = 3e-1
            @test coeffs(pol2)[end] ≈ 2 atol = 3e-1
        end
    end

    if true #test bayesian predictive distribution (pdf, cdf, pdf_deriv) for btg 4
        #rangeθ = [2.0 5; 4 7; 5 10]  #number of length scales is height of rangeθ
        (pdf, cdf, dpdf) = solve(btg4)
        a = y0 -> dpdf(x0, Fx0, y0)
        b = y0 -> pdf(x0, Fx0, y0)
        c = y0 -> cdf(x0, Fx0, y0)
        #println("Checking derivative of btg_pdf...")
        (_, _, plt1, pol1) = checkDerivative(b, a, 0.5, nothing, 4, 12, 10) #first arg is function, second arg is derivative
        (_, _, plt2, pol2) = checkDerivative(c, b, 0.5, nothing, 4, 12, 10) #first arg is function, second arg is derivative
        @testset "comp_btg_dist4" begin
            @test coeffs(pol1)[end] ≈ 2 atol = 3e-1
            @test coeffs(pol2)[end] ≈ 2 atol = 3e-1
        end
        #display(plt1)
        #println(pol)
        #println("Plotting dpdf, pdf, and cdf...")
        if false
            plt(a, .1, 1, 100)
            plt(b, .1, 1, 100)
            plt!(c, .1, 1, 100, title = "PDF and CDF of Bayesian Predictive Distribution")
        end
    end
    end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
function info(eventlog::EventLog)
    n_traces = length(eventlog.traces)
    if n_traces == 1
        infotext = "EventLog with " * string(n_traces) * " trace."
    else
        infotext = "EventLog with " * string(n_traces) * " traces."
    end
    print(infotext)
    return infotext
end

Base.length(eventlog::EventLog) = length(eventlog.traces)

Base.length(trace::Trace) = length(trace.events)

function Base.filter(eventlog::EventLog, filter_function::Function, field::Symbol)
    filtered_eventlog = deepcopy(eventlog)
    trace_filter = map(filter_function, getproperty(filtered_eventlog, field))
    filtered_field = getproperty(filtered_eventlog, field)[trace_filter]
    setproperty!(filtered_eventlog, field, filtered_field)
    return filtered_eventlog
end

function Base.filter!(eventlog::EventLog, filter_function::Function, field::Symbol)
    trace_filter = map(filter_function, getproperty(eventlog, field))
    filtered_field = getproperty(eventlog, field)[trace_filter]
    setproperty!(eventlog, field, filtered_field)
    return eventlog
end


text:
include("to_include.jl")

println(f_aux(2,3))


text:
for i in results[3]["result"]
    print(i["sourceUrl"], "\n")
end

protocols = []
for i in results
    tmp = [j["sourceUrl"] for j in i["result"]]
    append!(protocols, tmp)
end

unique(protocols)

l = unique(deepcopy(results[!, :domain]))

for i in l
    print(i, "\n")
end

t = unique(
    [
        i
        for i in l
            if !startswith(i, "www")
                & !startswith(i, "de")
                & !startswith(i, "en")
                & !startswith(i, "de-de")
                & !startswith(i, "nl")
                & !startswith(i, "uk")
                & !startswith(i, "m.")
                & !startswith(i, "it")
                & !startswith(i, "it-it")
                & !startswith(i, "pt")
                & !startswith(i, "pl")
    ]
)

for i in t
    print(i, "\n")
end


text:
### Triangle ###
function asMatOnePad(v1::SVector{3,T}, v2::SVector{3,T}, v3::SVector{3,T}) where {T}
    A = @SMatrix [
    v1[1] v2[1] v3[1];
    v1[2] v2[2] v3[2];
    v1[3] v2[3] v3[3];
    one(T)  one(T)  one(T)
    ]
    return A
end

function asMat(v1::SVector{3,T}, v2::SVector{3,T}, v3::SVector{3,T}) where {T}
    A = @SMatrix [
    v1[1] v2[1] v3[1];
    v1[2] v2[2] v3[2];
    v1[3] v2[3] v3[3]
    ]
    return A
end

for funName in (:asMatOnePad, :asMat)
    @eval begin
        function $funName(sv::SVector{3,SVector{3,T}}) where {T}
            return $funName(sv[1], sv[2], sv[3])
        end
        function $funName(sv::NTuple{3,SVector{3,T}}) where {T}
            return $funName(sv[1], sv[2], sv[3])
        end
    end
end

### Tetrahedron ###
function asMatOnePad(v1::SVector{3,T}, v2::SVector{3,T}, v3::SVector{3,T}, v4::SVector{3,T}) where {T}
    A = @SMatrix [
        v1[1]   v2[1]   v3[1]   v4[1];
        v1[2]   v2[2]   v3[2]   v4[2];
        v1[3]   v2[3]   v3[3]   v4[3];
        one(T)  one(T)  one(T)  one(T)
    ]
    return A
end

function asMat(v1::SVector{3,T}, v2::SVector{3,T}, v3::SVector{3,T}, v4::SVector{3,T}) where {T}
    A = @SMatrix [
        v1[1] v2[1] v3[1] v4[1];
        v1[2] v2[2] v3[2] v4[2];
        v1[3] v2[3] v3[3] v4[3]
    ]
    return A
end

for funName in (:asMatOnePad, :asMat)
    @eval begin
        function $funName(sv::SVector{4,SVector{3,T}}) where {T}
            return $funName(sv[1], sv[2], sv[3], sv[4])
        end
        function $funName(sv::NTuple{4,SVector{3,T}}) where {T}
            return $funName(sv[1], sv[2], sv[3], sv[4])
        end
    end
end


text:
"""
Run MTM problem and plot results
"""

using Term
# install_stacktrace()
install_term_logger()

import jcontrol.control: ControlOptions, Bounds
import jcontrol.bicycle: State
import jcontrol: run_mtm, Track

print("\n\n" * hLine("start"; style="bold green"))

# --------------------------------- kinematic -------------------------------- #
# track, bike, control_model, solution = run_mtm(
#     :kinematics,  # model type
#     3;  # supports density
#     control_options=:default,
#     showtrials=50,
#     n_iter=5000,
#     timed=false,
#     showplots=true,
# )

# ---------------------------------- dynamic---------------------------------- #

track = Track(;start_waypoint=2, keep_n_waypoints=-1)

# icond = State(; u=10)
fcond = State(; u=30, ω=0, δ=0)

track, bike, control_model, solution = run_mtm(
    :dynamics,  # model type
    1;  # supports density
    track=track,
    control_options=:default,
    icond=nothing,
    fcond=fcond,
    showtrials=30,
    n_iter=5000,
    timed=false,
    showplots=true,
)

print("\n", hLine("done"; style="bold blue") * "\n\n")

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
module Subsequences

export longest_common_subsequence, longest_contiguous_subsequence

# -------

function longest_common_subsequence(a, b; result_base = "", join_fn = string)
    lengths = zeros(length(a) + 1, length(b) + 1)

    for (i, x) in enumerate(a)
        for (j, y) in enumerate(b)
            if x == y
                lengths[i+1, j+1] = lengths[i, j] + 1
            else
                lengths[i+1, j+1] = max(lengths[i+1, j], lengths[i, j+1])
            end
        end
    end

    x, y = length(a) + 1, length(b) + 1
    result = result_base
    a_start, a_end = 0, 0
    b_start, b_end = 0, 0

    while x != 1 && y != 1
        if lengths[x, y] == lengths[x-1, y]
            x -= 1
        elseif lengths[x, y] == lengths[x, y-1]
            y -= 1
        else
            result = join_fn(a[x-1], result)

            a_start = x - 1
            b_start = y - 1

            if (a_end == 0)
                a_end = x - 1
                b_end = y - 1
            end

            x -= 1
            y -= 1
        end
    end

    result, a_start:a_end, b_start:b_end
end

longest_common_subsequence(a::Array, b::Array) = longest_common_subsequence(a, b, result_base = [], join_fn = vcat)

# -------

function longest_contiguous_subsequence(a, b)
    m = zeros(Int, length(a) + 1, length(b) + 1)
    longest, x_longest, y_longest = 0, 0, 0

    for x in 2:(length(a) + 1)
        for y in 2:(length(b) + 1)

            if a[x - 1] == b[y - 1]
                m[x, y] = m[x - 1, y - 1] + 1

                if m[x, y] > longest
                    longest = m[x, y]
                    x_longest = x
                    y_longest = y
                end
            else
                m[x, y] = 0
            end

        end
    end

    a_range = (x_longest - longest):(x_longest - 1)
    b_range = (y_longest - longest):(y_longest - 1)
    a[a_range], a_range, b_range
end

end


text:
### Concrete Types

type ODECompositeSolution{uType,tType,rateType,P,A,IType} <: AbstractODESolution
  u::uType
  t::tType
  k::rateType
  prob::P
  alg::A
  interp::IType
  alg_choice::Vector{Int}
  dense::Bool
  tslocation::Int
  retcode::Symbol
end
(sol::ODECompositeSolution)(t,deriv::Type=Val{0};idxs=nothing) = sol.interp(t,idxs,deriv)
(sol::ODECompositeSolution)(v,t,deriv::Type=Val{0};idxs=nothing) = sol.interp(v,t,idxs,deriv)

type ODECompositeTestSolution{uType,uType2,uEltype,tType,rateType,P,A,IType} <: AbstractODETestSolution
  u::uType
  u_analytic::uType2
  errors::Dict{Symbol,uEltype}
  t::tType
  k::rateType
  prob::P
  alg::A
  interp::IType
  alg_choice::Vector{Int}
  dense::Bool
  tslocation::Int
  retcode::Symbol
end
(sol::ODECompositeTestSolution)(t,deriv::Type=Val{0};idxs=nothing) = sol.interp(t,idxs,deriv)
(sol::ODECompositeTestSolution)(v,t,deriv::Type=Val{0};idxs=nothing) = sol.interp(v,t,idxs,deriv)

function build_solution{uType,tType,isinplace}(
        prob::AbstractODEProblem{uType,tType,isinplace},
        alg::OrdinaryDiffEqCompositeAlgorithm,t,u;dense=false,alg_choice=[1],
        k=[],interp = (tvals) -> nothing, retcode = :Default, kwargs...)
  ODECompositeSolution(u,t,k,prob,alg,interp,alg_choice,dense,0,retcode)
end

function build_solution{uType,tType,isinplace}(
        prob::AbstractODETestProblem{uType,tType,isinplace},
        alg::OrdinaryDiffEqCompositeAlgorithm,t,u,;dense=false,
        k=[],interp = (tvals) -> nothing,alg_choice=[1],
        timeseries_errors=true,dense_errors=true,
        calculate_error = true,retcode = :Default, kwargs...)
  u_analytic = Vector{uType}(0)
  errors = Dict{Symbol,eltype(u[1])}()
  sol = ODECompositeTestSolution(u,u_analytic,errors,t,k,prob,alg,interp,alg_choice,dense,0,retcode)
  if calculate_error
    calculate_solution_errors!(sol;timeseries_errors=timeseries_errors,dense_errors=dense_errors)
  end
  sol
end


text:
const DEFAULT_ETA = 1e-3
const DEFAULT_STEPS = 3e3

abstract type Rescaler end

struct OrderGuess
    name::String
    testloss::Real
    trainloss::Real
    x::Array
    y::Array
    preds::Array
end

struct OrderEstimator
    name::String
    initparams::Array{Float64, 1}
    funform::Base.Callable
    scaler::Base.Callable
end

function OrderEstimator(name::String, initparams::Array{Float64, 1}, funform::Base.Callable)
    OrderEstimator(name, initparams, funform, PosScaler)
end

function getscale(a::Array)
    m = mean(a)
    s = std(a) |> (x -> x > 0 ? x : 1)
    return [m, s, minimum((a .- m) ./ s) - 1]
end

struct PosScaler <: Rescaler
    m::Real
    s::Real
    m2::Real

    PosScaler(a::Array) = new(getscale(a)...)
end

struct LogScaler <: Rescaler
    m::Real
    s::Real
    m2::Real

    LogScaler(a::Array) = new(getscale(log.(a))...)
end

transform(scaler::PosScaler, a::Array) = (a .- scaler.m) ./ scaler.s .- scaler.m2
function transform(scaler::LogScaler, a::Array)
    (log.(a) .- scaler.m) ./ scaler.s .- scaler.m2
end

invert(scaler::PosScaler, a::Array) = (a .+ scaler.m2) * scaler.s .+ scaler.m
invert(scaler::LogScaler, a::Array) = exp.((a .+ scaler.m2) * scaler.s .+ scaler.m)

function orders()
    return [
    OrderEstimator(
        "O(1)",
        [1.0],
        c -> (x -> c)
        ),
    OrderEstimator(
        "O(log n)",
        [0.0, 1.0, 2.0],
        (c, b, l) -> (x -> max(0.0, c) + max(1e-2, b) * log(max(1.01, l), x))
        ),
    OrderEstimator(
        "O(n)",
        [0.0, 1.0],
        (c, b) -> (x -> max(0.0,c) + max(0.001, b) * x)
        ),
    OrderEstimator(
        "O(n log n)",
        [0.0, 0.0, 1.0, 2.0],
        (c, b1, b2, l) -> (x -> max(0.0, c) + max(0, b1) * x + max(0.001, b2) * x * log(max(1.01, l), x))
        ),
    OrderEstimator(
        "O(n^p)",
        [0.0, 0.0, 1.0, 2.0],
        (c, b1, b2, p) -> (x -> max(0.0, c) + max(0, b1) * x + max(0.001, b2) * x ^ max(1.01, p))
        ),
    OrderEstimator(
        "O(k^n)",
        [0.0, 0.0, 1.0, 2.0],
        (c, b1, b2, k) -> (x -> max(0.0, c) + max(0, b1) * x + max(0.001, b2) * max(1.001, k) ^ x),
        LogScaler
        ),
    ]
end

function plot(guess::OrderGuess)
    plot(
        guess.x,
        hcat(guess.preds, guess.y),
        labels=["preds" "test"],
        legend=:topleft,
        title=guess.name
    )
end

function _optimize!(loss, parameters, eta, steps=1e5, breaklimit=1e-15)
    for _ in 1:steps
        l = loss(parameters...)
        l < breaklimit && break
        g = gradient(loss, parameters...) ./ l
        parameters .= parameters .- (g .* eta ./ sum(abs.(g)))
    end
end

function fit(order::OrderEstimator, x, y, eta=1e-3, steps=5e4)::Base.Callable
    scaler = order.scaler(y)

    prepped_y = transform(scaler, y)
    ps = order.initparams |> copy

    predscaler = order.scaler(order.funform(ps...).(x))

    model(ps) = x -> transform(predscaler, order.funform(ps...).(x))
    loss(ps...) = sum((model(ps)(x) .- prepped_y) .^ 2)

    ps[1] += mean(prepped_y .- model(ps)(x))
    _optimize!(loss, ps, eta, steps)

    return x -> invert(scaler, model(ps)(x))
end

function guessbigo(x, y, eta=DEFAULT_ETA, steps=DEFAULT_STEPS, trainrate=0.7)::Array{OrderGuess}

    te = (length(x) * trainrate) |> round |> Int
    train_x, train_y = (a -> a[1:te]).([x, y])
    test_x, test_y = (a -> a[te+1:end]).([x, y])

    recs = []
    for order in orders()
        estimator = fit(order, train_x, train_y, eta, steps)

        testerror = sum(abs.(estimator(test_x) .- test_y))
        trainerror = sum(abs.(estimator(train_x) .- train_y))
        push!(
            recs,
            OrderGuess(order.name, testerror, trainerror, x, y, estimator(x))
        )
    end
    return sort(recs, by=x-> x.testloss)
end
guessbigo(fun::Base.Callable, x=50:25:500, eta=DEFAULT_ETA, steps=DEFAULT_STEPS, trainrate=0.7) = guessbigo(x, fun.(x), eta, steps, trainrate)

best(guesses::Array{OrderGuess}) = sort(guesses, by=x-> x.testloss)[1].name

function bigos(report::RunReport, eta=DEFAULT_ETA, steps=DEFAULT_STEPS, trainrate=0.7)
    return Dict(
        [
            (
                report.names[i],
                guessbigo(
                    report.inputsizes,
                    report.times[:, i],
                    eta,
                    steps,
                    trainrate
                    ) |> best
            )
            for i in axes(report.names, 1)
        ]
    )
end


text:
function PLacExample_VS_TestPLac_M2Utility(ins)

    # Definition of the inputs for the ODEs
    if length(ins) == 1
        IPTG2,aTc1 = ins[1];
    else
        IPTG2,aTc1 = ins;
    end
    inputsM1 = [0,IPTG2,0,IPTG2,0,IPTG2];
    inputsM2 = [0,aTc1,IPTG2,0,0,aTc1,IPTG2,0,0,aTc1,IPTG2,0];

    # Solve ODEs

    solM1 = PLacExample_SolveAll(tsMS, pD1MS, spMS, inputsM1, ivss1MS, sampsMS, pre1MS);
    solM2 = TestPLac_M2_SolveAll(tsMS, pD2MS, spMS, inputsM2, ivss2MS, sampsMS, pre2MS);

    # Extracte wanted vectors (observables) with time reduction
    Obs1_M1 = 3 .*solM1[:,4,:];

    Obs1_M2 = 3 .*solM2[:,4,:];

    mu1_M1 = mean(Obs1_M1, dims=2);

    sd1_M1 = cov(Obs1_M1, dims=2).+(0.1*Matrix((Diagonal(ones(size(Obs1_M1)[1])))));

    mu1_M2 = mean(Obs1_M2, dims=2);

    sd1_M2 = cov(Obs1_M2, dims=2).+(0.1*Matrix((Diagonal(ones(size(Obs1_M2)[1])))));

    BHD1 = BhattacharyyaDist(mu1_M1[:,1], mu1_M2[:,1], sd1_M1, sd1_M2);

    util = mean(BHD1);

    return(util)
end

max_stars_repo_name: DavidGomezC/BOMBSbeta

max_issues_repo_name: DavidGomezC/BOMBSbeta

max_forks_repo_name: DavidGomezC/BOMBSbeta


text:
abstract abstract_newton_direction # see src/newton_direction/

type class_settings
	max_it::Int64

	kappa_tau_tol::Float64
	gap_tol::Float64

	primal_tol::Float64
	dual_tol::Float64

	primal_infeas_tol::Float64
	dual_infeas_tol::Float64
	unbounded_value::Float64

	diagonal_modification::Float64

	linear_system_solver::abstract_linear_system_solver
	newton_solver::abstract_newton_direction

	delta_min::Float64
	delta_max::Float64
	delta_start::Float64
	delta_increase::Float64
	delta_decrease::Float64

	verbose::Bool

	function class_settings()
		this = new();

		# termination criteron
		this.max_it = 600;

		this.kappa_tau_tol = 1e-8
		this.gap_tol = 1e-8;
		this.primal_tol = 1e-8;
		this.dual_tol = 1e-8;

		this.primal_infeas_tol = 1e-8;
		this.dual_infeas_tol = 1e-8;
		this.unbounded_value = 1e4;

		# algorithm parameters
		this.diagonal_modification = 1e-8;

		# algorithm strategies
		#this.linear_system_solver = linear_solver_MATLAB();
		this.linear_system_solver = linear_solver_MUMPS();
		#this.newton_solver = class_newton_hsd();
		this.newton_solver = class_homogeneous_newton();

		# delta parameters
		this.delta_min = 1e-8;
		this.delta_max = 1e20;
		this.delta_start = 1.1;
		this.delta_increase = 8;
		this.delta_decrease = 1.0/3.0;

		this.verbose = true;

		return this
	end
end

max_forks_repo_name: ohinder/Everest


text:
# Spectra and time series plots from virtual satellite data.
#
# t rho vx vy p bz ex ey
#
# Hongyang Zhou, hyzhou@umich.edu

using SignalAnalysis, DelimitedFiles, Glob
using SignalAnalysis.Units
using Plots, Measures

files = glob("satellite*.csv")

samplerate = 2Hz
color = :dense
left_margin = 5mm

const ax = fill(plot(),10,1)

for file in files
   println("file: ", file)
   data = readdlm(file, ','; header=true)
   t0 = data[1][1]
   p = signal(data[1][:,5] .* 1e9, samplerate) # [nPa]
   bz = signal(data[1][:,6] .* 1e9, samplerate)# [nT]
   ex = signal(data[1][:,7] .* 1e3, samplerate)# [mV/m]
   ey = signal(data[1][:,8] .* 1e3, samplerate)# [mV/m]

   l = @layout[[a1;a2;a3;a4;a5] grid(2,2){0.7w}]
   #ax[1] = plot(title="Power Spectra at [$(file[end-5:end-4]), 0, 0]Re",
   #   framestyle=nothing,showaxis=false,xticks=false,yticks=false,margin=0mm)

   x = parse(Int, file[end-5:end-4])
   ax[2] = plot(data[1][:,1], data[1][:,2] ./ 1e6; label="density", title="[$x, 0, 0]Re",
      ylabel="[amu/cc]", left_margin)
   ax[3] = plot(data[1][:,1], data[1][:,3:4] ./ 1e3; label=["Vx" "Vy"],
      ylabel="[km/s]", left_margin)
   ax[4] = plot(data[1][:,1], samples(p); label="p", ylabel="[nPa]", left_margin)
   ax[5] = plot(data[1][:,1], samples(bz); label="bz", ylabel="[nT]", left_margin)
   ax[6] = plot(data[1][:,1], data[1][:,7:8] .* 1e3; label=["Ex" "Ey"],
      ylabel="[mV/m]", xlabel="Time [s]", left_margin)

   # magnitude units in [db], or 10log10(power)
   ax[7] = specgram(p; nfft=128, t0, crange=100, color ,title="P")
   ax[8] = specgram(bz; nfft=128, t0, crange=100, color, title="Bz")
   ax[9] = specgram(ex; nfft=128, t0, crange=100, color, title="Ex")
   ax[10] = specgram(ey; nfft=128, t0, crange=100, color, title="Ey")
   plot(ax[2:10]..., layout=l, size=(1500,1000), legend=true)
   savefig("spectrum_"*file[end-5:end-4]*"Re.png")
end


text:
# For this PS3, I primary used Prof.Ransom's code as a referecne point to write down my code
# Also, I tried to write explantation with each line of code to better understand what is going on
# Also Prof. Ransom mentioned that its extremely imp to visualize what is going on before writng a loop - This has helped me quite a bit

cd("D:\\Semester 3_Fall 2020\\ECON 6453 - Econometrics III\\ProblemSets\\PS3-gev")
using Random, Distributions, JLD2, DataFrames, CSV, LinearAlgebra, Statistics, FreqTables, TexTables, Optim, HTTP, GLM, ForwardDiff

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 4: All Wrapping Function
#:::::::::::::::::::::::::::::::::::::::::::::::::::
function problemset_3()

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 1: Multinomial logit with Alternate Specific Covariates
#:::::::::::::::::::::::::::::::::::::::::::::::::::
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS3-gev/nlsw88w.csv"

df = CSV.read(HTTP.get(url).body)

X = [df.age df.white df.collgrad]
Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4,df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
y = df.occupation

function mlogit(alpha, X, Z, y)
    beta = alpha[1:end-1] # Indexing the parameter so coefficient on X is the first element
    gamma = alpha[end] # Indexing the parameter so coefficient on Z is the last element
    K = size(X,2) # K is the no. of covariates in X, which in this case are 3
    J = length(unique(y)) # J is the no. of coefficients, which are 8 in this case (So J-1 = 7)
    N = length(y) # N is a vector with the same length as y (not sure why we are defining N)
    bigY = zeros(N,J) # Simply defining the matrix of dependent variables Y, which are going to be dummies in case of Multinomial Logit

    for j=1:J # Simply defining j to be something, next step, j represents column
        bigY[:,j] = y.==j # Making a dummies for each of the corresponding occupation in y
    end
    bigBeta = [reshape(beta,K,J-1) zeros(K)] # Reshaping into a matrix with K rows and J-1 columns
    # Adding a column of zeros in the end to make it conform with bigY and also to normalize the parameters of the last occupation to be zero

    T = promote_type(eltype(X),eltype(alpha)) # To get the standard errors we promote the type
    num   = zeros(T,N,J) # Define numerator T-type , N-rows, J-columns
    dem   = zeros(T,N) # Define denominator T-type , N-rows
    for j=1:J
        num[:,j] = exp.(X*bigBeta[:,j] .+ (Z[:,j] .- Z[:,J])*gamma) # Numerator Likelihood function
        dem .+= num[:,j] # Denominator Likelihood function: Simply the Sum of all the numerators
    end

    P = num./repeat(dem,1,J) # Putting together numerator and denominator

    loglike = -sum( bigY.*log.(P) ) # Since Optim minimizes the function, we take its negative to maximize

    return loglike # This is used for storing the function parameters
end

startvals = [2*rand(7*size(X,2)).-1; .1]
td = TwiceDifferentiable(alpha -> mlogit(alpha, X, Z, y), startvals; autodiff = :forward)

# Use Optim #
α_hat_multinomial = optimize(td, startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
α_hat_multinomial_ad = α_hat_multinomial.minimizer

# Find Hessian #
H  = Optim.hessian!(td, α_hat_multinomial_ad) # Finding the Hessian Matrix
α_hat_multinomial_ad_se = sqrt.(diag(inv(H))) # Standard Errors of alpha_hat

println([α_hat_multinomial_ad α_hat_multinomial_ad_se])

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 2: Interpretation of gamma_hat
#:::::::::::::::::::::::::::::::::::::::::::::::::::
#  Since it is a linear-log model, The coefficient gamma/100 represents the change in utility with a 1% increase in wage.

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3: Nested Logit Model
#:::::::::::::::::::::::::::::::::::::::::::::::::::
function nested_logit(alpha, X, Z, y, nesting_structure)

    beta = alpha[1:end-3]
    lambda = alpha[end-2:end-1]
    gamma = alpha[end]
    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==j
    end
    bigBeta = [repeat(beta[1:K],1,length(nesting_structure[1])) repeat(beta[K+1:2K],1,length(nesting_structure[2])) zeros(K)]
    # Betas here vary with the nest structure

    T = promote_type(eltype(X),eltype(alpha))
    num   = zeros(T,N,J)
    lidx  = zeros(T,N,J) # Linear Index: Numerator is divided by lambda in choice probabilities
    # lidx is a part of the numerator without the exponent
    dem   = zeros(T,N)

    for j=1:J
        if j in nesting_structure[1]
            lidx[:,j] = exp.( (X*bigBeta[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)./lambda[1] )
        elseif j in nesting_structure[2]
            lidx[:,j] = exp.( (X*bigBeta[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)./lambda[2] )
        else
            lidx[:,j] = exp.(zeros(N))
        end
    end

    for j=1:J
        if j in nesting_structure[1]
            num[:,j] = lidx[:,j].*sum(lidx[:,nesting_structure[1][:]];dims=2).^(lambda[1]-1)
        elseif j in nesting_structure[2]
            num[:,j] = lidx[:,j].*sum(lidx[:,nesting_structure[2][:]];dims=2).^(lambda[2]-1)
            # only sum up the colums of the nesting structure
        else
            num[:,j] = lidx[:,j]
        end

        dem .+= num[:,j]
    end

    P = num./repeat(dem,1,J)

    loglike = -sum( bigY.*log.(P) )

    return loglike
end

nesting_structure = [[1 2 3], [4 5 6 7]]
startvals = [2*rand(2*size(X,2)).-1; 1; 1; .1]

td2 = TwiceDifferentiable(alpha -> nested_logit(alpha, X, Z, y, nesting_structure), startvals; autodiff = :forward)

# Using Optim
α_hat_nlogit = optimize(td2, startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
α_hat_nlogit_ad = α_hat_nlogit.minimizer

# Estimating Hessian
H  = Optim.hessian!(td2, α_hat_nlogit_ad)
α_hat_nlogit_ad_se = sqrt.(diag(inv(H)))

println([α_hat_nlogit_ad α_hat_nlogit_ad_se])

return nothing
end

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 4: All Wrapping Function
#:::::::::::::::::::::::::::::::::::::::::::::::::::
problemset_3()


text:
# SAM File Format
# ===============

module SAM

import Automa
import Automa.RegExp: @re_str
import Automa.Stream: @mark, @markpos, @relpos, @abspos
import BioAlignments
import BioGenerics: BioGenerics, isfilled, header
import BioGenerics.Exceptions: missingerror
import BioGenerics.Automa: State
import BioSequences
import TranscodingStreams: TranscodingStreams, TranscodingStream

using Printf: @sprintf

#TODO: update import BioCore.RecordHelper: unsafe_parse_decimal
# r"[0-9]+" must match `data[range]`.
function unsafe_parse_decimal(::Type{T}, data::Vector{UInt8}, range::UnitRange{Int}) where {T<:Unsigned}
    x = zero(T)
    @inbounds for i in range
        x = Base.Checked.checked_mul(x, 10 % T)
        x = Base.Checked.checked_add(x, (data[i] - UInt8('0')) % T)
    end
    return x
end

# r"[-+]?[0-9]+" must match `data[range]`.
function unsafe_parse_decimal(::Type{T}, data::Vector{UInt8}, range::UnitRange{Int}) where {T<:Signed}
    lo = first(range)
    if data[lo] == UInt8('-')
        sign = T(-1)
        lo += 1
    elseif data[lo] == UInt8('+')
        sign = T(+1)
        lo += 1
    else
        sign = T(+1)
    end
    x = zero(T)
    @inbounds for i in lo:last(range)
        x = Base.Checked.checked_mul(x, 10 % T)
        x = Base.Checked.checked_add(x, (data[i] - UInt8('0')) % T)
    end
    return sign * x
end

# #TODO: update BioCore.Ragel.State (will likely change with TrnscodingStreams).
# import BufferedStreams: BufferedStreams, BufferedInputStream
# # A type keeping track of a ragel-based parser's state.
# mutable struct State{T<:BufferedInputStream}
#     stream::T      # input stream
#     cs::Int        # current DFA state of Ragel
#     linenum::Int   # line number: parser is responsible for updating this
#     finished::Bool # true if finished (regardless of where in the stream we are)
# end

# function State(initstate::Int, input::BufferedInputStream)
#     return State(input, initstate, 1, false)
# end

include("flags.jl")
include("metainfo.jl")
include("record.jl")
include("header.jl")
include("reader.jl")
include("readrecord.jl")
include("writer.jl")

end


text:
export multiple_shooting_objective

function generate_loss_func(loss,t,i)
  new_loss = nothing
  if typeof(loss)<:L2Loss && i==1
    new_loss = L2Loss(t,loss.data[:,i:i+length(t)])
  elseif typeof(loss)<:L2Loss
    i = (i-1)*length(t)
    if i+length(t) < size(loss.data)[2]
      new_loss = L2Loss(t,loss.data[:,i:i+length(t)])
    else
      new_loss = L2Loss(t,loss.data[:,i:size(loss.data)[2]])
    end
  end
  new_loss
end

struct Merged_Solution{T1,T2,T3}
  u::T1
  t::T2
  sol::T3
end

function multiple_shooting_objective(prob::DiffEqBase.DEProblem,alg,loss,
                              regularization=nothing;prior=nothing,
                              mpg_autodiff = false,discontinuity_weight=1.0,
                              verbose_opt = false,
                              prob_generator = STANDARD_PROB_GENERATOR,
                              autodiff_prototype = mpg_autodiff ? zeros(init_N_params) : nothing,
                              autodiff_chunk = mpg_autodiff ? ForwardDiff.Chunk(autodiff_prototype) : nothing,
                              kwargs...)
  cost_function = function (p)
    N = length(p)-length(prob.p)
    time_len = Int(floor((length(loss.t)*length(prob.u0))/N))
    time_dur = loss.t[1:time_len+1]
    sol = []
    loss_val = 0
    j = 1
    for i in 1:length(prob.u0):N
        tmp_prob = remake(prob;u0=p[i:i+length(prob.u0)-1],p=p[N+1:N+length(prob.p)],tspan=(time_dur[1],time_dur[end]))
        if typeof(loss) <: Union{L2Loss,LogLikeLoss}
          push!(sol,solve(tmp_prob,alg;saveat=time_dur,save_everystep=false,dense=false,kwargs...))
          if (j+1)*time_len < length(loss.t)
            time_dur = loss.t[j*time_len+1:(j+1)*time_len+1]
          else
            time_dur = loss.t[j*time_len+1:Int(length(loss.t))]
          end
        else
          push!(sol,solve(tmp_prob,alg;kwargs...))
        end
        j = j+1
    end
    if any((s.retcode != :Success for s in sol))
        return Inf
    end
    u = [i for j in 1:length(sol) for i in sol[j].u[1:end-1]]
    t = [i for j in 1:length(sol) for i in sol[j].t[1:end-1]]
    push!(u,sol[end].u[end])
    push!(t,sol[end].t[end])
    sol_loss = Merged_Solution(u,t,sol)
    sol_new = DiffEqBase.build_solution(prob,alg,loss.t,sol_loss.u,
                                        retcode = :Success)
    loss_val = loss(sol_new)
    if prior != nothing
      loss_val += prior_loss(prior,p[end-length(prior):end])
    end
    if regularization != nothing
      loss_val += regularization(p)
    end
    for i in 2:length(sol)
      loss_val += discontinuity_weight*sum(sol[i][1] - sol[i-1][end])^2
    end
    if verbose_opt
      count::Int += 1
      if mod(count,verbose_steps) == 0
        println("Iteration: $count")
        println("Current Cost: $loss_val")
        println("Parameters: $p")
      end
    end
    loss_val
  end
    if mpg_autodiff
      gcfg = ForwardDiff.GradientConfig(cost_function, autodiff_prototype, autodiff_chunk)
      g! = (x, out) -> ForwardDiff.gradient!(out, cost_function, x, gcfg)
    else
      g! = (x, out) -> Calculus.finite_difference!(cost_function,x,out,:central)
    end
    cost_function2 = function (p,grad)
      if length(grad)>0
          g!(p,grad)
      end
      cost_function(p)
    end
  DiffEqObjective(cost_function,cost_function2)
end


text:
"Create a function call from a list of arguments"
function rlang_p(f, args...; kwargs...)
    argn = length(args)+length(kwargs)
    s = l = protect(allocArray(LangSxp,argn+1))
    try
        setcar!(s, sexp(f))
        for argv in args
            s = cdr(s)
            setcar!(s, sexp(argv))
        end
        for (key,argv) in kwargs
            s = cdr(s)
            settag!(s, sexp(key))
            setcar!(s, sexp(argv))
        end
    finally
        unprotect(1)
    end
    l
end

"Create a function call from a function pointer and a list of arguments and return it as an RObject, which can then be evaulated"
rlang(f, args...; kwargs...) = RObject(rlang_p(f,args...; kwargs...))

"""
Evaluate a function in the global environment. The first argument corresponds
to the function to be called. It can be either a FunctionSxp type, a SymSxp or
a Symbol."""
rcall_p(f,args...;kwargs...) = reval_p(rlang_p(f,args...;kwargs...))

"""
Evaluate a function in the global environment. The first argument corresponds
to the function to be called. It can be either a FunctionSxp type, a SymSxp or
a Symbol.
"""
rcall(f,args...;kwargs...) = RObject(rcall_p(f,args...;kwargs...))

(f::RObject{S})(args...;kwargs...) where S<:FunctionSxp = rcall(f,args...;kwargs...)


text:
## Day 20: Trench Map ##
########################

day20_input = readlines(open("input-20"))
day20_test = readlines(open("test-20"))

function parse_rules(input)
    rules = zeros(Bool, 512)
    for i = 1:512
        rules[i] = input[1][i] == '#'
    end
    rules
end

@assert parse_rules(day20_test)[1] == false
@assert parse_rules(day20_input)[1] == true

function parse_image(input)
    w = length(input[3])
    h = length(input) - 2
    image = zeros(Bool, (w, h))
    for x = 1:w, y = 1:h
        image[x, y] = input[y+2][x] == '#'
    end
    image
end

# An image is actually a Matrix of booleans + a boolean at infinity.
# This is because the rules can force us to update all the infinite values,
# if rules[1] == true.

function update_image(image, inf_pixel, rules)
    # Turns out we need to extend by 2 pixles on each side to also get the
    w, h = size(image) .+ (4, 4)
    output = zeros(Bool, (w, h))
    for x = 1:w, y = 1:h
        lookup = get_lookup_code(image, inf_pixel, x - 2, y - 2)
        output[x, y] = rules[lookup+1]
    end
    (output, inf_pixel ? output[512] : output[1])
end

function get_lookup_code(image, inf_pixel, cx, cy)
    w, h = size(image)
    result = 0
    for dx = -1:1, dy = -1:1
        x = cx + dx
        y = cy + dy
        if x <= 0 || x > w || y <= 0 || y > h
            result = 2 * result + inf_pixel
        else
            result = 2 * result + image[x, y]
        end
    end
    result
end

function enhance(input, n)
    rules = parse_rules(input)
    image = transpose(parse_image(input))
    inf_pixel = false
    for _ = 1:n
        image, inf_pixel = update_image(image, inf_pixel, rules)
    end
    image
end

@assert 35 == @show sum(enhance(day20_test, 2))
@assert 5819 == @show sum(enhance(day20_input, 2))

@assert 3351 == @show sum(enhance(day20_test, 50))
@assert 18516 == @show sum(enhance(day20_input, 50))


text:
using IncrementalInference
using Test

@testset "Test delete tree clique" begin

fg = generateCanonicalFG_lineStep(3;
                                   poseEvery=1,
                                   landmarkEvery=3,
                                   posePriorsAt=[],
                                   landmarkPriorsAt=[0],
                                   sightDistance=2,
                                   solverParams=SolverParams(algorithms=[:default, :parametric]))

# getSolverParams(fg).graphinit = false
# getSolverParams(fg).treeinit = true
# getSolverParams(fg).dbg = true
# getSolverParams(fg).useMsgLikelihoods = true
# getSolverParams(fg).dbg = true
# getSolverParams(fg).drawtree = true
# getSolverParams(fg).showtree = true

smtasks = Task[]
oldtree = solveTree!(fg; smtasks=smtasks, verbose=false, recordcliqs=ls(fg));

@test IIF.isRoot(oldtree, IIF.CliqueId(1))
@test IIF.isRoot(oldtree, IIF.getClique(oldtree,1))
@test !IIF.isRoot(oldtree, IIF.CliqueId(2))
@test !IIF.isRoot(oldtree, IIF.CliqueId(3))

IIF.deleteClique!(oldtree, IIF.CliqueId(1))

@test IIF.isRoot(oldtree, IIF.CliqueId(2))
@test IIF.isRoot(oldtree, IIF.CliqueId(3))

@test IIF.getClique(oldtree, :x0) == IIF.getClique(oldtree, IIF.CliqueId(3)) == IIF.getClique(oldtree, 1)
@test IIF.getClique(oldtree, :x3) == IIF.getClique(oldtree, IIF.CliqueId(2)) == IIF.getClique(oldtree, 2)

# drawTree(oldtree, show=true)

tree = solveTree!(fg, oldtree; smtasks=smtasks, verbose=false, recordcliqs=ls(fg));

# csmAnimate(tree, hists, frames=1)

end

@testset "test tree show and list functions" begin

##

fg = generateCanonicalFG_Kaess(graphinit=false)

eo = [:l2; :l1; :x1; :x2; :x3]

tree = buildTreeReset!(fg, eo)

show(tree)

show(tree[:x3])
show(tree[:x1])
show(tree[:l2])

rtl = ls(tree)

# test the return is same structure and relatively sorted
@test 3 == length(rtl)
@test rtl[1][1] == 1
@test intersect(rtl[1][2], [:x3; :x2]) |> length == 2
@test rtl[2][1] == 2
@test intersect(rtl[2][2], [:x1; :l1]) |> length == 2
@test rtl[3][1] == 3
@test intersect(rtl[3][2], [:l2;]) |> length == 1

# test individual tree list

# the root clique with two children
rtl1 = ls(tree, 1)
rtl1 = ls(tree, :x3)
@test rtl1.parent |> length == 0
@test rtl1.children |> length == 2
@test intersect( (x->x[1]).(rtl1.children), [2,3]) |> length == 2
@test intersect(rtl1.children[1][2], [:x1; :l1]) |> length == 2
@test intersect(rtl1.children[2][2], [:l2]) |> length == 1

# first leaf clique with a parent
rtl2 = ls(tree, 2)
rtl2 = ls(tree, :x1)
rtl2 = ls(tree, :l1)

@test rtl2.parent |> length == 1
@test rtl2.children |> length == 0
@test intersect( (x->x[1]).(rtl2.parent), [1]) |> length == 1
@test length(rtl2.children) == 0

# second leaf clique with a parent
rtl3 = ls(tree, 3)
rtl3 = ls(tree, :l2)

@test rtl3.parent |> length == 1
@test rtl3.children |> length == 0
@test intersect( (x->x[1]).(rtl3.parent), [1]) |> length == 1
@test length(rtl3.children) == 0

##

end

##


text:
"""
    NumPyArrays extends PyCall to provide for additional conversion of Julia
    arrays into NumPy arrays without copying.

```jldoctest
julia> using NumPyArrays, PyCall

julia> rA = reinterpret(UInt8, zeros(Int8, 4,4));

julia> pytypeof(PyObject(rA))
PyObject <class 'list'>

julia> pytypeof(NumPyArray(rA))
PyObject <class 'numpy.ndarray'>

julia> pytypeof(PyObject(NumPyArray(rA)))
PyObject <class 'numpy.ndarray'>

julia> sA = @view collect(1:16)[5:9];

julia> pytypeof(PyObject(sA))
PyObject <class 'list'>

julia> pytypeof(NumPyArray(sA))
PyObject <class 'numpy.ndarray'>
```
"""
module NumPyArrays

# Portions of this code are derived directly from PyCall.jl

@static if isdefined(Base, :Experimental) &&
           isdefined(Base.Experimental, Symbol("@optlevel"))
    Base.Experimental.@optlevel 1
end

export NumPyArray, pytypeof

# Imports for _NumPyArray
import PyCall: NpyArray, PYARR_TYPES, @npyinitialize, npy_api, npy_type
import PyCall: @pycheck, NPY_ARRAY_ALIGNED, NPY_ARRAY_WRITEABLE, pyembed
import PyCall: PyObject, PyPtr

# General imports
import PyCall: pyimport, pytype_query, pytypeof, PyArray

"""
    KnownImmutableArraysWithParent{T} where T <: PyCall.PYARR_TYPES

Immutable `AbstractArray`s in `Base` that have a non-immutable parent that can be embedded in the PyCall GC
"""
const KnownImmutableArraysWithParent{T} = Union{SubArray{T}, Base.ReinterpretArray{T}, Base.ReshapedArray{T}, Base.PermutedDimsArray{T}} where T

"""
    KnownStridedArrays{T} where T <: PyCall.PYARR_TYPES

`AbstractArray`s in `Base` where the method `strides` is applicable
"""
const KnownStridedArrays{T} = StridedArray{T} where T

"""
    NumPyArray{T,N}(po::PyObject)

NumPyArray is a wrapper around a PyCall.PyArray. It is an AbstractArray.
The main purpose of a NumPyArray is so to provide a constructor to generalize
the conversion of Julia arrays into NumPyArrays. `T` is the element type of the array.
`N` is the number of dimensions.

For other uses, such as wrapping an existing array from NumPy, use `PyCall.PyArray`.

Use `PyObject` and `PyArray` methods to convert `NumPyArray` into those types.
"""
mutable struct NumPyArray{T,N} <: AbstractArray{T,N}
    pa::PyArray{T,N}
end
NumPyArray{T,N}(po::PyObject) where {T,N} = NumPyArray{T,N}(PyArray(po))
NumPyArray(po::PyObject) = NumPyArray(PyArray(po))

"""
    NumPyArray(a::AbstractArray, [revdims::Bool])

Convert an AbstractArray where `isapplicable(strides, a)` is `true` to a NumPy array.
The AbstractArray must either be mutable or have a mutable parent. Optionally,
transpose the dimensions of the array if `revdims` is `true`.
"""
NumPyArray(a::AbstractArray{T}) where T <: PYARR_TYPES = NumPyArray(a, false)
function NumPyArray(a::KnownStridedArrays{T}, revdims::Bool) where T <: PYARR_TYPES
    _NumPyArray(a, revdims)
end
function NumPyArray(a::AbstractArray{T}, revdims::Bool) where T <: PYARR_TYPES
    # For a general AbstractArray, we do not know if strides applies
    if applicable(strides, a)
        _NumPyArray(a, revdims)
    else
        error("Only AbstractArrays where strides is applicable can be converted to NumPyArrays.")
    end
end
NumPyArray(o::PyPtr) = NumPyArray(PyObject(o))

# Modified PyCall.NpyArray to accept AbstractArray{T}, assumes strides is applicable
function _NumPyArray(a::AbstractArray{T}, revdims::Bool) where T <: PYARR_TYPES
    @npyinitialize
    size_a = revdims ? reverse(size(a)) : size(a)
    strides_a = revdims ? reverse(strides(a)) : strides(a)
    p = @pycheck ccall(npy_api[:PyArray_New], PyPtr,
        (PyPtr,Cint,Ptr{Int},Cint, Ptr{Int},Ptr{T}, Cint,Cint,PyPtr),
        npy_api[:PyArray_Type],
        ndims(a), Int[size_a...], npy_type(T),
        Int[strides_a...] * sizeof(eltype(a)), a, sizeof(eltype(a)),
        NPY_ARRAY_ALIGNED | NPY_ARRAY_WRITEABLE,
        C_NULL)
    return NumPyArray{T,ndims(a)}(p, a)
end

# Make a NumPyArray that embeds a reference to keep, to prevent Julia
# from garbage-collecting keep until o is finalized.
# See also PyObject(o::PyPtr, keep::Any) from which this is derived
NumPyArray{T,N}(o::PyPtr, keep::Any) where {T,N} = numpyembed(NumPyArray{T,N}(PyObject(o)), keep)

# PyCall already has convert(::Type{PyObject}, o) = PyObject(o)
#Base.convert(::Type{PyObject}, a::NumPyArray) = a.po
PyObject(a::NumPyArray) = PyObject(PyArray(a))
Base.convert(::Type{PyArray}, a::NumPyArray) = PyArray(a)
PyArray(a::NumPyArray) = a.pa
Base.convert(::Type{Array}, a::NumPyArray{T}) where T = convert(Array{T}, PyArray(a))
Base.convert(T::Type{<:Array}, a::NumPyArray) = convert(T, PyArray(a))

# See PyCall.pyembed(po::PyObject, jo::Any)
function numpyembed(a::NumPyArray{T,N}, jo::Any) where {T,N}
    if isimmutable(jo)
        if applicable(parent, jo)
            return NumPyArray{T,N}(pyembed(PyObject(a), parent(jo)))
        else
            throw(ArgumentError("numpyembed: immutable argument without a parent is not allowed"))
        end
    else
        return NumPyArray{T,N}(pyembed(PyObject(a), jo))
    end
end
numpyembed(a::NumPyArray, jo::KnownImmutableArraysWithParent) = numpyembed(a, jo.parent)

# AbstractArray interface, provided as a convenience. Conversion to PyArray is recommended
Base.size(a::NumPyArray) = size(PyArray(a))
Base.length(a::NumPyArray) = length(PyArray(a))
Base.getindex(a::NumPyArray, args...) = getindex(PyArray(a), args...)
Base.setindex!(a::NumPyArray, args...) = setindex!(PyArray(a), args...)
Base.axes(a::NumPyArray) = axes(PyArray(a))
# Strides should be added to PyArray
Base.strides(a::NumPyArray{T}) where T = PyArray(a).st
Base.stride(a::NumPyArray{T}) where T = stride(PyArray(a))

Base.pointer(a::NumPyArray, args...) = pointer(PyArray(a), args...)
Base.unsafe_convert(t::Type{Ptr{T}}, a::NumPyArray{T}) where T = Base.unsafe_convert(t, PyArray(a))

Base.similar(a::NumPyArray, ::Type{T}, dims::Dims) where {T} = similar(PyArray(a), T, dims)

# Aliasing some PyCall functions. Conversion to PyObject or PyArray is recommended
pytypeof(a::NumPyArray) = pytypeof(PyObject(PyArray(a)))

end # module NumPyArrays