text:
int_rules_3_5 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 3.5*Miscellaneous*logarithms =#
    @apply_utils Antiderivative(~u * log(~v), ~x) => With([w = DerivativeDivides(~v, ~u * (1 - ~v), ~x)], w * PolyLog(2, 1 - ~v) <-- Not(FalseQ(w)))
    @apply_utils Antiderivative((~(a') + ~(b') * log(~u)) * log(~v) * ~w, ~x) => With([z = DerivativeDivides(~v, ~w * (1 - ~v), ~x)], z * (~a + ~b * log(~u)) * PolyLog(2, 1 - ~v) - ~b * Antiderivative(SimplifyIntegrand((z * PolyLog(2, 1 - ~v) * D(~u, ~x)) / ~u, ~x), ~x) <-- Not(FalseQ(z))) <-- FreeQ([~a, ~b], ~x) && InverseFunctionFreeQ(~u, ~x)
    @apply_utils Antiderivative(log(~(c') * log(~(d') * (~x) ^ ~(n')) ^ ~(p')), ~x) => ~x * log(~c * log(~d * (~x) ^ ~n) ^ ~p) - ~n * ~p * Antiderivative(1 / log(~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * log(~(c') * log(~(d') * (~x) ^ ~(n')) ^ ~(p'))) / ~x, ~x) => (log(~d * (~x) ^ ~n) * (~a + ~b * log(~c * log(~d * (~x) ^ ~n) ^ ~p))) / ~n - ~b * ~p * log(~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * log(~(c') * log(~(d') * (~x) ^ ~(n')) ^ ~(p'))), ~x) => ((~e * ~x) ^ (~m + 1) * (~a + ~b * log(~c * log(~d * (~x) ^ ~n) ^ ~p))) / (~e * (~m + 1)) - ((~b * ~n * ~p) / (~m + 1)) * Antiderivative((~e * ~x) ^ ~m / log(~d * (~x) ^ ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(a') + ~(b') * log(~(c') * (~RFx) ^ ~(p'))) ^ ~(n'), ~x) => ~x * (~a + ~b * log(~c * (~RFx) ^ ~p)) ^ ~n - ~b * ~n * ~p * Antiderivative(SimplifyIntegrand((~x * (~a + ~b * log(~c * (~RFx) ^ ~p)) ^ (~n - 1) * D(~RFx, ~x)) / ~RFx, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (RationalFunctionQ(~RFx, ~x) && IGtQ(~n, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * log(~(c') * (~RFx) ^ ~(p'))) ^ ~(n') / (~(d') + ~(e') * ~x), ~x) => (log(~d + ~e * ~x) * (~a + ~b * log(~c * (~RFx) ^ ~p)) ^ ~n) / ~e - ((~b * ~n * ~p) / ~e) * Antiderivative((log(~d + ~e * ~x) * (~a + ~b * log(~c * (~RFx) ^ ~p)) ^ (~n - 1) * D(~RFx, ~x)) / ~RFx, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && (RationalFunctionQ(~RFx, ~x) && IGtQ(~n, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * log(~(c') * (~RFx) ^ ~(p'))) ^ ~(n'), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * log(~c * (~RFx) ^ ~p)) ^ ~n) / (~e * (~m + 1)) - ((~b * ~n * ~p) / (~e * (~m + 1))) * Antiderivative(SimplifyIntegrand(((~d + ~e * ~x) ^ (~m + 1) * (~a + ~b * log(~c * (~RFx) ^ ~p)) ^ (~n - 1) * D(~RFx, ~x)) / ~RFx, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~p], ~x) && (RationalFunctionQ(~RFx, ~x) && (IGtQ(~n, 0) && ((EqQ(~n, 1) || IntegerQ(~m)) && NeQ(~m, -1))))
    @apply_utils Antiderivative(log(~(c') * (~RFx) ^ ~(n')) / (~d + ~(e') * (~x) ^ 2), ~x) => With([u = IntHide(1 / (~d + ~e * (~x) ^ 2), ~x)], u * log(~c * (~RFx) ^ ~n) - ~n * Antiderivative(SimplifyIntegrand((u * D(~RFx, ~x)) / ~RFx, ~x), ~x)) <-- FreeQ([~c, ~d, ~e, ~n], ~x) && (RationalFunctionQ(~RFx, ~x) && Not(PolynomialQ(~RFx, ~x)))
    @apply_utils Antiderivative(log(~(c') * (~Px) ^ ~(n')) / ~Qx, ~x) => With([u = IntHide(1 / ~Qx, ~x)], u * log(~c * (~Px) ^ ~n) - ~n * Antiderivative(SimplifyIntegrand((u * D(~Px, ~x)) / ~Px, ~x), ~x)) <-- FreeQ([~c, ~n], ~x) && (QuadraticQ([~Qx, ~Px], ~x) && EqQ(D(~Px / ~Qx, ~x), 0))
    @apply_utils Antiderivative(~RGx * (~(a') + ~(b') * log(~(c') * (~RFx) ^ ~(p'))) ^ ~(n'), ~x) => With([u = ExpandIntegrand((~a + ~b * log(~c * (~RFx) ^ ~p)) ^ ~n, ~RGx, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (RationalFunctionQ(~RFx, ~x) && (RationalFunctionQ(~RGx, ~x) && IGtQ(~n, 0)))
    @apply_utils Antiderivative(~RGx * (~(a') + ~(b') * log(~(c') * (~RFx) ^ ~(p'))) ^ ~(n'), ~x) => With([u = ExpandIntegrand(~RGx * (~a + ~b * log(~c * (~RFx) ^ ~p)) ^ ~n, ~x)], Antiderivative(u, ~x) <-- SumQ(u)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (RationalFunctionQ(~RFx, ~x) && (RationalFunctionQ(~RGx, ~x) && IGtQ(~n, 0)))
    @apply_utils Antiderivative(~RFx * (~(a') + ~(b') * log(~u)), ~x) => With([lst = SubstForFractionalPowerOfLinear(~RFx * (~a + ~b * log(~u)), ~x)], lst(2) * lst(4) * Subst(Antiderivative(lst(1), ~x), ~x, lst(3) ^ (1 / lst(2))) <-- Not(FalseQ(lst))) <-- FreeQ([~a, ~b], ~x) && RationalFunctionQ(~RFx, ~x)
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * log(1 + ~(e') * ((~F) ^ (~(c') * (~(a') + ~(b') * ~x))) ^ ~(n')), ~x) => (-((~f + ~g * ~x) ^ ~m) * PolyLog(2, -(~e) * ((~F) ^ (~c * (~a + ~b * ~x))) ^ ~n)) / (~b * ~c * ~n * log(~F)) + ((~g * ~m) / (~b * ~c * ~n * log(~F))) * Antiderivative((~f + ~g * ~x) ^ (~m - 1) * PolyLog(2, -(~e) * ((~F) ^ (~c * (~a + ~b * ~x))) ^ ~n), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~e, ~f, ~g, ~n], ~x) && GtQ(~m, 0)
    @apply_utils Antiderivative((~(f') + ~(g') * ~x) ^ ~(m') * log(~d + ~(e') * ((~F) ^ (~(c') * (~(a') + ~(b') * ~x))) ^ ~(n')), ~x) => (((~f + ~g * ~x) ^ (~m + 1) * log(~d + ~e * ((~F) ^ (~c * (~a + ~b * ~x))) ^ ~n)) / (~g * (~m + 1)) - ((~f + ~g * ~x) ^ (~m + 1) * log(1 + (~e / ~d) * ((~F) ^ (~c * (~a + ~b * ~x))) ^ ~n)) / (~g * (~m + 1))) + Antiderivative((~f + ~g * ~x) ^ ~m * log(1 + (~e / ~d) * ((~F) ^ (~c * (~a + ~b * ~x))) ^ ~n), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (GtQ(~m, 0) && NeQ(~d, 1))
    @apply_utils Antiderivative(log(~(d') + ~(e') * ~x + ~(f') * sqrt(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2)), ~x) => ~x * log(~d + ~e * ~x + ~f * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)) + (((~f) ^ 2 * ((~b) ^ 2 - 4 * ~a * ~c)) / 2) * Antiderivative(~x / ((2 * ~d * ~e - ~b * (~f) ^ 2) * (~a + ~b * ~x + ~c * (~x) ^ 2) - ~f * ((~b * ~d - 2 * ~a * ~e) + (2 * ~c * ~d - ~b * ~e) * ~x) * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && EqQ((~e) ^ 2 - ~c * (~f) ^ 2, 0)
    @apply_utils Antiderivative(log(~(d') + ~(e') * ~x + ~(f') * sqrt(~(a') + ~(c') * (~x) ^ 2)), ~x) => ~x * log(~d + ~e * ~x + ~f * sqrt(~a + ~c * (~x) ^ 2)) - ~a * ~c * (~f) ^ 2 * Antiderivative(~x / (~d * ~e * (~a + ~c * (~x) ^ 2) + ~f * (~a * ~e - ~c * ~d * ~x) * sqrt(~a + ~c * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f], ~x) && EqQ((~e) ^ 2 - ~c * (~f) ^ 2, 0)
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * log(~(d') + ~(e') * ~x + ~(f') * sqrt(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2)), ~x) => ((~g * ~x) ^ (~m + 1) * log(~d + ~e * ~x + ~f * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2))) / (~g * (~m + 1)) + (((~f) ^ 2 * ((~b) ^ 2 - 4 * ~a * ~c)) / (2 * ~g * (~m + 1))) * Antiderivative((~g * ~x) ^ (~m + 1) / ((2 * ~d * ~e - ~b * (~f) ^ 2) * (~a + ~b * ~x + ~c * (~x) ^ 2) - ~f * ((~b * ~d - 2 * ~a * ~e) + (2 * ~c * ~d - ~b * ~e) * ~x) * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m], ~x) && (EqQ((~e) ^ 2 - ~c * (~f) ^ 2, 0) && (NeQ(~m, -1) && IntegerQ(2 * ~m)))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * log(~(d') + ~(e') * ~x + ~(f') * sqrt(~(a') + ~(c') * (~x) ^ 2)), ~x) => ((~g * ~x) ^ (~m + 1) * log(~d + ~e * ~x + ~f * sqrt(~a + ~c * (~x) ^ 2))) / (~g * (~m + 1)) - ((~a * ~c * (~f) ^ 2) / (~g * (~m + 1))) * Antiderivative((~g * ~x) ^ (~m + 1) / (~d * ~e * (~a + ~c * (~x) ^ 2) + ~f * (~a * ~e - ~c * ~d * ~x) * sqrt(~a + ~c * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~m], ~x) && (EqQ((~e) ^ 2 - ~c * (~f) ^ 2, 0) && (NeQ(~m, -1) && IntegerQ(2 * ~m)))
    @apply_utils Antiderivative(~(v') * log(~(d') + ~(e') * ~x + ~(f') * sqrt(~u)), ~x) => Antiderivative(~v * log(~d + ~e * ~x + ~f * sqrt(ExpandToSum(~u, ~x))), ~x) <-- FreeQ([~d, ~e, ~f], ~x) && (QuadraticQ(~u, ~x) && (Not(QuadraticMatchQ(~u, ~x)) && (EqQ(~v, 1) || MatchQ(~v, (~(g') * ~x) ^ ~(m') <-- FreeQ([g, m], ~x)))))
    @apply_utils Antiderivative(log(~(c') * (~x) ^ ~(n')) ^ ~(r') / (~x * (~(a') * (~x) ^ ~(m') + ~(b') * log(~(c') * (~x) ^ ~(n')) ^ ~q)), ~x) => log(~a * (~x) ^ ~m + ~b * log(~c * (~x) ^ ~n) ^ ~q) / (~b * ~n * ~q) - ((~a * ~m) / (~b * ~n * ~q)) * Antiderivative((~x) ^ (~m - 1) / (~a * (~x) ^ ~m + ~b * log(~c * (~x) ^ ~n) ^ ~q), ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~q, ~r], ~x) && EqQ(~r, ~q - 1)
    @apply_utils Antiderivative((log(~(c') * (~x) ^ ~(n')) ^ ~(r') * (~(a') * (~x) ^ ~(m') + ~(b') * log(~(c') * (~x) ^ ~(n')) ^ ~q) ^ ~(p')) / ~x, ~x) => Antiderivative(ExpandIntegrand(log(~c * (~x) ^ ~n) ^ ~r / ~x, (~a * (~x) ^ ~m + ~b * log(~c * (~x) ^ ~n) ^ ~q) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p, ~q, ~r], ~x) && (EqQ(~r, ~q - 1) && IGtQ(~p, 0))
    @apply_utils Antiderivative((log(~(c') * (~x) ^ ~(n')) ^ ~(r') * (~(a') * (~x) ^ ~(m') + ~(b') * log(~(c') * (~x) ^ ~(n')) ^ ~q) ^ ~(p')) / ~x, ~x) => (~a * (~x) ^ ~m + ~b * log(~c * (~x) ^ ~n) ^ ~q) ^ (~p + 1) / (~b * ~n * ~q * (~p + 1)) - ((~a * ~m) / (~b * ~n * ~q)) * Antiderivative((~x) ^ (~m - 1) * (~a * (~x) ^ ~m + ~b * log(~c * (~x) ^ ~n) ^ ~q) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p, ~q, ~r], ~x) && (EqQ(~r, ~q - 1) && NeQ(~p, -1))
    @apply_utils Antiderivative((~(d') * (~x) ^ ~(m') + ~(e') * log(~(c') * (~x) ^ ~(n')) ^ ~(r')) / (~x * (~(a') * (~x) ^ ~(m') + ~(b') * log(~(c') * (~x) ^ ~(n')) ^ ~q)), ~x) => (~e * log(~a * (~x) ^ ~m + ~b * log(~c * (~x) ^ ~n) ^ ~q)) / (~b * ~n * ~q) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~q, ~r], ~x) && (EqQ(~r, ~q - 1) && EqQ(~a * ~e * ~m - ~b * ~d * ~n * ~q, 0))
    @apply_utils Antiderivative((~u + ~(d') * (~x) ^ ~(m') + ~(e') * log(~(c') * (~x) ^ ~(n')) ^ ~(r')) / (~x * (~(a') * (~x) ^ ~(m') + ~(b') * log(~(c') * (~x) ^ ~(n')) ^ ~q)), ~x) => (~e * log(~a * (~x) ^ ~m + ~b * log(~c * (~x) ^ ~n) ^ ~q)) / (~b * ~n * ~q) + Antiderivative(~u / (~x * (~a * (~x) ^ ~m + ~b * log(~c * (~x) ^ ~n) ^ ~q)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~q, ~r], ~x) && (EqQ(~r, ~q - 1) && EqQ(~a * ~e * ~m - ~b * ~d * ~n * ~q, 0))
    @apply_utils Antiderivative((~(d') * (~x) ^ ~(m') + ~(e') * log(~(c') * (~x) ^ ~(n')) ^ ~(r')) / (~x * (~(a') * (~x) ^ ~(m') + ~(b') * log(~(c') * (~x) ^ ~(n')) ^ ~q)), ~x) => (~e * log(~a * (~x) ^ ~m + ~b * log(~c * (~x) ^ ~n) ^ ~q)) / (~b * ~n * ~q) - ((~a * ~e * ~m - ~b * ~d * ~n * ~q) / (~b * ~n * ~q)) * Antiderivative((~x) ^ (~m - 1) / (~a * (~x) ^ ~m + ~b * log(~c * (~x) ^ ~n) ^ ~q), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~q, ~r], ~x) && (EqQ(~r, ~q - 1) && NeQ(~a * ~e * ~m - ~b * ~d * ~n * ~q, 0))
    @apply_utils Antiderivative(((~(d') * (~x) ^ ~(m') + ~(e') * log(~(c') * (~x) ^ ~(n')) ^ ~(r')) * (~(a') * (~x) ^ ~(m') + ~(b') * log(~(c') * (~x) ^ ~(n')) ^ ~q) ^ ~(p')) / ~x, ~x) => (~e * (~a * (~x) ^ ~m + ~b * log(~c * (~x) ^ ~n) ^ ~q) ^ (~p + 1)) / (~b * ~n * ~q * (~p + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p, ~q, ~r], ~x) && (EqQ(~r, ~q - 1) && (NeQ(~p, -1) && EqQ(~a * ~e * ~m - ~b * ~d * ~n * ~q, 0)))
    @apply_utils Antiderivative(((~(d') * (~x) ^ ~(m') + ~(e') * log(~(c') * (~x) ^ ~(n')) ^ ~(r')) * (~(a') * (~x) ^ ~(m') + ~(b') * log(~(c') * (~x) ^ ~(n')) ^ ~q) ^ ~(p')) / ~x, ~x) => (~e * (~a * (~x) ^ ~m + ~b * log(~c * (~x) ^ ~n) ^ ~q) ^ (~p + 1)) / (~b * ~n * ~q * (~p + 1)) - ((~a * ~e * ~m - ~b * ~d * ~n * ~q) / (~b * ~n * ~q)) * Antiderivative((~x) ^ (~m - 1) * (~a * (~x) ^ ~m + ~b * log(~c * (~x) ^ ~n) ^ ~q) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p, ~q, ~r], ~x) && (EqQ(~r, ~q - 1) && (NeQ(~p, -1) && NeQ(~a * ~e * ~m - ~b * ~d * ~n * ~q, 0)))
    @apply_utils Antiderivative((~(d') * (~x) ^ ~(m') + ~(e') * (~x) ^ ~(m') * log(~(c') * (~x) ^ ~(n')) + ~(f') * log(~(c') * (~x) ^ ~(n')) ^ ~(q')) / (~x * (~(a') * (~x) ^ ~(m') + ~(b') * log(~(c') * (~x) ^ ~(n')) ^ ~q) ^ 2), ~x) => (~d * log(~c * (~x) ^ ~n)) / (~a * ~n * (~a * (~x) ^ ~m + ~b * log(~c * (~x) ^ ~n) ^ ~q)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~q], ~x) && (EqQ(~e * ~n + ~d * ~m, 0) && EqQ(~a * ~f + ~b * ~d * (~q - 1), 0))
    @apply_utils Antiderivative((~d + ~(e') * log(~(c') * (~x) ^ ~(n'))) / (~(a') * ~x + ~(b') * log(~(c') * (~x) ^ ~(n')) ^ ~q) ^ 2, ~x) => (-(~e) * log(~c * (~x) ^ ~n)) / (~a * (~a * ~x + ~b * log(~c * (~x) ^ ~n) ^ ~q)) + ((~d + ~e * ~n) / ~a) * Antiderivative(1 / (~x * (~a * ~x + ~b * log(~c * (~x) ^ ~n) ^ ~q)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~q], ~x) && EqQ(~d + ~e * ~n * ~q, 0)
    @apply_utils Antiderivative(log(~u), ~x) => ~x * log(~u) - Antiderivative(SimplifyIntegrand((~x * D(~u, ~x)) / ~u, ~x), ~x) <-- InverseFunctionFreeQ(~u, ~x)
    @apply_utils Antiderivative(log(~u), ~x) => ~x * log(~u) - Antiderivative(SimplifyIntegrand(~x * Simplify(D(~u, ~x) / ~u), ~x), ~x) <-- ProductQ(~u)
    @apply_utils Antiderivative(log(~u) / (~(a') + ~(b') * ~x), ~x) => (log(~a + ~b * ~x) * log(~u)) / ~b - (1 / ~b) * Antiderivative(SimplifyIntegrand((log(~a + ~b * ~x) * D(~u, ~x)) / ~u, ~x), ~x) <-- FreeQ([~a, ~b], ~x) && (RationalFunctionQ(D(~u, ~x) / ~u, ~x) && (NeQ(~a, 0) || Not(BinomialQ(~u, ~x) && EqQ(BinomialDegree(~u, ~x) ^ 2, 1))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~(m') * log(~u), ~x) => ((~a + ~b * ~x) ^ (~m + 1) * log(~u)) / (~b * (~m + 1)) - (1 / (~b * (~m + 1))) * Antiderivative(SimplifyIntegrand(((~a + ~b * ~x) ^ (~m + 1) * D(~u, ~x)) / ~u, ~x), ~x) <-- FreeQ([~a, ~b, ~m], ~x) && (InverseFunctionFreeQ(~u, ~x) && NeQ(~m, -1))
    @apply_utils Antiderivative(log(~u) / ~Qx, ~x) => With([v = IntHide(1 / ~Qx, ~x)], v * log(~u) - Antiderivative(SimplifyIntegrand((v * D(~u, ~x)) / ~u, ~x), ~x)) <-- QuadraticQ(~Qx, ~x) && InverseFunctionFreeQ(~u, ~x)
    @apply_utils Antiderivative((~u) ^ (~(a') * ~x) * log(~u), ~x) => (~u) ^ (~a * ~x) / ~a - Antiderivative(SimplifyIntegrand(~x * (~u) ^ (~a * ~x - 1) * D(~u, ~x), ~x), ~x) <-- FreeQ(~a, ~x) && InverseFunctionFreeQ(~u, ~x)
    @apply_utils Antiderivative(~v * log(~u), ~x) => With([w = IntHide(~v, ~x)], Dist(log(~u), w, ~x) - Antiderivative(SimplifyIntegrand((w * D(~u, ~x)) / ~u, ~x), ~x) <-- InverseFunctionFreeQ(w, ~x)) <-- InverseFunctionFreeQ(~u, ~x)
    @apply_utils Antiderivative(~v * log(~u), ~x) => With([w = IntHide(~v, ~x)], Dist(log(~u), w, ~x) - Antiderivative(SimplifyIntegrand(w * Simplify(D(~u, ~x) / ~u), ~x), ~x) <-- InverseFunctionFreeQ(w, ~x)) <-- ProductQ(~u)
    @apply_utils Antiderivative(log(~v) * log(~w), ~x) => (~x * log(~v) * log(~w) - Antiderivative(SimplifyIntegrand((~x * log(~w) * D(~v, ~x)) / ~v, ~x), ~x)) - Antiderivative(SimplifyIntegrand((~x * log(~v) * D(~w, ~x)) / ~w, ~x), ~x) <-- InverseFunctionFreeQ(~v, ~x) && InverseFunctionFreeQ(~w, ~x)
    @apply_utils Antiderivative(~u * log(~v) * log(~w), ~x) => With([z = IntHide(~u, ~x)], (Dist(log(~v) * log(~w), z, ~x) - Antiderivative(SimplifyIntegrand((z * log(~w) * D(~v, ~x)) / ~v, ~x), ~x)) - Antiderivative(SimplifyIntegrand((z * log(~v) * D(~w, ~x)) / ~w, ~x), ~x) <-- InverseFunctionFreeQ(z, ~x)) <-- InverseFunctionFreeQ(~v, ~x) && InverseFunctionFreeQ(~w, ~x)
    @apply_utils Antiderivative((~f) ^ (~(a') * log(~u)), ~x) => Antiderivative((~u) ^ (~a * log(~f)), ~x) <-- FreeQ([~a, ~f], ~x)
    #= If(TrueQ($LoadShowSteps), Antiderivative((~u)/(~x),~x) := With([lst=FunctionOfLog(u,x)], ShowStep("","Antiderivative(F(log(a*x^n))/x,x)","Subst(Antiderivative(F(x),x),x,log(a*x^n) )/n",Hold( 1/lst((3))*Subst(Antiderivative(lst((1)),x),x,log(lst((2)))))) <--  Not(FalseQ(lst))) <-- SimplifyFlag && NonsumQ(u), Antiderivative((~u)/(~x),~x) := With([lst=FunctionOfLog(u,x)], 1/lst((3))*Subst(Antiderivative(lst((1)),x),x,log(lst((2)))) <--  Not(FalseQ(lst))) <-- NonsumQ(u)) =#
    If(TrueQ($LoadShowSteps), Antiderivative((~u), ~x) := With([lst = FunctionOfLog(Cancel(x*u), x)], ShowStep("", "Antiderivative(F(log(a*x^n))/x,x)", "Subst(Antiderivative(F(x),x),x,log(a*x^n))/n", Hold( 1/lst((3))* Subst(Antiderivative(lst((1)), x), x, log(lst((2)))))) <-- Not(FalseQ(lst))) <-- SimplifyFlag && NonsumQ(u), Antiderivative((~u), ~x) := With([lst = FunctionOfLog(Cancel(x*u), x)], 1/lst((3))*Subst(Antiderivative(lst((1)), x), x, log(lst((2)))) <-- Not(FalseQ(lst))) <-- NonsumQ(u))
    @apply_utils Antiderivative(~(u') * log(Gamma(~v)), ~x) => (log(Gamma(~v)) - LogGamma(~v)) * Antiderivative(~u, ~x) + Antiderivative(~u * LogGamma(~v), ~x)
    @apply_utils Antiderivative(~(u') * (~(a') * (~x) ^ ~(m') + ~(b') * (~x) ^ ~(r') * log(~(c') * (~x) ^ ~(n')) ^ ~(q')) ^ ~(p'), ~x) => Antiderivative(~u * (~x) ^ (~p * ~r) * (~a * (~x) ^ (~m - ~r) + ~b * log(~c * (~x) ^ ~n) ^ ~q) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p, ~q, ~r], ~x) && IntegerQ(~p)
end

max_stars_repo_path: src/rules/3 Logarithms/3.5 Miscellaneous logarithms.jl


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using ACE
using Test, Printf, LinearAlgebra, StaticArrays, BenchmarkTools, Test
using JuLIP, JuLIP.Testing
using JuLIP: evaluate!, evaluate_d!, evaluate, evaluate_d
using ACE: PolyTransform, PolyCutoff1s, PolyCutoff2s

##

@testset "ACE.jl" begin
    include("pairpots/test_jacobi.jl")
    include("pairpots/test_transforms.jl")
    include("pairpots/test_basics.jl")
    include("pairpots/test_repulsion.jl")
    include("pairpots/test_orthpolys.jl")
    # ----------------------
    include("test_ylm.jl")
    include("test_rylm.jl")
    include("test_cg.jl")
    # ----------------------
    include("test_basis.jl")
    include("test_filter_rpi_alg.jl")
    include("test_fast.jl")
    include("test_real.jl")
    include("test_multispecies.jl")
    include("test_orth.jl")
    # include("test_descriptor.jl")
    # ----------------------
    include("test_compressA.jl")
    # ---------------------- bonds stuff
    include("bonds/test_cylindrical.jl")
    include("bonds/test_fourier.jl")
    include("bonds/test_envpairbasis.jl")
end

ext: jl


text:
typealias HANDLE Ptr{Void}
typealias DWORD UInt32
typealias WORD UInt16
typealias LPTSTR Cwchar_t
typealias LPBYTE Ptr{Char}

immutable STARTUPINFO
	cb::DWORD
	lpReserved::LPTSTR
	lpDesktop::LPTSTR
	lpTitle::LPTSTR
	dwX::DWORD
	dwY::DWORD
	dwXSize::DWORD
	dwYSize::DWORD
	dwXCountChars::DWORD
	dwYCountChars::DWORD
	dwFillAttribute::DWORD
	dwFlags::DWORD
	wShowWindow::WORD
	cbReserved2::WORD
	lpReserved2::LPBYTE
	hStdInput::HANDLE
	hStdOutput::HANDLE
	hStdError::HANDLE
	STARTUPINFO() = new()
end

immutable PROCESS_INFORMATION
	hProcess::HANDLE
	hThread::HANDLE
	dwProcessId::DWORD
	dwThreadId::DWORD
	PROCESS_INFORMATION() = new(C_NULL,C_NULL,C_NULL,C_NULL)
end

CreateProcess(cmd) = begin
si = [STARTUPINFO()]
	pi = [PROCESS_INFORMATION()]
	ccall(:CreateProcessW, Cchar,
	(Ptr{Cwchar_t}, Ptr{Cwchar_t}, Ptr{Int}, Ptr{Int}, Cchar, Int64,
	Ptr{UInt8}, Ptr{UInt8}, Ptr{PROCESS_INFORMATION}, Ptr{STARTUPINFO}),
	C_NULL,
	utf16(cmd),
	C_NULL,
	C_NULL,
	0,
	0,
	C_NULL, C_NULL,
	convert(Ptr{PROCESS_INFORMATION}, pointer(pi)),
	convert(Ptr{STARTUPINFO}, pointer(pi)))
end


text:
@doc raw"""
    Game(s, t)

A two-player, two-strategy, symmetric game. The games are characterized by the diagonal elements of
a payoff matrix

```math
\begin{bmatrix}
    1 & S \\
    T & 0
\end{bmatrix}
```

with ``-0.5 \leq S \leq 0.5`` and ``0.5 \leq T \leq 1.5``,
"""
struct Game
    payoffs::SMatrix{2, 2, Float64}
    function Game(s::Number, t::Number)
        if -0.5 ≤ s ≤ 0.5
            if 0.5 ≤ t ≤ 1.5
                return new([1.0 s; t 0.0])
            end
            throw(DomainError(t, "expected 0.5 ≤ t ≤ 1.5"))
        end
        throw(DomainError(s, "expected -0.5 ≤ s ≤ 0.5"))
    end
end

"""
    sparam(game)

Get the S-parameter of the game.
"""
sparam(g::Game) = g.payoffs[1,2]

"""
    tparam(game)

Get the T-parameter of the game.
"""
tparam(g::Game) = g.payoffs[2,1]

"""
    quadrant(game)

Get the quadrant of S-T space in which the game resides.
"""
function quadrant(g::Game)
    S, T = sparam(g), tparam(g)
    if one(T) ≤ T
        zero(S) ≤ S ? 1 : 4
    else
        zero(S) ≤ S ? 2 : 3
    end
end

"""
    name(quad::Int)

Get the classical name of the games in a quadrant, e.g. Hawk-Dove, Harmony, Stag Hunt or Prisoner's
Dilemma.
"""
function name(quad::Int)
    if quad === 1
        "Hawk-Dove"
    elseif quad === 2
        "Harmony"
    elseif quad === 3
        "Stag Hunt"
    elseif quad === 4
        "Prisoner's Dilemma"
    else
        error("invalid quadrant: $q")
    end
end

"""
    name(game)

Get the classical name of the game, e.g. Hawk-Dove, Harmony, Stag Hunt or Prisoner's Dilemma.
"""
name(g::Game) = name(quadrant(g))

"""
    play(g::Game, a, b)

Play the game with the strategies `a` and `b`, returning the payoff recieved by agent playing `a`.
"""
play(g::Game, a::Int, b::Int) = g.payoffs[a, b]

Base.show(io::IO, g::Game) = print(io, "Game($(sparam(g)), $(tparam(g)))")

const GamesIterator = let ParameterRange = typeof(-0.5:0.1:0.5)
    Iterators.ProductIterator{Tuple{ParameterRange, ParameterRange}}
end

"""
    Games(sstep, tstep)

An iterator over the games with ``S`` and ``T`` parameters varied over their ranges with the
provided step sizes.
"""
struct Games
    iterator::GamesIterator
    Games(sstep, tstep) = let ss = -0.5:sstep:0.5, ts = 0.5:tstep:1.5
        new(product(ss, ts))
    end
end

Base.iterate(g::Games) = let iteree = iterate(g.iterator)
    if iteree === nothing
        iteree
    else
        element, state = iteree
        Game(element...), state
    end
end

Base.iterate(g::Games, state) = let iteree = iterate(g.iterator, state)
    if iteree === nothing
        iteree
    else
        element, state = iteree
        Game(element...), state
    end
end

Base.length(g::Games) = length(g.iterator)

@doc raw"""
    Sigmoid(β)

A rule specifying that the agent will switch its strategy with probability

```math
P(dp) = \frac{1}{1 + e^{-\beta dp}}
```

where ``dp`` is the difference in the two payoffs.
"""
struct Sigmoid <: AbstractRule
    β::Float64
    function Sigmoid(β::Float64 = 1.0)
        if β < zero(β)
            throw(DomainError(β, "must be non-negative"))
        end
        new(β)
    end
end

apply(s::Sigmoid, dp::Float64) = 1.0 / (1.0 + exp(-s.β * dp))

@doc raw"""
    Heaviside(ϵ)

A rule specifying that the agent will switch its strategy with probability

```math
P(dp) = \begin{cases}
    0.0, & dp \leq -\epsilon \\
    0.5, & -\epsilon < dp < \epsilon \\
    1.0, & \epsilon < dp
\end{cases}
```

where ``dp`` is the difference in the two payoffs.
"""
struct Heaviside <: AbstractRule
    ϵ::Float64
    function Heaviside(ϵ::Float64 = 1e-3)
        if ϵ < zero(ϵ)
            throw(DomainError(ϵ, "must be non-negative"))
        end
        new(ϵ)
    end
end

apply(h::Heaviside, dp::Float64) = (abs(dp) < h.ϵ) ? 0.5 : float(dp > 0)

lattice_graph(m::Int, n::Int) = if m < 1
    throw(DomainError(m, "lattice must have at least 1 row"))
elseif n < 1
    throw(DomainError(m, "lattice must have at least 1 column"))
else
    crosspath(m, path_graph(n))
end

"""
    RandomScheme(p=0.5)

The "random" update scheme. Under this scheme, the agent will randomly choose strategy 2 with
probability `p`. This is primarily for validation purposes.
"""
struct RandomScheme <: AbstractScheme
    p::Float64
end
RandomScheme() = RandomScheme(0.5)

function decide(r::RandomScheme,
                arena::AbstractArena{G, RandomScheme},
                ss::AbstractVector{Int},
                ps::AbstractArray{Float64,2}) where {G <: SimpleGraph}
    Int.(rand(length(ss)) .> r.p) .+ 1
end

function decide(::RandomScheme, ::AbstractArena, ::AbstractVector{Int}, ::AbstractVector{Float64})
    error("this method is ill-defined")
end
decide(::RandomScheme, ::AbstractArena, ::Int, ::Float64) = error("this method is ill-defined")

"""
    CounterFactural([rule = Sigmoid()])

The "counter-factual" update scheme. Under this scheme, the agent will switch their strategy
according to some probability. The probability the agent will switch is determined by applying the
scheme's rule to the payoff difference.
"""
struct CounterFactual{Rule <: AbstractRule} <: AbstractScheme
    rule::Rule
end
CounterFactual() = CounterFactual{Sigmoid}(Sigmoid())

function decide(cf::S,
                arena::AbstractArena{G, S},
                s::Int,
                dp::Float64) where {G, S <: CounterFactual}
    (rand() <= apply(cf.rule, dp)) ? (3 - s) : s
end

"""
    MimicRandom()

The "mimic random" update scheme. Under this scheme, the agent will switch their strategy at
uniformly at random to the strategy of one of their neighbors or their own.
"""
struct MimicRandom <: AbstractScheme end

function decide(m::MimicRandom,
                arena::AbstractArena{G, MimicRandom},
                ss::AbstractVector{Int},
                dps::AbstractArray{Float64,2}) where {G <: SimpleGraph}
    map(i -> ss[rand([i; neighbors(arena, i)])], 1:length(ss))
end

function decide(::MimicRandom, ::AbstractArena, ::AbstractVector{Int}, ::AbstractVector{Float64})
    error("this method is ill-defined")
end
decide(::MimicRandom, ::AbstractArena, ::Int, ::Float64) = error("this method is ill-defined")

"""
    MimicBest()

The "mimic best" update scheme. Under this scheme, the agent will switch their strategy to that
of their neighbors who recieved the greatest payout, possibly retaining their own.
"""
struct MimicBest <: AbstractScheme end

function decide(m::MimicBest,
                arena::AbstractArena{G, MimicBest},
                ss::AbstractVector{Int},
                ps::AbstractArray{Float64,2}) where {G <: SimpleGraph}
    map(1:length(ss)) do i
        ns = neighbors(arena, i)
        m, j = findmax(ps[ns])
        m > ps[i] ? ss[ns[j]] : ss[i]
    end
end

function decide(::MimicBest, ::AbstractArena, ::AbstractVector{Int}, ::AbstractVector{Float64})
    error("this method is ill-defined")
end
decide(::MimicBest, ::AbstractArena, ::Int, ::Float64) = error("this method is ill-defined")

"""
    MimicBiased()

The "mimic biased" update scheme. Under this scheme, the agent will switch their strategy
randomly to that of their neighbors with a log-probability proportional to their payoffs.
"""
struct MimicBiased <: AbstractScheme end

function decide(m::MimicBiased,
                arena::AbstractArena{G, MimicBiased},
                ss::AbstractVector{Int},
                ps::AbstractArray{Float64,2}) where {G <: SimpleGraph}
    map(1:length(ss)) do i
        idx = [i; neighbors(arena, i)]
        p = exp.(ps[idx])

        perm = sortperm(p)
        p = p[perm]
        idx = idx[perm]

        j = findfirst(cumsum(p) .>= sum(p)*rand())
        ss[idx[j]]
    end
end

function decide(::MimicBiased, ::AbstractArena, ::AbstractVector{Int}, ::AbstractVector{Float64})
    error("this method is ill-defined")
end
decide(::MimicBiased, ::AbstractArena, ::Int, ::Float64) = error("this method is ill-defined")

"""
    Arena(game, graph, scheme) <: AbstractArena

An arena in which the agents, situated on the vertices of a graph, play a game with their neighbors.
Every pair of agents plays the same game. With each round of play, the agents use the scheme to
decide what their next strategy will be.

In this arena, the agents accumulate payoff by playing with every neighbor. The decision of what
strategy to use in the next time step is based on the difference between the total payoffs for each
of the agents' strategies.
"""
struct Arena{Graph, Scheme} <: AbstractArena{Graph, Scheme}
    game::Game
    graph::Graph
    scheme::Scheme
    function Arena(ga::Game, gr::Graph, s::Scheme) where {Graph, Scheme}
        if nv(gr) < 2
            throw(DomainError(gr, "graph cannot be empty"))
        elseif has_self_loops(gr)
            throw(DomainError(gr, "graph cannot have self-loops"))
        elseif any(map(c -> length(c) < 2, connected_components(gr)))
            throw(DomainError(gr, "each component of the graph must have at least 2 elements"))
        end
        new{Graph, Scheme}(ga, gr, s)
    end
end


text:
using Tk
w = Toplevel("Component Interaction Example")
fr = Frame(w)
pack(fr, {:expand=>true, :fill => "both"})
## The task: For a minimal "application", write a program that
## presents a form with three components to the user: A numeric input
## field ("Value") and two buttons ("increment" and "random").

value = Entry(fr, "")
increment = Button(fr, "Increment")
random = Button(fr, "Random")

formlayout(value, "Value:")
formlayout(increment, " ")
formlayout(random, " ")

set_value(value, "0") ## The field is initialized to zero.

tk_bind(increment, "command") do path  ## increment its value with the "increment" button.
  val = get_value(value) | float
  set_value(value, string(val + 1))
end

function validate_command(path, P)
    try
        if length(P) > 0 parsefloat(P) end
        tcl("expr", "TRUE")
    catch e
        tcl("expr", "FALSE")
    end
end
function invalid_command(path, W)
    println("Invalid value")
    tcl(W, "delete", "@0", "end")
end

tk_configure(value, {:validate=>"key", :validatecommand=>validate_command, :invalidcommand=>invalid_command })

## Pressing the "random" button presents a confirmation dialog, and resets the field's value to a random value if the answer is "Yes".
tk_bind(random, "command") do path
   out = Messagebox(w, "Randomize input", "Select a new random number?")
   if out == "ok"
       new_value = floor(100*rand(1))[1]
       set_value(value, string(new_value))
   end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# MobilityTheories.jl

"""
    polaronmobility(Trange, ε_Inf, ε_S, freq, effectivemass; verbose::Bool=false)

Solves the Feynman polaron problem variationally with finite temperature Osaka energies. From the resulting v, and w parameters, calculates polaron structure (wave function size, etc.).  Uses FHIP, Kadanoff (Boltzmann relaxation time) and Hellwarth direct contour integration to predict a temperature-dependent mobility for the material system.

# Arguments
- `Trange::range`: temperature range.
- `ε_Inf`: reduced optical dielectric constant.
- `ε_S`: reduced static dielectric constant.
- `freq`: characteristic dielectric phonon frequency (THz).
- `effectivemass`: bare-band effective mass (mₑ).

Returns a structure of type `Polaron`, containing arrays of useful
information.  Also prints a lot of information to the standard out - which
may be more useful if you're just inquiring as to a particular data point,
rather than plotting a temperature-dependent parameter.

As an example, to calculate the electron polaron in MAPI at 300 K:
# Examples
```jldoctest
polaronmobility(300, 4.5, 24.1, 2.25, 0.12)
```
"""
function polaronmobility(Trange, ε_Inf, ε_S, freq, effectivemass; verbose::Bool=false)
    println("\n\nPolaron mobility for system ε_Inf=$ε_Inf, ε_S=$ε_S, freq=$freq,
                 effectivemass=$effectivemass; with Trange $Trange ...")

    # Internally we we use 'mb' for the 'band mass' in SI units, of the effecitve-mass of the electron
    mb=effectivemass*MassElectron
    ω = (2*pi)*freq*1e12 # angular-frequency, of the phonon mode

    α=frohlichalpha(ε_Inf, ε_S,  freq,    effectivemass)
    #α=2.395939683378253 # Hard coded; from MAPI params, 4.5, 24.1, 2.25THz, 0.12me

    v=7.1 # starting guess for v,w variational parameters
    w=6.5

    @printf("Polaron mobility input parameters: ε_Inf=%f ε_S=%f freq=%g α=%f \n",ε_Inf, ε_S, freq, α )
    @printf("Derived params in SI: ω =%g mb=%g \n",ω ,mb)

    # Empty struct for storing data
    # A slightly better way of doing this ф_ф ...
    p=Polaron()
    # populate data structure with (athermal) parameters supplied...
    append!(p.α,α) # appending so as not to mess with type immutability
    append!(p.mb,mb)
    append!(p.ω,ω)

    # We define βred as the subsuming the energy of the phonon; i.e. kbT c.f. ħω
    for T in Trange
        β=1/(kB*T)
        βred=ħ*ω*β
        append!(p.βred,βred)
        @printf("T: %f β: %.3g βred: %.3g ħω  = %.3g meV\t",T,β,βred, 1E3* ħ*ω  / q)

        if T==0
            v,w=feynmanvw(α, v=v,w=w)
        else
            v,w=feynmanvw(α,βred, v=v,w=w)
        end

        @printf("\n Polaron Parameters:  v= %.4f  w= %.4f",v,w)

        # From 1962 Feynman, definition of v and w in terms of the coupled Mass and spring-constant
        # See Page 1007, just after equation (18)
        # Units of M appear to be 'electron masses'
        # Unsure of units for k, spring coupling constant
        k=(v^2-w^2)
        M=(v^2-w^2)/w^2
        append!(p.k,k)
        append!(p.M,M)

        @printf("  ||   M=%f  k=%f\t",M,k)
        @printf("\n Bare-band effective mass: %f Polaron effective mass: %f Polaron mass enhancement: %f%%",effectivemass,effectivemass*(1+M),M*100)

        @printf("\n Polaron frequency (SI) v= %.2g Hz  w= %.2g Hz",
                v*ω /(2*pi),w*ω /(2*pi))

		# (46) in Feynman1955
		meSmallAlpha(α )=α /6 + 0.025*α ^2
		# (47) In Feynman1955
		meLargeAlpha(α )=16*α ^4 / (81*π ^4)
		#meLargeAlpha(α )=202*(α /10)^4
        if (verbose) # asymptotic solutions - not that interesting when you have the actual ones!
		    @printf("\n Feynman1955(46,47): meSmallAlpha(α)= %.3f meLargeAlpha(α)= %.3f",
                    meSmallAlpha(α),meLargeAlpha(α))
            @printf("\n Feynman1962: Approximate ~ Large alpha limit, v/w = %.2f  =~approx~= alpha^2 = %.2f ",
                    v/w,α ^2)
        end

        # POLARON SIZE
        @printf("\n Polaron size (rf), following Schultz1959. (s.d. of Gaussian polaron ψ )")
        # Schultz1959 - rather nicely he actually specifies everything down into units!
        # just before (2.4) in Schultz1959
        mu=((v^2-w^2)/v^2)
        # (2.4)
        rf=sqrt(3/(2*mu*v))
        # (2.4) SI scaling inferred from units in (2.5a) and Table II
        rfsi=rf*sqrt(2*me*ω )
        @printf("\n\t Schultz1959(2.4): rf= %g (int units) = %g m [SI]",rf,rfsi )
        append!(p.rfsi,rfsi)

        if (verbose)
        scale=sqrt(2*mb*ω) # Note we're using mb;
        #band effective-mass in SI units (i.e. meff*melectron)

        rfa=(3/(0.44*α ))^0.5 # As given in Schultz1959(2.5a), but that 0.44 is actually 4/9
        @printf("\n\t Schultz1959(2.5a) with 0.44: Feynman α→0 expansion: rfa= %g (int units) = %g m [SI]",rfa,scale*rfa )
        rfa=(3/((4/9)*α ))^0.5 # Rederived from Feynman1955, 8-8-2017; Yellow 2017.B Notebook pp.33-34
        @printf("\n\t Schultz1959(2.5a) with 4/9 re-derivation: Feynman α→0 expansion: rfa= %g (int units) = %g m [SI]",rfa,scale*rfa )
        append!(p.rfsmallalpha,scale*rfa)

        rfb=3*(pi/2)^0.5 * α
        @printf("\n\t Schultz1959(2.5b): Feynman α→∞ expansion: rf= %g (int units) = %g m [SI]",rfb,scale*rfb )

        # Schultz1959 - Between (5.7) and (5.8) - resonance of Feynman SHM system
        phononfreq=sqrt(k/M)
        @printf("\n\t Schultz1959(5.7-5.8): fixed-e: phononfreq= %g (int units) = %g [SI, Hz] = %g [meV]",
            phononfreq,phononfreq*ω /(2*pi), phononfreq*hbar*ω *1000/q)

        phononfreq=sqrt(k/mu) # reduced mass
        @printf("\n\t Schultz1959(5.7-5.8): reducd mass: phononfreq= %g (int units) = %g [SI, Hz] = %g [meV]",
            phononfreq,phononfreq*ω /(2*pi), phononfreq*hbar*ω *1000/q)

        @printf("\n\t Schultz1959: electronfreq= %g (int units) = %g [SI, Hz] = %g [meV]",
            sqrt(k/1),sqrt(k/1)*ω /(2*pi), sqrt(k/1)*hbar*ω *1000/q)
        @printf("\n\t Schultz1959: combinedfreq= %g (int units) = %g [SI, Hz] = %g [meV]",
            sqrt(k/(1+M)),sqrt(k/(1+M))*ω /(2*pi), sqrt(k/(1+M))*hbar*ω *1000/q)

        # Devreese1972: 10.1103/PhysRevB.5.2367
        # p.2371, RHS.
        @printf("\n Devreese1972: (Large Alpha) Franck-Condon frequency = %.2f", 4*α ^2/(9*pi))
        end

        # F(v,w,β,α)=-(A(v,w,β)+B(v,w,β,α)+C(v,w,β)) #(62a) - Hellwarth 1999
        if (T>0)
            @printf("\n Polaron Free Energy: A= %f B= %f C= %f F= %f",A(v,w,βred),B(v,w,βred,α),C(v,w,βred),F(v,w,βred,α))
            @printf("\t = %f meV",1000.0 * F(v,w,βred,α) * ħ*ω  / q)
            append!(p.A,A(v,w,βred))
            append!(p.B,B(v,w,βred,α))
            append!(p.C,C(v,w,βred))
            append!(p.F,F(v,w,βred,α))
        else # Athermal case; Enthalpy
            @printf("\n Polaron Enthalpy: F= %f = %f meV \n",F(v,w,α), 1_000*F(v,w,α) * ħ*ω/q)

            return # return early, as if T=0, all mobility theories = infinite / fall over
        end

        # FHIP
        #    - low-T mobility, final result of Feynman1962
        # [1.60] in Devreese2016 page 36; 6th Edition of Frohlich polaron notes (ArXiv)
        # I believe here β is in SI (expanded) units
        @printf("\n Polaron Mobility theories:")
        μ=(w/v)^3 * (3*q)/(4*mb*ħ*ω^2*α*β) * exp(ħ*ω*β) * exp((v^2-w^2)/(w^2*v))
        @printf("\n\tμ(FHIP)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)
        append!(p.T,T)
        append!(p.FHIPμ,μ*100^2)

        # Kadanoff
        #     - low-T mobility, constructed around Boltzmann equation.
        #     - Adds factor of 3/(2*beta) c.f. FHIP, correcting phonon emission behaviour
        # [1.61] in Devreese2016 - Kadanoff's Boltzmann eqn derived mob
        μ=(w/v)^3 * (q)/(2*mb*ω*α) * exp(ħ*ω*β) * exp((v^2-w^2)/(w^2*v))
        @printf("\n\tμ(Kadanoff,via Devreese2016)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)

        append!(p.Kμ,μ*100^2)

        ######
        # OK, now deep-diving into Kadanoff1963 itself to extract Boltzmann equation components
        # Particularly right-hand-side of page 1367
        #
        # From 1963 Kadanoff, (9), define eqm. number of phonons (just from T and phonon omega)
        Nbar=(exp(βred) - 1)^-1
        @printf("\n\t\tEqm. Phonon. pop. Nbar: %f ",Nbar)
        if (verbose)
            @printf("\n\texp(Bred): %f exp(-Bred): %f exp(Bred)-1: %f",exp(βred),exp(-βred),exp(βred)-1)
        end
        Nbar=exp(-βred)
        #Note - this is only way to get Kadanoff1963 to be self-consistent with
        #FHIP, and later statements (Devreese) of the Kadanoff mobility.
        #It suggests that Kadanoff used the wrong identy for Nbar in 23(b) for
        #the Gamma0 function, and should have used a version with the -1 to
        #account for Bose / phonon statistics

        #myv=sqrt(k*(1+M)/M) # cross-check maths between different papers
        #@printf("\nv: %f myv: %f\n",v,myv)

        # Between 23 and 24 in Kadanoff 1963, for small momenta skip intergration --> Gamma0
        Gamma0=2*α * Nbar * (M+1)^(1/2) * exp(-M/v)
        Gamma0*=ω  #* (ω *hbar) # Kadanoff 1963 uses hbar=omega=mb=1 units
            # Factor of omega to get it as a rate relative to phonon frequency
            # Factor of omega*hbar to get it as a rate per energy window
        μ=q/( mb*(M+1) * Gamma0 ) #(25) Kadanoff 1963, with SI effective mass

        if (verbose) # these are cross-checks
            @printf("\n\tμ(Kadanoff1963 [Eqn. 25]) = %f m^2/Vs \t = %.2f cm^2/Vs",μ,μ*100^2)
            @printf("\n\t\t Eqm. Phonon. pop. Nbar: %f ",Nbar)
        end

        @printf("\n\t\tGamma0 = %g rad/s = %g /s ",
                Gamma0, Gamma0/(2*pi))
        @printf(" \n\t\tTau=1/Gamma0 = %g s = %f ps",
            2*pi/Gamma0, 2*pi*1E12/Gamma0)
        Eloss=hbar*ω * Gamma0/(2*pi) # Simply Energy * Rate
        @printf("\n\t\tEnergy Loss = %g J/s = %g meV/ps",Eloss,Eloss * 1E3 / (q*1E12) )
        append!(p.Tau, 2*pi*1E12/Gamma0) # Boosted into ps ?

        # Hellwarth1999 - directly do contour integration in Feynman1962, for
        # finite temperature DC mobility
        # Hellwarth1999 Eqn (2) and (1) - These are going back to the general
        # (pre low-T limit) formulas in Feynman1962.  to evaluate these, you
        # need to do the explicit contour integration to get the polaron
        # self-energy
        R=(v^2-w^2)/(w^2*v) # inline, page 300 just after Eqn (2)

        #b=R*βred/sinh(b*βred*v/2) # This self-references b! What on Earth?
        # OK! I now understand that there is a typo in Hellwarth1999 and
        # Biaggio1997. They've introduced a spurious b on the R.H.S. compared to
        # the original, Feynman1962:
        b=R*βred/sinh(βred*v/2) # Feynman1962 version; page 1010, Eqn (47b)

        a=sqrt( (βred/2)^2 + R*βred*coth(βred*v/2))
        k(u,a,b,v) = (u^2+a^2-b*cos(v*u))^(-3/2)*cos(u) # integrand in (2)
        K=quadgk(u->k(u,a,b,v),0,Inf)[1] # numerical quadrature integration of (2)

        #Right-hand-side of Eqn 1 in Hellwarth 1999 // Eqn (4) in Baggio1997
        RHS= α/(3*sqrt(π)) * βred^(5/2) / sinh(βred/2) * (v^3/w^3) * K
        μ=RHS^-1 * (q)/(ω*mb)
        @printf("\n\tμ(Hellwarth1999)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)
        append!(p.Hμ,μ*100^2)

        if (verbose)
        # Hellwarth1999/Biaggio1997, b=0 version... 'Setting b=0 makes less than 0.1% error'
        # So let's test this
        R=(v^2-w^2)/(w^2*v) # inline, page 300 just after Eqn (2)
        b=0
        a=sqrt( (βred/2)^2 + R*βred*coth(βred*v/2))
        #k(u,a,b,v) = (u^2+a^2-b*cos(v*u))^(-3/2)*cos(u) # integrand in (2)
        K=quadgk(u->k(u,a,b,v),0,Inf)[1] # numerical quadrature integration of (2)

        #Right-hand-side of Eqn 1 in Hellwarth 1999 // Eqn (4) in Baggio1997
        RHS= α/(3*sqrt(π)) * βred^(5/2) / sinh(βred/2) * (v^3/w^3) * K
        μ=RHS^-1 * (q)/(ω*mb)
        @printf("\n\tμ(Hellwarth1999,b=0)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)
        @printf("\n\tError due to b=0; %f",(100^2*μ-p.Hμ[length(p.Hμ)])/(100^2*μ))
        #append!(Hμs,μ*100^2)
        end
        @printf("\n") # blank line at end of spiel.

        # Recycle previous variation results (v,w) as next guess
        initial=[v,w] # Caution! Might cause weird sticking in local minima

        append!(p.v,v)
        append!(p.w,w)

    end

    return(p)
end

"""
    make_polaron(ϵ_optic, ϵ_static, phonon_freq, m_eff; temp = 300.0, efield_freq = 0.0, volume = nothing, ir_activity = nothing, N_params = 1, rtol = 1e-3, verbose = false, threads = false)

Solves the Feynman polaron problem variationally with finite temperature
Osaka energies. From the resulting v and w parameters, calculates polaron
structure (wave function size, etc.).

Uses FHIP to predict a temperature- and frequency-dependent mobility, complex conductivity and impedance for the material system.

Can evaluate polaron properties for materials with multiple phonon branches using infrared activities and phonon branch frequencies.

# Arguments
- `ϵ_optic`: reduced optical dielectric constant.
- `ϵ_static`: reduced static dielectric constant.
- `phonon_freq`: vector of characteristic dielectric phonon frequencies (THz).
- `m_eff`: bare-band effective-mass (mₑ).
- `Trange`: temperature value or range.
- `Ω`: electric field frequency value or range (THz).
- `volume`: unit cell volume for the material (m³). `nothing` for one phonon mode.
- `ir_activity`: vector of infrared activities. `nothing` for one phonon  mode.
- `N_params`: number of variational parameters to minimise the polaron energy.
- `rtol`: relative tolerance for the accuracy of any quadrature integrations.
- 'verbose': `true` to print progress. `false` to print nothing.
- 'threads'; `true` activates multithreading through use of @tullio macro. Tullio utilises einsum notation.

Returns a structure of type NewPolaron, containing arrays of useful
information.  Also prints a lot of information to the standard out - which
may be more useful if you're just inquiring as to a particular data point,
rather than plotting a temperature-dependent parameter.

As an example, to calculate the electron polaron in MAPI, at temperatures 0:100:400 K and electric field frequencies 0.0:0.1:5.0 THz, and inclusive of 15 optical phonon modes:

# Examples
```jldoctest
make_polaron(
    4.5,
    24.1,
    [4.016471586720514, 3.887605410774121, 3.5313112232401513, 2.755392921480459, 2.4380741812443247, 2.2490917637719408, 2.079632190634424, 2.0336707697261187, 1.5673011873879714, 1.0188379384951798, 1.0022960504442775, 0.9970130778462072, 0.9201781906386209, 0.800604081794174, 0.5738689505255512],
    0.12;
    temp = 0.0:100.0:400.0,
    efield_freq = 0.0:0.1:5.0,
    volume = (6.29e-10)^3,
    ir_activity = [0.08168931020200264, 0.006311654262282101, 0.05353548710183397, 0.021303020776321225, 0.23162784335484837, 0.2622203718355982, 0.23382298607799906, 0.0623239656843172, 0.0367465760261409, 0.0126328938653956, 0.006817361620021601, 0.0103757951973341, 0.01095811116040592, 0.0016830270365341532, 0.00646428491253749],
    N_params = 1)
```

"""
function make_polaron(ϵ_optic, ϵ_static, phonon_freq, m_eff, Trange, Ω; volume = nothing, ir_activity = nothing, rtol = 1e-4, N_params = 1, verbose = false, threads = false)

    # Number of phonon modes.
    N_modes = length(phonon_freq)

    if verbose
        println("Calculating Fröhlich alpha...")
    end

    # Convert THz phonon frequencies to 2π THz.
    ω = 2π .* phonon_freq

    if N_modes == 1
        # One phonon mode.

        # Calculate Frohlich alpha coupling parameters.
        α = frohlichalpha(ϵ_optic, ϵ_static, phonon_freq, m_eff)
    else
        # Multiple phonon modes

        # Calculate contribution to the ionic delectric constant for each phonon mode.
        ϵ_ionic = [ϵ_ionic_mode(i, j, volume) for (i, j) in zip(phonon_freq, ir_activity)]

        # Calculate contribution to Frohlich alpha for each phonon mode.
        α = [multi_frohlichalpha(ϵ_optic, i, sum(ϵ_ionic), j, m_eff) for (i, j) in zip(ϵ_ionic, phonon_freq)]
    end

    if verbose
        show(IOContext(stdout, :limit => true), round.(α, digits = 3))
        print("\n\n")
        println("Calculating thermodynamic betas...")
    end

    # Calculate reduced thermodynamic betas for each phonon mode.
    # If temperature is zero, return Inf.
    @tullio threads = threads betas[m, i] := Trange[i] == 0.0 ? Inf64 : ħ * ω[m] / (kB * Trange[i]) * 1e12 (i in eachindex(Trange))

    if verbose
        show(IOContext(stdout, :limit => true), round.(betas, digits = 3))
        print("\n\n")
        println("Calculating variational parameters...")
        global count = 0
        global processes = length(Trange)
    end

    # Calculate variational parameters for each temperature from multiple phonon frequencies.
    @tullio threads = threads params[i] := Trange[i] == 0.0 ? var_params(α; v = 5.6, w = 3.4, ω = ω, rtol = rtol, N = N_params, verbose = verbose) : var_params(α, betas[:, i]; v = 5.6, w = 3.4, ω = ω, rtol = rtol, T = Trange[i], N = N_params, verbose = verbose) (i in eachindex(Trange))

    # Separate tuples of variational parameters into a list of 'v' and 'w' parameters for each temperature.
    @tullio threads = threads v_params[i] := params[i][1] (i in eachindex(Trange))
    @tullio threads = threads w_params[i] := params[i][2] (i in eachindex(Trange))

    if verbose
        println("\nv: ")
        show(IOContext(stdout, :limit => true), round.(v_params, digits = 3))
        print("\n\n")
        println("w: ")
        show(IOContext(stdout, :limit => true), round.(w_params, digits = 3))
        print("\n\n")
        println("Calculating spring constants...")
    end

    # Calculate fictitious spring constants for each temperature.
    @tullio threads = threads spring_constants[i] := v_params[i]^2 - w_params[i]^2

    if verbose
        show(IOContext(stdout, :limit => true), round.(spring_constants, digits = 3))
        print("\n\n")
        println("Calculating fictitious masses...")
    end

    # Calculate fictitious masses for each temperature.
    @tullio threads = threads masses[i] := spring_constants[i] / w_params[i]^2

    if verbose
        show(IOContext(stdout, :limit => true), round.(masses, digits = 3))
        print("\n\n")
        println("Calculating free energies...")
        global count = 0
        global processes = length(Trange)
    end

    # Calculate ground-state free energies for each temperature.
    @tullio threads = threads energies[i] := Trange[i] == 0.0 ? multi_F(v_params[i], w_params[i], α; ω = ω, rtol = rtol, verbose = verbose) * 1000 * ħ / eV * 1e12 : multi_F(v_params[i], w_params[i], α, betas[:, i]; ω = ω, rtol = rtol, T = Trange[i],  verbose = verbose) * 1000 * ħ / eV * 1e12 (i in eachindex(Trange))

    if verbose
        show(IOContext(stdout, :limit => true), round.(energies, digits = 3))
        print("\n\n")
        println("Calculating complex impedances...")
        global count = 0
        global processes = length(Trange) * length(Ω)
    end

    # Calculate complex impedances for each frequency and temperature. Returns a matrix.
    @tullio threads = threads impedances[f, i] := Ω[f] == Trange[i] == 0.0 ? 0.0 + 1im * 0.0 : polaron_complex_impedence(Ω[f], betas[:, i], α, v_params[i], w_params[i]; ω = ω, rtol = rtol, T = Trange[i], verbose = verbose) / eV^2 * 1e12 * me * m_eff * volume * 100^3 (i in eachindex(Trange), f in eachindex(Ω))

    if verbose
        show(IOContext(stdout, :limit => true), round.(impedances, digits = 3))
        print("\n\n")
        println("Calculating complex conductivities...")
    end

    # Calculate complex conductivities for each frequency and temperature. Returns a matrix.
    @tullio threads = threads conductivities[f, i] := Ω[f] == Trange[i] == 0.0 ? Inf64 + 1im * 0.0 : 1 / impedances[f, i] (i in eachindex(Trange), f in eachindex(Ω))

    if verbose
        show(IOContext(stdout, :limit => true), round.(conductivities, digits = 3))
        print("\n\n")
        println("Calculating mobilities...")
        global count = 0
        global processes = length(Trange)
    end

    # Calculates the dc mobility for each temperature.
    @tullio threads = threads mobilities[i] := Trange[i] == 0.0 ? Inf64 : polaron_mobility(betas[:, i], α, v_params[i], w_params[i]; ω = ω, rtol = rtol, T = Trange[i], verbose = verbose) * eV / (1e12 * me * m_eff) * 100^2 (i in eachindex(Trange))

    if verbose
        show(IOContext(stdout, :limit => true), round.(mobilities, digits = 3))
        print("\n")
    end

    # Return Polaron mutable struct with evaluated data.
    return NewPolaron(α, Trange, betas, phonon_freq, v_params, w_params, spring_constants, masses, energies, Ω, impedances, conductivities, mobilities)
end

"""
    make_polaron(α, Trange, Ω; ω = 1.0, rtol = 1e-4, verbose = false, threads = false)

Same as above but from a model system with specified alpha values rather than from material properties. Here we only have one phonon mode with a normalised frequency `ω = 1.0`.
"""
function make_polaron(α, Trange, Ω; ω = 1.0, rtol = 1e-4, verbose = false, threads = false)

    if verbose
        println("\nCalculating thermodynamic betas...")
    end

    # Calculate reduced thermodynamic betas for each phonon mode.
    # If temperature is zero, return Inf.
    @tullio threads = threads betas[i] := Trange[i] == 0.0 ? Inf64 : ω / Trange[i]

    if verbose
        show(IOContext(stdout, :limit => true), round.(betas, digits = 3))
        print("\n\n")
        println("Calculating variational parameters...")
        global count = 0
        global processes = length(Trange) * length(α)
    end

    # Calculate variational parameters for each alpha parameter and temperature. Returns a Matrix of tuples.
    @tullio threads = threads params[j, i] := Trange[i] == 0.0 ? var_params(α[j]; v = 5.6, w = 3.4, ω = ω, rtol = rtol, verbose = verbose) : var_params(α[j], betas[i]; v = 5.6, w = 3.4, ω = ω, rtol = rtol, T = Trange[i], verbose = verbose) (i in eachindex(Trange), j in eachindex(α))

    # Separate tuples of variational parameters into a Matrices of 'v' and 'w' parameters for each alpha parameter and temperature.
    @tullio threads = threads v_params[j, i] := params[j, i][1]
    @tullio threads = threads w_params[j, i] := params[j, i][2]

    if verbose
        println("\nv: ")
        show(IOContext(stdout, :limit => true), round.(v_params, digits = 3))
        print("\n\n")
        println("w: ")
        show(IOContext(stdout, :limit => true), round.(w_params, digits = 3))
        print("\n\n")
        println("Calculating spring constants...")
    end

    # Calculate fictitious spring constants for each alpha parameter and temperature. Returns a Matrix.
    @tullio threads = threads spring_constants[j, i] := v_params[j, i]^2 - w_params[j, i]^2

    if verbose
        show(IOContext(stdout, :limit => true), round.(spring_constants, digits = 3))
        print("\n\n")
        println("Calculating fictitious masses...")
    end

    # Calculate fictitious masses for each alpha parameter and temperature. Returns a Matrix.
    @tullio threads = threads masses[j, i] := spring_constants[j, i] / w_params[j, i]^2

    if verbose
        show(IOContext(stdout, :limit => true), round.(masses, digits = 3))
        print("\n\n")
        println("Calculating free energies...")
        global count = 0
        global processes = length(T) * length(α)
    end

    # Calculate ground-state free energies for each alpha parameter and temperature. Returns a Matrix.
    @tullio threads = threads energies[j, i] := Trange[i] == 0.0 ? multi_F(v_params[j, i], w_params[j, i], α[j]; ω = ω, rtol = rtol, verbose = verbose) : multi_F(v_params[j, i], w_params[j, i], α[j], betas[i]; ω = ω, rtol = rtol, T = Trange[i],  verbose = verbose) (i in eachindex(Trange), j in eachindex(α))

    if verbose
        show(IOContext(stdout, :limit => true), round.(energies, digits = 3))
        print("\n\n")
        println("Calculating complex impedances...")
        global count = 0
        global processes = length(T) * length(Ω) * length(α)
    end

    # Calculate complex impedances for each alpha parameter, frequency and temperature. Returns a 3D Array.
    @tullio threads = threads impedances[k, j, i] := Ω[j] == Trange[i] == 0.0 ? 0.0 + 1im * 0.0 : polaron_complex_impedence(Ω[j], [betas[i]], α[k], v_params[k, i], w_params[k, i]; ω = ω, rtol = rtol, T = Trange[i],  verbose = verbose) (i in eachindex(Trange), j in eachindex(Ω), k in eachindex(α))

    if verbose
        show(IOContext(stdout, :limit => true), round.(impedances, digits = 3))
        print("\n\n")
        println("Calculating complex conductivities...")
    end

    # Calculate complex conductivities for each alpha parameter, frequency and temperature. Returns a 3D array.
    @tullio threads = threads conductivities[k, j, i] := Ω[j] == Trange[i] == 0.0 ? Inf64 + 1im * 0.0 : 1 / impedances[k, j, i] (i in eachindex(Trange), j in eachindex(Ω), k in eachindex(α))

    if verbose
        show(IOContext(stdout, :limit => true), round.(conductivities, digits = 3))
        print("\n\n")
        println("Calculating mobilities...")
        global count = 0
        global processes = length(T) * length(α)
    end

    # Calculates the dc mobility for each alpha parameter and each temperature.
    @tullio threads = threads mobilities[j, i] := Trange[i] == 0.0 ? Inf64 : polaron_mobility(betas[i], α[j], v_params[j, i], w_params[j, i]; ω = ω, rtol = rtol, T = Trange[i], verbose = verbose) (i in eachindex(Trange), j in eachindex(α))

    if verbose
        show(IOContext(stdout, :limit => true), round.(mobilities, digits = 3))
        print("\n")
    end

    # Return Polaron mutable struct with evaluated data.
    return NewPolaron(α, Trange, betas, ω, v_params, w_params, spring_constants, masses, energies, Ω, impedances, conductivities, mobilities)
end

"""
    save_polaron(p::NewPolaron, prefix)

Saves data from 'polaron' into file "prefix".
This is a .jdl file for storing the polaron data whilst preserving types. Allows for saving multidimensional arrays that sometimes arise in the polaron data.
Each parameter in the NewPolaron type is saved as a dictionary entry. E.g. NewPolaron.α is saved under JLD.load("prefix.jld")["alpha"].
"""
function save_polaron(polaron::NewPolaron, prefix)

    println("Saving polaron data to $prefix.jld ...")

    JLD.save("$prefix.jld",
    "alpha", polaron.α,
    "temperature", polaron.T,
    "beta", polaron.β,
    "phonon freq", polaron.ω,
    "v", polaron.v,
    "w", polaron.w,
    "spring", polaron.κ,
    "mass", polaron.M,
    "energy", polaron.F,
    "efield freq", polaron.Ω,
    "impedance", polaron.Z,
    "conductivity", polaron.σ,
    "mobility", polaron.μ
    )

    println("... Polaron data saved.")
end

"""
    load_polaron(p::NewPolaron, prefix)

Loads data from file "polaron_file_path" into a NewPolaron type.
"""
function load_polaron(polaron_file_path)

    println("Loading polaron data from $polaron_file_path ...")

    data = JLD.load("$polaron_file_path")

    polaron = NewPolaron(
        data["alpha"],
        data["temperature"],
        data["beta"],
        data["phonon freq"],
        data["v"],
        data["w"],
        data["spring"],
        data["mass"],
        data["energy"],
        data["efield freq"],
        data["impedance"],
        data["conductivity"],
        data["mobility"]
    )

    println("... Polaron loaded.")

    return polaron
end

"""
    savepolaron(fileprefix, p::Polaron)

Saves data from polaron 'p' into file "fileprefix".
This is a simple space-delimited text file, with each entry a separate temperature, for plotting with Gnuplot or similar.

Structure of file is written to the header:
# Ts, βreds, Kμs, Hμs, FHIPμs, vs, ws, ks, Ms, As, Bs, Cs, Fs, Taus, rfsis
# 1    2     3    4     5      6   7   8   9  10  11  12  13    14     15
"""
function savepolaron(fileprefix, p::Polaron)
    println("Saving data to $fileprefix.dat ...")
    f=open("$fileprefix.dat","w")

    @printf(f,"# %s\n",fileprefix) # put name / material at header
    @printf(f,"# Params in SI: ω =%g mb=%g \n",p.ω[1] ,p.mb[1])
    @printf(f,"# Alpha parameter: α = %f  \n",p.α[1] )

    @printf(f,"# Ts, βreds, Kμs, Hμs, FHIPμs, vs, ws, ks, Ms, As, Bs, Cs, Fs, Taus, rfsis\n")
    @printf(f,"#  1    2     3    4     5      6   7   8   9  10  11  12  13    14     15\n") # columns for GNUPLOT etc.

    for i in 1:length(p.T)
        @printf(f,"%d %03f %g %g %g %g %g %g %g %g %g %g %g %g %g \n",
        p.T[i], p.βred[i], p.Kμ[i], p.Hμ[i], p.FHIPμ[i],
        p.v[i], p.w[i],
        p.k[i], p.M[i], p.A[i], p.B[i], p.C[i], p.F[i],
        p.Tau[i], p.rfsi[i])
    end
    close(f)
end


text:
<element name="post-item">

        {{ post.title }}
        ({{ post.hostname }})

        {{ post.karma }} points
        by {{ post.user.username }}

            {{ post.fuzzyTime }}

        |

            0 comments

</element>


text:
module TestSIMDIndexing
using Test

@testset "$file" for file in [
        "test_mappedarrays.jl",
        "test_fillarrays.jl",
        "test_ambiguity.jl",
        ]
    include(file)
end

end  # module


text:
# __precompile__()

module MPSWriter

using SparseArrays

# order, indices, weights
const SOS = Tuple{Int, Vector{Int}, Vector{Float64}}

function getrowsense(rowlb::Vector{T1}, rowub::Vector{T2}) where {T1 <: Real, T2 <:Real}
    @assert length(rowlb) == length(rowub)
    row_sense = Array{Symbol}(undef,length(rowub))
    hasranged = false
    for r=1:length(rowlb)
        @assert rowlb[r] <= rowub[r]
    	if (rowlb[r] == -Inf && rowub[r] != Inf) || (rowlb[r] == typemin(eltype(rowlb)) && rowub[r] != typemax(eltype(rowub)))
    		row_sense[r] = :(<=) # LE constraint
    	elseif (rowlb[r] != -Inf && rowub[r] == Inf)  || (rowlb[r] != typemin(eltype(rowlb)) && rowub[r] == typemax(eltype(rowub)))
    		row_sense[r] = :(>=) # GE constraint
    	elseif rowlb[r] == rowub[r]
    		row_sense[r] = :(==) # Eq constraint
        elseif (rowlb[r] == -Inf && rowub[r] == Inf)
            error("Cannot have a constraint with no bounds")
    	else
            row_sense[r] = :ranged
            hasranged = true
    	end
    end
    row_sense, hasranged
end

function writerows!(io::IO, row_sense::Vector{Symbol}, rownames::Vector{String})
    # Objective and constraint names
    println(io, "ROWS\n N  OBJ")
    sensechar = ' '
    for i in 1:length(row_sense)
        if row_sense[i] == :(<=)
            sensechar = 'L'
        elseif row_sense[i] == :(==)
            sensechar = 'E'
        elseif row_sense[i] == :(>=)
            sensechar = 'G'
        elseif row_sense[i] == :ranged
            sensechar = 'E'
        else
            error("Unknown row sense $(row_sense[i])")
        end
        println(io, " $sensechar  $(rownames[i])")
    end
end

const SUPPORTEDVARIABLETYPES = [:Bin, :Int, :Cont, :Fixed]

function writecolumns!(io::IO, A::AbstractMatrix, colcat, c::Vector, sense::Symbol, colnames::Vector{String}, rownames::Vector{String})
    @assert sense == :Min || sense == :Max
    @assert length(colcat) == length(c) == size(A)[2]

    integer_group = false

    println(io, "COLUMNS")

    for col in 1:length(c)
        colty = colcat[col]
        if !(colty in SUPPORTEDVARIABLETYPES)
            error("The MPS file writer does not currently support variables of the type $colty")
        end
        if (colty == :Bin || colty == :Int) && !integer_group
            println(io, "    MARKER    'MARKER'                 'INTORG'")
            integer_group = true
        elseif (colty == :Cont || colty == :Fixed) && integer_group
            println(io, "    MARKER    'MARKER'                 'INTEND'")
            integer_group = false
        end
        inconstraint = writecolumn!(io, A, col, colnames, rownames)
    	if abs(c[col]) > 1e-10 || !inconstraint # Non-zeros
    		# Flip signs for maximisation
            _println(io, "    $(rpad(colnames[col], 8))  $(rpad("OBJ", 8))  ", (sense==:Max ? -1 : 1)*c[col])
	    end
    end
    if integer_group
        println(io, "    MARKER    'MARKER'                 'INTEND'")
    end
end

function _println(io::IO, s::String, x::Number)
    print(io, s)
    print_shortest(io, x)
    println(io)
end

function writecolumn!(io::IO, A::AbstractSparseMatrix{Tv,Ti}, col::Int, colnames::Vector{String}, rownames::Vector{String}) where {Tv,Ti}
    inconstraint = false
    rows, vals = rowvals(A), nonzeros(A)
    for j in nzrange(A, col)
        _println(io, "    $(rpad(colnames[col], 8))  $(rpad(rownames[rows[j]], 8))  ", vals[j])
        inconstraint = true
    end
    inconstraint
end

function writecolumn!(io::IO, A::Array{T, 2}, col::Int, colnames::Vector{String}, rownames::Vector{String}) where {T}
    inconstraint = false
    for row in 1:size(A)[1]
        if abs(A[row, col]) > 1e-10 # Non-zero
            _println(io, "    $(rpad(colnames[col], 8))  $(rpad(rownames[row], 8))  ", A[row, col])
            inconstraint = true
        end
    end
    inconstraint
end

function writerhs!(io::IO, rowlb, rowub, row_sense::Vector{Symbol}, rownames::Vector{String})
    @assert length(rowlb) == length(rowub) == length(row_sense)
    println(io, "RHS")
    for row in 1:length(rowlb)
        _println(io, "    rhs       $(rpad(rownames[row], 8))  ", row_sense[row] == :(<=) ? rowub[row] : rowlb[row])
    end
end

function writeranges!(io::IO, rowlb, rowub, row_sense::Vector{Symbol}, rownames::Vector{String})
    @assert length(rowlb) == length(rowub) == length(row_sense)
    println(io, "RANGES")
    for row=1:length(row_sense)
        if row_sense[row] == :ranged
            _println(io, "    rhs       $(rpad(rownames[row], 8))  ", rowub[row] - rowlb[row])
        end
    end
end

function writebounds!(io::IO, collb, colub, colnames::Vector{String})
    @assert length(collb) == length(colub)
    println(io, "BOUNDS")
    for col in 1:length(collb)
        if colub[col] == Inf
            if collb[col] == -Inf
                writebound!(io, "FR", col, colnames[col])
                continue
            else
                writebound!(io, "PL", col, colnames[col])
            end
        else
            writebound!(io, "UP", col, colub[col], colnames[col])
        end
        if collb[col] == -Inf
            writebound!(io, "MI", col, colnames[col])
        elseif collb[col] != 0
            writebound!(io, "LO", col, collb[col], colnames[col])
        end
    end
end

function writebound!(io::IO, ty::String, vidx::Int, colname::String)
    @assert ty in ["FR", "MI", "PL"]
    println(io, " $ty BOUNDS    $colname")
end

function writebound!(io::IO, ty::String, vidx::Int, val::Real, colname::String)
    @assert ty in ["LO", "UP"]
    _println(io, " $ty BOUNDS    $(rpad(colname, 8))  ", val)
end

function writesos!(io::IO, sos::Vector{SOS}, maxvarindex::Int, colnames::Vector{String})
    println(io, "SOS")
    for i=1:length(sos)
        order   = sos[i][1]
        indices = sos[i][2]
        weights = sos[i][3]
        @assert length(indices) == length(weights)
        println(io, " S$(order) SOS$i")
        for (idx, weight) in zip(indices, weights)
            @assert idx > 0 && idx <= maxvarindex
            _println(io, "    $(rpad(colnames[idx], 8))  ", weight)
        end
    end
end

function writequad!(io::IO, Q::AbstractSparseMatrix{Tv,Ti}, sense::Symbol, colnames::Vector{String}) where {Tv,Ti}
    @assert sense == :Min || sense == :Max
    println(io, "QMATRIX")
    rows = rowvals(Q)
    vals = nonzeros(Q)
    if sense == :Max
        vals *= -1
    end
    for i = 1:size(Q)[2]
        for j in nzrange(Q, i)
            _println(io, "    $(rpad(colnames[rows[j]],8)) $(rpad(colnames[i],8))  ", vals[j])
        end
    end
end

function writequad!(io::IO, Q::AbstractMatrix, sense::Symbol, colnames::Vector{String})
    @assert sense == :Min || sense == :Max
    sgn = (sense == :Max ? -1 : 1)
    println(io, "QMATRIX")
    for i = 1:size(Q)[2]
        for j in 1:size(Q)[1]
            if abs(Q[j, i]) > 1e-10
                _println(io, "    $(rpad(colnames[j],8)) $(rpad(colnames[i],8))  ", sgn*Q[j, i])
            end
        end
    end
end

function write(io::IO,
    A::AbstractMatrix,       # the constraint matrix
    collb::Vector,  # vector of variable lower bounds
    colub::Vector,  # vector of variable upper bounds
    c::Vector,      # vector containing variable objective coefficients
    rowlb::Vector,  # constraint lower bounds
    rowub::Vector,  # constraint upper bounds
    sense::Symbol,           # model sense
    colcat::Vector{Symbol},  # constraint types
    sos::Vector{SOS},        # SOS information
    Q::AbstractMatrix,      #  Quadratic objectives 0.5 * x' Q x
    modelname::AbstractString ="MPSWriter_jl",  # MPS model name
    colnames::Vector{String}  = ["V$i" for i in 1:length(c)],
    rownames::Vector{String}  = ["C$i" for i in 1:length(rowub)]
)
    # Max width (8) for names in fixed MPS format
    # TODO: replace names by alpha-numeric
    #   to enable more variables and constraints.
    #   Although to be honest no one should be using
    #   this with that many constraints/variables
    # @assert length(c) <= 9999999
    # @assert length(rowub) <= 9999999

    # Sanity checks
    @assert length(rowlb) == length(rowub)
    @assert length(collb) == length(colub) == length(c) == length(colcat)
    @assert sense == :Min || sense == :Max
    @assert size(Q) == (0,0) || size(Q) == (length(collb), length(collb))

    println(io, "NAME          $modelname")
    row_sense, hasranged = getrowsense(rowlb, rowub)
    writerows!(io, row_sense, rownames)
    writecolumns!(io, A, colcat, c, sense, colnames, rownames)
    writerhs!(io, rowlb, rowub, row_sense, rownames)
    if hasranged
        writeranges!(io, rowlb, rowub, row_sense, rownames)
    end
    writebounds!(io, collb, colub, colnames)
    if length(sos) > 0
        writesos!(io, sos, length(collb), colnames)
    end
    if length(Q) > 0
        writequad!(io, Q, sense, colnames)
    end
    println(io, "ENDATA")
end

end # module


text:
using BandedMatrices, FillArrays, Test
import LinearAlgebra: axpy!
import LazyArrays: DenseColumnMajor
import BandedMatrices: BandedColumns, bandeddata

@testset "BandedMatrix SubArray" begin
    @testset "BandedMatrix SubArray interface" begin
        A = brand(10,10,1,2)
        V = view(A,2:4,3:6)
        @test isbanded(V)
        @test bandwidths(V) == (2,1)
        @test BandedMatrices.MemoryLayout(V) == BandedColumns(DenseColumnMajor())
        @test all(Matrix(BandedMatrices._BandedMatrix(bandeddata(V), size(V,1), bandwidths(V)...)) .===
                Matrix(V))
    end

    @testset "BandedMatrix SubArray arithmetic" begin
        let A = brand(10,10,1,2), B = brand(20,20,1,2)
            @test isa(view(B,1:10,1:10),BandedMatrices.BandedSubBandedMatrix{Float64})

            B2 = copy(B)
            @test (2.0A+B[1:10,1:10]) ≈ axpy!(2.0,A,view(B2,1:10,1:10))
            @test (2.0A+B[1:10,1:10]) ≈ B2[1:10,1:10]

            A2 = copy(A)
            @test (2.0B[1:10,1:10]+A) ≈ axpy!(2.0,view(B,1:10,1:10),A2)
            @test (2.0B[1:10,1:10]+A) ≈ A2
        end

        let A = brand(20,20,1,2), B = brand(20,20,1,2)
            @test isa(view(B,1:10,1:10),BandedMatrices.BandedSubBandedMatrix{Float64})

            B2 = copy(B)
            @test (2.0A[1:10,1:10]+B[1:10,1:10]) ≈ axpy!(2.0,view(A,1:10,1:10),view(B2,1:10,1:10))
            @test (2.0A[1:10,1:10]+B[1:10,1:10]) ≈ B2[1:10,1:10]

            B2 = copy(B)
            @test (2.0A[1:10,:]+B[1:10,:]) ≈ axpy!(2.0,view(A,1:10,:),view(B2,1:10,:))
            @test (2.0A[1:10,:]+B[1:10,:]) ≈ B2[1:10,:]

            B2 = copy(B)
            @test (2.0A[:,1:10]+B[:,1:10]) ≈ axpy!(2.0,view(A,:,1:10),view(B2,:,1:10))
            @test (2.0A[:,1:10]+B[:,1:10]) ≈ B2[:,1:10]

            B2 = copy(B)
            @test (2.0A[:,:]+B[:,:]) ≈ axpy!(2.0,view(A,:,:),view(B2,:,:))
            @test (2.0A[:,:]+B[:,:]) ≈ B2[:,:]
        end

        let A = brand(10,10,1,2), B = brand(20,10,1,2)
            B2 = copy(B)
            @test (2.0A+B[1:10,1:10]) ≈ axpy!(2.0,A,view(B2,1:10,:))
            @test (2.0A+B[1:10,1:10]) ≈ B2[1:10,1:10]

            A2 = copy(A)
            @test (2.0B[1:10,1:10]+A) ≈ axpy!(2.0,view(B,1:10,:),A2)
            @test (2.0B[1:10,1:10]+A) ≈ A2
        end

        let A = brand(10,10,1,2), B = brand(10,20,1,2)
            B2 = copy(B)
            @test (2.0A+B[1:10,1:10]) ≈ axpy!(2.0,A,view(B2,:,1:10))
            @test (2.0A+B[1:10,1:10]) ≈ B2[1:10,1:10]

            A2 = copy(A)
            @test (2.0B[1:10,1:10]+A) ≈ axpy!(2.0,view(B,:,1:10),A2)
            @test (2.0B[1:10,1:10]+A) ≈ A2
        end
    end

    @testset "BandedMatrix SubArray *" begin
        let A=brand(10,10,1,2), v=rand(20)
            @test A*view(v,1:10) ≈ Matrix(A)*v[1:10]
            @test A*view(v,1:2:20) ≈ Matrix(A)*v[1:2:20]

            M=rand(20,20)
            @test A*view(M,1:10,1:10) ≈ Matrix(A)*M[1:10,1:10]
            @test A*view(M,1:2:20,1:2:20) ≈ Matrix(A)*M[1:2:20,1:2:20]
        end

        let A=brand(10,10,1,2), B=brand(20,10,1,2), C=brand(10,20,1,2),
            D=brand(20,20,1,2), M=rand(10,10), V=view(rand(20,20),1:2:20,1:2:20)

            for S in (view(A,:,:),view(B,1:10,:),view(C,:,1:10),view(D,1:10,1:10),
                        view(D,2:11,1:10),view(D,1:10,2:11),view(D,11:20,11:20))
                @test A*S ≈ Matrix(A)*Matrix(S)
                @test S*A ≈ Matrix(S)*Matrix(A)
                @test S*S ≈ Matrix(S)*Matrix(S)

                @test M*S ≈ M*Matrix(S)
                @test S*M ≈ Matrix(S)*M

                @test V*S ≈ Matrix(V)*Matrix(S)
                @test S*V ≈ Matrix(S)*Matrix(V)
            end
        end
    end

    @testset "BandedMatrix SubArray conversion" begin
        @testset "Container $C" for C in (Matrix{Int32}, MyMatrix{Int32})
            matrix = convert(C, rand(Int32, 10, 12))
            T = eltype(matrix)

            banded = @inferred BandedMatrix{T, C}(matrix, (1, 1))
            indices = 2:8, 3:9, 5:9
            @testset "view $kr, $jr" for kr in indices, jr in indices
                bview = @inferred view(banded, kr, jr)
                @test bview isa BandedMatrices.BandedSubBandedMatrix{T, C}
                @test @inferred(convert(BandedMatrix, bview)) isa BandedMatrix{T, C}
                @test convert(BandedMatrix, bview) == banded[kr, jr]
            end
        end
    end
end


text:
module ArraysBenchs

include("KernelsBenchs.jl")

include("ApplyBenchs.jl")

end # module


text:
################################################################
"""
    @snoopi_bot config::BotConfig snoopScript

macro that generates precompile files and includes them in the package. Calls other bot functions.

# Examples

`@snoopi_bot` the examples that call the package functions.

```julia
using SnoopCompile

@snoopi_bot "MatLang" begin
  using MatLang
  examplePath = joinpath(dirname(dirname(pathof(MatLang))), "examples")
  include(joinpath(examplePath,"Language_Fundamentals", "usage_Entering_Commands.jl"))
  include(joinpath(examplePath,"Language_Fundamentals", "usage_Matrices_and_Arrays.jl"))
  include(joinpath(examplePath,"Language_Fundamentals", "Data_Types", "usage_Numeric_Types.jl"))
end
```
"""
macro snoopi_bot(config::BotConfig, snoopScript)

    packageName = config.packageName
    blacklist = config.blacklist
    subst = config.subst
    ################################################################
    packagePath = joinpath(pwd(),"src","$packageName.jl")
    precompilePath, precompileFolder = precompile_pather(packageName)

    quote
        packageSym = Symbol($packageName)
        ################################################################
        using SnoopCompile
        ################################################################
        precompile_deactivator($packagePath, $precompilePath);
        ################################################################

        ### Log the compiles
        data = @snoopi begin
            $(esc(snoopScript))
        end

        ################################################################
        ### Parse the compiles and generate precompilation scripts
        pc = SnoopCompile.parcel(data, subst = $subst, blacklist = $blacklist)
        onlypackage = Dict( packageSym => sort(pc[packageSym]) )
        SnoopCompile.write($precompileFolder,onlypackage)
        ################################################################
        precompile_activator($packagePath, $precompilePath)
    end

end

macro snoopi_bot(packageName::String, snoopScript)
    config = BotConfig(packageName)
    return quote
        @snoopi_bot $config $(esc(snoopScript))
    end
end
macro snoopi_bot(configExpr, snoopScript)
    config = eval(configExpr)
    return quote
        @snoopi_bot $config $(esc(snoopScript))
    end
end

"""
    @snoopi_bot config::BotConfig

If you do not have additional examples, you can use your runtests.jl file. To do that use:

```julia
using SnoopCompile

# using runtests:
@snoopi_bot "MatLang"
```
"""
macro snoopi_bot(config::BotConfig)

    packageName = config.packageName

    package = Symbol(packageName)
    snoopScript = esc(quote
        using $(package)
        runtestpath = joinpath(dirname(dirname(pathof( $package ))), "test", "runtests.jl")
        include(runtestpath)
    end)
    return quote
        @snoopi_bot $config $(esc(snoopScript))
    end
end

macro snoopi_bot(packageName::String)
    config = BotConfig(packageName)
    return quote
        @snoopi_bot $config
    end
end

macro snoopi_bot(configExpr)
    config = eval(configExpr)
    return quote
        @snoopi_bot $config
    end
end


text:
function image_widget(fn)
    img = Gtk.GtkImageLeaf(fn)
    vbox = Gtk.GtkBoxLeaf(:v)
    push!(vbox, Gtk.GtkLabelLeaf(fn))
    push!(vbox, img)
    show(img)
    vbox
end

function replace_refimg(tmpfn, reffn)
    try
        mkpath(reffn)
    end
    cp(tmpfn, reffn, remove_destination = true)
    info("Replaced reference image $reffn with $tmpfn")
end

"Show a Gtk popup with both images and a confirmation whether we should replace the new image with the old one"
function replace_refimg_dialog(tmpfn, reffn_old; reffn_new = reffn_old)
    # add the images
    imgbox = Gtk.GtkBoxLeaf(:h)
    push!(imgbox, image_widget(tmpfn))
    push!(imgbox, image_widget(reffn_old))

    win = Gtk.GtkWindowLeaf("Should we make this the new reference image?")
    push!(win, Gtk.GtkFrameLeaf(imgbox))

    showall(win)

    # now ask the question
    if Gtk.ask_dialog("Should we make this the new reference image?", "No", "Yes")
        replace_refimg(tmpfn, reffn_new)
    end

    destroy(win)
end


text:
function gersten_alphabet(n::Integer; commutative::Bool = true)
    indexing = [(i, j) for i in 1:n for j in 1:n if i ≠ j]
    S = [ϱ(i, j) for (i, j) in indexing]

    if !commutative
        append!(S, [λ(i, j) for (i, j) in indexing])
    end

    return Alphabet(S)
end

function _commutation_rule(
    ::Type{W},
    A::Alphabet,
    x::S,
    y::S,
) where {S,T,W<:AbstractWord{T}}
    return W(T[A[x], A[y]]) => W(T[A[y], A[x]])
end

function _pentagonal_rule(
    ::Type{W},
    A::Alphabet,
    x::S,
    y::S,
    z::S,
) where {S,T,W<:AbstractWord{T}}
    # x·y·x⁻¹·y⁻¹ => z, i.e. z·y·x => x·y
    return W(T[A[z], A[y], A[x]]) => W(T[A[x], A[y]])
end
function _hexagonal_rule(
    ::Type{W},
    A::Alphabet,
    x::S,
    y::S,
    z::S,
    w::S,
) where {S,T,W<:AbstractWord{T}}
    # x·y⁻¹·z => z·w⁻¹·x
    return W(T[A[x], A[inv(y)], A[z]]) => W(T[A[z], A[w^-1], A[x]])
end

gersten_relations(n::Integer; commutative) =
    gersten_relations(Word{UInt8}, n, commutative = commutative)

function gersten_relations(::Type{W}, n::Integer; commutative) where {W<:AbstractWord}
    @assert n > 1 "Gersten relations are defined only for n>1, got n=$n"
    A = gersten_alphabet(n, commutative = commutative)
    @assert length(A) <= KnuthBendix._max_alphabet_length(W) "Type $W can not represent words over alphabet with $(length(A)) letters."

    rels = Pair{W,W}[]

    for (i, j, k, l) in Iterators.product(1:n, 1:n, 1:n, 1:n)
        if i ≠ j && k ≠ l && k ≠ i && k ≠ j && l ≠ i
            push!(rels, _commutation_rule(W, A, ϱ(i, j), ϱ(k, l)))

            commutative && continue

            push!(rels, _commutation_rule(W, A, λ(i, j), λ(k, l)))
        end
    end

    if !commutative
        for (i, j, k, l) in Iterators.product(1:n, 1:n, 1:n, 1:n)
            if (i ≠ j && k ≠ l && k ≠ j && l ≠ i)
                push!(rels, _commutation_rule(W, A, ϱ(i, j), λ(k, l)))
                push!(rels, _commutation_rule(W, A, λ(i, j), ϱ(k, l)))
            end
        end
    end

    # pentagonal rule:
    # x*y*inv(x)*inv(y)=>z

    for (i, j, k) in Iterators.product(1:n, 1:n, 1:n)
        if (i ≠ j && k ≠ i && k ≠ j)
            push!(rels, _pentagonal_rule(W, A, ϱ(i, j)^-1, ϱ(j, k)^-1, ϱ(i, k)^-1))
            push!(rels, _pentagonal_rule(W, A, ϱ(i, j)^-1, ϱ(j, k), ϱ(i, k)))

            commutative && continue

            push!(rels, _pentagonal_rule(W, A, ϱ(i, j), λ(j, k), ϱ(i, k)^-1))
            push!(rels, _pentagonal_rule(W, A, ϱ(i, j), λ(j, k)^-1, ϱ(i, k)))

            # the same as above, but with ϱ ↔ λ:
            push!(rels, _pentagonal_rule(W, A, λ(i, j)^-1, λ(j, k)^-1, λ(i, k)^-1))
            push!(rels, _pentagonal_rule(W, A, λ(i, j)^-1, λ(j, k), λ(i, k)))

            push!(rels, _pentagonal_rule(W, A, λ(i, j), ϱ(j, k), λ(i, k)^-1))
            push!(rels, _pentagonal_rule(W, A, λ(i, j), ϱ(j, k)^-1, λ(i, k)))
        end
    end

    if !commutative
        for (i, j) in Iterators.product(1:n, 1:n)
            if i ≠ j
                push!(rels, _hexagonal_rule(W, A, ϱ(i, j), ϱ(j, i), λ(i, j), λ(j, i)))
                w = W([A[ϱ(i, j)], A[ϱ(j, i)^-1], A[λ(i, j)]])
                push!(rels, w^2 => inv(A, w)^2)
            end
        end
    end

    return A, rels
end


text:
# using Pipelines

function julia_program_warn(p::JuliaProgram)
    if nthreads() == 1
        @warn "Submitting a JuliaProgram with 1-threaded Julia session is not recommended because it might block schedulers. Starting Julia with multi-threads is suggested. Help: https://docs.julialang.org/en/v1/manual/multi-threading/#Starting-Julia-with-multiple-threads" maxlog=1
    end
end
julia_program_warn(p::CmdProgram) = nothing

"""
    Job(p::Program; kwargs...)
    Job(p::Program, inputs; kwargs...)
    Job(p::Program, inputs, outputs; kwargs...)

Create `Job` by using `Program` from Pipelines.jl package.
The 3 methods are a wrapper around `run(::Program, ...)`.

`kwargs...` include keyword arguments of `Job(::BaseAbstractCmd, ...)` and `run(::Program, ...)`.
"""
function Job(p::Program;
    name::AbstractString = p.name,
    user::AbstractString = "",
    ncpu::Int64 = 1,
    mem::Int64 = 0,
    schedule_time::Union{DateTime,Period} = DateTime(0),
    wall_time::Period = Week(1),
    priority::Int = 20,
    dependency = Vector{Pair{Symbol,Int64}}(),
    stdout=nothing,
    stderr=nothing,
    kwargs...
)
    julia_program_warn(p)
    task = @task run(p; stdout=stdout, stderr=stderr, kwargs...)
    stdout_file = format_stdxxx_file(stdout)
    stderr_file = format_stdxxx_file(stderr)

    Job(generate_id(), name, user, ncpu, mem, schedule_time, DateTime(0), DateTime(0), DateTime(0), wall_time, QUEUING, priority, dependency, task, stdout_file, stderr_file)
end

function Job(p::Program, inputs;
    name::AbstractString = p.name,
    user::AbstractString = "",
    ncpu::Int64 = 1,
    mem::Int64 = 0,
    schedule_time::Union{DateTime,Period} = DateTime(0),
    wall_time::Period = Week(1),
    priority::Int = 20,
    dependency = Vector{Pair{Symbol,Int64}}(),
    stdout=nothing,
    stderr=nothing,
    kwargs...
)
    julia_program_warn(p)
    task = @task run(p, inputs; stdout=stdout, stderr=stderr, kwargs...)
    stdout_file = format_stdxxx_file(stdout)
    stderr_file = format_stdxxx_file(stderr)

    Job(generate_id(), name, user, ncpu, mem, schedule_time, DateTime(0), DateTime(0), DateTime(0), wall_time, QUEUING, priority, dependency, task, stdout_file, stderr_file)
end

function Job(p::Program, inputs, outputs;
    name::AbstractString = p.name,
    user::AbstractString = "",
    ncpu::Int64 = 1,
    mem::Int64 = 0,
    schedule_time::Union{DateTime,Period} = DateTime(0),
    wall_time::Period = Week(1),
    priority::Int = 20,
    dependency = Vector{Pair{Symbol,Int64}}(),
    stdout=nothing,
    stderr=nothing,
    kwargs...
)
    julia_program_warn(p)
    task = @task run(p, inputs, outputs; stdout=stdout, stderr=stderr, kwargs...)
    stdout_file = format_stdxxx_file(stdout)
    stderr_file = format_stdxxx_file(stderr)

    Job(generate_id(), name, user, ncpu, mem, schedule_time, DateTime(0), DateTime(0), DateTime(0), wall_time, QUEUING, priority, dependency, task, stdout_file, stderr_file)
end

program_close_io = JuliaProgram(
    name = "Close Julia IO",
    id_file = ".close-julia-io",
    inputs = "io" => IO,
    main = (inputs, outputs) -> begin
        close(inputs["io"])
        Dict{String,Any}()
    end
)

"""
    close_in_future(io::IO, job::Job)
    close_in_future(io::IO, jobs::Vector{Job})

Close `io` after `job` is in PAST status (either DONE/FAILED/CANCELLED). It is userful if jobs uses `::IO` as `stdout` or `stderr`, because the program will not close `::IO` manually!
"""
function close_in_future(io::IO, job::Job)
    close_job = Job(program_close_io, "io" => io, touch_run_id_file=false, dependency = [PAST => job.id])
    submit!(close_job)
    close_job
end
function close_in_future(io::IO, jobs::Vector{Job})
    deps = [PAST => x.id for x in jobs]
    close_job = Job(program_close_io, "io" => io, touch_run_id_file=false, dependency = deps)
    submit!(close_job)
    close_job
end


text:
module SearchBenchmarks

using Plots
import Random

include("base.jl") # basic types and utilities for benchmarking
export all_search_functions, Benchmark

# it follows the list of all search implementations.
# each file defines a search function with two methods:
# search initialization and search kernel (using initialization result).
include("naivesearch.jl")
export naivesearch
include("bloom_v0.jl")
include("bloom_v0i.jl")
include("bloom_v1.jl")
include("bloom_v2.jl")
include("bloom_v3.jl")
include("bloom_v4.jl")
export bloom_v0, bloom_v0i, bloom_v1, bloom_v2, bloom_v3, bloom_v4
include("bloom_best.jl")
include("bloom_best2.jl")
include("bloom_best3.jl")
include("bloom_bestc.jl")
export bloom_best, bloom_best2, bloom_best3,bloom_bestc
include("bmsearch.jl")
export bmsearch
#add more search implementations...

include("benchmark.jl") # do and visualize benchmark
export benchmark

end # module


text:
using Quiqbox
using Test
using Quiqbox
using Suppressor: @suppress_out

@testset "Optimization.jl" begin

errorThreshold1 = 1e-10
errorThreshold2 = 1e-4

# Floating basis set
nucCoords = [[-0.7,0.0,0.0], [0.7,0.0,0.0]]
nuc = ["H", "H"]

ecMethods = [Quiqbox.defaultECmethod(:RHF, nuc, nucCoords),
             Quiqbox.defaultECmethod(:UHF, nuc, nucCoords)]

E_t1s = [-1.775682556228, -1.775682514824]
# X₁, X₂, Y₁, Y₂, Z₁, Z₂, α₁, α₂, d₁, d₂
par_t1s = [[-0.686629461207,  0.686629461207,  0.0,             0.0,
             0.0,             0.0,             1.331636672728,  0.311858635704,
             0.454798446619,  0.662643931788],
           [-0.686629463257,  0.686629463257,  0.0,             0.0,
             0.0,             0.0,             1.331636672181,  0.311858633503,
             0.454798445732,  0.662643932396]]

grad_t1s = [[-0.056475936092,  0.056475936092,  0.0,             0.0,
              0.0,             0.0,            -0.125188275118,  0.017527948867,
             -0.107797228786,  0.073985454102],
            [-0.056475926135,  0.056475926135,  0.0,             0.0,
              0.0,             0.0,            -0.125188271319,  0.017527947101,
             -0.107797224221,  0.073985450757]]

for (ecMethod, E_t1, par_t1, grad_t1) in zip(ecMethods, E_t1s, par_t1s, grad_t1s)

    bfSource1 = genBasisFunc(missing, ("STO-2G", "H"))[]
    gfs1 = bfSource1.gauss |> collect
    cens = makeCenter.(nucCoords)
    bs1 = genBasisFunc.(cens, Ref(gfs1))
    pars1 = uniqueParams!(bs1, filterMapping=true)

    local Es1L, pars1L, grads1L
    @suppress_out begin
        Es1L, pars1L, grads1L = optimizeParams!(bs1, pars1, nuc, nucCoords, ecMethod,
                                                maxSteps = 200)
    end

    @test isapprox(Es1L[end], E_t1, atol=errorThreshold1)
    @test isapprox(pars1L[end, :], par_t1, atol=errorThreshold1)
    @test isapprox(grads1L[end, :], grad_t1, atol=errorThreshold1)

end

# Grid-based basis set
grid = GridBox(1, 3.0)
gf2 = GaussFunc(0.7,1)
bs2 = genBasisFunc.(grid.box, Ref([gf2]))

pars2 = uniqueParams!(bs2, filterMapping=true)[[1,4]]

local Es2L, pars2L, grads2L
@suppress_out begin
    Es2L, pars2L, grads2L = optimizeParams!(bs2, pars2, nuc, nucCoords, maxSteps = 100)
end

E_t2 = -1.16666630
# L, α
par_t2  = [2.90239973, 0.17642659]
grad_t2 = [0.48936670, 0.10233850]

@test isapprox(Es2L[end], E_t2, atol=errorThreshold2)
@test isapprox(pars2L[end, :], par_t2, atol=errorThreshold2)
@test isapprox(grads2L[end, :], grad_t2, atol=errorThreshold2)

end


text:
export VoronoiDiagram2

export bounded_faces,
       bounded_halfedges,
       ccb,
       ccb_halfedges,
       degree,
       down,
       dual,
       face,
       faces,
       finite_edges,
       halfedge,
       halfedges,
       has_source,
       has_target,
       incident_halfedges,
       incident_halfedges,
       is_bisector,
       is_halfedge_on_ccb,
       is_incident_edge,
       is_incident_face,
       is_ray,
       is_segment,
       is_unbounded,
       is_valid,
       is_valid,
       left,
       locate,
       next,
       number_of_connected_components,
       number_of_faces,
       number_of_halfedges,
       number_of_vertices,
       opposite,
       previous,
       right,
       site,
       sites,
       source,
       swap,
       target,
       twin,
       unbounded_faces,
       unbounded_halfedges,
       up,
       vertices

VoronoiDiagram2{T}(ps::AbstractVector) where T =
    VoronoiDiagram2(collect(CxxRef{_triangulation_point(T)}, CxxRef.(ps)))
VoronoiDiagram2{T}(ps...) where T = VoronoiDiagram2(collect(ps))

@cxxdereference Base.insert!(vd::VoronoiDiagram2, ps::AbstractVector) =
    insert!(vd, collect(CxxRef{_triangulation_point(typeof(vd).parameters[1])},
                        CxxRef.(ps)))
@cxxdereference finite_edges(vd::VoronoiDiagram2) =
    let dg = dual(vd)
        dual.(Ref(dg), edges(dg))
    end