text:
int_rules_1_3_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.3.3*Miscellaneous*algebraic*functions =#
    @apply_utils Antiderivative(~u / (~(e') * sqrt(~(a') + ~(b') * ~x) + ~(f') * sqrt(~(c') + ~(d') * ~x)), ~x) => (~c / (~e * (~b * ~c - ~a * ~d))) * Antiderivative((~u * sqrt(~a + ~b * ~x)) / ~x, ~x) - (~a / (~f * (~b * ~c - ~a * ~d))) * Antiderivative((~u * sqrt(~c + ~d * ~x)) / ~x, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ(~a * (~e) ^ 2 - ~c * (~f) ^ 2, 0))
    @apply_utils Antiderivative(~u / (~(e') * sqrt(~(a') + ~(b') * ~x) + ~(f') * sqrt(~(c') + ~(d') * ~x)), ~x) => (-(~d) / (~e * (~b * ~c - ~a * ~d))) * Antiderivative(~u * sqrt(~a + ~b * ~x), ~x) + (~b / (~f * (~b * ~c - ~a * ~d))) * Antiderivative(~u * sqrt(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ(~b * (~e) ^ 2 - ~d * (~f) ^ 2, 0))
    @apply_utils Antiderivative(~u / (~(e') * sqrt(~(a') + ~(b') * ~x) + ~(f') * sqrt(~(c') + ~(d') * ~x)), ~x) => ~e * Antiderivative((~u * sqrt(~a + ~b * ~x)) / ((~a * (~e) ^ 2 - ~c * (~f) ^ 2) + (~b * (~e) ^ 2 - ~d * (~f) ^ 2) * ~x), ~x) - ~f * Antiderivative((~u * sqrt(~c + ~d * ~x)) / ((~a * (~e) ^ 2 - ~c * (~f) ^ 2) + (~b * (~e) ^ 2 - ~d * (~f) ^ 2) * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~a * (~e) ^ 2 - ~c * (~f) ^ 2, 0) && NeQ(~b * (~e) ^ 2 - ~d * (~f) ^ 2, 0))
    @apply_utils Antiderivative(~(u') / (~(d') * (~x) ^ ~(n') + ~(c') * sqrt(~(a') + ~(b') * (~x) ^ ~(p'))), ~x) => (-(~b) / (~a * ~d)) * Antiderivative(~u * (~x) ^ ~n, ~x) + (1 / (~a * ~c)) * Antiderivative(~u * sqrt(~a + ~b * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (EqQ(~p, 2 * ~n) && EqQ(~b * (~c) ^ 2 - (~d) ^ 2, 0))
    @apply_utils Antiderivative((~x) ^ ~(m') / (~(d') * (~x) ^ ~(n') + ~(c') * sqrt(~(a') + ~(b') * (~x) ^ ~(p'))), ~x) => -(~d) * Antiderivative((~x) ^ (~m + ~n) / (~a * (~c) ^ 2 + (~b * (~c) ^ 2 - (~d) ^ 2) * (~x) ^ (2 * ~n)), ~x) + ~c * Antiderivative(((~x) ^ ~m * sqrt(~a + ~b * (~x) ^ (2 * ~n))) / (~a * (~c) ^ 2 + (~b * (~c) ^ 2 - (~d) ^ 2) * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (EqQ(~p, 2 * ~n) && NeQ(~b * (~c) ^ 2 - (~d) ^ 2, 0))
    @apply_utils Antiderivative(1 / ((~a + ~(b') * (~x) ^ 3) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([r = Numerator(Rt(~a / ~b, 3)), s = Denominator(Rt(~a / ~b, 3))], (r / (3 * ~a)) * Antiderivative(1 / ((r + s * ~x) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) + (r / (3 * ~a)) * Antiderivative((2r - s * ~x) / (((r ^ 2 - r * s * ~x) + s ^ 2 * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~d, ~e, ~f], ~x) && PosQ(~a / ~b)
    @apply_utils Antiderivative(1 / ((~a + ~(b') * (~x) ^ 3) * sqrt(~(d') + ~(f') * (~x) ^ 2)), ~x) => With([r = Numerator(Rt(~a / ~b, 3)), s = Denominator(Rt(~a / ~b, 3))], (r / (3 * ~a)) * Antiderivative(1 / ((r + s * ~x) * sqrt(~d + ~f * (~x) ^ 2)), ~x) + (r / (3 * ~a)) * Antiderivative((2r - s * ~x) / (((r ^ 2 - r * s * ~x) + s ^ 2 * (~x) ^ 2) * sqrt(~d + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~d, ~f], ~x) && PosQ(~a / ~b)
    @apply_utils Antiderivative(1 / ((~a + ~(b') * (~x) ^ 3) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([r = Numerator(Rt(-(~a) / ~b, 3)), s = Denominator(Rt(-(~a) / ~b, 3))], (r / (3 * ~a)) * Antiderivative(1 / ((r - s * ~x) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) + (r / (3 * ~a)) * Antiderivative((2r + s * ~x) / ((r ^ 2 + r * s * ~x + s ^ 2 * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~d, ~e, ~f], ~x) && NegQ(~a / ~b)
    @apply_utils Antiderivative(1 / ((~a + ~(b') * (~x) ^ 3) * sqrt(~(d') + ~(f') * (~x) ^ 2)), ~x) => With([r = Numerator(Rt(-(~a) / ~b, 3)), s = Denominator(Rt(-(~a) / ~b, 3))], (r / (3 * ~a)) * Antiderivative(1 / ((r - s * ~x) * sqrt(~d + ~f * (~x) ^ 2)), ~x) + (r / (3 * ~a)) * Antiderivative((2r + s * ~x) / ((r ^ 2 + r * s * ~x + s ^ 2 * (~x) ^ 2) * sqrt(~d + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~d, ~f], ~x) && NegQ(~a / ~b)
    @apply_utils Antiderivative((~u * (~A + ~(B') * (~x) ^ 4)) / sqrt(~v), ~x) => With([a = Coeff(~v, ~x, 0), b = Coeff(~v, ~x, 2), c = Coeff(~v, ~x, 4), d = Coeff(1 / ~u, ~x, 0), e = Coeff(1 / ~u, ~x, 2), f = Coeff(1 / ~u, ~x, 4)], ~A * Subst(Antiderivative(1 / (d - (b * d - a * e) * (~x) ^ 2), ~x), ~x, ~x / sqrt(~v)) <-- EqQ(a * ~B + ~A * c, 0) && EqQ(c * d - a * f, 0)) <-- FreeQ([~A, ~B], ~x) && (PolyQ(~v, (~x) ^ 2, 2) && PolyQ(1 / ~u, (~x) ^ 2, 2))
    @apply_utils Antiderivative(1 / ((~a + ~(b') * ~x) * sqrt(~c + ~(d') * (~x) ^ 2) * sqrt(~e + ~(f') * (~x) ^ 2)), ~x) => ~a * Antiderivative(1 / (((~a) ^ 2 - (~b) ^ 2 * (~x) ^ 2) * sqrt(~c + ~d * (~x) ^ 2) * sqrt(~e + ~f * (~x) ^ 2)), ~x) - ~b * Antiderivative(~x / (((~a) ^ 2 - (~b) ^ 2 * (~x) ^ 2) * sqrt(~c + ~d * (~x) ^ 2) * sqrt(~e + ~f * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x)
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) * sqrt(~(d') + ~(e') * ~x + ~(f') * sqrt(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2)), ~x) => ((2 * ((~f * (((5 * ~b * ~c * (~g) ^ 2 - 2 * (~b) ^ 2 * ~g * ~h) - 3 * ~a * ~c * ~g * ~h) + 2 * ~a * ~b * (~h) ^ 2) + ~c * ~f * ((10 * ~c * (~g) ^ 2 - ~b * ~g * ~h) + ~a * (~h) ^ 2) * ~x + 9 * (~c) ^ 2 * ~f * ~g * ~h * (~x) ^ 2 + 3 * (~c) ^ 2 * ~f * (~h) ^ 2 * (~x) ^ 3) - (~e * ~g - ~d * ~h) * ((5 * ~c * ~g - 2 * ~b * ~h) + ~c * ~h * ~x) * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2))) / (15 * (~c) ^ 2 * ~f * (~g + ~h * ~x))) * sqrt(~d + ~e * ~x + ~f * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (EqQ((~e * ~g - ~d * ~h) ^ 2 - (~f) ^ 2 * ((~c * (~g) ^ 2 - ~b * ~g * ~h) + ~a * (~h) ^ 2), 0) && EqQ((2 * (~e) ^ 2 * ~g - 2 * ~d * ~e * ~h) - (~f) ^ 2 * (2 * ~c * ~g - ~b * ~h), 0))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * (~u + ~(f') * (~(j') + ~(k') * sqrt(~v))) ^ ~(n'), ~x) => Antiderivative((~g + ~h * ~x) ^ ~m * (ExpandToSum(~u + ~f * ~j, ~x) + ~f * ~k * sqrt(ExpandToSum(~v, ~x))) ^ ~n, ~x) <-- FreeQ([~f, ~g, ~h, ~j, ~k, ~m, ~n], ~x) && (LinearQ(~u, ~x) && (QuadraticQ(~v, ~x) && (Not(LinearMatchQ(~u, ~x) && (QuadraticMatchQ(~v, ~x) && (EqQ(~j, 0) || EqQ(~f, 1)))) && EqQ((Coefficient(~u, ~x, 1) * ~g - ~h * (Coefficient(~u, ~x, 0) + ~f * ~j)) ^ 2 - (~f) ^ 2 * (~k) ^ 2 * ((Coefficient(~v, ~x, 2) * (~g) ^ 2 - Coefficient(~v, ~x, 1) * ~g * ~h) + Coefficient(~v, ~x, 0) * (~h) ^ 2), 0))))
    #= Antiderivative(1/((~d')+(~e')*(~x)+(~f')*sqrt((~a')+(~b')*(~x)+(~c')*(~x)^2)),~x) :=  Antiderivative((d+e*x)/(d^2-a*f^2+(2*d*e-b*f^2)*x),x) - f*Antiderivative(sqrt(a+b*x+c*x^2)/(d^2-a*f^2+(2*d*e-b*f^2)*x),x) <--  FreeQ([a,b,c,d,e,f],x) && EqQ(e^2-c*f^2,0) =#
    #= Antiderivative(1/((~d')+(~e')*(~x)+(~f')*sqrt((~a')+(~c')*(~x)^2)),~x) :=  Antiderivative((d+e*x)/(d^2-a*f^2+2*d*e*x),x) - f*Antiderivative(sqrt(a+c*x^2)/(d^2-a*f^2+2*d*e*x),x) <-- FreeQ([a,c,d,e,f],x) &&  EqQ(e^2-c*f^2,0) =#
    @apply_utils Antiderivative((~(g') + ~(h') * (~(d') + ~(e') * ~x + ~(f') * sqrt(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2)) ^ ~n) ^ ~(p'), ~x) => 2 * Subst(Antiderivative(((~g + ~h * (~x) ^ ~n) ^ ~p * ((((~d) ^ 2 * ~e - (~b * ~d - ~a * ~e) * (~f) ^ 2) - (2 * ~d * ~e - ~b * (~f) ^ 2) * ~x) + ~e * (~x) ^ 2)) / (-2 * ~d * ~e + ~b * (~f) ^ 2 + 2 * ~e * ~x) ^ 2, ~x), ~x, ~d + ~e * ~x + ~f * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~n], ~x) && (EqQ((~e) ^ 2 - ~c * (~f) ^ 2, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~(g') + ~(h') * (~(d') + ~(e') * ~x + ~(f') * sqrt(~a + ~(c') * (~x) ^ 2)) ^ ~n) ^ ~(p'), ~x) => (1 / (2 * ~e)) * Subst(Antiderivative(((~g + ~h * (~x) ^ ~n) ^ ~p * ((((~d) ^ 2 + ~a * (~f) ^ 2) - 2 * ~d * ~x) + (~x) ^ 2)) / (~d - ~x) ^ 2, ~x), ~x, ~d + ~e * ~x + ~f * sqrt(~a + ~c * (~x) ^ 2)) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~h, ~n], ~x) && (EqQ((~e) ^ 2 - ~c * (~f) ^ 2, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~(g') + ~(h') * (~u + ~(f') * sqrt(~v)) ^ ~n) ^ ~(p'), ~x) => Antiderivative((~g + ~h * (ExpandToSum(~u, ~x) + ~f * sqrt(ExpandToSum(~v, ~x))) ^ ~n) ^ ~p, ~x) <-- FreeQ([~f, ~g, ~h, ~n], ~x) && (LinearQ(~u, ~x) && (QuadraticQ(~v, ~x) && (Not(LinearMatchQ(~u, ~x) && QuadraticMatchQ(~v, ~x)) && (EqQ(Coefficient(~u, ~x, 1) ^ 2 - Coefficient(~v, ~x, 2) * (~f) ^ 2, 0) && IntegerQ(~p)))))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * (~(e') * ~x + ~(f') * sqrt(~(a') + ~(c') * (~x) ^ 2)) ^ ~(n'), ~x) => (1 / (2 ^ (~m + 1) * (~e) ^ (~m + 1))) * Subst(Antiderivative((~x) ^ ((~n - ~m) - 2) * (~a * (~f) ^ 2 + (~x) ^ 2) * (-(~a) * (~f) ^ 2 * ~h + 2 * ~e * ~g * ~x + ~h * (~x) ^ 2) ^ ~m, ~x), ~x, ~e * ~x + ~f * sqrt(~a + ~c * (~x) ^ 2)) <-- FreeQ([~a, ~c, ~e, ~f, ~g, ~h, ~n], ~x) && (EqQ((~e) ^ 2 - ~c * (~f) ^ 2, 0) && IntegerQ(~m))
    @apply_utils Antiderivative((~x) ^ ~(p') * (~g + ~(i') * (~x) ^ 2) ^ ~(m') * (~(e') * ~x + ~(f') * sqrt(~a + ~(c') * (~x) ^ 2)) ^ ~(n'), ~x) => (1 / (2 ^ (2 * ~m + ~p + 1) * (~e) ^ (~p + 1) * (~f) ^ (2 * ~m))) * (~i / ~c) ^ ~m * Subst(Antiderivative((~x) ^ (((~n - 2 * ~m) - ~p) - 2) * (-(~a) * (~f) ^ 2 + (~x) ^ 2) ^ ~p * (~a * (~f) ^ 2 + (~x) ^ 2) ^ (2 * ~m + 1), ~x), ~x, ~e * ~x + ~f * sqrt(~a + ~c * (~x) ^ 2)) <-- FreeQ([~a, ~c, ~e, ~f, ~g, ~i, ~n], ~x) && (EqQ((~e) ^ 2 - ~c * (~f) ^ 2, 0) && (EqQ(~c * ~g - ~a * ~i, 0) && (IntegersQ(~p, 2 * ~m) && (IntegerQ(~m) || GtQ(~i / ~c, 0)))))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x + ~(i') * (~x) ^ 2) ^ ~(m') * (~(d') + ~(e') * ~x + ~(f') * sqrt(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2)) ^ ~(n'), ~x) => (2 / (~f) ^ (2 * ~m)) * (~i / ~c) ^ ~m * Subst(Antiderivative(((~x) ^ ~n * ((((~d) ^ 2 * ~e - (~b * ~d - ~a * ~e) * (~f) ^ 2) - (2 * ~d * ~e - ~b * (~f) ^ 2) * ~x) + ~e * (~x) ^ 2) ^ (2 * ~m + 1)) / (-2 * ~d * ~e + ~b * (~f) ^ 2 + 2 * ~e * ~x) ^ (2 * (~m + 1)), ~x), ~x, ~d + ~e * ~x + ~f * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~n], ~x) && (EqQ((~e) ^ 2 - ~c * (~f) ^ 2, 0) && (EqQ(~c * ~g - ~a * ~i, 0) && (EqQ(~c * ~h - ~b * ~i, 0) && (IntegerQ(2 * ~m) && (IntegerQ(~m) || GtQ(~i / ~c, 0))))))
    @apply_utils Antiderivative((~g + ~(i') * (~x) ^ 2) ^ ~(m') * (~(d') + ~(e') * ~x + ~(f') * sqrt(~a + ~(c') * (~x) ^ 2)) ^ ~(n'), ~x) => (1 / (2 ^ (2 * ~m + 1) * ~e * (~f) ^ (2 * ~m))) * (~i / ~c) ^ ~m * Subst(Antiderivative(((~x) ^ ~n * ((((~d) ^ 2 + ~a * (~f) ^ 2) - 2 * ~d * ~x) + (~x) ^ 2) ^ (2 * ~m + 1)) / (-(~d) + ~x) ^ (2 * (~m + 1)), ~x), ~x, ~d + ~e * ~x + ~f * sqrt(~a + ~c * (~x) ^ 2)) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~i, ~n], ~x) && (EqQ((~e) ^ 2 - ~c * (~f) ^ 2, 0) && (EqQ(~c * ~g - ~a * ~i, 0) && (IntegerQ(2 * ~m) && (IntegerQ(~m) || GtQ(~i / ~c, 0)))))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x + ~(i') * (~x) ^ 2) ^ ~(m') * (~(d') + ~(e') * ~x + ~(f') * sqrt(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2)) ^ ~(n'), ~x) => (((~i / ~c) ^ (~m - 1 / 2) * sqrt(~g + ~h * ~x + ~i * (~x) ^ 2)) / sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~m * (~d + ~e * ~x + ~f * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~n], ~x) && (EqQ((~e) ^ 2 - ~c * (~f) ^ 2, 0) && (EqQ(~c * ~g - ~a * ~i, 0) && (EqQ(~c * ~h - ~b * ~i, 0) && (IGtQ(~m + 1 / 2, 0) && Not(GtQ(~i / ~c, 0))))))
    @apply_utils Antiderivative((~g + ~(i') * (~x) ^ 2) ^ ~(m') * (~(d') + ~(e') * ~x + ~(f') * sqrt(~a + ~(c') * (~x) ^ 2)) ^ ~(n'), ~x) => (((~i / ~c) ^ (~m - 1 / 2) * sqrt(~g + ~i * (~x) ^ 2)) / sqrt(~a + ~c * (~x) ^ 2)) * Antiderivative((~a + ~c * (~x) ^ 2) ^ ~m * (~d + ~e * ~x + ~f * sqrt(~a + ~c * (~x) ^ 2)) ^ ~n, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~i, ~n], ~x) && (EqQ((~e) ^ 2 - ~c * (~f) ^ 2, 0) && (EqQ(~c * ~g - ~a * ~i, 0) && (IGtQ(~m + 1 / 2, 0) && Not(GtQ(~i / ~c, 0)))))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x + ~(i') * (~x) ^ 2) ^ ~(m') * (~(d') + ~(e') * ~x + ~(f') * sqrt(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2)) ^ ~(n'), ~x) => (((~i / ~c) ^ (~m + 1 / 2) * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)) / sqrt(~g + ~h * ~x + ~i * (~x) ^ 2)) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~m * (~d + ~e * ~x + ~f * sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~i, ~n], ~x) && (EqQ((~e) ^ 2 - ~c * (~f) ^ 2, 0) && (EqQ(~c * ~g - ~a * ~i, 0) && (EqQ(~c * ~h - ~b * ~i, 0) && (ILtQ(~m - 1 / 2, 0) && Not(GtQ(~i / ~c, 0))))))
    @apply_utils Antiderivative((~g + ~(i') * (~x) ^ 2) ^ ~(m') * (~(d') + ~(e') * ~x + ~(f') * sqrt(~a + ~(c') * (~x) ^ 2)) ^ ~(n'), ~x) => (((~i / ~c) ^ (~m + 1 / 2) * sqrt(~a + ~c * (~x) ^ 2)) / sqrt(~g + ~i * (~x) ^ 2)) * Antiderivative((~a + ~c * (~x) ^ 2) ^ ~m * (~d + ~e * ~x + ~f * sqrt(~a + ~c * (~x) ^ 2)) ^ ~n, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~i, ~n], ~x) && (EqQ((~e) ^ 2 - ~c * (~f) ^ 2, 0) && (EqQ(~c * ~g - ~a * ~i, 0) && (ILtQ(~m - 1 / 2, 0) && Not(GtQ(~i / ~c, 0)))))
    @apply_utils Antiderivative((~w) ^ ~(m') * (~u + ~(f') * (~(j') + ~(k') * sqrt(~v))) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~w, ~x) ^ ~m * (ExpandToSum(~u + ~f * ~j, ~x) + ~f * ~k * sqrt(ExpandToSum(~v, ~x))) ^ ~n, ~x) <-- FreeQ([~f, ~j, ~k, ~m, ~n], ~x) && (LinearQ(~u, ~x) && (QuadraticQ([~v, ~w], ~x) && (Not(LinearMatchQ(~u, ~x) && (QuadraticMatchQ([~v, ~w], ~x) && (EqQ(~j, 0) || EqQ(~f, 1)))) && EqQ(Coefficient(~u, ~x, 1) ^ 2 - Coefficient(~v, ~x, 2) * (~f) ^ 2 * (~k) ^ 2, 0))))
    @apply_utils Antiderivative(1 / ((~a + ~(b') * (~x) ^ ~(n')) * sqrt(~(c') * (~x) ^ 2 + ~(d') * (~a + ~(b') * (~x) ^ ~(n')) ^ ~(p'))), ~x) => (1 / ~a) * Subst(Antiderivative(1 / (1 - ~c * (~x) ^ 2), ~x), ~x, ~x / sqrt(~c * (~x) ^ 2 + ~d * (~a + ~b * (~x) ^ ~n) ^ (2 / ~n))) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && EqQ(~p, 2 / ~n)
    @apply_utils Antiderivative(sqrt(~a + ~(b') * sqrt(~c + ~(d') * (~x) ^ 2)), ~x) => (2 * (~b) ^ 2 * ~d * (~x) ^ 3) / (3 * (~a + ~b * sqrt(~c + ~d * (~x) ^ 2)) ^ (3 / 2)) + (2 * ~a * ~x) / sqrt(~a + ~b * sqrt(~c + ~d * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~a) ^ 2 - (~b) ^ 2 * ~c, 0)
    @apply_utils Antiderivative(sqrt(~(a') * (~x) ^ 2 + ~(b') * ~x * sqrt(~c + ~(d') * (~x) ^ 2)) / (~x * sqrt(~c + ~(d') * (~x) ^ 2)), ~x) => ((sqrt(2) * ~b) / ~a) * Subst(Antiderivative(1 / sqrt(1 + (~x) ^ 2 / ~a), ~x), ~x, ~a * ~x + ~b * sqrt(~c + ~d * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2 * ~d, 0) && EqQ((~b) ^ 2 * ~c + ~a, 0))
    @apply_utils Antiderivative(sqrt(~(e') * ~x * (~(a') * ~x + ~(b') * sqrt(~c + ~(d') * (~x) ^ 2))) / (~x * sqrt(~c + ~(d') * (~x) ^ 2)), ~x) => Antiderivative(sqrt(~a * ~e * (~x) ^ 2 + ~b * ~e * ~x * sqrt(~c + ~d * (~x) ^ 2)) / (~x * sqrt(~c + ~d * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2 * ~d, 0) && EqQ((~b) ^ 2 * ~c * ~e + ~a, 0))
    @apply_utils Antiderivative(sqrt(~(c') * (~x) ^ 2 + ~(d') * sqrt(~a + ~(b') * (~x) ^ 4)) / sqrt(~a + ~(b') * (~x) ^ 4), ~x) => ~d * Subst(Antiderivative(1 / (1 - 2 * ~c * (~x) ^ 2), ~x), ~x, ~x / sqrt(~c * (~x) ^ 2 + ~d * sqrt(~a + ~b * (~x) ^ 4))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c) ^ 2 - ~b * (~d) ^ 2, 0)
    @apply_utils Antiderivative(((~(c') + ~(d') * ~x) ^ ~(m') * sqrt(~(b') * (~x) ^ 2 + sqrt(~a + ~(e') * (~x) ^ 4))) / sqrt(~a + ~(e') * (~x) ^ 4), ~x) => ((1 - im) / 2) * Antiderivative((~c + ~d * ~x) ^ ~m / sqrt(sqrt(~a) - im * ~b * (~x) ^ 2), ~x) + ((1 + im) / 2) * Antiderivative((~c + ~d * ~x) ^ ~m / sqrt(sqrt(~a) + im * ~b * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (EqQ(~e, (~b) ^ 2) && GtQ(~a, 0))
    @apply_utils Antiderivative(1 / ((~c + ~(d') * ~x) * sqrt(~a + ~(b') * (~x) ^ 3)), ~x) => (2 / (3 * ~c)) * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 3), ~x) + (1 / (3 * ~c)) * Antiderivative((~c - 2 * ~d * ~x) / ((~c + ~d * ~x) * sqrt(~a + ~b * (~x) ^ 3)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ(~b * (~c) ^ 3 - 4 * ~a * (~d) ^ 3, 0)
    @apply_utils Antiderivative(1 / ((~c + ~(d') * ~x) * sqrt(~a + ~(b') * (~x) ^ 3)), ~x) => ((-6 * ~a * (~d) ^ 3) / (~c * (~b * (~c) ^ 3 - 28 * ~a * (~d) ^ 3))) * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 3), ~x) + (1 / (~c * (~b * (~c) ^ 3 - 28 * ~a * (~d) ^ 3))) * Antiderivative(Simp(~c * (~b * (~c) ^ 3 - 22 * ~a * (~d) ^ 3) + 6 * ~a * (~d) ^ 4 * ~x, ~x) / ((~c + ~d * ~x) * sqrt(~a + ~b * (~x) ^ 3)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ(((~b) ^ 2 * (~c) ^ 6 - 20 * ~a * ~b * (~c) ^ 3 * (~d) ^ 3) - 8 * (~a) ^ 2 * (~d) ^ 6, 0)
    @apply_utils Antiderivative(1 / ((~c + ~(d') * ~x) * sqrt(~a + ~(b') * (~x) ^ 3)), ~x) => With([q = Rt(~b / ~a, 3)], (-q / ((1 + sqrt(3)) * ~d - ~c * q)) * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 3), ~x) + (~d / ((1 + sqrt(3)) * ~d - ~c * q)) * Antiderivative((1 + sqrt(3) + q * ~x) / ((~c + ~d * ~x) * sqrt(~a + ~b * (~x) ^ 3)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ(((~b) ^ 2 * (~c) ^ 6 - 20 * ~a * ~b * (~c) ^ 3 * (~d) ^ 3) - 8 * (~a) ^ 2 * (~d) ^ 6, 0)
    @apply_utils Antiderivative((~e + ~(f') * ~x) / ((~c + ~(d') * ~x) * sqrt(~a + ~(b') * (~x) ^ 3)), ~x) => ((2 * ~e) / ~d) * Subst(Antiderivative(1 / (1 + 3 * ~a * (~x) ^ 2), ~x), ~x, (1 + (2 * ~d * ~x) / ~c) / sqrt(~a + ~b * (~x) ^ 3)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~d * ~e - ~c * ~f, 0) && (EqQ(~b * (~c) ^ 3 - 4 * ~a * (~d) ^ 3, 0) && EqQ(2 * ~d * ~e + ~c * ~f, 0)))
    @apply_utils Antiderivative((~e + ~(f') * ~x) / ((~c + ~(d') * ~x) * sqrt(~a + ~(b') * (~x) ^ 3)), ~x) => ((-2 * ~e) / ~d) * Subst(Antiderivative(1 / (9 - ~a * (~x) ^ 2), ~x), ~x, (1 + (~f * ~x) / ~e) ^ 2 / sqrt(~a + ~b * (~x) ^ 3)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~d * ~e - ~c * ~f, 0) && (EqQ(~b * (~c) ^ 3 + 8 * ~a * (~d) ^ 3, 0) && EqQ(2 * ~d * ~e + ~c * ~f, 0)))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) / ((~c + ~(d') * ~x) * sqrt(~a + ~(b') * (~x) ^ 3)), ~x) => ((2 * ~d * ~e + ~c * ~f) / (3 * ~c * ~d)) * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 3), ~x) + ((~d * ~e - ~c * ~f) / (3 * ~c * ~d)) * Antiderivative((~c - 2 * ~d * ~x) / ((~c + ~d * ~x) * sqrt(~a + ~b * (~x) ^ 3)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~d * ~e - ~c * ~f, 0) && ((EqQ(~b * (~c) ^ 3 - 4 * ~a * (~d) ^ 3, 0) || EqQ(~b * (~c) ^ 3 + 8 * ~a * (~d) ^ 3, 0)) && NeQ(2 * ~d * ~e + ~c * ~f, 0)))
    @apply_utils Antiderivative((~e + ~(f') * ~x) / ((~c + ~(d') * ~x) * sqrt(~a + ~(b') * (~x) ^ 3)), ~x) => With([k = Simplify((~d * ~e + 2 * ~c * ~f) / (~c * ~f))], (((1 + k) * ~e) / ~d) * Subst(Antiderivative(1 / (1 + (3 + 2k) * ~a * (~x) ^ 2), ~x), ~x, (1 + ((1 + k) * ~d * ~x) / ~c) / sqrt(~a + ~b * (~x) ^ 3))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~d * ~e - ~c * ~f, 0) && (EqQ(((~b) ^ 2 * (~c) ^ 6 - 20 * ~a * ~b * (~c) ^ 3 * (~d) ^ 3) - 8 * (~a) ^ 2 * (~d) ^ 6, 0) && EqQ(6 * ~a * (~d) ^ 4 * ~e - ~c * ~f * (~b * (~c) ^ 3 - 22 * ~a * (~d) ^ 3), 0)))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) / ((~c + ~(d') * ~x) * sqrt(~a + ~(b') * (~x) ^ 3)), ~x) => (-((6 * ~a * (~d) ^ 4 * ~e - ~c * ~f * (~b * (~c) ^ 3 - 22 * ~a * (~d) ^ 3))) / (~c * ~d * (~b * (~c) ^ 3 - 28 * ~a * (~d) ^ 3))) * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 3), ~x) + ((~d * ~e - ~c * ~f) / (~c * ~d * (~b * (~c) ^ 3 - 28 * ~a * (~d) ^ 3))) * Antiderivative((~c * (~b * (~c) ^ 3 - 22 * ~a * (~d) ^ 3) + 6 * ~a * (~d) ^ 4 * ~x) / ((~c + ~d * ~x) * sqrt(~a + ~b * (~x) ^ 3)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~d * ~e - ~c * ~f, 0) && (EqQ(((~b) ^ 2 * (~c) ^ 6 - 20 * ~a * ~b * (~c) ^ 3 * (~d) ^ 3) - 8 * (~a) ^ 2 * (~d) ^ 6, 0) && NeQ(6 * ~a * (~d) ^ 4 * ~e - ~c * ~f * (~b * (~c) ^ 3 - 22 * ~a * (~d) ^ 3), 0)))
    #= Antiderivative(((~e)+(~f')*(~x))/(((~c)+(~d')*(~x))*sqrt((~a)+(~b')*(~x)^3)),~x) :=  With([q=(1+sqrt(3))*f/e], 4*3^(1/4)*sqrt(2-sqrt(3))*f*(1+sqrt(3)+q*x)^2*sqrt((1+q^3*x^3)/(1+ sqrt(3)+q*x)^4)/(q*sqrt(a+b*x^3))* Subst(Antiderivative(1/(((1-sqrt(3))*d-c*q+((1+sqrt(3))*d-c*q)*x)* sqrt(7-4*sqrt(3)-2*(3-2*sqrt(3))*x^2-x^4)),x),x,(-1+sqrt(3)-q*x) /(1+sqrt(3)+q*x))) <-- FreeQ([a,b,c,d,e,f],x) && NeQ(d*e-c*f,0) &&  EqQ(b*e^3-2*(5+3*sqrt(3))*a*f^3,0) &&  NeQ(b*c^3-2*(5-3*sqrt(3))*a*d^3,0) =#
    @apply_utils Antiderivative((~e + ~(f') * ~x) / ((~c + ~(d') * ~x) * sqrt(~a + ~(b') * (~x) ^ 3)), ~x) => With([q = Simplify(((1 + sqrt(3)) * ~f) / ~e)], ((4 * 3 ^ (1 / 4) * sqrt(2 - sqrt(3)) * ~f * (1 + q * ~x) * sqrt(((1 - q * ~x) + q ^ 2 * (~x) ^ 2) / (1 + sqrt(3) + q * ~x) ^ 2)) / (q * sqrt(~a + ~b * (~x) ^ 3) * sqrt((1 + q * ~x) / (1 + sqrt(3) + q * ~x) ^ 2))) * Subst(Antiderivative(1 / ((((1 - sqrt(3)) * ~d - ~c * q) + ((1 + sqrt(3)) * ~d - ~c * q) * ~x) * sqrt(1 - (~x) ^ 2) * sqrt((7 - 4 * sqrt(3)) + (~x) ^ 2)), ~x), ~x, ((-1 + sqrt(3)) - q * ~x) / (1 + sqrt(3) + q * ~x))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~d * ~e - ~c * ~f, 0) && (EqQ(~b * (~e) ^ 3 - 2 * (5 + 3 * sqrt(3)) * ~a * (~f) ^ 3, 0) && NeQ(~b * (~c) ^ 3 - 2 * (5 - 3 * sqrt(3)) * ~a * (~d) ^ 3, 0)))
    @apply_utils Antiderivative((~e + ~(f') * ~x) / ((~c + ~(d') * ~x) * sqrt(~a + ~(b') * (~x) ^ 3)), ~x) => With([q = Simplify(((-1 + sqrt(3)) * ~f) / ~e)], ((4 * 3 ^ (1 / 4) * sqrt(2 + sqrt(3)) * ~f * (1 - q * ~x) * sqrt((1 + q * ~x + q ^ 2 * (~x) ^ 2) / ((1 - sqrt(3)) - q * ~x) ^ 2)) / (q * sqrt(~a + ~b * (~x) ^ 3) * sqrt(-((1 - q * ~x)) / ((1 - sqrt(3)) - q * ~x) ^ 2))) * Subst(Antiderivative(1 / (((1 + sqrt(3)) * ~d + ~c * q + ((1 - sqrt(3)) * ~d + ~c * q) * ~x) * sqrt(1 - (~x) ^ 2) * sqrt(7 + 4 * sqrt(3) + (~x) ^ 2)), ~x), ~x, ((1 + sqrt(3)) - q * ~x) / (-1 + sqrt(3) + q * ~x))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~d * ~e - ~c * ~f, 0) && (EqQ(~b * (~e) ^ 3 - 2 * (5 - 3 * sqrt(3)) * ~a * (~f) ^ 3, 0) && NeQ(~b * (~c) ^ 3 - 2 * (5 + 3 * sqrt(3)) * ~a * (~d) ^ 3, 0)))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) / ((~c + ~(d') * ~x) * sqrt(~a + ~(b') * (~x) ^ 3)), ~x) => With([q = Rt(~b / ~a, 3)], (((1 + sqrt(3)) * ~f - ~e * q) / ((1 + sqrt(3)) * ~d - ~c * q)) * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 3), ~x) + ((~d * ~e - ~c * ~f) / ((1 + sqrt(3)) * ~d - ~c * q)) * Antiderivative((1 + sqrt(3) + q * ~x) / ((~c + ~d * ~x) * sqrt(~a + ~b * (~x) ^ 3)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~d * ~e - ~c * ~f, 0) && (NeQ(((~b) ^ 2 * (~c) ^ 6 - 20 * ~a * ~b * (~c) ^ 3 * (~d) ^ 3) - 8 * (~a) ^ 2 * (~d) ^ 6, 0) && NeQ(((~b) ^ 2 * (~e) ^ 6 - 20 * ~a * ~b * (~e) ^ 3 * (~f) ^ 3) - 8 * (~a) ^ 2 * (~f) ^ 6, 0)))
    @apply_utils Antiderivative((~f + ~(g') * ~x + ~(h') * (~x) ^ 2) / ((~c + ~(d') * ~x + ~(e') * (~x) ^ 2) * sqrt(~a + ~(b') * (~x) ^ 3)), ~x) => -2 * ~g * ~h * Subst(Antiderivative(1 / (2 * ~e * ~h - (~b * ~d * ~f - 2 * ~a * ~e * ~h) * (~x) ^ 2), ~x), ~x, (1 + (2 * ~h * ~x) / ~g) / sqrt(~a + ~b * (~x) ^ 3)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ(~b * ~d * ~f - 2 * ~a * ~e * ~h, 0) && (EqQ(~b * (~g) ^ 3 - 8 * ~a * (~h) ^ 3, 0) && (EqQ((~g) ^ 2 + 2 * ~f * ~h, 0) && EqQ((~b * ~d * ~f + ~b * ~c * ~g) - 4 * ~a * ~e * ~h, 0))))
    @apply_utils Antiderivative((~f + ~(g') * ~x + ~(h') * (~x) ^ 2) / ((~c + ~(e') * (~x) ^ 2) * sqrt(~a + ~(b') * (~x) ^ 3)), ~x) => (-(~g) / ~e) * Subst(Antiderivative(1 / (1 + ~a * (~x) ^ 2), ~x), ~x, (1 + (2 * ~h * ~x) / ~g) / sqrt(~a + ~b * (~x) ^ 3)) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~g, ~h], ~x) && (EqQ(~b * (~g) ^ 3 - 8 * ~a * (~h) ^ 3, 0) && (EqQ((~g) ^ 2 + 2 * ~f * ~h, 0) && EqQ(~b * ~c * ~g - 4 * ~a * ~e * ~h, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * (~x) ^ 3) / (~c + ~(d') * ~x), ~x) => (~b / ~d) * Antiderivative((~x) ^ 2 / sqrt(~a + ~b * (~x) ^ 3), ~x) + ((~b * ~c) / (~d) ^ 3) * Antiderivative((~c - ~d * ~x) / sqrt(~a + ~b * (~x) ^ 3), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ(~b * (~c) ^ 3 - ~a * (~d) ^ 3, 0)
    @apply_utils Antiderivative(sqrt(~a + ~(b') * (~x) ^ 3) / (~c + ~(d') * ~x), ~x) => ((~b / ~d) * Antiderivative((~x) ^ 2 / sqrt(~a + ~b * (~x) ^ 3), ~x) + ((~b * ~c) / (~d) ^ 3) * Antiderivative((~c - ~d * ~x) / sqrt(~a + ~b * (~x) ^ 3), ~x)) - ((~b * (~c) ^ 3 - ~a * (~d) ^ 3) / (~d) ^ 3) * Antiderivative(1 / ((~c + ~d * ~x) * sqrt(~a + ~b * (~x) ^ 3)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ(~b * (~c) ^ 3 - ~a * (~d) ^ 3, 0)
    @apply_utils Antiderivative(1 / ((~c + ~(d') * ~x) * (~a + ~(b') * (~x) ^ 3) ^ (1 / 3)), ~x) => ((sqrt(3) * atan((1 - (2 ^ (1 / 3) * Rt(~b, 3) * (~c - ~d * ~x)) / (~d * (~a + ~b * (~x) ^ 3) ^ (1 / 3))) / sqrt(3))) / (2 ^ (4 / 3) * Rt(~b, 3) * ~c) + log((~c + ~d * ~x) ^ 2 * (~c - ~d * ~x)) / (2 ^ (7 / 3) * Rt(~b, 3) * ~c)) - (3 * log(Rt(~b, 3) * (~c - ~d * ~x) + 2 ^ (2 / 3) * ~d * (~a + ~b * (~x) ^ 3) ^ (1 / 3))) / (2 ^ (7 / 3) * Rt(~b, 3) * ~c) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ(~b * (~c) ^ 3 + ~a * (~d) ^ 3, 0)
    @apply_utils Antiderivative(1 / ((~c + ~(d') * ~x) * (~a + ~(b') * (~x) ^ 3) ^ (1 / 3)), ~x) => (1 / (2 * ~c)) * Antiderivative(1 / (~a + ~b * (~x) ^ 3) ^ (1 / 3), ~x) + (1 / (2 * ~c)) * Antiderivative((~c - ~d * ~x) / ((~c + ~d * ~x) * (~a + ~b * (~x) ^ 3) ^ (1 / 3)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ(2 * ~b * (~c) ^ 3 - ~a * (~d) ^ 3, 0)
    @apply_utils Antiderivative((~e + ~(f') * ~x) / ((~c + ~(d') * ~x) * (~a + ~(b') * (~x) ^ 3) ^ (1 / 3)), ~x) => ((sqrt(3) * ~f * atan((1 + (2 * Rt(~b, 3) * (2 * ~c + ~d * ~x)) / (~d * (~a + ~b * (~x) ^ 3) ^ (1 / 3))) / sqrt(3))) / (Rt(~b, 3) * ~d) + (~f * log(~c + ~d * ~x)) / (Rt(~b, 3) * ~d)) - (3 * ~f * log(Rt(~b, 3) * (2 * ~c + ~d * ~x) - ~d * (~a + ~b * (~x) ^ 3) ^ (1 / 3))) / (2 * Rt(~b, 3) * ~d) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ(~d * ~e + ~c * ~f, 0) && EqQ(2 * ~b * (~c) ^ 3 - ~a * (~d) ^ 3, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) / ((~(c') + ~(d') * ~x) * (~a + ~(b') * (~x) ^ 3) ^ (1 / 3)), ~x) => (~f / ~d) * Antiderivative(1 / (~a + ~b * (~x) ^ 3) ^ (1 / 3), ~x) + ((~d * ~e - ~c * ~f) / ~d) * Antiderivative(1 / ((~c + ~d * ~x) * (~a + ~b * (~x) ^ 3) ^ (1 / 3)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x)
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ 3) ^ (2 / 3) / (~c + ~(d') * ~x), ~x) => ((~a + ~b * (~x) ^ 3) ^ (2 / 3) / (2 * ~d) - ((~b * ~c) / (~d) ^ 2) * Antiderivative(~x / (~a + ~b * (~x) ^ 3) ^ (1 / 3), ~x)) + (1 / (~d) ^ 2) * Antiderivative((~a * (~d) ^ 2 + ~b * (~c) ^ 2 * ~x) / ((~c + ~d * ~x) * (~a + ~b * (~x) ^ 3) ^ (1 / 3)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x)
    @apply_utils Antiderivative(1 / ((~c + ~(d') * ~x) * (~a + ~(b') * (~x) ^ 3) ^ (2 / 3)), ~x) => With([q = Rt(~b, 3)], ((((-(~d) * atan((1 + (2 * q * ~x) / (~a + ~b * (~x) ^ 3) ^ (1 / 3)) / sqrt(3))) / (2 * sqrt(3) * q ^ 2 * (~c) ^ 2) + (sqrt(3) * ~d * atan((1 + (2 * q * (2 * ~c + ~d * ~x)) / (~d * (~a + ~b * (~x) ^ 3) ^ (1 / 3))) / sqrt(3))) / (2 * q ^ 2 * (~c) ^ 2)) - (~d * log(~c + ~d * ~x)) / (2 * q ^ 2 * (~c) ^ 2)) - (~d * log(q * ~x - (~a + ~b * (~x) ^ 3) ^ (1 / 3))) / (4 * q ^ 2 * (~c) ^ 2)) + (3 * ~d * log(q * (2 * ~c + ~d * ~x) - ~d * (~a + ~b * (~x) ^ 3) ^ (1 / 3))) / (4 * q ^ 2 * (~c) ^ 2)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ(2 * ~b * (~c) ^ 3 - ~a * (~d) ^ 3, 0)
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Px * (~c + ~(d') * ~x) ^ ~q * (~a + ~(b') * (~x) ^ 3) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand(((~c) ^ 3 + (~d) ^ 3 * (~x) ^ 3) ^ ~q * (~a + ~b * (~x) ^ 3) ^ ~p, ((~x) ^ ~m * ~Px) / (((~c) ^ 2 - ~c * ~d * ~x) + (~d) ^ 2 * (~x) ^ 2) ^ ~q, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~p], ~x) && (PolyQ(~Px, ~x) && (ILtQ(~q, 0) && (IntegerQ(~m) && (RationalQ(~p) && EqQ(Denominator(~p), 3)))))
    @apply_utils Antiderivative(~(Px') * (~c + ~(d') * ~x) ^ ~q * (~a + ~(b') * (~x) ^ 3) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand(((~c) ^ 3 + (~d) ^ 3 * (~x) ^ 3) ^ ~q * (~a + ~b * (~x) ^ 3) ^ ~p, ~Px / (((~c) ^ 2 - ~c * ~d * ~x) + (~d) ^ 2 * (~x) ^ 2) ^ ~q, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && (PolyQ(~Px, ~x) && (ILtQ(~q, 0) && (RationalQ(~p) && EqQ(Denominator(~p), 3))))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Px * (~c + ~(d') * ~x + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(b') * (~x) ^ 3) ^ ~(p'), ~x) => (1 / (~c) ^ ~q) * Antiderivative(ExpandIntegrand(((~c) ^ 3 - (~d) ^ 3 * (~x) ^ 3) ^ ~q * (~a + ~b * (~x) ^ 3) ^ ~p, ((~x) ^ ~m * ~Px) / (~c - ~d * ~x) ^ ~q, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~p], ~x) && (PolyQ(~Px, ~x) && (EqQ((~d) ^ 2 - ~c * ~e, 0) && (ILtQ(~q, 0) && (IntegerQ(~m) && (RationalQ(~p) && EqQ(Denominator(~p), 3))))))
    @apply_utils Antiderivative(~(Px') * (~c + ~(d') * ~x + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(b') * (~x) ^ 3) ^ ~(p'), ~x) => (1 / (~c) ^ ~q) * Antiderivative(ExpandIntegrand(((~c) ^ 3 - (~d) ^ 3 * (~x) ^ 3) ^ ~q * (~a + ~b * (~x) ^ 3) ^ ~p, ~Px / (~c - ~d * ~x) ^ ~q, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && (PolyQ(~Px, ~x) && (EqQ((~d) ^ 2 - ~c * ~e, 0) && (ILtQ(~q, 0) && (RationalQ(~p) && EqQ(Denominator(~p), 3)))))
    @apply_utils Antiderivative((~c + ~(d') * (~x) ^ ~(n')) ^ ~q * (~a + ~(b') * (~x) ^ ~(nn')) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~a + ~b * (~x) ^ ~nn) ^ ~p, (~c / ((~c) ^ 2 - (~d) ^ 2 * (~x) ^ (2 * ~n)) - (~d * (~x) ^ ~n) / ((~c) ^ 2 - (~d) ^ 2 * (~x) ^ (2 * ~n))) ^ -(~q), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~nn, ~p], ~x) && (Not(IntegerQ(~p)) && (ILtQ(~q, 0) && IGtQ(log(2, ~nn / ~n), 0)))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~c + ~(d') * (~x) ^ ~(n')) ^ ~q * (~a + ~(b') * (~x) ^ ~(nn')) ^ ~p, ~x) => ((~e * ~x) ^ ~m / (~x) ^ ~m) * Antiderivative(ExpandIntegrand((~x) ^ ~m * (~a + ~b * (~x) ^ ~nn) ^ ~p, (~c / ((~c) ^ 2 - (~d) ^ 2 * (~x) ^ (2 * ~n)) - (~d * (~x) ^ ~n) / ((~c) ^ 2 - (~d) ^ 2 * (~x) ^ (2 * ~n))) ^ -(~q), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~nn, ~p], ~x) && (Not(IntegerQ(~p)) && (ILtQ(~q, 0) && IGtQ(log(2, ~nn / ~n), 0)))
    @apply_utils Antiderivative((~x) ^ ~(m') / (~c + ~(d') * (~x) ^ ~n + ~(e') * sqrt(~a + ~(b') * (~x) ^ ~n)), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ ((~m + 1) / ~n - 1) / (~c + ~d * ~x + ~e * sqrt(~a + ~b * ~x)), ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && IntegerQ((~m + 1) / ~n))
    @apply_utils Antiderivative(~(u') / (~c + ~(d') * (~x) ^ ~n + ~(e') * sqrt(~a + ~(b') * (~x) ^ ~n)), ~x) => ~c * Antiderivative(~u / (((~c) ^ 2 - ~a * (~e) ^ 2) + ~c * ~d * (~x) ^ ~n), ~x) - ~a * ~e * Antiderivative(~u / ((((~c) ^ 2 - ~a * (~e) ^ 2) + ~c * ~d * (~x) ^ ~n) * sqrt(~a + ~b * (~x) ^ ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && EqQ(~b * ~c - ~a * ~d, 0)
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.3 Miscellaneous/1.3.3 Miscellaneous algebraic functions.jl


text:
# This file needs `src/brusselator_produce.jl` to have run with appropriate parameters.
# The `saveat` parameter of that file is the time output of the frames.
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))
using Random, OrdinaryDiffEq
import GLMakie

a = 9.0
b = 10.2
L = 50.0 # size of domain
N = 250  # no. of grid points
d = 1.7

prefix = "brusselator2D_dense"
fname = savename(prefix, @strdict(a, b, d, L, N, ic), "jld2")
data = load(datadir("Brusselator", fname))
uout = data["u"]
tvec = data["t"]

heatobs = GLMakie.Observable(uout[1][1,:,:])
tobs = GLMakie.Observable(0.0)
titobs = GLMakie.lift(t -> "t = $(t)", tobs)

fig = GLMakie.Figure(resolution = (600, 550))
ax = fig[1,1] = GLMakie.Axis(fig; title = titobs)
hmap = GLMakie.heatmap!(ax, heatobs; colormap = :inferno, colorrange = (0, 1.5))
cb = GLMakie.Colorbar(fig[1, 2], hmap; width = 20)
display(fig)

GLMakie.record(
        fig, projectdir("animations", "11", "brusselator_d=$(d)_b=$(b).mp4"),
        1:length(tvec); framerate = 15
    ) do i

    tobs[] = tvec[i]
    heatobs[] = uout[i][1,:,:]
end

ext: jl

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using LibExpat

xdoc = raw"""<inventory title="OmniCorp Store #45x10^3">

    <item upc="123456789" stock="12">
      <name>Invisibility Cream</name>
      <price>14.50</price>
      <description>Makes you invisible</description>
    </item>
    <item upc="445322344" stock="18">
      <name>Levitation Salve</name>
      <price>23.99</price>
      <description>Levitate yourself for up to 3 hours per application</description>
    </item>

    <item upc="485672034" stock="653">
      <name>Blork and Freen Instameal</name>
      <price>4.95</price>
      <description>A tasty meal in a tablet; just add water</description>
    </item>
    <item upc="132957764" stock="44">
      <name>Grob winglets</name>
      <price>3.56</price>
      <description>Tender winglets of Grob. Just add water</description>
    </item>

</inventory>
"""

debracket(s) = replace(s, r".+\>(.+)\<.+" => s"\1")

etree = xp_parse(xdoc)
firstshow = LibExpat.find(etree, "//item")[1]
println("The first item's node XML entry is:\n", firstshow, "\n\n")

prices = LibExpat.find(etree, "//price")
println("Prices:")
for p in prices
    println("\t", debracket(string(p)))
end
println("\n")

namearray = LibExpat.find(etree, "//name")
println("Array of names of items:\n\t", map(s -> debracket(string(s)), namearray))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Documenter, JuliaDB, IndexedTables

@info "makedocs"
makedocs(
   clean = true,
   debug = true,
   format = Documenter.HTML(),
   sitename = "JuliaDB.jl",
   pages = [
        "index.md",
        "basics.md",
        "operations.md",
        "joins.md",
        "onlinestats.md",
        "plotting.md",
        "missing_values.md",
        "out_of_core.md",
        "ml.md",
        "tutorial.md",
        "api.md",
   ]
)

@info "deploydocs"
deploydocs(
    repo = "github.com/JuliaComputing/JuliaDB.jl.git"
)


text:
##############################################################################
##
## Default values for unspecified objects
##
## Sometimes needed when dealing with NA's for which some value must exist in
## the underlying data vector
##
##############################################################################

baseval{T <: String}(s::Type{T}) = ""
baseval(x::Any) = zero(x)

# slow, but maintains order and seems to work:
function _setdiff(a::Vector, b::Vector)
    idx = Int[]
    for i in 1:length(a)
        if !contains(b, a[i])
            push!(idx, i)
        end
    end
    a[idx]
end

## Issue: this doesn't maintain the order in a:
## setdiff(a::Vector, b::Vector) = elements(Set(a...) - Set(b...))

# MAYBE try to get in base?
function fill(x::Vector, lengths::Vector{Int})
    if length(x) != length(lengths)
        error("vector lengths must match")
    end
    res = similar(x, sum(lengths))
    i = 1
    for idx in 1:length(x)
        tmp = x[idx]
        for kdx in 1:lengths[idx]
            res[i] = tmp
            i += 1
        end
    end
    res
end

function _uniqueofsorted(x::Vector)
    idx = fill(true, length(x))
    lastx = x[1]
    for i = 2:length(x)
        if lastx == x[i]
            idx[i] = false
        else
            lastx = x[i]
        end
    end
    x[idx]
end

function make_unique{S<:ByteString}(names::Vector{S})
    x = Index()
    names = copy(names)
    dups = Int[]
    for i in 1:length(names)
        if has(x, names[i])
            push!(dups, i)
        else
            push!(x, names[i])
        end
    end
    for i in dups
        nm = names[i]
        newnm = nm
        k = 1
        while true
            newnm = "$(nm)_$k"
            if !has(x, newnm)
                push!(x, newnm)
                break
            end
            k += 1
        end
        names[i] = newnm
    end
    names
end

# TODO: move to set.jl? call nointer or nodupes?
# reasonably fast approach: foreach argument, iterate over
# the (presumed) set, checking for duplicates, adding to a hash table as we go
function nointer(ss...)
    d = Dict{Any,Int}(0)
    for s in ss
        for item in s
            ct = get(d, item, 0)
            if ct == 0 # we're good, add it
                d[item] = 1
            else
                return false
            end
        end
    end
    return true
end

function concat{T1,T2}(v1::Vector{T1}, v2::Vector{T2})
    # concatenate vectors, converting to type Any if needed.
    if T1 == T2 && T1 != Any
        [v1, v2]
    else
        res = Array(Any, length(v1) + length(v2))
        res[1:length(v1)] = v1
        res[length(v1)+1 : length(v1)+length(v2)] = v2
        res
    end
end

function _same_set(a, b)
    # there are definitely MUCH faster ways of doing this
    length(a) == length(b) && all(sort(a) == sort(b))
end

function generate_column_names(n::Int)
    convert(Vector{ByteString}, map(i -> "x" * string(i), 1:n))
end

function repeat(dv::AbstractVector, n::Int)
    m = length(dv)
    if rem(n, m) != 0
        error("New size must be a multiple of input vector's length")
    end
    new_dv = similar(dv, n)
    for i in 0:(fld(n, m) - 1)
        for j in 1:m
            new_j = i * m + j
            new_dv[new_j] = dv[j]
        end
    end
    return new_dv
end

function ancestors(t::Type)
    a = {t}
    while t != Any
        t = super(t)
        push!(a, t)
    end
    return a
end

function common_ancestors(s::Type, t::Type)
    return filter(e -> contains(ancestors(s), e), ancestors(t))
end

earliest_common_ancestor(s::Type, t::Type) = first(common_ancestors(s, t))

# Need to do type inference
# Like earliest_common_ancestor, but ignores NA
function _dv_most_generic_type(vals)
    # iterate over vals tuple to find the most generic non-NA type
    toptype = None
    for i = 1:length(vals)
        if !isna(vals[i])
            toptype = promote_type(toptype, typeof(vals[i]))
        end
    end
    if !method_exists(baseval, (toptype, ))
        error("No baseval exists for type: $(toptype)")
    end
    return toptype
end


text:
module inference

using ergm.sampler, ergm.models
using Statistics
import StatsBase
export mcmc_mle

function mcmc_estimate_means(
    sampler_initial_state,
    model :: ExponentialFamily,
    estimation_steps,
    burn_in,
    sample_interval
    )
    # seed MC with one of our samples (may not be the best choice...)
    s = GibbsSampler(
        sampler_initial_state,
        model,
        burn_in,
        sample_interval
    )
    _, stats = sample(s, estimation_steps)
    mean_stats = mean(stats, dims=1)
    mean_stats
end

function mcmc_mle(
    observations, model :: ExponentialFamily,
    gradient_descent_steps, estimation_steps, learning_rate
    )

    α = learning_rate
    θ = initial_guess
    θs = [θ]
    Ls = []

    observation_stats = [get_stats(model.stats, o) for o ∈ observations]
    target_Es = mean(reduce(hcat, observation_stats), dims=2)

    # need better heuristic
    burn_in = 10
    sample_interval = 10

    for i ∈ 1:gradient_descent_steps
        # seed sampler with random observation
        G0 = StatsBase.sample(observations)
        Es = mcmc_estimate_means(
            G0,
            model,
            estimation_steps,
            burn_in,
            sample_interval
        )
        dθ = target_Es - Es
        θ .+= α * dθ
        L = undef
        push!(Ls, L)
    end

    Ls
end

end


text:
using Bombe
using Test

import RData

let og_array = RData.load("test/Example_Log_Likelihood_Array.RData")["x"]
    global ll_arr = copy(permutedims(og_array, [3, 1, 2]))
end
let og_weights = RData.load("test/weightMatrix.RData")["weightMatrix"]
    global r_weights = exp.(permutedims(reshape(og_weights, 500, 2, 32), [3, 1, 2]))
end
rel_eff = RData.load("test/Rel_Eff.RData")["rel_eff"]
r_loo = RData.load("test/Example_Loo.RData")["example_loo"]
with_rel_eff = psis_loo(ll_arr; rel_eff = rel_eff)

@testset "Bombe.jl" begin

end

max_stars_repo_name: ParadaCarleton/Bombe

max_issues_repo_name: ParadaCarleton/Bombe

max_forks_repo_name: ParadaCarleton/Bombe


text:
using Test

using Distances
using Clustering

@testset "kmedoids() (k-medoids)" begin

Random.seed!(34568)

d = 3
n = 200
k = 10

X = rand(d, n)
costs = pairwise(SqEuclidean(), X, dims=2)
@assert size(costs) == (n, n)

R = kmedoids(costs, k)
@test isa(R, KmedoidsResult)
@test length(R.medoids) == length(unique(R.medoids))
@test all(a -> 1 <= a <= k, R.assignments)
@test R.assignments[R.medoids] == collect(1:k) # Every medoid should belong to its own cluster
@test sum(R.counts) == n
@test R.acosts == costs[LinearIndices((n, n))[CartesianIndex.(R.medoids[R.assignments], 1:n)]]
@test isapprox(sum(R.acosts), R.totalcost)
@test R.converged

# this data set has three obvious groups:
# group 1: [1, 3, 4], values: [1, 2, 3]
# group 2: [2, 5, 7], values: [6, 7, 8]
# group 3: [6, 8, 9], values: [21, 20, 22]
#

X = reshape(map(Float64, [1, 6, 2, 3, 7, 21, 8, 20, 22]), 1, 9)
costs = pairwise(SqEuclidean(), X, dims=2)

R = kmedoids!(costs, [1, 2, 6])
@test isa(R, KmedoidsResult)
@test R.medoids == [3, 5, 6]
@test R.assignments == [1, 2, 1, 1, 2, 3, 2, 3, 3]
@test R.counts == [3, 3, 3]
@test R.acosts ≈ [1, 1, 0, 1, 0, 0, 1, 1, 1]
@test R.totalcost ≈ 6.0
@test R.converged

end


text:
include("utils.jl")

include("byrow_test.jl")
include("dare_test.jl")
include("dlyap_test.jl")
include("simula01_test.jl")


text:
# TODO: choose method based on problem size
defaultmethod() = :Logarithmic

mutable struct PWLData
    counter::Int
    PWLData() = new(0)
end

function initPWL!(m::JuMP.Model)
    if !haskey(m.ext, :PWL)
        m.ext[:PWL] = PWLData()
    end
    return nothing
end

const VarOrAff = Union{JuMP.Variable,JuMP.AffExpr}

function piecewiselinear(m::JuMP.Model, x::VarOrAff, d, f::Function; method=defaultmethod())
    initPWL!(m)
    fd = [f(xx) for xx in d]
    return piecewiselinear(m, x, d, fd; method=method)
end

piecewiselinear(m::JuMP.Model, x::VarOrAff, d, fd; method=defaultmethod()) =
    piecewiselinear(m, x, UnivariatePWLFunction(d, fd); method=method)

function piecewiselinear(m::JuMP.Model, x::VarOrAff, pwl::UnivariatePWLFunction; method=defaultmethod())
    initPWL!(m)
    counter = m.ext[:PWL].counter
    counter += 1
    m.ext[:PWL].counter = counter
    d = [_x[1] for _x in pwl.x]
    fd = pwl.z
    n = length(d)

    if n != length(fd)
        error("You provided a different number of breakpoints ($n) and function values at the breakpoints ($(length(fd)))")
    end
    if n == 0
        error("I don't know how to handle a piecewise linear function with no breakpoints")
    end

    z = JuMP.@variable(m, lowerbound=minimum(fd), upperbound=maximum(fd), basename="z_$counter")

    if n == 1
        JuMP.@constraint(m, x ==  d[1])
        JuMP.@constraint(m, z == fd[1])
        return z
    end

    if method == :Incremental
        δ = JuMP.@variable(m, [1:n], lowerbound=0, upperbound=1, basename="δ_$counter")
        y = JuMP.@variable(m, [1:n-1], Bin, basename="y_$counter")
        JuMP.@constraint(m, x ==  d[1] + sum(δ[i]*( d[i+1]- d[i]) for i in 1:n-1))
        JuMP.@constraint(m, z == fd[1] + sum(δ[i]*(fd[i+1]-fd[i]) for i in 1:n-1))
        for i in 1:n-1
            JuMP.@constraint(m, δ[i+1] ≤ y[i])
            JuMP.@constraint(m, y[i] ≤ δ[i])
        end
    elseif method == :MC
        x̂ = JuMP.@variable(m, [1:n-1],      basename="x̂_$counter")
        ẑ = JuMP.@variable(m, [1:n-1],      basename="ẑ_$counter")
        y = JuMP.@variable(m, [1:n-1], Bin, basename="y_$counter")
        JuMP.@constraint(m, sum(y) == 1)
        JuMP.@constraint(m, sum(x̂) == x)
        JuMP.@constraint(m, sum(ẑ) == z)
        Δ = [(fd[i+1]-fd[i])/(d[i+1]-d[i]) for i in 1:n-1]
        for i in 1:n-1
            JuMP.@constraints(m, begin
                x̂[i] ≥ d[i]  *y[i]
                x̂[i] ≤ d[i+1]*y[i]
                ẑ[i] == fd[i]*y[i] + Δ[i]*(x̂[i]-d[i]*y[i])
            end)
        end
    elseif method in (:DisaggLogarithmic,:DLog)
        γ = JuMP.@variable(m, [i=1:n,j=max(1,i-1):min(n-1,i)], lowerbound=0, upperbound=1)
        JuMP.@constraint(m, sum(γ) == 1)
        JuMP.@constraint(m, γ[1,1]* d[1] + sum((γ[i,i-1]+γ[i,i])* d[i] for i in 2:n-1) + γ[n,n-1]* d[n] == x)
        JuMP.@constraint(m, γ[1,1]*fd[1] + sum((γ[i,i-1]+γ[i,i])*fd[i] for i in 2:n-1) + γ[n,n-1]*fd[n] == z)
        r = ceil(Int, log2(n-1))
        H = reflected_gray_codes(r)
        y = JuMP.@variable(m, [1:r], Bin)
        for j in 1:r
            JuMP.@constraint(m, sum((γ[i,i]+γ[i+1,i])*H[i][j] for i in 1:(n-1)) == y[j])
        end
    else
        # V-formulation methods
        λ = JuMP.@variable(m, [1:n], lowerbound=0, upperbound=1, basename="λ_$counter")
        JuMP.@constraint(m, sum(λ) == 1)
        JuMP.@constraint(m, sum(λ[i]* d[i] for i in 1:n) == x)
        JuMP.@constraint(m, sum(λ[i]*fd[i] for i in 1:n) == z)

        if method in (:Logarithmic,:Log)
            sos2_logarithmic_formulation!(m, λ)
        elseif method in (:LogarithmicIB,:LogIB)
            sos2_logarithmic_IB_formulation!(m, λ)
        elseif method == :CC
            sos2_cc_formulation!(m, λ)
        elseif method in (:ZigZag,:ZZB)
            sos2_zigzag_formulation!(m, λ)
        elseif method in (:ZigZagInteger,:ZZI)
            sos2_zigzag_general_integer_formulation!(m, λ)
        elseif method == :GeneralizedCelaya
            sos2_generalized_celaya_formulation!(m, λ)
        elseif method == :SymmetricCelaya
            sos2_symmetric_celaya_formulation!(m, λ)
        elseif method == :SOS2
            JuMP.addSOS2(m, [i*λ[i] for i in 1:n])
        else
            error("Unrecognized method $method")
        end
    end

    return z
end

function sos2_cc_formulation!(m::JuMP.Model, λ)
    counter = m.ext[:PWL].counter
    n = length(λ)
    y = JuMP.@variable(m, [1:n-1], Bin, basename="y_$counter")
    JuMP.@constraint(m, sum(y) == 1)
    JuMP.@constraint(m, λ[1] ≤ y[1])
    for i in 2:n-1
        JuMP.@constraint(m, λ[i] ≤ y[i-1] + y[i])
    end
    JuMP.@constraint(m, λ[n] ≤ y[n-1])
    return nothing
end

# function sos2_mc_formulation!(m::JuMP.Model, λ) # not currently used
#     counter = m.ext[:PWL].counter
#     n = length(λ)
#     γ = JuMP.@variable(m, [1:n-1, 1:n],     basename="γ_$counter")
#     y = JuMP.@variable(m, [1:n-1],     Bin, basename="y_$counter")
#     JuMP.@constraint(m, sum(y) == 1)
#     JuMP.@constraint(m, sum(γ[i,:] for i in 1:n-1) .== λ)
#     for i in 1:n-1
#         JuMP.@constraint(m, γ[i,i] + γ[i,i+1] ≥ y[i])
#     end
#     return nothing
# end

function sos2_logarithmic_formulation!(m::JuMP.Model, λ)
    counter = m.ext[:PWL].counter
    n = length(λ)-1
    k = ceil(Int,log2(n))
    y = JuMP.@variable(m, [1:k], Bin, basename="y_$counter")
    sos2_encoding_constraints!(m, λ, y, reflected_gray_codes(k), unit_vector_hyperplanes(k))
    return nothing
end

function sos2_logarithmic_IB_formulation!(m::JuMP.Model, λ) # IB: independent branching
    counter = m.ext[:PWL].counter
    n = length(λ)-1
    k = ceil(Int,log2(n))
    y = JuMP.@variable(m, [1:k], Bin, basename="y_$counter")
    _H = reflected_gray_codes(k)
    d = length(_H)
    H = Dict(i => _H[i] for i in 1:d)
    H[0] = H[1]
    H[d+1] = H[d]
    for j in 1:k
        JuMP.@constraints(m, begin
            sum(λ[i] for i in 1:(n+1) if H[i-1][j] == H[i][j] == 1) ≤     y[j]
            sum(λ[i] for i in 1:(n+1) if H[i-1][j] == H[i][j] == 0) ≤ 1 - y[j]
        end)
    end
    return nothing
end

function sos2_zigzag_formulation!(m::JuMP.Model, λ)
    counter = m.ext[:PWL].counter
    n = length(λ)-1
    k = ceil(Int,log2(n))
    y = JuMP.@variable(m, [1:k], Bin, basename="y_$counter")
    sos2_encoding_constraints!(m, λ, y, zigzag_codes(k), zigzag_hyperplanes(k))
    return nothing
end

function sos2_zigzag_general_integer_formulation!(m::JuMP.Model, λ)
    counter = m.ext[:PWL].counter
    n = length(λ)-1
    k = ceil(Int,log2(n))
    # TODO: tighter upperbounds
    y = JuMP.@variable(m, [i=1:k], Int, lowerbound=0, upperbound=2^(k-i), basename="y_$counter")
    sos2_encoding_constraints!(m, λ, y, integer_zigzag_codes(k), unit_vector_hyperplanes(k))
    return nothing
end

function sos2_generalized_celaya_formulation!(m::JuMP.Model, λ)
    counter = m.ext[:PWL].counter
    n = length(λ)-1
    k = ceil(Int,log2(n))
    codes = generalized_celaya_codes(k)
    lb = [minimum(t[i] for t in codes) for i in 1:k]
    ub = [maximum(t[i] for t in codes) for i in 1:k]
    y = JuMP.@variable(m, [i=1:k], Int, lowerbound=lb[i], upperbound=ub[i], basename="y_$counter")
    sos2_encoding_constraints!(m, λ, y, codes, generalized_celaya_hyperplanes(k))
    return nothing
end

function sos2_symmetric_celaya_formulation!(m::JuMP.Model, λ)
    counter = m.ext[:PWL].counter
    n = length(λ)-1
    k = ceil(Int,log2(n))
    codes = symmetric_celaya_codes(k)
    lb = [minimum(t[i] for t in codes) for i in 1:k]
    ub = [maximum(t[i] for t in codes) for i in 1:k]
    y = JuMP.@variable(m, [i=1:k], Int, lowerbound=lb[i], upperbound=ub[i], basename="y_$counter")
    sos2_encoding_constraints!(m, λ, y, codes, symmetric_celaya_hyperplanes(k))
    return nothing
end

function sos2_encoding_constraints!(m, λ, y, h, B)
    n = length(λ)-1
    for b in B
        JuMP.@constraints(m, begin
            dot(b,h[1])*λ[1] + sum(min(dot(b,h[v]),dot(b,h[v-1]))*λ[v] for v in 2:n) + dot(b,h[n])*λ[n+1] ≤ dot(b,y)
            dot(b,h[1])*λ[1] + sum(max(dot(b,h[v]),dot(b,h[v-1]))*λ[v] for v in 2:n) + dot(b,h[n])*λ[n+1] ≥ dot(b,y)
        end)
    end
    return nothing
end

function reflected_gray_codes(k::Int)
    if k == 0
        return Vector{Int}[]
    elseif k == 1
        return [[0],[1]]
    else
        codes′ = reflected_gray_codes(k-1)
        return vcat([vcat(code,0) for code in codes′],
                    [vcat(code,1) for code in reverse(codes′)])
    end
end

function zigzag_codes(k::Int)
    if k <= 0
        error("Invalid code length $k")
    elseif k == 1
        return [[0],[1]]
    else
        codes′ = zigzag_codes(k-1)
        return vcat([vcat(code,0) for code in codes′],
                    [vcat(code,1) for code in codes′])
    end
end

function integer_zigzag_codes(k::Int)
    if k <= 0
        error("Invalid code length $k")
    elseif k == 1
        return [[0],[1]]
    else
        codes′ = integer_zigzag_codes(k-1)
        offset = [2^(j-2) for j in k:-1:2]
        return vcat([vcat(code,        0) for code in codes′],
                    [vcat(code.+offset,1) for code in codes′])
    end
end

function generalized_celaya_codes(k::Int)
    if k <= 0
        error("Invalid code length $k")
    elseif k == 1
        return [[0],[1]]
    elseif k == 2
        return [[0,0],[1,1],[1,0],[0,1]]
    else
        codes′ = generalized_celaya_codes(k-1)
        n = length(codes′)
        hp = Int(n/2)
        firstcodes  = [codes′[i] for i in 1:hp]
        secondcodes = [codes′[i] for i in (hp+1):n]
        return vcat([vcat(codes′[i], 0) for i in     1 :  n],
                    [vcat(codes′[i], 1) for i in     1 : hp],
                    [vcat(codes′[i],-1) for i in (hp+1):  n])
    end
end

function symmetric_celaya_codes(k::Int)
    if k <= 0
        error("Invalid code length $k")
    elseif k == 1
        return [[0],[1]]
    else
        codes′ = generalized_celaya_codes(k-1)
        n = length(codes′)
        hp = Int(n/2)
        firstcodes  = [codes′[i] for i in 1:hp]
        secondcodes = [codes′[i] for i in (hp+1):n]
        return vcat([vcat(codes′[i], 0) for i in     1 :  n],
                    [vcat(codes′[i], 1) for i in     1 : hp],
                    [vcat(codes′[i],-1) for i in (hp+1):  n])
    end
end

function zigzag_hyperplanes(k::Int)
    hps = Vector{Int}[]
    for i in 1:k
        hp = zeros(Int, k)
        hp[i] = 1
        for j in (i+1):k
            hp[j] = 2^(j-i-1)
        end
        push!(hps, hp)
    end
    return hps
end

function unit_vector_hyperplanes(k::Int)
    hps = Vector{Int}[]
    for i in 1:k
        hp = zeros(Int,k)
        hp[i] = 1
        push!(hps, hp)
    end
    return hps
end

generalized_celaya_hyperplanes(k::Int) = compute_hyperplanes(generalized_celaya_codes(k))

symmetric_celaya_hyperplanes(k::Int) = compute_hyperplanes(symmetric_celaya_codes(k))

function compute_hyperplanes(C::Vector{Vector{T}}) where T <: Number
    n = length(C)
    k = length(C[1])

    d = zeros(n-1,k)
    for i in 1:n-1
        d[i,:] = C[i+1] - C[i]
    end

    if k <= 1
        error("Cannot process codes of length $k")
    elseif k == 2
        @assert n == 4
        spanners = Vector{Float64}[]
        for i in 1:n-1
            v = canonical!(d[i,:])
            v = [v[2], -v[1]]
            push!(spanners, v)
        end
        indices = [1]
        approx12 = isapprox(spanners[1],spanners[2])
        approx13 = isapprox(spanners[1],spanners[3])
        approx23 = isapprox(spanners[2],spanners[3])
        if !approx12
            push!(indices,2)
        end
        if !approx13 && !(!approx12 && approx23)
            push!(indices,3)
        end
        return spanners[indices]
    end

    indices = [1,2]
    spanners = Vector{Float64}[]
    while !isempty(indices)
        if indices == [n-1]
            break
        end
        if rank(d[indices,:]) == length(indices) && length(indices) <= k-1
            if length(indices) == k-1
                nullsp = nullspace(d[indices,:])
                @assert size(nullsp,2) == 1
                v = vec(nullsp)
                push!(spanners, canonical!(v))
            end
            if indices[end] != n-1
                push!(indices, indices[end]+1)
            else
                pop!(indices)
                indices[end] += 1
            end
        else
            if indices[end] != n-1
                indices[end] += 1
            else
                pop!(indices)
                indices[end] += 1
            end
        end
    end

    keepers = [1]
    for i in 2:length(spanners)
        alreadyin = false
        for j in keepers
            if isapprox(spanners[i], spanners[j])
                alreadyin = true
                break
            end
        end
        if !alreadyin
            push!(keepers, i)
        end
    end

    return spanners[keepers]
end

function canonical!(v::Vector{Float64})
    normalize!(v)
    for j in 1:length(v)
        if abs(v[j]) < 1e-8
            v[j] = 0
        end
    end
    sgn = sign(v[findfirst([x!=0 for x in v])])
    for j in 1:length(v)
        if abs(v[j]) < 1e-8
            v[j] = 0
        else
            v[j] *= sgn
        end
    end
    return v
end

function optimal_IB_scheme!(m::JuMP.Model, λ, pwl, subsolver) # IB: independent branching
    m.ext[:OptimalIB] = Int[]

    if !haskey(m.ext, :OptimalIBCache)
        m.ext[:OptimalIBCache] = Dict()
    end

    T = pwl.T
    n = maximum(maximum(t) for t in T)
    J = 1:n
    E = fill(false, n, n)
    for t in T
        for i in t, j in t
            E[i,j] = true
        end
    end

    if haskey(m.ext[:OptimalIBCache], pwl.T)
        xx, yy = m.ext[:OptimalIBCache][pwl.T]
        t = size(xx,1)
    else
        t = ceil(Int, log2(length(T)))
        xx = Array{Float64,2}(undef,0,0)
        yy = Array{Float64,2}(undef,0,0)

        while true
            model = JuMP.Model(solver=subsolver)
            JuMP.@variable(model, x[1:t,1:n], Bin)
            JuMP.@variable(model, y[1:t,1:n], Bin)
            JuMP.@variable(model, z[1:t,1:n,1:n], Bin)

            for j in 1:t
                for r in J, s in J
                    if r >= s
                        continue
                    end
                    JuMP.@constraints(model, begin
                        z[j,r,s] <= x[j,r] + x[j,s]
                        z[j,r,s] <= x[j,r] + y[j,r]
                        z[j,r,s] <= x[j,s] + y[j,s]
                        z[j,r,s] <= y[j,r] + y[j,s]
                        z[j,r,s] >= x[j,r] + y[j,s] - 1
                        z[j,r,s] >= x[j,s] + y[j,r] - 1
                    end)
                end
                for r in J
                    JuMP.@constraint(model, x[j,r] + y[j,r] <= 1)
                end
            end

            for r in J, s in J
                if r >= s
                    continue
                end
                if E[r,s]
                    JuMP.@constraint(model, sum(z[j,r,s] for j in 1:t) == 0)
                else
                    JuMP.@constraint(model, sum(z[j,r,s] for j in 1:t) >= 1)
                end
            end

            JuMP.@objective(model, Min, sum(x) + sum(y))
            stat = JuMP.solve(model)
            xx = JuMP.getvalue(x)
            yy = JuMP.getvalue(y)
            if any(isnan, xx) || any(isnan, yy)
                t += 1
            else
                break
            end
        end
        m.ext[:OptimalIBCache][pwl.T] = (xx,yy)
    end

    y = JuMP.@variable(m, [1:t], Bin)

    dˣ = [_x[1] for _x in pwl.x]
    dʸ = [_x[2] for _x in pwl.x]
    uˣ, uʸ = unique(dˣ), unique(dʸ)
    @assert issorted(uˣ)
    @assert issorted(uʸ)
    nˣ, nʸ = length(uˣ), length(uʸ)
    ˣtoⁱ = Dict(uˣ[i] => i for i in 1:nˣ)
    ʸtoʲ = Dict(uʸ[i] => i for i in 1:nʸ)

    for i in 1:t
        JuMP.@constraints(m, begin
            sum(λ[ˣtoⁱ[pwl.x[j][1]],ʸtoʲ[pwl.x[j][2]]] for j in J if xx[i,j] == 1) ≤     y[i]
            sum(λ[ˣtoⁱ[pwl.x[j][1]],ʸtoʲ[pwl.x[j][2]]] for j in J if yy[i,j] == 1) ≤ 1 - y[i]
        end)
    end
    push!(m.ext[:OptimalIB], t)

    return nothing
end

piecewiselinear(m::JuMP.Model, x::VarOrAff, y::VarOrAff, dˣ, dʸ, f::Function; method=defaultmethod()) =
    piecewiselinear(m, x, y, BivariatePWLFunction(dˣ, dʸ, f); method=method)

function piecewiselinear(m::JuMP.Model, x₁::VarOrAff, x₂::VarOrAff, pwl::BivariatePWLFunction; method=defaultmethod(), subsolver=nothing)
    if (method == :OptimalIB) && (subsolver === nothing)
        error("No MIP solver provided to construct optimal IB scheme. Pass a solver object to the piecewiselinear function, e.g. piecewiselinear(m, x₁, x₂, bivariatefunc, method=:OptimalIB, subsolver=GurobiSolver())")
    end

    initPWL!(m)
    counter = m.ext[:PWL].counter
    counter += 1
    m.ext[:PWL].counter = counter
    dˣ = [_x[1] for _x in pwl.x]
    dʸ = [_x[2] for _x in pwl.x]
    uˣ, uʸ = unique(dˣ), unique(dʸ)
    @assert issorted(uˣ)
    @assert issorted(uʸ)

    T = pwl.T

    nˣ, nʸ = length(uˣ), length(uʸ)

    if nˣ == 0 || nʸ == 0
        error("I don't know how to handle a piecewise linear function with zero breakpoints")
    elseif nˣ == 1
        @assert length(dʸ) == nʸ == length(pwl.z)
        return piecewiselinear(m, x₂, UnivariatePWLFunction(dʸ, [pwl.z[i] for i in 1:nʸ]))
    elseif nʸ == 1
        @assert length(dˣ) == nˣ == length(pwl.z)
        return piecewiselinear(m, x₁, UnivariatePWLFunction(dˣ, [pwl.z[i] for i in 1:nˣ]))
    end

    ˣtoⁱ = Dict(uˣ[i] => i for i in 1:nˣ)
    ʸtoʲ = Dict(uʸ[i] => i for i in 1:nʸ)

    fd = Array{Float64}(undef, nˣ, nʸ)
    for (v,fv) in zip(pwl.x, pwl.z)
        # i is the linear index into pwl.x...really want (i,j) pair
        fd[ˣtoⁱ[v[1]],ʸtoʲ[v[2]]] = fv
    end

    z = JuMP.@variable(m, lowerbound=minimum(fd), upperbound=maximum(fd), basename="z_$counter")

    if method == :MC
        x̂₁ = JuMP.@variable(m, [T],      basename="x̂₁_$counter")
        x̂₂ = JuMP.@variable(m, [T],      basename="x̂₂_$counter")
        ẑ  = JuMP.@variable(m, [T],      basename="ẑ_$counter")
        y  = JuMP.@variable(m, [T], Bin, basename="y_$counter")
        JuMP.@constraint(m, sum(y)  == 1)
        JuMP.@constraint(m, sum(x̂₁) == x₁)
        JuMP.@constraint(m, sum(x̂₂) == x₂)
        JuMP.@constraint(m, sum(ẑ)  == z)
        for t in T
            @assert length(t) == 3
            r¹,  r²,  r³  = pwl.x[t[1]], pwl.x[t[2]], pwl.x[t[3]]
            fz¹, fz², fz³ = pwl.z[t[1]], pwl.z[t[2]], pwl.z[t[3]]
            for P in [[1,2,3],[2,3,1],[3,1,2]]
                p¹, p², p³ = [r¹, r², r³][P]
                p¹, p², p³ = pwl.x[t[P[1]]], pwl.x[t[P[2]]], pwl.x[t[P[3]]]

                A = [p¹[1] p¹[2] 1
                     p²[1] p²[2] 1
                     p³[1] p³[2] 1]
                b = [0,0,1]
                q = A \ b
                @assert isapprox(q[1]*p¹[1] + q[2]*p¹[2] + q[3], 0, atol=1e-4)
                @assert isapprox(q[1]*p²[1] + q[2]*p²[2] + q[3], 0, atol=1e-4)
                @assert isapprox(q[1]*p³[1] + q[2]*p³[2] + q[3], 1, atol=1e-4)
                JuMP.@constraint(m, q[1]*x̂₁[t] + q[2]*x̂₂[t] + q[3]*y[t] ≥ 0)
            end
            A = [r¹[1] r¹[2] 1
                 r²[1] r²[2] 1
                 r³[1] r³[2] 1]
            b = [fz¹, fz², fz³]
            q = A \ b
            @assert isapprox(q[1]*r¹[1] + q[2]*r¹[2] + q[3], fz¹, atol=1e-4)
            @assert isapprox(q[1]*r²[1] + q[2]*r²[2] + q[3], fz², atol=1e-4)
            @assert isapprox(q[1]*r³[1] + q[2]*r³[2] + q[3], fz³, atol=1e-4)
            JuMP.@constraint(m, ẑ[t] == q[1]*x̂₁[t] + q[2]*x̂₂[t] + q[3]*y[t])
        end
    elseif method in (:DisaggLogarithmic,:DLog)
        T = pwl.T
        X = pwl.x
        Z = pwl.z
        n = length(X)
        γ = JuMP.@variable(m, [t=T,v=t], lowerbound=0, upperbound=1, basename="γ_$counter")
        Tv = Dict(v => Any[] for v in 1:n)
        for t in T, v in t
            push!(Tv[v], t)
        end
        JuMP.@constraint(m, sum(γ) == 1)
        JuMP.@constraint(m, sum(sum(γ[t,i] for t in Tv[i]) * X[i][1] for i in 1:n) == x₁)
        JuMP.@constraint(m, sum(sum(γ[t,i] for t in Tv[i]) * X[i][2] for i in 1:n) == x₂)
        JuMP.@constraint(m, sum(sum(γ[t,i] for t in Tv[i]) * Z[i]    for i in 1:n) == z)
        r = ceil(Int, log2(length(T)))
        H = reflected_gray_codes(r)
        y = JuMP.@variable(m, [1:r], Bin, basename="y_$counter")
        for j in 1:r
            JuMP.@constraint(m, sum(sum(γ[T[i],v] for v in T[i])*H[i][j] for i in 1:length(T)) == y[j])
        end
    else
        # V-formulation methods
        λ = JuMP.@variable(m, [1:nˣ,1:nʸ], lowerbound=0, upperbound=1, basename="λ_$counter")
        JuMP.@constraint(m, sum(λ) == 1)
        JuMP.@constraint(m, sum(λ[i,j]*uˣ[i]   for i in 1:nˣ, j in 1:nʸ) == x₁)
        JuMP.@constraint(m, sum(λ[i,j]*uʸ[j]   for i in 1:nˣ, j in 1:nʸ) == x₂)
        JuMP.@constraint(m, sum(λ[i,j]*fd[i,j] for i in 1:nˣ, j in 1:nʸ) == z)

        if method == :CC
            T = pwl.T
            y = JuMP.@variable(m, [T], Bin, basename="y_$counter")
            JuMP.@constraint(m, sum(y) == 1)
            Ts = Dict((i,j) => Vector{Int}[] for i in 1:nˣ, j in 1:nʸ)
            for t in T, ind in t
                rˣ, rʸ = pwl.x[ind]
                i, j = ˣtoⁱ[rˣ], ʸtoʲ[rʸ]
                push!(Ts[(i,j)], t)
            end
            for i in 1:nˣ, j in 1:nʸ
                JuMP.@constraint(m, λ[i,j] ≤ sum(y[t] for t in Ts[(i,j)]))
            end
        elseif method == :Incremental
            error("Incremental formulation for bivariate functions is not currently implemented.")
            # TODO: implement algorithm of Geissler et al. (2012)
        elseif method == :OptimalIB
            optimal_IB_scheme!(m, λ, pwl, subsolver)
        else
            # formulations with SOS2 along each dimension
            Tx = [sum(λ[tˣ,tʸ] for tˣ in 1:nˣ) for tʸ in 1:nʸ]
            Ty = [sum(λ[tˣ,tʸ] for tʸ in 1:nʸ) for tˣ in 1:nˣ]

            if method in (:Logarithmic,:Log)
                sos2_logarithmic_formulation!(m, Tx)
                sos2_logarithmic_formulation!(m, Ty)
            elseif method in (:LogarithmicIB,:LogIB)
                sos2_logarithmic_IB_formulation!(m, Tx)
                sos2_logarithmic_IB_formulation!(m, Ty)
            elseif method in (:ZigZag,:ZZB)
                sos2_zigzag_formulation!(m, Tx)
                sos2_zigzag_formulation!(m, Ty)
            elseif method in (:ZigZagInteger,:ZZI)
                sos2_zigzag_general_integer_formulation!(m, Tx)
                sos2_zigzag_general_integer_formulation!(m, Ty)
            elseif method == :GeneralizedCelaya
                sos2_generalized_celaya_formulation!(m, Tx)
                sos2_generalized_celaya_formulation!(m, Ty)
            elseif method == :SymmetricCelaya
                sos2_symmetric_celaya_formulation!(m, Tx)
                sos2_symmetric_celaya_formulation!(m, Ty)
            elseif method == :SOS2
                γˣ = JuMP.@variable(m, [1:nˣ], lowerbound=0, upperbound=1, basename="γˣ_$counter")
                γʸ = JuMP.@variable(m, [1:nʸ], lowerbound=0, upperbound=1, basename="γʸ_$counter")
                JuMP.@constraint(m, [tˣ in 1:nˣ], γˣ[tˣ] == sum(λ[tˣ,tʸ] for tʸ in 1:nʸ))
                JuMP.@constraint(m, [tʸ in 1:nʸ], γʸ[tʸ] == sum(λ[tˣ,tʸ] for tˣ in 1:nˣ))
                JuMP.addSOS2(m, [i*γˣ[i] for i in 1:nˣ])
                JuMP.addSOS2(m, [i*γʸ[i] for i in 1:nʸ])
            else
                error("Unrecognized method $method")
            end

            pattern = pwl.meta[:structure]
            if pattern == :UnionJack
                numT = 0
                minˣ, minʸ = uˣ[1], uʸ[1]
                # find the index of the bottom-left point
                idx = findfirst(pwl.x) do w; w == (minˣ, minʸ); end
                for t in pwl.T
                    if idx in t
                        numT += 1
                    end
                end
                @assert 1 <= numT <= 2

                w = JuMP.@variable(m, category=:Bin, basename="w_$counter")
                # if numT==1, then bottom-left is contained in only one point, and so needs separating; otherwise numT==2, and need to offset by one
                JuMP.@constraints(m, begin
                    sum(λ[tx,ty] for tx in 1:2:nˣ, ty in    numT :2:nʸ) ≤     w
                    sum(λ[tx,ty] for tx in 2:2:nˣ, ty in (3-numT):2:nʸ) ≤ 1 - w
                end)
            elseif pattern == :K1
                # assumption is that the triangulation is uniform, as defined in types.jl
                w = JuMP.@variable(m, [1:2], Bin, basename="w_$counter")
                JuMP.@constraints(m, begin
                    sum(λ[tx,ty] for tx in 1:nˣ, ty in 1:nʸ if mod(tx,2) == mod(ty,2) && (tx+ty) in 2:4:(nˣ+nʸ)) ≤ w[1]
                    sum(λ[tx,ty] for tx in 1:nˣ, ty in 1:nʸ if mod(tx,2) == mod(ty,2) && (tx+ty) in 4:4:(nˣ+nʸ)) ≤ 1 - w[1]
                    sum(λ[tx,ty] for tx in 1:nˣ, ty in 1:nʸ if mod(tx,2) != mod(ty,2) && (tx+ty) in 3:4:(nˣ+nʸ)) ≤ w[2]
                    sum(λ[tx,ty] for tx in 1:nˣ, ty in 1:nʸ if mod(tx,2) != mod(ty,2) && (tx+ty) in 5:4:(nˣ+nʸ)) ≤ 1 - w[2]
                end)
            elseif pattern == :OptimalTriangleSelection
                m.ext[:OptimalTriSelect] = Int[]

                if !haskey(m.ext, :OptimalTriSelectCache)
                    m.ext[:OptimalTriSelectCache] = Dict()
                end

                J = [(i,j) for i in 1:nˣ, j in 1:nʸ]
                E = Set{Tuple{Tuple{Int,Int},Tuple{Int,Int}}}()
                for t in T
                    @assert length(t) == 3
                    IJ = [(ˣtoⁱ[pwl.x[i][1]],ʸtoʲ[pwl.x[i][2]]) for i in t]
                    im = minimum(ij[1] for ij in IJ)
                    iM = maximum(ij[1] for ij in IJ)
                    jm = minimum(ij[2] for ij in IJ)
                    jM = maximum(ij[2] for ij in IJ)
                    @assert im < iM
                    @assert im < iM
                    if ((im,jm) in IJ) && ((iM,jM) in IJ)
                        push!(E, ((im,jM),(iM,jm)))
                    elseif ((im,jM) in IJ) && ((iM,jm) in IJ)
                        push!(E, ((im,jm),(iM,jM)))
                    else
                        error()
                    end
                end

                if haskey(m.ext[:OptimalTriSelectCache], E)
                    xx, yy = m.ext[:OptimalTriSelectCache][E]
                    t = JuMP.size(xx,1)
                else
                    if subsolver === nothing
                        error("No MIP solver provided to construct optimal triangle selection. Pass a solver object to the piecewiselinear function, e.g. piecewiselinear(m, x₁, x₂, bivariatefunc, method=:Logarithmic, subsolver=GurobiSolver())")
                    end
                    t = 1
                    xx, yy = Array(Float64,0,0), Array(Float64,0,0)
                    while true
                        subm = JuMP.Model(solver=subsolver)
                        JuMP.@variable(subm, xˢᵘᵇ[1:t,J], Bin)
                        JuMP.@variable(subm, yˢᵘᵇ[1:t,J], Bin)
                        JuMP.@variable(subm, zˢᵘᵇ[1:t,J,J], Bin)

                        for j in 1:t
                            for r in J, s in J
                                # lexicographic ordering on points on grid
                                if r[1] > s[1] || (r[1] == s[1] && r[2] >= s[2])
                                    continue
                                end
                                JuMP.@constraints(subm, begin
                                    zˢᵘᵇ[j,r,s] <= xˢᵘᵇ[j,r] + xˢᵘᵇ[j,s]
                                    zˢᵘᵇ[j,r,s] <= xˢᵘᵇ[j,r] + yˢᵘᵇ[j,r]
                                    zˢᵘᵇ[j,r,s] <= xˢᵘᵇ[j,s] + yˢᵘᵇ[j,s]
                                    zˢᵘᵇ[j,r,s] <= yˢᵘᵇ[j,r] + yˢᵘᵇ[j,s]
                                    zˢᵘᵇ[j,r,s] >= xˢᵘᵇ[j,r] + yˢᵘᵇ[j,s] - 1
                                    zˢᵘᵇ[j,r,s] >= xˢᵘᵇ[j,s] + yˢᵘᵇ[j,r] - 1
                                end)
                            end
                            for r in J
                                JuMP.@constraint(subm, xˢᵘᵇ[j,r] + yˢᵘᵇ[j,r] <= 1)
                            end
                        end

                        for r in J, s in J
                            # lexicographic ordering on points on grid
                            (r[1] > s[1] || (r[1] == s[1] && r[2] >= s[2])) && continue
                            if (r,s) in E
                                JuMP.@constraint(subm, sum(zˢᵘᵇ[j,r,s] for j in 1:t) >= 1)
                            elseif max(abs(r[1]-s[1]), abs(r[2]-s[2])) == 1
                                JuMP.@constraint(subm, sum(zˢᵘᵇ[j,r,s] for j in 1:t) == 0)
                            end
                        end

                        JuMP.@objective(subm, Min, sum(xˢᵘᵇ) + sum(yˢᵘᵇ))
                        stat = JuMP.solve(subm)
                        if any(isnan, subm.colVal)
                            t += 1
                        else
                            xx = JuMP.getvalue(xˢᵘᵇ)
                            yy = JuMP.getvalue(yˢᵘᵇ)
                            m.ext[:OptimalTriSelectCache][E] = (xx,yy)
                            break
                        end
                    end
                end
                y = JuMP.@variable(m, [1:t], Bin, basename="Δselect_$counter")

                for i in 1:t
                    JuMP.@constraints(m, begin
                        sum(λ[v[1],v[2]] for v in J if xx[i,v] ≈ 1) ≤     y[i]
                        sum(λ[v[1],v[2]] for v in J if yy[i,v] ≈ 1) ≤ 1 - y[i]
                    end)
                end
                push!(m.ext[:OptimalTriSelect], t)
            elseif pattern in (:Stencil,:Stencil9)
                w = JuMP.@variable(m, [1:3,1:3], Bin, basename="w_$counter")
                for oˣ in 1:3, oʸ in 1:3
                    innoT = fill(true, nˣ, nʸ)
                    for (i,j,k) in pwl.T
                        xⁱ, xʲ, xᵏ = pwl.x[i], pwl.x[j], pwl.x[k]
                        iiˣ, iiʸ = ˣtoⁱ[xⁱ[1]], ʸtoʲ[xⁱ[2]]
                        jjˣ, jjʸ = ˣtoⁱ[xʲ[1]], ʸtoʲ[xʲ[2]]
                        kkˣ, kkʸ = ˣtoⁱ[xᵏ[1]], ʸtoʲ[xᵏ[2]]
                        # check to see if one of the points in the triangle falls on the grid
                        if (mod1(iiˣ,3) == oˣ && mod1(iiʸ,3) == oʸ) || (mod1(jjˣ,3) == oˣ && mod1(jjʸ,3) == oʸ) || (mod1(kkˣ,3) == oˣ && mod1(kkʸ,3) == oʸ)
                            innoT[iiˣ,iiʸ] = false
                            innoT[jjˣ,jjʸ] = false
                            innoT[kkˣ,kkʸ] = false
                        end
                    end
                    JuMP.@constraints(m, begin
                        sum(λ[i,j] for i in oˣ:3:nˣ, j in oʸ:3:nʸ) ≤  1 - w[oˣ,oʸ]
                        sum(λ[i,j] for i in 1:nˣ, j in 1:nʸ if innoT[i,j]) ≤ w[oˣ,oʸ]
                    end)
                end
            else
                @assert pattern in (:Upper, :Lower, :BestFit, :Random)
                # Eⁿᵉ[i,j] = true means that we must cover the edge {(i,j),(i+1,j+1)}
                Eⁿᵉ = fill(false, nˣ-1, nʸ-1)
                for (i,j,k) in pwl.T
                    xⁱ, xʲ, xᵏ = pwl.x[i], pwl.x[j], pwl.x[k]
                    iiˣ, iiʸ = ˣtoⁱ[xⁱ[1]], ʸtoʲ[xⁱ[2]]
                    jjˣ, jjʸ = ˣtoⁱ[xʲ[1]], ʸtoʲ[xʲ[2]]
                    kkˣ, kkʸ = ˣtoⁱ[xᵏ[1]], ʸtoʲ[xᵏ[2]]
                    IJ = [(iiˣ,iiʸ), (jjˣ,jjʸ), (kkˣ,kkʸ)]
                    im = min(iiˣ, jjˣ, kkˣ)
                    iM = max(iiˣ, jjˣ, kkˣ)
                    jm = min(iiʸ, jjʸ, kkʸ)
                    jM = max(iiʸ, jjʸ, kkʸ)
                    if ((im,jM) in IJ) && ((iM,jm) in IJ)
                        Eⁿᵉ[im,jm] = true
                    else
                        @assert (im,jm) in IJ && (iM,jM) in IJ
                    end
                end

                # diagonal lines running from SW to NE. Grouped with an offset of 3.
                wⁿᵉ = JuMP.@variable(m, [0:2], Bin, basename="wⁿᵉ_$counter")
                for o in 0:2
                    Aᵒ = Set{Tuple{Int,Int}}()
                    Bᵒ = Set{Tuple{Int,Int}}()
                    for offˣ in o:3:(nˣ-2)
                        SWinA = true # whether we put the SW corner of the next triangle to cover in set A
                        for i in (1+offˣ):(nˣ-1)
                            j = i - offˣ
                            if !(1 ≤ i ≤ nˣ-1)
                                continue
                            end
                            if !(1 ≤ j ≤ nʸ-1)
                                continue # should never happen
                            end
                            if Eⁿᵉ[i,j] # if we need to cover the edge...
                                if SWinA # figure out which set we need to put it in; this depends on previous triangle in our current line
                                    push!(Aᵒ, (i  ,j  ))
                                    push!(Bᵒ, (i+1,j+1))
                                else
                                    push!(Aᵒ, (i+1,j+1))
                                    push!(Bᵒ, (i  ,j  ))
                                end
                                SWinA = !SWinA
                            end
                        end
                    end
                    for offʸ in (3-o):3:(nʸ-1)
                        SWinA = true
                        for j in (offʸ+1):(nʸ-1)
                            i = j - offʸ
                            if !(1 ≤ i ≤ nˣ-1)
                                continue
                            end
                            if Eⁿᵉ[i,j]
                                if SWinA
                                    push!(Aᵒ, (i  ,j  ))
                                    push!(Bᵒ, (i+1,j+1))
                                else
                                    push!(Aᵒ, (i+1,j+1))
                                    push!(Bᵒ, (i  ,j  ))
                                end
                                SWinA = !SWinA
                            end
                        end
                    end
                    JuMP.@constraints(m, begin
                        sum(λ[i,j] for (i,j) in Aᵒ) ≤     wⁿᵉ[o]
                        sum(λ[i,j] for (i,j) in Bᵒ) ≤ 1 - wⁿᵉ[o]
                    end)
                end

                wˢᵉ = JuMP.@variable(m, [0:2], Bin, basename="wˢᵉ_$counter")
                for o in 0:2
                    Aᵒ = Set{Tuple{Int,Int}}()
                    Bᵒ = Set{Tuple{Int,Int}}()
                    for offˣ in o:3:(nˣ-2)
                        SEinA = true
                        # for i in (1+offˣ):-1:1
                            # j = offˣ - i + 2
                        for j in 1:(nʸ-1)
                            i = nˣ - j - offˣ
                            if !(1 ≤ i ≤ nˣ-1)
                                continue
                            end
                            if !Eⁿᵉ[i,j]
                                if SEinA
                                    push!(Aᵒ, (i+1,j  ))
                                    push!(Bᵒ, (i  ,j+1))
                                else
                                    push!(Aᵒ, (i  ,j+1))
                                    push!(Bᵒ, (i+1,j  ))
                                end
                                SEinA = !SEinA
                            end
                        end
                    end
                    for offʸ in (3-o):3:(nʸ-1)
                        SEinA = true
                        for j in (offʸ+1):(nʸ-1)
                            i = nˣ - j + offʸ
                            if !(1 ≤ i ≤ nˣ-1)
                                continue
                            end
                            if !Eⁿᵉ[i,j]
                                if SEinA
                                    push!(Aᵒ, (i+1,j  ))
                                    push!(Bᵒ, (i  ,j+1))
                                else
                                    push!(Aᵒ, (i  ,j+1))
                                    push!(Bᵒ, (i+1,j  ))
                                end
                                SEinA = !SEinA
                            end
                        end
                    end
                    JuMP.@constraints(m, begin
                        sum(λ[i,j] for (i,j) in Aᵒ) ≤     wˢᵉ[o]
                        sum(λ[i,j] for (i,j) in Bᵒ) ≤ 1 - wˢᵉ[o]
                    end)
                end
            end
        end
    end

    return z
end


text:
# ---
# title: 1458. Max Dot Product of Two Subsequences
# id: problem1458
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/max-dot-product-of-two-subsequences/description/>
# hidden: true
# ---
#
# Given two arrays `nums1` and `nums2`.
#
# Return the maximum dot product between  **non-empty** subsequences of nums1
# and nums2 with the same length.
#
# A subsequence of a array is a new array which is formed from the original
# array by deleting some (can be none) of the characters without disturbing the
# relative positions of the remaining characters. (ie, `[2,3,5]` is a
# subsequence of `[1,2,3,4,5]` while `[1,5,3]` is not).
#
#
#
# **Example 1:**
#
#
#
#     Input: nums1 = [2,1,-2,5], nums2 = [3,0,-6]
#     Output: 18
#     Explanation: Take subsequence [2,-2] from nums1 and subsequence [3,-6] from nums2.
#     Their dot product is (2*3 + (-2)*(-6)) = 18.
#
# **Example 2:**
#
#
#
#     Input: nums1 = [3,-2], nums2 = [2,-6,7]
#     Output: 21
#     Explanation: Take subsequence [3] from nums1 and subsequence [7] from nums2.
#     Their dot product is (3*7) = 21.
#
# **Example 3:**
#
#
#
#     Input: nums1 = [-1,-1], nums2 = [1,1]
#     Output: -1
#     Explanation: Take subsequence [-1] from nums1 and subsequence [1] from nums2.
#     Their dot product is -1.
#
#
#
# **Constraints:**
#
#   * `1 <= nums1.length, nums2.length <= 500`
#   * `-1000 <= nums1[i], nums2[i] <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
###############################################################################
#
# Nelder-Mead method to find peak
#
###############################################################################
function find_peak(
            f::F,
            ms::NelderMeadMethod{FT},
            tol::SolutionToleranceND{FT}
) where {F<:Function, FT<:AbstractFloat}
    # unpack values
    @unpack cen_x, con_x, exp_x, ref_x, x_inis = ms;
    nX = ms.N;
    nS = ms.N + 1;

    # initialize the simplex
    for irow in 1:nS
        ms.simplex[irow] .= x_inis
        if irow < nS
            ms.simplex[irow][irow] = x_inis[irow] / 2;
        end
        ms.simplex[irow][end] = f(ms.simplex[irow])
    end

    # find the optimal simplex
    count_all = 0;
    judges = [false for i in 1:nX];
    while true
        # 1. sort the simplex from high to low
        sort!(ms.simplex, by=x->x[end], rev=true);
        for icol in 1:nX
            judges[icol] = ( abs(ms.simplex[1][icol] - ms.simplex[nS][icol]) <= tol.tol[icol]);
        end
        all(judges) ? break : nothing;

        # 2. calculate the centroid of the N best
        cen_x .= FT(0)
        for irow in 1:nX
            cen_x .+= ms.simplex[irow]
        end
        cen_x ./= nX;

        # 3. reflection of the worst
        ref_x  .= cen_x;
        ref_x .*= 2;
        ref_x .-= ms.simplex[nS];
        ref_y   = f(ref_x);

        if ms.simplex[1][end] >= ref_y > ms.simplex[nX][end]
            ms.simplex[end]     .= ref_x;
            ms.simplex[end][end] = ref_y;

        # 4. expansion of the reflection
        elseif ref_y > ms.simplex[1][end]
            exp_x  .= ref_x;
            exp_x .*= 2;
            exp_x .-= cen_x;
            exp_y   = f(exp_x);

            if exp_y > ref_y
                ms.simplex[end]     .= exp_x;
                ms.simplex[end][end] = exp_y;
            else
                ms.simplex[end]     .= ref_x;
                ms.simplex[end][end] = ref_y;
            end

        # 5. contraction of the worst
        else
            con_x  .= cen_x;
            con_x .+= ms.simplex[nS];
            con_x ./= 2;
            con_y   = f(con_x);

            if con_y > ms.simplex[nS][end]
                ms.simplex[end]     .= con_x;
                ms.simplex[end][end] = con_y;

            # 6. shrink
            else
                for irow in 2:nS
                    ms.simplex[irow]    .+= ms.simplex[1];
                    ms.simplex[irow][end] = f(ms.simplex[irow]);
                end
            end
        end

        # 7. iteration ++
        count_all += 1;

        # 8. determine whether to break
        count_all > tol.n_limit ? break : nothing;
    end

    return view(ms.simplex[1], 1:nX)
end

function find_peak(
            f::F,
            ms::NelderMeadMethod{FT},
            tol::ResidualTolerance{FT}
) where {F<:Function, FT<:AbstractFloat}
    # unpack values
    @unpack cen_x, con_x, exp_x, ref_x, x_inis = ms;
    nX = ms.N;
    nS = ms.N + 1;

    # initialize the simplex
    for irow in 1:nS
        ms.simplex[irow] .= x_inis
        if irow < nS
            ms.simplex[irow][irow] = x_inis[irow] / 2;
        end
        ms.simplex[irow][end] = f(ms.simplex[irow])
    end

    # find the optimal simplex
    count_all = 0;
    while true
        # 1. sort the simplex from high to low
        sort!(ms.simplex, by=x->x[end], rev=true);
        if (ms.simplex[1][end] - ms.simplex[end][end]) < tol.tol
            break
        end

        # 2. calculate the centroid of the N best
        cen_x .= FT(0)
        for irow in 1:nX
            cen_x .+= ms.simplex[irow]
        end
        cen_x ./= nX;

        # 3. reflection of the worst
        ref_x  .= cen_x;
        ref_x .*= 2;
        ref_x .-= ms.simplex[nS];
        ref_y   = f(ref_x);

        if ms.simplex[1][end] >= ref_y > ms.simplex[nX][end]
            ms.simplex[end]     .= ref_x;
            ms.simplex[end][end] = ref_y;

        # 4. expansion of the reflection
        elseif ref_y > ms.simplex[1][end]
            exp_x  .= ref_x;
            exp_x .*= 2;
            exp_x .-= cen_x;
            exp_y   = f(exp_x);

            if exp_y > ref_y
                ms.simplex[end]     .= exp_x;
                ms.simplex[end][end] = exp_y;
            else
                ms.simplex[end]     .= ref_x;
                ms.simplex[end][end] = ref_y;
            end

        # 5. contraction of the worst
        else
            con_x  .= cen_x;
            con_x .+= ms.simplex[nS];
            con_x ./= 2;
            con_y   = f(con_x);

            if con_y > ms.simplex[nS][end]
                ms.simplex[end]     .= con_x;
                ms.simplex[end][end] = con_y;

            # 6. shrink
            else
                for irow in 2:nS
                    ms.simplex[irow]    .+= ms.simplex[1];
                    ms.simplex[irow][end] = f(ms.simplex[irow]);
                end
            end
        end

        # 7. iteration ++
        count_all += 1;

        # 8. determine whether to break
        count_all > tol.n_limit ? break : nothing;
    end

    return view(ms.simplex[1], 1:nX)
end


text:
include("Granular_Module.jl")
include("Fast_Hertz.jl")
include("Min_Energy.jl")

using JuLIP
using NeighbourLists
using Main.Fast_Hertz
using Main.Minimise_Energy
using LinearAlgebra
using DelimitedFiles
################################
# get an initial configuration
################################
Grain_Num = 1024
α = 2.0
ϕ = 0.842

Jamming_Info = readdlm("jamming_config.txt")

x = Jamming_Info[:,1]'
y = Jamming_Info[:,2]'
z = zeros(Grain_Num)'
X = [x; y; z]
X = vecs(X)

grain_type = Jamming_Info[:,4] .+ 1
r = Jamming_Info[:,3]

Config = gran.setup_cell(Grain_Num, ϕ, α)
Config.X .= X
Config.M = r
Config.Z = grain_type

ϕ = π*sum(r.^2)
set_data!(Config, :Phi, ϕ)
set_data!(Config, :Rad, r)

##########################################
# minimise the potential energy
##########################################
@time Config1 = Main.Minimise_Energy.min_Energy(Config)


text:
module PropertyGraphs

using LightGraphs
using Bijections

export
    NoEdgeProperties,
    NoVertexProperties,
    SimplePropertyGraph,
    SimplePropertyDiGraph,
    PropertyGraph,
    label,
    labels,
    write_dot

include("property_graph.jl")
include("dot.jl")

end


text:
using Test, StorageColorTypes

@testset "rgb types" begin
    include("rgbtypes.jl")
end

@testset "HSV/HSL constructors" begin
    @testset "$C constructor" for C in (HSV24, HSL24, HLS24, HLS)
        @test C(100, 1.0, 0.5) isa Color{Float32, 3}
        @test C(1N0f8, 0.6N0f8, 0.0N0f8) === C(1, 0.6, 0)
    end
end

@testset "Lab/Luv constructors" begin
    @testset "$C constructor" for C in (Lab24, Luv24)
        @test C(50, 40, 30.0) isa Color{Float32, 3}
        @test C(1N0f8, 0.6N0f8, 0.0N0f8) === C(1, 0.6, 0)
    end
end


text:
import Distributions
# @testset "$(format("Model Fitting: choose[Gaussian v.s. Poisson]"))" begin
#     let
#         total_passed = 0
#         for i = 1:1000
#             sample = rand(Distributions.Gamma(rand(1:100),rand(1:100)),5000)
#             if choose(Gaussian(),Gamma(),data=sample) == :Gamma
#                 total_passed = total_passed + 1
#             end
#         end
#         @test total_passed/1000.0 > 0.9
#         @info @sprintf("[Success rate][Gamma[α∈(1,100),θ∈(1,100)](true) v.s. Gaussian]: %f\n",total_passed/1000.0)
#     end

#     let
#         total_passed = 0
#         for i = 1:1000
#             sample = rand(Distributions.Gamma(rand(1:10),rand(1:10)),5000)
#             if choose(Gaussian(),Gamma(),data=sample,
#                       comp=Comparator{MGF}(MGF(),eval_at = collect(0.01:0.001:0.02))) == :Gamma
#                 total_passed = total_passed + 1
#             end
#         end
#         @test total_passed/1000 > 0.9
#         @info @sprintf("[Success rate][Gamma[α∈(1,10),θ∈(1,10)](true) v.s. Gaussian]: %f\n",total_passed/1000.0)
#     end

#     let
#         total_passed = 0
#         for i = 1:1000
#             sample = rand(Distributions.Gamma(rand(1:10),rand(1:10)),5000)
#             if choose(Gaussian(),Gamma(),data=sample,
#                       comp=Comparator{MGF}(MGF())) == :Gamma
#                 total_passed = total_passed + 1
#             end
#         end
#         @test total_passed/1000 > 0.9
#         @info @sprintf("[Success rate][Gamma[α∈(1,10),θ∈(1,10)](true) v.s. Gaussian]: %f\n",total_passed/1000.0)
#     end

#     let
#         total_passed = 0
#         for i = 1:1000
#             sample = rand(Distributions.Gamma(rand(1:100),rand(1:100)),5000)
#             if choose(Gaussian(),Gamma(),data=sample,
#                       comp=Comparator{MGF}(MGF())) == :Gamma
#                 total_passed = total_passed + 1
#             end
#         end
#         @test total_passed/1000 > 0.9
#         @info @sprintf("[Success rate][Gamma[α∈(1,100),θ∈(1,100)](true) v.s. Gaussian]: %f\n",total_passed/1000.0)
#     end
#     let
#         total_passed = 0
#         for i = 1:1000
#             sample = rand(Distributions.Gaussian(rand(50:100),rand(1:5)),5000)
#             if choose(Gaussian(),Gamma(),data=sample,
#                       comp=Comparator{MGF}(MGF())) == :Gaussian
#                 total_passed = total_passed + 1
#             end
#         end
#         @test total_passed/1000 > 0.9
#         @info @sprintf("[Success rate][Gamma[α∈(1,100),θ∈(1,100)] v.s. Gaussian(true)]: %f\n",total_passed/1000.0)
#     end
# end

@testset "$(format("Model Fitting: check[continuous/integral]"))" begin
    # test the continuous case

    [let
     tc1 = rand(100) * 100
     @test check(:continuous,tc1) == true
     @test check(:integral,tc1)   == false
     tc2 = rand(100,100) * 100
     @test check(:continuous,tc2) == true
     @test check(:integral,tc2)   == false
     end for i = 1:5]
    # test the integral case
    [let
     tc1 = rand(1:10,100)
     @test check(:continuous,tc1) == false
     @test check(:integral,tc1)   == true
     tc2 = rand(1:10,100,100)
     @test check(:continuous,tc2) == false
     @test check(:integral,tc2)   == true
     end for i = 1:5]
end

@testset "$(format("Model Fitting: check[Bernoulli]"))"  begin
    [let
     tc1 = rand(0:1,100)
     @test check(:Bernoulli,tc1) ==  true
     @test check(:Bernoulli,tc1 .* 1.0) == true
     @test check(:Bernoulli,tc1 .* 1.5) == false
     tc2 = rand(0:1,100,100)
     @test check(:Bernoulli,tc2) == true
     @test check(:Bernoulli,tc2 .* 1.0) == true
     @test check(:Bernoulli,tc2 .* 1.5) == false
    end for i = 1:5]
end

@testset "$(format("Model Fitting: check[Poisson/NB]"))" begin
    let
        tc1 = rand(Distributions.Poisson(10),100)
        @test choose(:Poisson,:NegativeBinomial;data = tc1) == :Poisson
        @test choose(Poisson,NegativeBinomial;data=tc1)     == :Poisson
        @test choose(Poisson(),NegativeBinomial();data = tc1)  == :Poisson
    end
end

@testset "$(format("Model Fitting: check[support]"))" begin
    let
        tc = collect(1:10)
        @test check(:support,tc,layout=:flatten) == (1,10)
        @test check(:support,collect(-5:5),layout=:flatten) == (-5,5)
    end
end

@testset "$(format("Model Fitting: check[Gamma/Gaussian]"))" begin
    [let
        tc_vec = rand(Distributions.Gaussian(rand()*100,rand()*10),10000)
        @test choose(Gaussian,Gamma,data= tc_vec) == :Gaussian
        @test choose(Gaussian(),Gamma(),data= tc_vec) == :Gaussian
        @test choose(:Gaussian,:Gamma,data= tc_vec) == :Gaussian
        tc_mat = rand(Distributions.Gaussian(rand()*100,rand()*10),1000,1000)
        @test choose(Gaussian,Gamma,data= tc_mat) == :Gaussian
        @test choose(Gaussian(),Gamma(),data= tc_mat) == :Gaussian
        @test choose(:Gaussian,:Gamma,data= tc_mat) == :Gaussian
    end for i =1:5]
    [let
        tc_vec = rand(Distributions.Gamma(rand()*10,rand()*10),10000)
        @test choose(Gaussian,Gamma,data= tc_vec)     == :Gamma
        @test choose(Gaussian(),Gamma(),data= tc_vec) == :Gamma
        @test choose(:Gaussian,:Gamma,data= tc_vec)   == :Gamma
        tc_mat = rand(Distributions.Gamma(rand()*10,rand()*10),1000,1000)
        @test choose(Gaussian,Gamma,data= tc_mat)     == :Gamma
        @test choose(Gaussian(),Gamma(),data= tc_mat) == :Gamma
        @test choose(:Gaussian,:Gamma,data= tc_mat)   == :Gamma
        tc_mat = rand(Distributions.Gamma(rand()*100,rand()*100),1000,1000)
        @test choose(Gaussian,Gamma,data= tc_mat)     == :Gamma
        @test choose(Gaussian(),Gamma(),data= tc_mat) == :Gamma
        @test choose(:Gaussian,:Gamma,data= tc_mat)   == :Gamma
    end for i=1:5]
end


text:
export MtlHeapDescriptor, MtlHeap

const MTLHeapDescriptor = Ptr{MtHeapDescriptor}

mutable struct MtlHeapDescriptor
    handle::MTLHeapDescriptor
end

function MtlHeapDescriptor()
    handle = mtNewHeapDescriptor()
    obj = MtlHeapDescriptor(handle)
    finalizer(unsafe_destroy!, obj)
    return obj
end

function unsafe_destroy!(desc::MtlHeapDescriptor)
    mtRelease(desc.handle)
end

Base.unsafe_convert(::Type{MTLHeapDescriptor}, d::MtlHeapDescriptor) = d.handle

Base.:(==)(a::MtlHeapDescriptor, b::MtlHeapDescriptor) = a.handle == b.handle
Base.hash(dev::MtlHeapDescriptor, h::UInt) = hash(dev.handle, h)

## properties

Base.propertynames(::MtlHeapDescriptor) =
    (:type, :storageMode, :cpuCacheMode, :hazardTrackingMode, :resourceOptions, :size)

function Base.getproperty(o::MtlHeapDescriptor, f::Symbol)
    if f === :type
        mtHeapDescriptorType(o)
    elseif f === :storageMode
        mtHeapDescriptorStorageMode(o)
    elseif f === :cpuCacheMode
        mtHeapDescriptorCPUCacheMode(o)
    elseif f === :hazardTrackingMode
        mtHeapDescriptorHazardTrackingMode(o)
    elseif f === :resourceOptions
        mtHeapDescriptorResourceOptions(o)
    elseif f === :size
        mtHeapDescriptorSize(o)
    else
        getfield(o, f)
    end
end

function Base.setproperty!(o::MtlHeapDescriptor, f::Symbol, val)
    if f === :type
        mtHeapDescriptorTypeSet(o, val)
    elseif f === :storageMode
        mtHeapDescriptorStorageModeSet(o, val)
    elseif f === :cpuCacheMode
        mtHeapDescriptorCpuCacheModeSet(o, val)
    elseif f === :hazardTrackingMode
        mtHeapDescriptorHazardTrackingModeSet(o, val)
    elseif f === :resourceOptions
        mtHeapDescriptorResourceOptionsSet(o, val)
    elseif f === :size
        mtHeapDescriptorSizeSet(o, val)
    else
        setfield!(o, f, val)
    end
end

## display

function Base.show(io::IO, ::MIME"text/plain", o::MtlHeapDescriptor)
    print(io, "MtlHeapDescriptor: ")
    for f = propertynames(o)
        print(io, "\n $f : $(getproperty(o, f))")
    end
end

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

const MTLHeap = Ptr{MtHeap}

mutable struct MtlHeap
    handle::MTLHeap
    device::MtlDevice
end

Base.unsafe_convert(::Type{MTLHeap}, d::MtlHeap) = d.handle

Base.:(==)(a::MtlHeap, b::MtlHeap) = a.handle == b.handle
Base.hash(dev::MtlHeap, h::UInt) = hash(dev.handle, h)

function MtlHeap(device::MtlDevice, desc::MtlHeapDescriptor)
    handle = mtDeviceNewHeapWithDescriptor(device, desc)
    obj = MtlHeap(handle, device)
    finalizer(unsafe_destroy!, obj)
    return obj
end

function unsafe_destroy!(desc::MtlHeap)
    mtRelease(desc.handle)
end

## properties

Base.propertynames(::MtlHeap) =
    (:device, :label, :type, :storageMode, :cpuCacheMode, :hazardTrackingMode,
      :resourceOptions, :size, :usedSize, :currentAllocatedSize, :maxAvailableSizeWithAlignment)

function Base.getproperty(o::MtlHeap, f::Symbol)
    if f === :device
        mtHeapDevice(heap)
    elseif f === :label
        ptr = mtHeapLabel(o)
        ptr == C_NULL ? nothing : unsafe_string(ptr)
    elseif f === :type
        mtHeapType(o)
    elseif f === :storageMode
        mtHeapStorageMode(o)
    elseif f === :cpuCacheMode
        mtHeapCPUCacheMode(o)
    elseif f === :hazardTrackingMode
        mtHeapHazardTrackingMode(o)
    elseif f === :resourceOptions
        mtHeapResourceOptions(o)
    elseif f === :size
        mtHeapSize(o)
    elseif f === :usedSize
        mtHeapUsedSize(o)
    elseif f === :currentAllocatedSize
        mtHeapCurrentAllocatedSize(o)
    elseif f === :maxAvailableSizeWithAlignment
        maxAvailableSizeWithAlignment(o)
    else
        getfield(o, f)
    end
end

## display

function Base.show(io::IO, ::MIME"text/plain", o::MtlHeap)
    print(io, "MtlHeap: ")
    for f = propertynames(o)
        print(io, "\n $f : $(getproperty(o, f))")
    end
end


text:
workspace();
using BoxChopper;

nbox = 100^2

boxes = [sort(randn(3,2),2) for i = 1:nbox]
r₀s = map(box->mean(box,2)[:,1], boxes)
nouts = [randn(3) for i = 1:nbox]

volfrac(boxes[1], nouts[1], r₀s[1])  # to force compilation

function benchmark()
    @time begin
        for i = 1:nbox
            volfrac(boxes[i], nouts[i], r₀s[i])
        end
    end
end

benchmark()  # run @time inside function to eliminate REPL effect


text:
using HMMs
using Brobdingnag
using Random

file = open( joinpath( dirname(dirname(pathof(HMMs))), "data", "thmm3_20040630" ), "r" )
hmm = read( file, HMMs.HMM{HMMs.GenTDist, Brob, Float64} )
close( file )
HMMs.reorder!( hmm )

Random.seed!(1)
y = rand( hmm, 1000 );
HMMs.setobservations( hmm, y );

initial = hmm.initialprobabilities
alpha = HMMs.forwardprobabilities( hmm );
penultimate = convert(Vector{Float64}, alpha[end,:]/sum(alpha[end,:]))

yp = rand( hmm, 1 )

HMMs.initialize( hmm )
@assert( hmm.initialprobabilities == convert( Vector{Float64}, alpha[end,:]/sum(alpha[end,:]) ) )

HMMs.update( hmm, yp[1] )
final = hmm.initialprobabilities

HMMs.clear( hmm )
HMMs.setobservations( hmm, hmm.y )
hmm.initialprobabilities = initial
alpha2 = HMMs.forwardprobabilities( hmm );

penultimate2 = convert(Vector{Float64}, alpha2[end-1,:]/sum(alpha2[end-1,:]))
@assert( penultimate == penultimate )

final2 = convert(Vector{Float64}, alpha2[end,:]/sum(alpha2[end,:]))
@assert( maximum(abs.(final .- final2)) < 1e-12 )


text:
const R = 8314.46261815324 #J/(kmol*K)
export R
const one_atm = 101325.0 #pascal
export one_atm