text:
int_rules_4_1_1_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.1.2*(g*cos)^p*(a+b*sin)^m =#
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~(p') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => (1 / ((~b) ^ ~p * ~f)) * Subst(Antiderivative((~a + ~x) ^ (~m + (~p - 1) / 2) * (~a - ~x) ^ ((~p - 1) / 2), ~x), ~x, ~b * sin(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (IntegerQ((~p - 1) / 2) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (GeQ(~p, -1) || Not(IntegerQ(~m + 1 / 2)))))
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~(p') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => (1 / ((~b) ^ ~p * ~f)) * Subst(Antiderivative((~a + ~x) ^ ~m * ((~b) ^ 2 - (~x) ^ 2) ^ ((~p - 1) / 2), ~x), ~x, ~b * sin(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (IntegerQ((~p - 1) / 2) && NeQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)), ~x) => (-(~b) * (~g * cos(~e + ~f * ~x)) ^ (~p + 1)) / (~f * ~g * (~p + 1)) + ~a * Antiderivative((~g * cos(~e + ~f * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~p], ~x) && (IntegerQ(2 * ~p) || NeQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (~a / ~g) ^ (2 * ~m) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ (2 * ~m + ~p) / (~a - ~b * sin(~e + ~f * ~x)) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegerQ(~m) && (LtQ(~p, -1) && GeQ(2 * ~m + ~p, 0))))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (~b * (~g * cos(~e + ~f * ~x)) ^ (~p + 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m) / (~a * ~f * ~g * ~m) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (EqQ(Simplify(~m + ~p + 1), 0) && Not(ILtQ(~p, 0))))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (~b * (~g * cos(~e + ~f * ~x)) ^ (~p + 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m) / (~a * ~f * ~g * Simplify(2 * ~m + ~p + 1)) + (Simplify(~m + ~p + 1) / (~a * Simplify(2 * ~m + ~p + 1))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ ~p * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (ILtQ(Simplify(~m + ~p + 1), 0) && (NeQ(2 * ~m + ~p + 1, 0) && Not(IGtQ(~m, 0)))))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (~b * (~g * cos(~e + ~f * ~x)) ^ (~p + 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1)) / (~f * ~g * (~m - 1)) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (EqQ((2 * ~m + ~p) - 1, 0) && NeQ(~m, 1)))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (-(~b) * (~g * cos(~e + ~f * ~x)) ^ (~p + 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1)) / (~f * ~g * (~m + ~p)) + ((~a * ((2 * ~m + ~p) - 1)) / (~m + ~p)) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ ~p * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IGtQ(Simplify(((2 * ~m + ~p) - 1) / 2), 0) && NeQ(~m + ~p, 0)))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (-(~b) * (~g * cos(~e + ~f * ~x)) ^ (~p + 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m) / (~a * ~f * ~g * (~p + 1)) + ((~a * (~m + ~p + 1)) / ((~g) ^ 2 * (~p + 1))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ (~p + 2) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (GtQ(~m, 0) && (LeQ(~p, -2 * ~m) && IntegersQ(~m + 1 / 2, 2 * ~p))))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (-2 * ~b * (~g * cos(~e + ~f * ~x)) ^ (~p + 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1)) / (~f * ~g * (~p + 1)) + (((~b) ^ 2 * ((2 * ~m + ~p) - 1)) / ((~g) ^ 2 * (~p + 1))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ (~p + 2) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (GtQ(~m, 1) && (LtQ(~p, -1) && IntegersQ(2 * ~m, 2 * ~p))))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)) / sqrt(~(g') * cos(~(e') + ~(f') * ~x)), ~x) => ((~a * sqrt(1 + cos(~e + ~f * ~x)) * sqrt(~a + ~b * sin(~e + ~f * ~x))) / (~a + ~a * cos(~e + ~f * ~x) + ~b * sin(~e + ~f * ~x))) * Antiderivative(sqrt(1 + cos(~e + ~f * ~x)) / sqrt(~g * cos(~e + ~f * ~x)), ~x) + ((~b * sqrt(1 + cos(~e + ~f * ~x)) * sqrt(~a + ~b * sin(~e + ~f * ~x))) / (~a + ~a * cos(~e + ~f * ~x) + ~b * sin(~e + ~f * ~x))) * Antiderivative(sin(~e + ~f * ~x) / (sqrt(~g * cos(~e + ~f * ~x)) * sqrt(1 + cos(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && EqQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (-(~b) * (~g * cos(~e + ~f * ~x)) ^ (~p + 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1)) / (~f * ~g * (~m + ~p)) + ((~a * ((2 * ~m + ~p) - 1)) / (~m + ~p)) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ ~p * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (GtQ(~m, 0) && (NeQ(~m + ~p, 0) && IntegersQ(2 * ~m, 2 * ~p))))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (~g * (~g * cos(~e + ~f * ~x)) ^ (~p - 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * (~m + ~p)) + (((~g) ^ 2 * (~p - 1)) / (~a * (~m + ~p))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ (~p - 2) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (LtQ(~m, -1) && (GtQ(~p, 1) && ((GtQ(~m, -2) || (EqQ(2 * ~m + ~p + 1, 0) || EqQ(~m, -2) && IntegerQ(~p))) && (NeQ(~m + ~p, 0) && IntegersQ(2 * ~m, 2 * ~p))))))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (2 * ~g * (~g * cos(~e + ~f * ~x)) ^ (~p - 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * (2 * ~m + ~p + 1)) + (((~g) ^ 2 * (~p - 1)) / ((~b) ^ 2 * (2 * ~m + ~p + 1))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ (~p - 2) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (LeQ(~m, -2) && (GtQ(~p, 1) && (NeQ(2 * ~m + ~p + 1, 0) && (Not(ILtQ(~m + ~p + 1, 0)) && IntegersQ(2 * ~m, 2 * ~p))))))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (~b * (~g * cos(~e + ~f * ~x)) ^ (~p + 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m) / (~a * ~f * ~g * (2 * ~m + ~p + 1)) + ((~m + ~p + 1) / (~a * (2 * ~m + ~p + 1))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ ~p * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (LtQ(~m, -1) && (NeQ(2 * ~m + ~p + 1, 0) && IntegersQ(2 * ~m, 2 * ~p))))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p / (~a + ~(b') * sin(~(e') + ~(f') * ~x)), ~x) => (~g * (~g * cos(~e + ~f * ~x)) ^ (~p - 1)) / (~b * ~f * (~p - 1)) + ((~g) ^ 2 / ~a) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (GtQ(~p, 1) && IntegerQ(2 * ~p)))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p / (~a + ~(b') * sin(~(e') + ~(f') * ~x)), ~x) => (~b * (~g * cos(~e + ~f * ~x)) ^ (~p + 1)) / (~a * ~f * ~g * (~p - 1) * (~a + ~b * sin(~e + ~f * ~x))) + (~p / (~a * (~p - 1))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~p], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (Not(GeQ(~p, 1)) && IntegerQ(2 * ~p)))
    @apply_utils Antiderivative(sqrt(~(g') * cos(~(e') + ~(f') * ~x)) / sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)), ~x) => ((~g * sqrt(1 + cos(~e + ~f * ~x)) * sqrt(~a + ~b * sin(~e + ~f * ~x))) / (~a + ~a * cos(~e + ~f * ~x) + ~b * sin(~e + ~f * ~x))) * Antiderivative(sqrt(1 + cos(~e + ~f * ~x)) / sqrt(~g * cos(~e + ~f * ~x)), ~x) - ((~g * sqrt(1 + cos(~e + ~f * ~x)) * sqrt(~a + ~b * sin(~e + ~f * ~x))) / (~b + ~b * cos(~e + ~f * ~x) + ~a * sin(~e + ~f * ~x))) * Antiderivative(sin(~e + ~f * ~x) / (sqrt(~g * cos(~e + ~f * ~x)) * sqrt(1 + cos(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && EqQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ (3 / 2) / sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)), ~x) => (~g * sqrt(~g * cos(~e + ~f * ~x)) * sqrt(~a + ~b * sin(~e + ~f * ~x))) / (~b * ~f) + ((~g) ^ 2 / (2 * ~a)) * Antiderivative(sqrt(~a + ~b * sin(~e + ~f * ~x)) / sqrt(~g * cos(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && EqQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p / sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)), ~x) => (-2 * ~b * (~g * cos(~e + ~f * ~x)) ^ (~p + 1)) / (~f * ~g * (2 * ~p - 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (3 / 2)) + ((2 * ~a * (~p - 2)) / (2 * ~p - 1)) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ ~p / (~a + ~b * sin(~e + ~f * ~x)) ^ (3 / 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (GtQ(~p, 2) && IntegerQ(2 * ~p)))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p / sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x)), ~x) => (-(~b) * (~g * cos(~e + ~f * ~x)) ^ (~p + 1)) / (~a * ~f * ~g * (~p + 1) * sqrt(~a + ~b * sin(~e + ~f * ~x))) + ((~a * (2 * ~p + 1)) / (2 * (~g) ^ 2 * (~p + 1))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ (~p + 2) / (~a + ~b * sin(~e + ~f * ~x)) ^ (3 / 2), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (LtQ(~p, -1) && IntegerQ(2 * ~p)))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => (((~a) ^ ~m * (~g * cos(~e + ~f * ~x)) ^ (~p + 1)) / (~f * ~g * (1 + sin(~e + ~f * ~x)) ^ ((~p + 1) / 2) * (1 - sin(~e + ~f * ~x)) ^ ((~p + 1) / 2))) * Subst(Antiderivative((1 + (~b / ~a) * ~x) ^ (~m + (~p - 1) / 2) * (1 - (~b / ~a) * ~x) ^ ((~p - 1) / 2), ~x), ~x, sin(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~p], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && IntegerQ(~m))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => (((~a) ^ 2 * (~g * cos(~e + ~f * ~x)) ^ (~p + 1)) / (~f * ~g * (~a + ~b * sin(~e + ~f * ~x)) ^ ((~p + 1) / 2) * (~a - ~b * sin(~e + ~f * ~x)) ^ ((~p + 1) / 2))) * Subst(Antiderivative((~a + ~b * ~x) ^ (~m + (~p - 1) / 2) * (~a - ~b * ~x) ^ ((~p - 1) / 2), ~x), ~x, sin(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && Not(IntegerQ(~m)))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (-((~g * cos(~e + ~f * ~x)) ^ (~p + 1)) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * sin(~e + ~f * ~x)) / (~f * ~g * (~p + 1)) + (1 / ((~g) ^ 2 * (~p + 1))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ (~p + 2) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1) * (~a * (~p + 2) + ~b * (~m + ~p + 2) * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (LtQ(0, ~m, 1) && (LtQ(~p, -1) && (IntegersQ(2 * ~m, 2 * ~p) || IntegerQ(~m)))))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (-((~g * cos(~e + ~f * ~x)) ^ (~p + 1)) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1) * (~b + ~a * sin(~e + ~f * ~x))) / (~f * ~g * (~p + 1)) + (1 / ((~g) ^ 2 * (~p + 1))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ (~p + 2) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 2) * ((~b) ^ 2 * (~m - 1) + (~a) ^ 2 * (~p + 2) + ~a * ~b * (~m + ~p + 1) * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (GtQ(~m, 1) && (LtQ(~p, -1) && (IntegersQ(2 * ~m, 2 * ~p) || IntegerQ(~m)))))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (-(~b) * (~g * cos(~e + ~f * ~x)) ^ (~p + 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 1)) / (~f * ~g * (~m + ~p)) + (1 / (~m + ~p)) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ ~p * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - 2) * ((~b) ^ 2 * (~m - 1) + (~a) ^ 2 * (~m + ~p) + ~a * ~b * ((2 * ~m + ~p) - 1) * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~p], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (GtQ(~m, 1) && (NeQ(~m + ~p, 0) && (IntegersQ(2 * ~m, 2 * ~p) || IntegerQ(~m)))))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (~g * (~g * cos(~e + ~f * ~x)) ^ (~p - 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * (~m + 1)) + (((~g) ^ 2 * (~p - 1)) / (~b * (~m + 1))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ (~p - 2) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * sin(~e + ~f * ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (LtQ(~m, -1) && (GtQ(~p, 1) && IntegersQ(2 * ~m, 2 * ~p))))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (-(~b) * (~g * cos(~e + ~f * ~x)) ^ (~p + 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~f * ~g * ((~a) ^ 2 - (~b) ^ 2) * (~m + 1)) + (1 / (((~a) ^ 2 - (~b) ^ 2) * (~m + 1))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ ~p * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * (~a * (~m + 1) - ~b * (~m + ~p + 2) * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~p], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (LtQ(~m, -1) && IntegersQ(2 * ~m, 2 * ~p)))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (~g * (~g * cos(~e + ~f * ~x)) ^ (~p - 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * (~m + ~p)) + (((~g) ^ 2 * (~p - 1)) / (~b * (~m + ~p))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ (~p - 2) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~b + ~a * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (GtQ(~p, 1) && (NeQ(~m + ~p, 0) && IntegersQ(2 * ~m, 2 * ~p))))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => ((~g * cos(~e + ~f * ~x)) ^ (~p + 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * (~b - ~a * sin(~e + ~f * ~x))) / (~f * ~g * ((~a) ^ 2 - (~b) ^ 2) * (~p + 1)) + (1 / ((~g) ^ 2 * ((~a) ^ 2 - (~b) ^ 2) * (~p + 1))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ (~p + 2) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (((~a) ^ 2 * (~p + 2) - (~b) ^ 2 * (~m + ~p + 2)) + ~a * ~b * (~m + ~p + 3) * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (LtQ(~p, -1) && IntegersQ(2 * ~m, 2 * ~p)))
    @apply_utils Antiderivative(1 / (sqrt(~(g') * cos(~(e') + ~(f') * ~x)) * sqrt(~a + ~(b') * sin(~(e') + ~(f') * ~x))), ~x) => ((2 * sqrt(2) * sqrt(~g * cos(~e + ~f * ~x)) * sqrt((~a + ~b * sin(~e + ~f * ~x)) / ((~a - ~b) * (1 - sin(~e + ~f * ~x))))) / (~f * ~g * sqrt(~a + ~b * sin(~e + ~f * ~x)) * sqrt((1 + cos(~e + ~f * ~x) + sin(~e + ~f * ~x)) / ((1 + cos(~e + ~f * ~x)) - sin(~e + ~f * ~x))))) * Subst(Antiderivative(1 / sqrt(1 + ((~a + ~b) * (~x) ^ 4) / (~a - ~b)), ~x), ~x, sqrt((1 + cos(~e + ~f * ~x) + sin(~e + ~f * ~x)) / ((1 + cos(~e + ~f * ~x)) - sin(~e + ~f * ~x)))) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && NeQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => ((~g * (~g * cos(~e + ~f * ~x)) ^ (~p - 1) * (1 - sin(~e + ~f * ~x)) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * ((-((~a - ~b)) * (1 - sin(~e + ~f * ~x))) / ((~a + ~b) * (1 + sin(~e + ~f * ~x)))) ^ (~m / 2)) / (~f * (~a + ~b) * (~m + 1))) * HypergeometricFunctions._₂F₁(~m + 1, ~m / 2 + 1, ~m + 2, (2 * (~a + ~b * sin(~e + ~f * ~x))) / ((~a + ~b) * (1 + sin(~e + ~f * ~x)))) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && EqQ(~m + ~p + 1, 0))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => ((~g * cos(~e + ~f * ~x)) ^ (~p + 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~f * ~g * (~a - ~b) * (~p + 1)) + (~a / ((~g) ^ 2 * (~a - ~b))) * Antiderivative(((~g * cos(~e + ~f * ~x)) ^ (~p + 2) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m) / (1 - sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && EqQ(~m + ~p + 2, 0))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (((~g * cos(~e + ~f * ~x)) ^ (~p + 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~f * ~g * (~a - ~b) * (~p + 1)) - ((~b * (~m + ~p + 2)) / ((~g) ^ 2 * (~a - ~b) * (~p + 1))) * Antiderivative((~g * cos(~e + ~f * ~x)) ^ (~p + 2) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m, ~x)) + (~a / ((~g) ^ 2 * (~a - ~b))) * Antiderivative(((~g * cos(~e + ~f * ~x)) ^ (~p + 2) * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m) / (1 - sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && ILtQ(~m + ~p + 2, 0))
    @apply_utils Antiderivative(sqrt(~(g') * cos(~(e') + ~(f') * ~x)) / (~a + ~(b') * sin(~(e') + ~(f') * ~x)), ~x) => With([q = Rt(-((~a) ^ 2) + (~b) ^ 2, 2)], (((~a * ~g) / (2 * ~b)) * Antiderivative(1 / (sqrt(~g * cos(~e + ~f * ~x)) * (q + ~b * cos(~e + ~f * ~x))), ~x) - ((~a * ~g) / (2 * ~b)) * Antiderivative(1 / (sqrt(~g * cos(~e + ~f * ~x)) * (q - ~b * cos(~e + ~f * ~x))), ~x)) + ((~b * ~g) / ~f) * Subst(Antiderivative(sqrt(~x) / ((~g) ^ 2 * ((~a) ^ 2 - (~b) ^ 2) + (~b) ^ 2 * (~x) ^ 2), ~x), ~x, ~g * cos(~e + ~f * ~x))) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && NeQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative(1 / (sqrt(~(g') * cos(~(e') + ~(f') * ~x)) * (~a + ~(b') * sin(~(e') + ~(f') * ~x))), ~x) => With([q = Rt(-((~a) ^ 2) + (~b) ^ 2, 2)], ((-(~a) / (2q)) * Antiderivative(1 / (sqrt(~g * cos(~e + ~f * ~x)) * (q + ~b * cos(~e + ~f * ~x))), ~x) - (~a / (2q)) * Antiderivative(1 / (sqrt(~g * cos(~e + ~f * ~x)) * (q - ~b * cos(~e + ~f * ~x))), ~x)) + ((~b * ~g) / ~f) * Subst(Antiderivative(1 / (sqrt(~x) * ((~g) ^ 2 * ((~a) ^ 2 - (~b) ^ 2) + (~b) ^ 2 * (~x) ^ 2)), ~x), ~x, ~g * cos(~e + ~f * ~x))) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && NeQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => ((~g * (~g * cos(~e + ~f * ~x)) ^ (~p - 1) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * (~m + ~p) * ((-(~b) * (1 - sin(~e + ~f * ~x))) / (~a + ~b * sin(~e + ~f * ~x))) ^ ((~p - 1) / 2) * ((~b * (1 + sin(~e + ~f * ~x))) / (~a + ~b * sin(~e + ~f * ~x))) ^ ((~p - 1) / 2))) * AppellF1(-(~p) - ~m, (1 - ~p) / 2, (1 - ~p) / 2, (1 - ~p) - ~m, (~a + ~b) / (~a + ~b * sin(~e + ~f * ~x)), (~a - ~b) / (~a + ~b * sin(~e + ~f * ~x))) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~p], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (ILtQ(~m, 0) && Not(IGtQ(~m + ~p + 1, 0))))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => ((~g * (~g * cos(~e + ~f * ~x)) ^ (~p - 1)) / (~f * (1 - (~a + ~b * sin(~e + ~f * ~x)) / (~a - ~b)) ^ ((~p - 1) / 2) * (1 - (~a + ~b * sin(~e + ~f * ~x)) / (~a + ~b)) ^ ((~p - 1) / 2))) * Subst(Antiderivative((-(~b) / (~a - ~b) - (~b * ~x) / (~a - ~b)) ^ ((~p - 1) / 2) * (~b / (~a + ~b) - (~b * ~x) / (~a + ~b)) ^ ((~p - 1) / 2) * (~a + ~b * ~x) ^ ~m, ~x), ~x, sin(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && Not(IGtQ(~m, 0)))
    @apply_utils Antiderivative((~(g') * sec(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => (~g) ^ (2 * IntPart(~p)) * (~g * cos(~e + ~f * ~x)) ^ FracPart(~p) * (~g * sec(~e + ~f * ~x)) ^ FracPart(~p) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m / (~g * cos(~e + ~f * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && Not(IntegerQ(~p))
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.1.2 (g cos)^p (a+b sin)^m.jl


text:
@deprecate relocatable!(ex::Expr) RelocatableExpr(ex)
@deprecate relocatable!(ex::RelocatableExpr) ex
@deprecate get_method whichtt
Base.@deprecate_binding FileModules ModuleExprsSigs
@deprecate firstlineno(ex) firstline(ex).line

ext: jl


text:
abstract type AbstractTrainableGaussianModel <: AbstractTrainablePotential end

### struct for defining a Gaussian force model
mutable struct TrainableGaussianModel <: AbstractTrainableGaussianModel
  dim::Int64
  nbasis::Int64
  shape::Tuple{Vararg{Int64}}
  mus::Vector{Vector{Float64}}
  sigs::Vector{Vector{Float64}}
  theta::Matrix{Float64}
  condition::Function
  f::Vector{Float64}
  em::Vector{Float64}
  function TrainableGaussianModel(shape::Tuple{Vararg{Int64}},mus::Vector{Vector{Float64}},sigs::Vector{Vector{Float64}},
                                  theta::Matrix{Float64},condition::Function)
    dim = length(shape)
    nbasis = reduce(*,shape,init=1)
    new(dim,nbasis,shape,mus,sigs,theta,condition,zeros(Float64,dim),zeros(Float64,nbasis))
  end
end

"""
Constructor with zeros for parameters and default true condition
"""
function TrainableGaussianModel(shape::Tuple{Vararg{Int64}},mus::Vector{Vector{Float64}},sigs::Vector{Vector{Float64}})
  dim::Int64 = length(shape)
  nbasis::Int64 = reduce(*,shape,init=1)
  theta::Matrix{Float64} = zeros(Float64, (nbasis,dim))
  condition::Fucntion = defaultcondition
  TrainableGaussianModel(shape,mus,sigs,theta,condition)
end

"""
Constructor with zeros for parameters
"""
function TrainableGaussianModel(shape::Tuple{Vararg{Int64}},mus::Vector{Vector{Float64}},sigs::Vector{Vector{Float64}},condition::Function)
  dim::Int64 = length(shape)
  nbasis::Int64 = reduce(*,shape,init=1)
  theta::Matrix{Float64} = zeros(Float64, (nbasis,dim))
  TrainableGaussianModel(shape,mus,sigs,theta,condition)
end

"""
Constructor with default true condition
"""
function TrainableGaussianModel(shape::Tuple{Vararg{Int64}},mus::Vector{Vector{Float64}},sigs::Vector{Vector{Float64}},theta::Matrix{Float64})
  dim::Int64 = length(shape)
  nbasis::Int64 = reduce(*,shape,init=1)
  condition::Function = defaultcondition
  TrainableGaussianModel(shape,mus,sigs,theta,condition)
end

"""
struct for defining a time-dependent Gaussian force model
make sure the last set of basis functions is the time-dependence
"""
mutable struct TrainableGaussianModelTD <: AbstractTrainableGaussianModel
  dim::Int64
  nbasis::Int64
  shape::Tuple{Vararg{Int64}}
  mus::Vector{Vector{Float64}}
  sigs::Vector{Vector{Float64}}
  theta::Matrix{Float64}
  condition::Function
  f::Vector{Float64}
  em::Vector{Float64}
  function TrainableGaussianModelTD(shape::Tuple{Vararg{Int64}},mus::Vector{Vector{Float64}},sigs::Vector{Vector{Float64}},
                         theta::Matrix{Float64},condition::Function)
    dim = length(shape)-1
    nbasis = reduce(*,shape,init=1)
    new(dim,nbasis,shape,mus,sigs,theta,condition,zeros(Float64,dim),zeros(Float64,nbasis))
  end
end

"""
Constructor with zeros for parameters and default true condition
"""
function TrainableGaussianModelTD(shape::Tuple{Vararg{Int64}},mus::Vector{Vector{Float64}},sigs::Vector{Vector{Float64}})
  dim::Int64 = length(shape)-1
  nbasis::Int64 = reduce(*,shape,init=1)
  theta::Matrix{Float64} = zeros(Float64, (nbasis,dim))
  condition::Function = defaultcondition
  TrainableGaussianModelTD(shape,mus,sigs,theta,condition)
end

"""
Constructor with zeros for parameters
"""
function TrainableGaussianModelTD(shape::Tuple{Vararg{Int64}},mus::Vector{Vector{Float64}},sigs::Vector{Vector{Float64}},condition::Function)
  dim::Int64 = length(shape)-1
  nbasis::Int64 = reduce(*,shape,init=1)
  theta::Matrix{Float64} = zeros(Float64, (nbasis,dim))
  TrainableGaussianModelTD(shape,mus,sigs,theta,condition)
end

"""
Constructor with default true condition
"""
function TrainableGaussianModelTD(shape::Tuple{Vararg{Int64}},mus::Vector{Vector{Float64}},sigs::Vector{Vector{Float64}},theta::Matrix{Float64})
  dim::Int64 = length(shape)-1
  nbasis::Int64 = reduce(*,shape,init=1)
  condition::Function = defaultcondition
  TrainableGaussianModelTD(shape,mus,sigs,theta,condition)
end

"""
Helper functions that get useful information about the model
"""
function getdimensionality(gm::AbstractTrainableGaussianModel)
  return gm.dim
end

# shape of the jacobian with respect to the parameters
function jacshape(gm::AbstractTrainableGaussianModel)
  return (gm.dim*gm.nbasis,gm.dim)
end

# TODO do i need this
function gradshape(gm::AbstractTrainableGaussianModel)
  return gm.dim*gm.nbasis
end

"""
Default condition on whether to use model in equation of motion
"""
function defaultcondition(sys::AbstractSystem)
  return true
end

"""
Functions that compute the value of the basis functions
"""
function computebasis!(sys::AbstractSystem, gm::TrainableGaussianModel)
  k::Vector{Float64} = exp.( -0.5 .* (sys.x[1] .- gm.mus[1]).^2 ./ gm.sigs[1] )
  for i in 2:gm.dim
    k = vec( k * (exp.( -0.5 .* (sys.x[i] .- gm.mus[i]).^2 ./ gm.sigs[i] ))' )
  end
  gm.em .= k
end

function computebasis!(sys::AbstractSystem, gm::TrainableGaussianModelTD)
  # forward
  k::Vector{Float64} = exp.( -0.5 .* (sys.x[1] .- gm.mus[1]).^2 ./ gm.sigs[1] )
  for i in 2:gm.dim
    k = vec( k * (exp.( -0.5 .* (sys.x[i] .- gm.mus[i]).^2 ./ gm.sigs[i] ))' )
  end
  k = vec( k * (exp.( -0.5 .* (sys.t .- gm.mus[end]).^2 ./ gm.sigs[end] ))' )
  gm.em .= k
end

"""
Force function with return value
"""
function force(sys::AbstractSystem, gm::AbstractTrainableGaussianModel; compbasis::Bool=true)
  if compbasis
    computebasis!(sys,gm)
  end
  fout::Vector{Float64} = zeros(Float64,gm.dim)
  for i in 1:gm.dim
    fout[i] = sum( @views gm.theta[:,i] .* gm.em )
  end
  return fout
end

"""
Force function with in-place changing
"""
function force!(sys::AbstractSystem, gm::AbstractTrainableGaussianModel; compbasis::Bool=true)
  if compbasis
    computebasis!(sys,gm)
  end
  for i in 1:gm.dim
    gm.f[i] = sum( @views gm.theta[:,i] .* gm.em )
  end
end

"""
Gradient wrt parameters for gaussian model
"""
function jacobian!(jac::Array{Float64},sys::AbstractSystem,gm::AbstractTrainableGaussianModel;
                   compbasis::Bool=false)
  if compbasis
    computebasis!(sys,gm)
  end
  ind0::Int64 = 1
  indf::Int64 = gm.nbasis
  for i in 1:gm.dim
    @. @views jac[ind0:indf,i] .= gm.em
    ind0 += gm.nbasis
    indf += gm.nbasis
  end
end

"""
Gradient wrt parameters for gaussian model
"""
function gradient!(grad::Array{Float64},sys::AbstractSystem,gm::AbstractTrainableGaussianModel;
                   compbasis::Bool=false)
  if compbasis
    computebasis!(sys,gm)
  end
  for i in 1:gm.dim
    grad[:,i] .= gm.em
  end
end

"""
Update parameters for the neural network model
"""
function updateparams!(lr::Float64,dtheta::Vector{Float64},gm::AbstractTrainableGaussianModel)
    # update coefficients for forces
    gm.theta .+= (lr .* reshape(dtheta, (gm.nbasis,gm.dim)))
end

"""
Set parameters for the neural network model
"""
function setparams!(theta::Vector{Float64},gm::AbstractTrainableGaussianModel)
    # update coefficients for forces
    gm.theta .= theta
end

function setparams!(gm::AbstractTrainableGaussianModel)
end


text:
#
# Draft ACEModel Interface:
#
# * A model incorporates both the parameterisation and the parameters.
# * E.g. a LinearACEModel would know about the basis and the coefficients
# * We can then perform the following operations:
#    - evaluate the model at a given configuration, with given parameters
#    - take gradient w.r.t. the configuration
#    - take gradient w.r.t. the parameters
#
# An advantage will be that we can do the parameter reduction trick and then
# use a fast evaluator to obtain the gradient w.r.t. configuration which is
# very expensive if we naively take derivatives w.r.t. the basis and then
# apply parameters. (this is not implemented below and will be the next step)
# Another advantage is that model, and full model construction are all stored
# together for future inspection.

"""
`struct LinearACEModel`: linear model for symmetric properties in terms of
a `SymmetricBasis`.

The typical way to construct a linear model is to first construct a basis
`basis`, some default coefficients `c` and then call
```julia
model = LinearACEModel(basis, c)
```

### Multiple properties

If `c::Vector{<: Number}` then the output of the model will be the property
encoded in the basis. But one can also use a single basis to produce
multiple properties (with different coefficients). This can be achieved by
simply supplying `c::Vector{SVector{N, T}}` where `N` will then be the
number of properties.
"""
struct LinearACEModel{TB, TP, TEV} <: AbstractACEModel
   basis::TB
   c::Vector{TP}
   evaluator::TEV
end

struct NaiveEvaluator end

function LinearACEModel(basis::SymmetricBasis, c = zeros(length(basis));
               evaluator = :standard)
   if evaluator == :naive
      ev = NaiveEvaluator()
   elseif evaluator == :standard
      ev = ProductEvaluator(basis, c)
   elseif evaluator == :recursive
      error("Recursive evaluator not yet implemented")
   else
      error("unknown evaluator")
   end
   return LinearACEModel(basis, c, ev)
end

# LinearACEModel(basis::SymmetricBasis, c::Vector, evaluator) =
#          LinearACEModel(basis, c, ev, VectorPool{eltype(c)})

# ------- parameter wrangling

nparams(m::LinearACEModel) = length(m.c)

params(m::LinearACEModel) = copy(m.c)

function set_params!(m::LinearACEModel, c)
   m.c[:] .= c
   set_params!(m.evaluator, m.basis, c)
   return m
end

set_params!(::NaiveEvaluator, args...) = nothing

# ------------------- FIO

==(V1::LinearACEModel, V2::LinearACEModel) =
      _allfieldsequal(V1, V2)

write_dict(V::LinearACEModel) =
      Dict( "__id__" => "ACE_LinearACEModel",
             "basis" => write_dict(V.basis),
                 "c" => write_dict(V.c),
         "evaluator" => write_dict(V.evaluator) )

function read_dict(::Val{:ACE_LinearACEModel}, D::Dict)
   basis = read_dict(D["basis"])
   c = read_dict(D["c"])
   # special evaluator version of the read_dict
   evaluator = read_dict(Val(Symbol(D["evaluator"]["__id__"])),
                         D["evaluator"], basis, c)
   return LinearACEModel(basis, c, evaluator)
end

write_dict(ev::NaiveEvaluator) =
      Dict("__id__" => "ACE_NaiveEvaluator" )

read_dict(::Val{:ACE_NaiveEvaluator}, D::Dict, args...) =
      NaiveEvaluator()

# ------------------- dispatching on the evaluators

evaluate(m::LinearACEModel, X::AbstractConfiguration) =
      evaluate(m, m.evaluator, X)

grad_config(m::LinearACEModel, X::AbstractConfiguration) =
      grad_config(m, m.evaluator, X)

# there is a canonical implementation for this, so no need to dispatch
function grad_params(m::LinearACEModel, cfg::AbstractConfiguration)
   _gi(Bi, c::Number) = Bi
   _gi(Bi, c::SVector{N, <: Number}) where {N} =
         SMatrix{N,N}( Diagonal([ Bi for _=1:N ]) )

   B = evaluate(m.basis, cfg)
   g = [ _gi(B[i], m.c[i]) for i = 1:length(B) ]
   release!(B)
   return g
end

# ∂_params ∂_config V
# currently doesn't work with multiple properties
grad_params_config(m::LinearACEModel, cfg::AbstractConfiguration) =
      evaluate_d(m.basis, cfg)

      # TODO: fix terminology, bring in linear with the _rrule_.... thing
adjoint_EVAL_D(m::LinearACEModel, cfg::AbstractConfiguration, w) =
      adjoint_EVAL_D(m, m.evaluator, cfg, w)

# ------------------- implementation of naive evaluator
#  this is only intended for testing, as it uses the naive evaluation of
#  the symmetric basis, rather than the conversion to the AA basis

function evaluate(m::LinearACEModel, ::NaiveEvaluator, cfg::AbstractConfiguration)
   B = evaluate(m.basis, cfg)
   val = contract(m.c, B)
   release!(B)
   return val
end

function grad_config(m::LinearACEModel, ::NaiveEvaluator,
                    cfg::AbstractConfiguration)
   dB = evaluate_d(m.basis, cfg)
   TG = promote_type(eltype(m.c), eltype(dB))
   g = zeros(TG, length(cfg))
   for ix = 1:length(cfg), ib = 1:length(m.basis)
      g[ix] += m.c[ib] * dB[ib, ix]
   end
   release!(dB)
   return g
end

function adjoint_EVAL_D(m::LinearACEModel, ::NaiveEvaluator,
                        X::AbstractConfiguration, w)
   dB = grad_params_config(m, X)
   g = zeros(size(dB, 1))
   for i = 1:length(g), j = 1:size(dB, 2)
      g[i] += dot(dB[i, j], w[j])
   end
   release_dB!(m.basis, dB)
   return g
end

# ------------------- an rrule for evaluating a linear model

import ChainRules
import ChainRules: rrule, @thunk, NoTangent, @not_implemented

function _adj_evaluate(dp, model::ACE.LinearACEModel, cfg)
   gp_ = ACE.grad_params(model, cfg)
   gp = [ a * dp for a in gp_ ]
   return NoTangent(), gp, _rrule_evaluate(dp, model, cfg)
end

# this is monkey-patching the rotten rrule inside of ACE
# ... and should replace that rrule. ALso introduce thunks to prevent
#     evaluating more than we need.
function ChainRules.rrule(::typeof(evaluate), model::ACE.LinearACEModel, cfg::AbstractConfiguration)
   return evaluate(model, cfg),
          dp -> _adj_evaluate(dp, model, cfg)
end

# rrule for the rrule ... this enables mixed second derivatives of the form
#   D^2 * / D p Dcfg.
# the code double-checks that indeed only those derivatives are needed!
function ChainRules.rrule(::typeof(_adj_evaluate), dp, model::ACE.LinearACEModel, cfg)
   # adj = (_, g_params, g_cfg)
   #   D(dq[1] * _ + dq[2] * g_params + dq[3] * g_cfg) / D(dp, model, cfg)
   #       0 = ^^^    ^^^ = 0
   #   D( dq[3] * g_cfg ) / D( dq, model, cfg )
   # but for simplicity ignore Dcfg for now (not yet implemented)
   # recall also that g_cfg = D (dp * eval(model, cfg)) / D cfg
   # dp should be a vector of the same length as the number of properties

   function _second_adj(dq_)
      # adj = (NoTangent(), gp1, g_cfg)
      # here we assume that only g_cfg was used, which means that
      # dq_[3] = force-like vector and dq_[2] == NoTangent()
      @assert dq_[1] == dq_[2] == ZeroTangent()
      @assert dq_[3] isa AbstractVector{<: ACE.DState}
      @assert length(dq_[3]) == length(cfg)
      dq = dq_[3]  # Vector of DStates

      # adj_n = ∑_j dq_j ⋅ ∂B_k / ∂r_j * θ_nk
      # dp ⋅ adj = ∑_n ∑_j dq_j ⋅ ∂B_k / ∂r_j * θ_nk * dp_n
      # grad[k] = ∑_j dq_j ⋅ ∂B_k / ∂r_j
      grad = ACE.adjoint_EVAL_D(model, model.evaluator, cfg, dq)

      # gradient w.r.t parameters:
      grad_params = [ gg * dp for gg in grad ]

      # gradient w.r.t. dp    # TODO: remove the |> Vector?
      grad_dp = sum( model.c[k] * grad[k] for k = 1:length(grad) )  |> Vector

      return NoTangent(), grad_dp, grad_params, NoTangent()
   end

   return _adj_evaluate(dp, model, cfg), _second_adj
end


text:
"""
Collision detection
"""
module Collision

export triangle_triangle, triangle_line, triangle_plane, plane_segment, triangle_edges, γpoint

using StaticArrays, SparseArrays, LinearAlgebra, Base.Iterators, IterTools

partial = (f::Function,y...)->(z...)->f(y...,z...)
second(x) = x[2]

# create a basis for Clifford Algebra using Gamma Matrices
# we will do some projective geometry with it.
# ∧ product works pretty well, but regressive product doesn't.  That's ok we just need join
const 𝑖 = 1.0im
γ0 = [1.0 0 0 0; 0 1 0 0; 0 0 -1 0; 0 0 0 -1]
γ1 = [0.0 0 0 1; 0 0 1 0; 0 -1 0 0; -1 0 0 0]
γ2 = [0.0 0 0 -𝑖; 0 0 𝑖 0; 0 𝑖 0 0; -𝑖 0 0 0]
γ3 = [0.0 0 1 0; 0 0 0 -1; -1 0 0 0; 0 1 0 0]

# (-+++) clifford algebra
# should do a performance comparison with StaticArrays, SparseArrays and Arrays
#const g1,g2,g3,g4 = sparse.(im .* (γ0, γ1, γ2, γ3))
#const g1,g2,g3,g4 = (im .* (γ0, γ1, γ2, γ3))
const g1,g2,g3,g4 = SMatrix{4,4,Complex{Float64}}.(im .* (γ0, γ1, γ2, γ3))

const 𝐼 = g1*g2*g3*g4
const ProjectivePoint = typeof(g1)
const ProjectiveLine = ProjectivePoint
const ProjectivePlane = ProjectivePoint

#∧(a::T, b::T) where T<:SparseMatrixCSC = a*b
∧(a::T, b::T) where T<:SArray = a*b
const join = ∧

# orientation of a pseudoscalar element.
import Base.sign
sign(ps) = sign(real(tr(ps*𝐼)))
volume(ps) = real(tr(ps*𝐼))*0.25

"""
    γpoint(x,y,z)

create a homogeneous point for collision detection using our specialized basis.
simply pass in the euclidean x,y,z coordinates
"""
γpoint(x,y,z) = x*g2 + y*g3 + z*g4 + 1.0g1
γpoint(xyz) = γpoint(xyz...)

""" true if triangle and line intersect """
function triangle_line((a,b,c)::Tuple{P,P,P}, (d,e)::Tuple{P,P}) where P <: ProjectivePoint
  de = d∧e
  crossings = map(sign, (de∧a∧b, de∧b∧c, de∧c∧a))

  # line passes through inside of triangle if line-crossing volumes are all same orientation
  # if product is zero then we are coplanar or exactly on an edge
  #!me can give false positive for miss where endpoint of segment is coplanar
  sum(crossings)^2 == 9.0 || prod(crossings) == 0.0
end

""" true if triangle and plane intersect """
triangle_plane((a,b,c)::Tuple{P,P,P}, p::V) where {V<:ProjectivePlane, P<:ProjectivePoint} =
  !(sign(p∧a) == sign(p∧b) == sign(p∧c))

""" true if plane and segment intersect """
plane_segment(abc::V, (d,e)::Tuple{P,P}) where {V<:ProjectivePlane, P<:ProjectivePoint} =
  sign(abc∧d) != sign(abc∧e)

"""
    triangle_triangle((a,b,c), (d,e,f))

return true if the triangle defined by points a,b,c is intersecting the triangle defined by points d,e,f
operates on projective points created by γpoint
"""
function triangle_triangle((a,b,c)::T, (d,e,f)::T) where {P<:ProjectivePoint, T<:Tuple{P,P,P}}
  abc = a∧b∧c
  def = d∧e∧f
  abc_edges = ((a,b), (b,c), (c,a))
  def_edges = ((d,e), (e,f), (f,d))

  check_plane(plane, edges) = ( plane_segment(plane, edge) for edge in edges )
  check_plane2(plane, edges) = imap(second, Iterators.filter(first, zip(check_plane(plane, edges), edges)))

  check_tri(tri, plane, edges) = ( triangle_line(tri, edge) for edge in check_plane2(plane, edges) )

  for ihit in check_tri((a,b,c), abc, def_edges)
    if ihit
      return true
    end
  end
  for ihit in check_tri((d,e,f), def, abc_edges)
    if ihit
      return true
    end
  end

  return false
end

function triangle_triangle((a,b,c)::T, (d,e,f)::T) where {P<:SVector, T<:Tuple{P,P,P}}
  triangle_triangle((γpoint(a), γpoint(b), γpoint(c)), (γpoint(d), γpoint(e), γpoint(f)))
end

triangle_triangle(a,b,c,d,e,f) = triangle_triangle((a,b,c), (d,e,f))

"""
    triangle_edges(a,b,c, d,e,f)

return true if the triangle defined by points a,b,c is intersecting the triangle defined by points d,e,f
the hit record will include all edges intersecting the interior of a triangle

the format of the hit record is:  ({true|false}, ({triangleID},({vertexID, vertexID})))
triangleID = 1 or 2 to indicate triangle abc or def and vertexID is in range 1:3 indexing segment in [a,b,c] or [d,e,f].

operates on projective points created by γpoint
"""
function triangle_edges(a::P, b::P, c::P, d::P, e::P, f::P) where {P<:ProjectivePoint}
  planes = (a∧b∧c, d∧e∧f)
  tris = ((a,b,c), (d,e,f))
  triedges = (((a,b), (b,c), (c,a)), ((d,e), (e,f), (f,d)))
  hitrecord = (((2,(1,2)), (2,(2,3)), (2,(3,1))), ((1,(1,2)), (1,(2,3)), (1,(3,1))))

  hits = Tuple{Int64,Tuple{Int64,Int64}}[]
  for itri in 1:2
    plane = planes[itri]
    tri = tris[itri]
    otheredges = triedges[(itri&1)+1]
    for iedge in 1:3
      edge = otheredges[iedge]
      if plane_segment(plane, edge) == 1
        if triangle_line(tri, edge)
          push!(hits, hitrecord[itri][iedge])
        end
      end
    end
  end

  (length(hits) > 0, hits)

  #==
  abc = a∧b∧c
  def = d∧e∧f
  abc_edges = ((a,b), (b,c), (c,a))
  iabc_edges = ((1,(1,2)), (1,(2,3)), (1,(3,1)))
  def_edges = ((d,e), (e,f), (f,d))
  idef_edges = ((2,(1,2)), (2,(2,3)), (2,(3,1)))

  check_plane(plane, edges) = ( plane_segment(plane, edge) for edge in edges )
  check_plane2(plane, edges) = imap(second, Iterators.filter(first, zip(check_plane(plane, edges), edges)))

  check_tri(tri, plane, edges) = ( triangle_line(tri, edge) for edge in check_plane2(plane, edges) )

  hits = vcat(Iterators.filter(first, zip(check_tri((a,b,c), abc, def_edges), idef_edges))|>collect,
              Iterators.filter(first, zip(check_tri((d,e,f), def, abc_edges), iabc_edges))|>collect )
  (length(hits) > 0, map(second,hits))
  ==#
end

function triangle_edges(a::P, b::P, c::P, d::P, e::P, f::P) where {P<:SVector}
  triangle_edges(γpoint(a), γpoint(b), γpoint(c), γpoint(d), γpoint(e), γpoint(f))
end

#==
function tri_edge(a,b,c, v,w)
  n = normalize((b-a)×(c-a))
end
==#

function meet_tri_edge(a, b, c, v, w)
  n = normalize((b-a)×(c-a))
  l = w-v
  l₀ = w
  d = ((a-l₀)⋅n)/(l⋅n)
  l₀ + l*d
end

function γmeet_tri_edge(a, b, c, v, w)
  ga,gb,gc,gv,gw = γpoint.((a,b,c,v,w))

  r = real(tr(ga*gb*gc*gw*𝐼))/real(tr(ga*gb*gc*gv*𝐼))
  p = (w - r*v)/(1.0-r)
end

"""
    segment_between(((a, b, c), (d, e, f)), [(itri, (v1, v2))])

return the longest segment traced on the face of a triangle, between the edges.
the indices are in the range [1,2,3] and refer to the triangle points
itri is in the range [1,2] and refers to either triangle abc or def
itri indicates the triangle the edge v1 and v2 refer to
"""
function segment_between(abc_def::T, edges) where {R<:Real,
                                                   P<:AbstractVector{R},
                                                   T<:Tuple{Tuple{P,P,P},
                                                            Tuple{P,P,P}}}
  @assert length(edges) > 1

  hitp = map(edges) do (s, (v1, v2))
    meet_tri_edge(abc_def[(s&1)+1]..., abc_def[s][[v1,v2]]...)
  end

  if length(hitp) == 2
    (hitp[2], hitp[1])
  else
    seg = (zero(typeof(hitp[1])), zero(typeof(hitp[1])))
    best = 0.0
    for i in 1:length(hitp)
      for j in (i+1):length(hitp)
        cand_seg = (hitp[j], hitp[i])
        cand = LinearAlgebra.norm_sqr(cand_seg[2]-cand_seg[1])
        if cand > best
          best = cand
          seg = cand_seg
        end
      end
    end
    seg
  end
end

function segment_between(abc_def::T, edges) where {R<:Real,
                                                   P<:AbstractVector{R},
                                                   T<:Tuple{AbstractVector{P},
                                                            AbstractVector{P}}}
  segment_between((Tuple(abc_def[1]), Tuple(abc_def[2])), edges)
end

end


text:
"""
    Graph{T}

A type representing an undirected graph
"""

"""
    Node{T<:Int}(index::Int, supply::T, demand::T)

Defines the `Node` struct
"""
mutable struct Node{T<:Int}
    index::Int
    supply::T
    demand::T
end

"""
    Edge{T<:Int}(index::Int, cost::T, capacity::T)

Defines the `Edge` struct
"""
mutable struct Edge{T<:Int}
    index::Int
    cost::T
    capacity::T
end

"""
    Graph{T<:Int}(nodes::Vector{Node{T}}, edges::Vector{Edge{T}}, first_out::Vector{Int}, end_out::Vector{Int})

Defines the `Graph` struct
"""
struct Graph{T<:Int}
    nodes::Vector{Node{T}}
    edges::Vector{Edge{T}}
    first_out::Vector{Int}
    end_node::Vector{Int}
end

"""
    adj(graph::Graph, node::Node)

Returns a `Vector` of pairs of (`Node`, `Edge`)
"""
function adj(graph::Graph, node::Node)::Vector{Tuple{Node, Edge}}
    node_ind = node.index
    edge_start_ind = graph.first_out[node_ind]
    edge_end_ind = graph.first_out[node_ind + 1]
    end_nodes = graph.end_node[edge_start_ind:edge_end_ind-1]
    nodes = graph.nodes[end_nodes]
    edges = graph.edges[edge_start_ind:edge_end_ind-1]
    return collect(zip(nodes, edges))
end

"""
    size(graph::Graph, dim::Integer)::Integer

Returns the size of the graph along the 'dims'
Where, 'dims' can be one of 1 == nodes and 2 == edges
"""
function Base.size(graph::Graph, dim::Integer)::Integer
    if dim == 1
        return size(graph.nodes, 1)
    elseif dim == 2
        return size(graph.edges, 1)
    else
        throw(ArgumentError("dims can only be either 1 (nodes) or 2 (edges)"))
    end
end

"""
    path(pred::Vector{Int64}, target::Int64, source::Int64)::Vector{Int64}

Returns the path from source to target given the `Vector` of predecessors
"""
function path(pred::Vector{Int64}, target::Int64, source::Int64)::Vector{Int64}
    path_arr = Vector{Int64}([target])
    x = target
    while x != source
        y = pred[x]
        push!(path_arr, y)
        x = y
    end
    return reverse(path_arr)
end


text:
using Revise
using GeodesicRendering
using AccretionGeometry
using CarterBoyerLindquist
using AccretionFormulae
using StaticArrays
using Plots

include("temperature_render.jl")
gr()

function flux(m, gp, max_time)
    g = AccretionFormulae.redshift(m, gp, max_time)
    r = gp.u[2]
    M = m.M
    a_star = m.a
    mdot = AccretionFormulae.mdot(M)
    flux = g^4 * AccretionFormulae.diss(mdot, r, a_star, M)
    return (flux, g)
end

"""
Calculating the index for each flux
"""
function en(pixel, e_max, n_bins)
    if typeof(pixel) == Tuple{Float64, Float64}
        flux = pixel[1]
        g = pixel[2]
        en = 6.4*g
        index = n_bins*(en/e_max)

        # rounding down to the nearest bin
        index = Int32(round(index-0.5))
        return (flux, index)
    else
        return NaN
    end
end

function iron_line_profile(;
                            mass = 1,
                            spin = 0.998,
                            obs_angle = 85.0,
                            tolerance = 1e-8,
                            dtmax = 1000.0,
                            size_multiplier::Int64 = 1,
                            resolution = 400,
                            fov = 3.0,
                            output = "plot",
                            normalised = true,
                            nbins = 100
                            )
    hmap, cache, title, new_img = temperature_render(;
                                            mass = mass,
                                            spin = spin,
                                            obs_angle = obs_angle,
                                            tolerance = tolerance,
                                            dtmax = dtmax,
                                            size_multiplier = size_multiplier,
                                            resolution = resolution,
                                            fov = fov
                                            )

    # define a common filter
    filter = ConstPointFunctions.filter_intersected

    # define the value function
    flux_redshift_vf = (
        PointFunction(flux) ∘ filter
    )

    # generate the images from the cache
    img = GeodesicRendering.apply(flux_redshift_vf, cache)

    # parameters
    e_min = 0
    e_max = 10
    n_bins = nbins

    bins = zeros(n_bins)
    x_vals = LinRange(0, e_max, n_bins)

    new = en.(img, e_max, n_bins)

    for tuple in new
        if typeof(tuple) == Tuple{Float64, Int32}
            bins[tuple[2]] += tuple[1]
        end
    end

    # normalising to peak at 1
    if normalised
        bins ./= maximum(bins)
    end
    plt = plot(x_vals, bins, xlims=(e_min, e_max), grid=false, framestyle=:box)
    xaxis!("Energy (keV)")
    yaxis!("Flux (Arbitrary Units)")

    if output == "plot"
        return plt, hmap, title, cache
    elseif output == "data"
        return x_vals, bins
    end
end

# plt, hmap, title, cache = iron_line_profile(tolerance=1e-12,
#                                             # size_multiplier=3,
#                                             # fov=12,
#                                             # dtmax=5,
#                                             obs_angle=30,
#                                             image_size_multiplier=1.5)
# display(plt)


text:
module JFVM

using PyCall
# I prefer not to use the following command for the issues that it has on windows machines
# pygui_start(:wx)
mayavis=0
try
  @pyimport mayavi.mlab as m
  mayavis=m
catch
  warn("Mayavi is not installed or could not be imported.")
end
using PyPlot

import Base: +, -, *, /, .*, ./
export MeshStructure, BoundaryCondition, CellValue, FaceValue, CellVector,
       arithmeticMean, geometricMean, harmonicMean, upwindMean, linearMean,
       tvdMean, createBC, boundaryConditionTerm, cellBoundary, solvePDE,
       divergenceTerm, gradientTerm, convectionUpwindTerm, createCellVector,
       convectionTerm, convectionTvdTerm, diffusionTerm, createCellVariable,
       createFaceVariable, copyCell, fluxLimiter, createMesh1D,
       createMesh2D, createMesh3D, createMeshRadial2D, createMeshCylindrical2D,
       createMeshCylindrical3D, createMeshCylindrical1D, solveLinearPDE,
       visualizeCells, linearSourceTerm, constantSourceTerm, transientTerm,
       solveMUMPSLinearPDE, faceEval, cellEval, permfieldlogrndg, permfieldlogrnde,
       plot, imshow, xlabel, ylabel, figure, legend, pcolor, contour, colorbar,
       visualizeCellVectors, JFVM_test, solveExplicitPDE, reshapeCell,
       cellVolume, reshapeInternalCell, internalCells, domainInt

include("fvmToolTypes.jl")
include("meshstructure.jl")
include("boundarycondition.jl")
include("domainVariables.jl")
include("diffusionterms.jl")
include("transientTerms.jl")
include("domainOperators.jl")
include("convectionTerms.jl")
include("averagingTerms.jl")
include("calculusTerms.jl")
include("sourceTerms.jl")
include("solveVisualizePDE.jl")
include("JFVMtools.jl")
include("jfvm_test.jl")

end # module


text:
#########
#File: c:\Users\digan\Dropbox\Dynamic_Networks\repos\ScoreDrivenExponentialRandomGraphs\_research\analysis_for_paper_revision\new_application_\load_wiki_talk_data.jl
#Project: c:\Users\digan\Dropbox\Dynamic_Networks\repos\ScoreDrivenExponentialRandomGraphs\_research\analysis_for_paper_revision\new_application_
#Created Date: Wednesday April 21st 2021
#Author: Domenico Di Gangi,  <digangidomenico@gmail.com>
#-----
#Last Modified: Friday May 7th 2021 5:48:20 pm
#Modified By:  Domenico Di Gangi
#-----
#Description:
#-----
########

#%%
using DrWatson
@quickactivate "ScoreDrivenExponentialRandomGraphs"
DrWatson.greet()
using Distributed
begin
using Logging
using ScoreDrivenERGM
import ScoreDrivenERGM:StaticNets,DynNets, Utilities, ErgmRcall
using LinearAlgebra
using CSV
using DataFrames
using Dates
using PyPlot
pygui(true)

nWorkers = 15
addprocs(nWorkers - nprocs())
@sync @everywhere begin
    using Pkg
    Pkg.activate(".")
    Pkg.instantiate()
    using ScoreDrivenERGM
    import ScoreDrivenERGM:StaticNets, DynNets
    import ScoreDrivenERGM.DynNets:SdErgm,SdErgmDirBin0Rec0, simulate_and_estimate_parallel
    using ScoreDrivenERGM.Utilities

    end
end

#%%
#load data
begin
data_path = "./data/wiki_tk/raw_data/"

df = CSV.read(data_path * "wiki_talk_daily_edges.csv", DataFrame)

df.date = Dates.Date.(df.date)

dates = unique(df.date)

# get sequence of edge lists at a given frequency
edge_list_T = [Matrix(gdf[[:source, :target]]) for gdf in groupby(df,[:date])]
edge_list_pres_T = [edge_list_pres(l) for l in edge_list_T]

n_links_T = [size(e)[1] for e in edge_list_T]
n_nodes_T = [ length(unique(reshape(e, :))) for e in edge_list_T]

end

# t_ref = argmax(n_nodes_T)
# n_links_sub = [100, 500, 1000, 2000, 3000, 5000, 10000]
# time = [@time DynNets.stats_from_mat(model1, edge_list_pres(edge_list_pres_T[t_ref][1:n_links, :])) for n_links in n_links_sub]
# plot(n_links_sub, time)
# time = [@time DynNets.stats_from_mat(model_rec_p_star, edge_list_pres(edge_list_pres_T[t_ref][1:n_links, :])) for n_links in n_links_sub]

begin
T = length(n_nodes_T)
c= Dict{String, Any}()
c["model"] =[model_rec_p_star, model_edge_gwd.([0, 0.5, 0.8, 1.2, 1.5, 2])..., model_edge_mutual_gwd.([0, 0.5, 0.8, 1.2, 1.5, 2])..., model_edge_gwesp]
c["Ti"] = [1]
c["Tf"] = [T]
c["f"] = ["d"]
list = sort(dict_list(c), by=x->(string(x["model"])))
d = list[1]
for d in list

    timeSim = @elapsed ch_stats = @sync @distributed vcat for edl in edge_list_pres_T[d["Ti"]:d["Tf"]]
        [DynNets.stats_from_mat(d["model"], edge_list_pres(edl))]
    end

    modelTag =  string(d["model"])
    ergmString = d["model"].staticModel.ergmTermsString

    res1 =  (;modelTag, ergmString, ch_stats) |>DrWatson.ntuple2dict |> DrWatson.tostringdict

    saveDict = merge(res1, d)

    d_save = deepcopy(d)
    d_save["ter"] = replace.( ergmString,  ", fixed = TRUE, cutoff=10" => "")

    saveName = replace.( savename(d_save, "jld2";allowedtypes = (Real, String, Symbol, NamedTuple, Tuple, ScoreDrivenERGM.DynNets.SdErgm) ), r"[\"]" => "")

    timeSave = @elapsed save( datadir("wiki_tk", "ch_stats", saveName), saveDict)

    Logging.@info("Time sim = $timeSim ,  time save = $timeSave ")

end

end

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
@testset "Space" begin

@testset "ContinuousSpace" begin
    s = ContinuousSpace(0., 1.)
    @test in(0.5, s) == true
    @test in(0.0, s) == true
    @test in(1.0, s) == true
    @test in(-1.0, s) == false
    @test in(-Inf, s) == false

    @test in(sample(s), s)
    @test eltype(s) == Float64
end

@testset "MultiContinuousSpace" begin
    s = MultiContinuousSpace([-1., -2.], [1., 2.])
    @test in([0, 0], s) == true
    @test in([0, 3], s) == false
    @test in(sample(s), s)
    @test eltype(s) == Array{Float64, 1}
end

@testset "DiscreteSpace" begin
    s = DiscreteSpace(10)

    @test in(0, s) == false
    @test in(1, s) == true
    @test in(5, s) == true
    @test in(10, s) == true

    @test in(sample(s), s)
    @test size(s) == (10,)
    @test length(s) == 10
end

@testset "MultiDiscreteSpace" begin
    s = MultiDiscreteSpace([2,3,2])

    @test in([0,0,0], s) == false
    @test in([1,1,1], s) == true
    @test in([3,3,3], s) == false

    @test size(s) == (2,3,2)
    @test length(s) == 2*3*2
    @test in(sample(s), s)
end

@testset "Space Tuple" begin
    @test in(([0.5], 5, [1, 1]),
             (MultiContinuousSpace(0,1, (1,)), DiscreteSpace(5), MultiDiscreteSpace([2,2]))) == true
    @test in(([0.5], 1, [1, 1]),
              (MultiContinuousSpace(0,1, (1,)), DiscreteSpace(5), MultiDiscreteSpace([2,2]))) == true
    @test in((),
             (MultiContinuousSpace(0,1, (1,)), DiscreteSpace(5), MultiDiscreteSpace([2,2]))) == false
end

@testset "Space Dict" begin
    @test in(
        Dict(
            "sensors" => Dict(
                "position" => [-10, 0, 10],
                "velocity" => [0.1, 0.2, 0.3],
                "front_cam" => (rand(10, 10, 3), rand(10, 10, 3)),
                "rear_cam" => rand(10,10,3)),
            "ext_controller" => [2, 1, 1],
            "inner_state" => Dict(
                "charge" => 35,
                "job_status" => Dict(
                    "task" => 3,
                    "progress" => [23]))),
        Dict(
            "sensors" =>  Dict(
                "position"=> MultiContinuousSpace(-100, 100, (3,)),
                "velocity"=> MultiContinuousSpace(-1, 1, (3,)),
                "front_cam"=> (MultiContinuousSpace(0, 1, (10, 10, 3)),
                               MultiContinuousSpace(0, 1, (10, 10, 3))),
                "rear_cam" => MultiContinuousSpace(0, 1, (10, 10, 3))),
            "ext_controller" => MultiDiscreteSpace([5, 2, 2]),
            "inner_state" => Dict(
                "charge" => DiscreteSpace(100),
                "job_status" => Dict(
                    "task" => DiscreteSpace(5),
                    "progress" => MultiContinuousSpace(0, 100, (1,)))))) == true
end

end


text:
#Bisection method for finding one root to f in the interval [a,b]
function bisection(f::Function, a::AbstractFloat, b::AbstractFloat;
    tol::AbstractFloat=1e-5, maxiter::Integer=200)
    fa = f(a)
    fa*f(b) <= 0 || error("No real root in [a,b]")
    i = 0
    local c
    while b-a > tol
    i += 1
    i != maxiter || error("Max iteration exceeded")
    c = (a+b)/2
    fc = f(c)
    if fc == 0
    break
    elseif fa*fc > 0
    a = c  # Root is in the right half of [a,b].
    fa = fc
    else
    b = c  # Root is in the left half of [a,b].
    end
    end
    return c
end

#Function for finding the real roots for the real polynomial f(x) = ax³+bx²+cx+d
function solvecubic(a::Real, b::Real, c::Real, d::Real)
    function one_real_root(a::Real, b::Real, g::Float64, h::Float64)
        R = -g/2 + sqrt(h)
        S = cbrt(R)
        T = -g/2 - sqrt(h)
        U = cbrt(T)
        x1 = (S + U) - (b/3a)
        x2 = -(S + U)/2 - (b/3a) - im*(S-U)*sqrt(3) / 2
        x3 = -(S + U)/2 - (b/3a) + im*(S-U)*sqrt(3) / 2
        return Union{Real,Complex}[x1, x2, x3]
    end

    function three_real_root(a::Real, b::Real, g::Float64, h::Float64)
        i = (g^2 / 4 - h)^(1/2)
        j = cbrt(i)
        K = acos(-(g/2i))
        L = -j
        M = cos(K/3)
        N = sqrt(3) * sin(K/3)
        P = -(b/3a)
        x1 = 2j * cos(K/3) - (b/3a)
        x2 = L * (M + N) + P
        x3 = L * (M - N) + P
        return sort([x1, x2, x3])
    end

    if a == zero(typeof(a))
        return Base.error("Can't solve when a=0 in a*x³ + b*x² + c*x + d =0")
    end
    f = (3c/a - (b/a)^2) / 3
    g = (2(b/a)^3 - (9*b*c/a^2) + 27d/a) / 27
    h = g^2 / 4 + f^3 / 27
    if f == zero(typeof(f)) && g == zero(typeof(g)) && h == zero(typeof(h))
        x1 = x2 = x3 = -cbrt(d/a)
        return [x1, x2, x3]
    elseif  h <= zero(typeof(h))
        return three_real_root(a, b, g, h)
    elseif h > zero(typeof(h))
        return one_real_root(a, b, g, h)
    end

    return Base.error("Can't solve when a=$(a) b=$(b) c=$(c) d=$(d)")
end

#Dynamic viscosity of water (Pa s) at T degree celsius
#Equation 5 in Korson et al. (1969)
ηₜ(T::AbstractFloat) = 1.0020*10^((1.1709*(20-T)-0.001827*(T-20)^2)/(T+89.93)-3)

#Saturation Vapor Pressure (Pa) at T degree celsius
#Equation 21 in Alduchov and Eskridge (1996)
SVPₜ(T::AbstractFloat) = 610.94*exp(17.625*T/(243.04+T))

#Vapour-pressure deficit (Pa) at T degree celsius and RH relative humidity (%)
VPDₜ(T::AbstractFloat,RH::AbstractFloat) = SVPₜ(T)*(1-RH/100)

#Calcualte day-light hours
function daylighthour(L::Real,d::Real)
    #Jenkins, 2013, The Sun's position in the sky (eq 17)
    #L = latitude (rad)
    #d = day of the year (d=1 is january 1)
    M = -0.0410+0.017202*d
    Φ = -1.3411+M+0.0334*sin(M)+0.0003*sin(2*M)
    ϵ = sin(23.4*pi/180)
    Φs = sin(Φ)
    return 24*(1-1/pi*acos(tan(L)*ϵ*Φs/sqrt(1-ϵ^2*Φs^2)))
end


text:
export CuDiag

const CuDiag{ElT,VecT}                           = Diag{ElT,VecT} where {VecT<:CuArray{ElT}}
const NonuniformCuDiagTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where {StoreT<:CuDiag}

CuArray(D::NonuniformCuDiagTensor) = CuArray(dense(D))

function NDTensors.dense(T::NonuniformCuDiagTensor{ElT}) where {ElT}
  R_data = CUDA.zeros(ElT, dim(inds(T)))
  diag_inds = diagind(reshape(R_data, dims(inds(T))), 0)
  R_data[diag_inds] .= data(store(T))
  return Tensor(Dense(R_data), inds(T))
end

function Base.promote_rule(::Type{Diag{ElT2,VecT2}}, ::Type{CuDense}) where {ElT2,VecT2<:Number}
    return promote_type(DenseT1, ElT2)
end

function Base.promote_rule(::Type{<:Tensor{ElT1,N1,StoreT1}},
                           ::Type{<:Tensor{ElT2,N2,StoreT2}}) where {ElT1,ElT2,
                                                                     N1,N2,
                                                                     StoreT1<:CuDense,StoreT2<:UniformDiag}
    ElT3 = promote_type(ElT1, ElT2)
    return Tensor{promote_type(ElT1,ElT2),N3,CuDense{ElT3, CuVector{ElT3}}} where {N3}
end
function Base.promote_rule(::Type{<:Tensor{ElT1,N1,StoreT1}}, ::Type{<:Tensor{ElT2,N2,StoreT2}}) where {ElT1,ElT2,
                                                                     N1,N2,
                                                                     StoreT1<:UniformDiag,StoreT2<:CuDense}
    ElT3 = promote_type(ElT1, ElT2)
    return Tensor{promote_type(ElT1,ElT2),N3,CuDense{ElT3, CuVector{ElT3}}} where {N3}
end

function Base.promote_rule(::Type{<:UniformDiag{ElT1}},
                           ::Type{<:NonuniformDiag{ElT2,VecT2}}) where {ElT1,ElT2,VecT2<:CuArray}
    ElT3 = promote_type(ElT1, ElT2)
    return NonuniformDiag{ElT3, CuVector{ElT3}}
end
function Base.promote_rule(::Type{<:NonuniformDiag{ElT2,VecT2}},
                           ::Type{<:UniformDiag{ElT1}}) where {ElT1,ElT2,VecT2<:CuArray}
    ElT3 = promote_type(ElT1, ElT2)
    return NonuniformDiag{ElT3, CuVector{ElT3}}
end
function Base.promote_rule(::Type{DenseT1}, ::Type{Diag{ElT2,VecT2}}) where {DenseT1<:CuDense, ElT2,VecT2<:Number}
  return promote_type(DenseT1, ElT2)
end

function contraction_output_type(TensorT1::Type{<:NonuniformCuDiagTensor},
                                 TensorT2::Type{<:CuDenseTensor},
                                 IndsR)
    return similar_type(promote_type(TensorT1, TensorT2), IndsR)
end
contraction_output_type(TensorT1::Type{<:CuDenseTensor},
                        TensorT2::Type{<:NonuniformCuDiagTensor},
                        IndsR) = contraction_output_type(TensorT2,TensorT1,IndsR)

function contraction_output_type(TensorT1::Type{<:DiagTensor{<:Number,<:CuDiag}},
                                 TensorT2::Type{<:DiagTensor{<:Number,<:CuDiag}},
                                 IndsR::Type)
  return similar_type(promote_type(TensorT1,TensorT2), IndsR)
end

function contraction_output_type(TensorT1::Type{<:UniformDiagTensor},
                                 TensorT2::Type{<:DiagTensor{<:Number,<:CuDiag}},
                                 IndsR::Type)
  return similar_type(promote_type(TensorT1,TensorT2), IndsR)
end
contraction_output_type(TensorT1::Type{<:DiagTensor{<:Number,<:CuDiag}},
                        TensorT2::Type{<:UniformDiagTensor},
                        IndsR::Type) = contraction_output_type(TensorT2,TensorT1,IndsR)

function zero_contraction_output(T1::UniformDiagTensor{ElT1, N1},
                                 T2::CuDenseTensor{ElT2,N2},
                                 indsR::IndsR) where {ElT1, N1, ElT2, N2, IndsR}
  ElT3 = promote_type(ElT1, ElT2)
  dat  = CUDA.zeros(ElT3, dim(indsR))
  return Tensor(Dense(dat), indsR)
end
zero_contraction_output(T2::CuDenseTensor{ElT2,N2}, T1::UniformDiagTensor{ElT1, N1}, indsR::IndsR) where {ElT1, N1, ElT2, N2, IndsR} = zero_contraction_output(T1, T2, indsR)

function zero_contraction_output(T1::TensorT1, T2::TensorT2, indsR) where {TensorT1<:NonuniformDiagTensor, TensorT2<:NonuniformCuDiagTensor}
  ElT3 = promote_type(eltype(TensorT1), eltype(TensorT2))
  dat  = CUDA.zeros(ElT3, length(data(store(T2))))
  return Tensor(Diag(dat), indsR)
end

function zero_contraction_output(T1::TensorT1, T2::TensorT2, indsR) where {TensorT1<:UniformDiagTensor, TensorT2<:NonuniformCuDiagTensor}
  ElT3 = promote_type(eltype(TensorT1), eltype(TensorT2))
  dat  = CUDA.zeros(ElT3, length(data(store(T2))))
  return Tensor(Diag(dat), indsR)
end
zero_contraction_output(T1::TensorT1, T2::TensorT2, indsR) where {TensorT2<:UniformDiagTensor, TensorT1<:NonuniformCuDiagTensor} = zero_contraction_output(T2, T1, indsR)

function zero_contraction_output(T1::TensorT1, T2::TensorT2, indsR) where {TensorT1<:DiagTensor, TensorT2<:CuDenseTensor}
  ElT3 = promote_type(eltype(TensorT1), eltype(TensorT2))
  dat  = CUDA.zeros(ElT3, length(data(store(T2))))
  return Tensor(Dense(dat), indsR)
end
zero_contraction_output(T1::TensorT1, T2::TensorT2, indsR) where {TensorT2<:DiagTensor, TensorT1<:CuDenseTensor} = zero_contraction_output(T2, T1, indsR)

function contraction_output(T1::UniformDiagTensor,T2::DiagTensor{Elt2,<:CuDiag},indsR) where{Elt2}
    zero_contraction_output(T1,T2,indsR)
end
contraction_output(T1::DiagTensor{Elt1, <:CuDiag},T2::UniformDiagTensor,indsR) where {Elt1} = contraction_output(T2,T1,indsR)

function contract!(C::CuDenseTensor{<:Number, NC},Clabels,
                   A::UniformDiagTensor{<:Number, NA},Alabels,
                   B::CuDenseTensor{<:Number, NB},Blabels) where {NC, NA, NB}
    Bstore = data(store(B))
    Astore = data(store(A))
    Cstore = data(store(C))
    copyto!(Cstore, Astore.*Bstore)
end

function contract!(C::CuDenseTensor{<:Number, NC},Clabels,
                   A::CuDenseTensor{<:Number, NA},Alabels,
                   B::UniformDiagTensor{<:Number, NB},Blabels) where {NC, NA, NB}
    contract!(C, Clabels, B, Blabels, A, Alabels)
end

function contract!(C::NonuniformCuDiagTensor{EltC,NC, IndsC},Clabels,
                   A::UniformDiagTensor{EltA, NA, IndsA},Alabels,
                   B::NonuniformCuDiagTensor{EltB, NB, IndsB},Blabels) where {EltC<:Number, EltB<:Number, EltA<:Number, NC, NB, NA, IndsA, IndsB, IndsC}
    Bstore = data(store(B))
    Astore = data(store(A))
    Cstore = data(store(C))
    copyto!(Cstore, Astore.*Bstore)
end

function contract!(C::NonuniformCuDiagTensor{EltC, NC, IndsC},Clabels,
                   B::NonuniformCuDiagTensor{EltB, NB, IndsB},Blabels,
                   A::UniformDiagTensor{EltA, NA, IndsA},Alabels) where {EltC<:Number, EltB<:Number, EltA<:Number, NC, NB, NA, IndsA, IndsB, IndsC}
    Bstore = data(store(B))
    Astore = data(store(A))
    Cstore = data(store(C))
    copyto!(Cstore, Astore.*Bstore)
end

function contract!(C::NonuniformCuDiagTensor{EltC, NC, IndsC},Clabels,
                   B::NonuniformCuDiagTensor{EltB, NB, IndsB},Blabels,
                   A::NonuniformCuDiagTensor{EltA, NA, IndsA},Alabels) where {EltC<:Number, EltB<:Number, EltA<:Number, NC, NB, NA, IndsA, IndsB, IndsC}
    Bstore = data(store(B))
    Astore = data(store(A))
    Cstore = data(store(C))
    copyto!(Cstore, Astore.*Bstore)
end

# Dense * NonuniformCuDiag
function contract!(C::CuDenseTensor,Clabels,
                   A::NonuniformCuDiagTensor,Alabels,
                   B::CuDenseTensor,Blabels)
    Astore = data(store(A))
    newAstore = CUDA.zeros(eltype(A), dims(inds(A))[1], dims(inds(A))[2])
    adi = diagind(newAstore, 0)
    newAstore[adi] .= Astore[:]
    newA = Tensor(Dense(vec(newAstore)), inds(A))
    contract!(C, Clabels, newA, Alabels, B, Blabels)
end

contract!(C::CuDenseTensor,Clabels,
                   A::CuDenseTensor,Alabels,
                   B::NonuniformCuDiagTensor,Blabels) = contract!(C, Clabels, B, Blabels, A, Alabels)


text:
module SimpleDiscretes

import Base: minimum, maximum
import Statistics

import Distributions
const Dst = Distributions

export SimpleDiscrete

struct SimpleDiscrete{T1<:Real, T2<:Real} <: Dst.DiscreteUnivariateDistribution
    xs::Vector{T1}
    ps::Vector{T2}
end

maximum(d::SimpleDiscrete) = isempty(d.xs) ? -Inf64 : maximum(d.xs)
minimum(d::SimpleDiscrete) = isempty(d.xs) ? +Inf64 : minimum(d.xs)

Statistics.mean(d::SimpleDiscrete) = sum(d.xs .* d.ps)
Statistics.var(d::SimpleDiscrete)  = Statistics.mean(d.xs.^2) - Statistics.mean(d.xs)^2

Dst.insupport(d::SimpleDiscrete, x::Real) = minimum(d) <= x && x <= maximum(d)

function Dst.cdf(d::SimpleDiscrete, x::Real)
    return sum(pᵢ for (xᵢ, pᵢ) in zip(md.mass_points, md.mass_probs) if xᵢ <= x)
end

function Dst.pdf(d::SimpleDiscrete{T1,T2}, x::Real) where {T1,T2}
    for (xᵢ , pᵢ) in zip(d.xs,d.ps)
        xᵢ ≈ x && return pᵢ
    end
    return zero(T2)
end

Dst.ncomponents(d::SimpleDiscrete) = length(d.xs)

Dst.params(d::SimpleDiscrete) = (d.xs, d.ps)

end


text:
include("findinstallexec.jl");

version = "Version0.1";

if Sys.isapple()
    brewstatus0 = Sys.which("brew");
    brewstatus1 = Sys.which("/usr/bin/brew");
    brewstatus2 = Sys.which("/usr/local/bin/brew");
    if brewstatus0 != nothing
        brew = "brew";
    elseif brewstatus1 != nothing
        brew = "/usr/bin/brew";
    elseif brewstatus2 != nothing
        brew = "/usr/local/bin/brew";
    else
        print("Homebrew is not installed on your computer.\n")
        print("Please visit https://brew.sh/ to install Homebrew.\n");
        print("Run this script again after installing Homebrew.\n");
        throw("Installation terminated because Homebrew is not installed yet");
    end
else
    brew = "";
end

print("Exasim find and install the required packages...\n");

gcc, appinstall = findinstallexec("g++", "gcc", brew, 0);
if appinstall==1
    gcc = findinstallexec("g++", "gcc", brew, 0);
end

mpi, appinstall = findinstallexec("mpicxx", "openmpi", brew, 0);
if appinstall==1
    mpi = findinstallexec("mpicxx", "openmpi", brew, 0);
end

nvcc, appinstall = findinstallexec("nvcc", "nvcc", brew, 10);
if appinstall==1
    nvcc = findinstallexec("nvcc", "nvcc", brew, 10);
end

metis, appinstall = findinstallexec("mpmetis", "metis", brew, 0);
if appinstall==1
    metis = findinstallexec("mpmetis", "metis", brew, 0);
end

#python = findinstallexec("python", "python", brew, 0);
gmsh, appinstall = findinstallexec("gmsh", "gmsh", brew, 0);
if appinstall==1
    gmsh = findinstallexec("gmsh", "gmsh", brew, 0);
end

if Sys.isapple()
    #julia = findinstallexec("julia", "julia", brew, 1);
    paraview,appinstall = findinstallexec("paraview", "paraview", brew, 1);
    if appinstall==1
        paraview = findinstallexec("paraview", "paraview", brew, 1);
    end
else
    #julia = findinstallexec("julia", "julia", brew, 0);
    paraview,appinstall = findinstallexec("paraview", "paraview", brew, 0);
    if appinstall==1
        paraview = findinstallexec("paraview", "paraview", brew, 0);
    end
end

if Sys.isapple()
    print("Installing " * "openblas" * " via brew.");
    run(string2cmd(brew * " install openblas"));
    print("Installing " * "lapack" * " via brew.");
    run(string2cmd(brew * " install lapack"));
else
    print("Installing " * "blas" * " via apt.");
    run(string2cmd("sudo apt install libblas-dev"));
    print("Installing " * "lapack" * " via apt.");
    run(string2cmd("sudo apt install liblapack-dev"));
end

cdir = pwd(); ii = findlast("Exasim", cdir);
versiondir = cdir[1:ii[end]] * "/" * version * "/Julia";

filem = versiondir * "/Preprocessing/initializepde.jl";
text = read(filem, String);

q = """ " """;
q = q[2];

oldgcc = "pde.cpucompiler = " * q * "g++" * q;
p = q * gcc * q;
newgcc = "pde.cpucompiler = " * p;
newtext = replace(text, oldgcc => newgcc);

oldmpi = "pde.mpicompiler = " * q * "mpicxx" * q;
p = q * mpi * q;
newmpi = "pde.mpicompiler = " * p;
newtext = replace(newtext, oldmpi => newmpi);

oldnvcc = "pde.gpucompiler = " * q * "nvcc" * q;
p = q * nvcc * q;
newnvcc = "pde.gpucompiler = " * p;
newtext = replace(newtext, oldnvcc => newnvcc);

oldmpirun = "pde.mpirun = " * q * "mpirun" * q;
p = q * mpi * q;
p = replace(p, "mpicxx" => "mpirun");
newmpirun = "pde.mpirun = " * p;
newtext = replace(newtext, oldmpirun => newmpirun);

oldmetis = "pde.metis = " * q * "mpmetis" * q;
p = q * metis * q;
newmetis = "pde.metis = " * p;
newtext = replace(newtext, oldmetis => newmetis);

oldgmsh = "pde.gmsh = " * q * "gmsh" * q;
p = q * gmsh * q;
newgmsh = "pde.gmsh = " * p;
newtext = replace(newtext, oldgmsh => newgmsh);

oldparaview = "pde.paraview = " * q * "paraview" * q;
p = q * paraview * q;
newparaview = "pde.paraview = " * p;
newtext = replace(newtext, oldparaview => newparaview);

fid = open(filem, "w");
write(fid, newtext);
close(fid);

max_stars_repo_name: cesmix/Exasim

max_issues_repo_name: cesmix/Exasim

max_forks_repo_name: cesmix/Exasim


text:
import eve
using DataFrames
o = get(eve.latestMarketLog())
t = eve.importMarketLog()
s = eve.summerizeItem(t)

df = eve.importMarketLog()


text:
const ⊗ = kron

function test_main()

  srand(1234)
  A = rand(3,3); A = 0.5*(A+A')
  B = rand(3,3); B = 0.5*(B+B')
  C = rand(3,3); C = 0.5*(C+C')

  term1 = A ⊗ eye(3) ⊗ eye(3)
  term2 = eye(3) ⊗ B ⊗ eye(3)
  term3 = eye(3) ⊗ eye(3) ⊗ C

  Lapl = term1 + term2 + term3

  invLapl = inv(Lapl)

  invterm1 = inv(term1)
  invterm2 = inv(term2)
  invterm3 = inv(term3)

  # these should be zeros
  diff1 = invterm1 - inv(A) ⊗ eye(3) ⊗ eye(3)
  diff2 = invterm2 - eye(3) ⊗ inv(B) ⊗ eye(3)
  diff3 = invterm3 - eye(3) ⊗ eye(3) ⊗ inv(C)

  # TODO: Find relationship between invterm1 and invLapl

end