text:
int_rules_1_3_4 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.3.4*Normalizing*algebraic*functions =#
    @apply_utils Antiderivative(~(u') * (~(c') * (~d * (~(a') + ~(b') * ~x)) ^ ~q) ^ ~p, ~x) => ((~c * (~d * (~a + ~b * ~x)) ^ ~q) ^ ~p / (~a + ~b * ~x) ^ (~p * ~q)) * Antiderivative(~u * (~a + ~b * ~x) ^ (~p * ~q), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~q, ~p], ~x) && (Not(IntegerQ(~q)) && Not(IntegerQ(~p)))
    @apply_utils Antiderivative(~(u') * (~(c') * (~(d') * (~(a') + ~(b') * ~x) ^ ~n) ^ ~q) ^ ~p, ~x) => ((~c * (~d * (~a + ~b * ~x) ^ ~n) ^ ~q) ^ ~p / (~a + ~b * ~x) ^ (~n * ~p * ~q)) * Antiderivative(~u * (~a + ~b * ~x) ^ (~n * ~p * ~q), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~q, ~p], ~x) && (Not(IntegerQ(~q)) && Not(IntegerQ(~p)))
    @apply_utils Antiderivative(~(u') * (~(c') * (~(a') + ~(b') * (~x) ^ ~(n')) ^ ~q) ^ ~p, ~x) => Simp((~c * (~a + ~b * (~x) ^ ~n) ^ ~q) ^ ~p / (~a + ~b * (~x) ^ ~n) ^ (~p * ~q)) * Antiderivative(~u * (~a + ~b * (~x) ^ ~n) ^ (~p * ~q), ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~p, ~q], ~x) && GeQ(~a, 0)
    @apply_utils Antiderivative(~(u') * (~(c') * (~a + ~(b') * (~x) ^ ~(n')) ^ ~q) ^ ~p, ~x) => Simp((~c * (~a + ~b * (~x) ^ ~n) ^ ~q) ^ ~p / (1 + (~b * (~x) ^ ~n) / ~a) ^ (~p * ~q)) * Antiderivative(~u * (1 + (~b * (~x) ^ ~n) / ~a) ^ (~p * ~q), ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~p, ~q], ~x) && Not(GeQ(~a, 0))
    @apply_utils Antiderivative(~(u') * (~(e') * (~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(q') * (~c + ~(d') * (~x) ^ ~(n')) ^ ~(q')) ^ ~p, ~x) => Antiderivative(~u * (~e * (~d / ~b) ^ ~q * (~a + ~b * (~x) ^ ~n) ^ (2 * ~q)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (IntegerQ(~q) && EqQ(~b * ~c - ~a * ~d, 0))
    @apply_utils Antiderivative(~(u') * (~(e') * (~(a') + ~(b') * (~x) ^ ~(n')) ^ ~q * (~c + ~(d') * (~x) ^ ~(n')) ^ ~q) ^ ~p, ~x) => Antiderivative(~u * (~e * ((-((~a) ^ 2) * ~d) / ~b + ~b * ~d * (~x) ^ (2 * ~n)) ^ ~q) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (IntegerQ(~q) && EqQ(~b * ~c + ~a * ~d, 0))
    #= Antiderivative((~u')*(((~a')+(~b')*(~x)^(~n'))*((~c)+(~d')*(~x)^(~n')))^(~p),~x) :=  Antiderivative(u*(a+b*x^n)^p*(c+d*x^n)^p,x) <-- FreeQ([a,b,c,d,n,p],x) &&  EqQ(b+d,0) && GtQ(a,0) && GtQ(c,0) =#
    @apply_utils Antiderivative(~(u') * (~(e') * (~(a') + ~(b') * (~x) ^ ~(n')) * (~c + ~(d') * (~x) ^ ~(n'))) ^ ~p, ~x) => Antiderivative(~u * (~a * ~c * ~e + (~b * ~c + ~a * ~d) * ~e * (~x) ^ ~n + ~b * ~d * ~e * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x)
    @apply_utils Antiderivative(~(u') * ((~(e') * (~(a') + ~(b') * (~x) ^ ~(n'))) / (~c + ~(d') * (~x) ^ ~(n'))) ^ ~p, ~x) => ((~b * ~e) / ~d) ^ ~p * Antiderivative(~u, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && EqQ(~b * ~c - ~a * ~d, 0)
    @apply_utils Antiderivative(~(u') * ((~(e') * (~(a') + ~(b') * (~x) ^ ~(n'))) / (~c + ~(d') * (~x) ^ ~(n'))) ^ ~p, ~x) => Antiderivative((~u * (~a * ~e + ~b * ~e * (~x) ^ ~n) ^ ~p) / (~c + ~d * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (GtQ(~b * ~d * ~e, 0) && GtQ(~c - (~a * ~d) / ~b, 0))
    #= Antiderivative((~u')*((~e')*((~a')+(~b')*(~x)^(~n'))/((~c)+(~d')*(~x)^(~n')))^(~p),~x) :=  Antiderivative(u*(a*e+b*e*x^n)^p/(c+d*x^n)^p,x) <-- FreeQ([a,b,c,d,e,n,p],x) &&  EqQ(b*c+a*d,0) && GtQ(b*e/d,0) && GtQ(c,0) =#
    #= Antiderivative((~u')*((~e')*((~a')+(~b')*(~x)^(~n'))/((~c)+(~d')*(~x)^(~n')))^(~p),~x) :=  Antiderivative(u*(-a*e-b*e*x^n)^p/(-c-d*x^n)^p,x) <-- FreeQ([a,b,c,d,e,n,p],x) &&  EqQ(b*c+a*d,0) && GtQ(b*e/d,0) && LtQ(c,0) =#
    @apply_utils Antiderivative(((~(e') * (~(a') + ~(b') * (~x) ^ ~(n'))) / (~c + ~(d') * (~x) ^ ~(n'))) ^ ~p, ~x) => With([q = Denominator(~p)], ((q * ~e * (~b * ~c - ~a * ~d)) / ~n) * Subst(Antiderivative(((~x) ^ (q * (~p + 1) - 1) * (-(~a) * ~e + ~c * (~x) ^ q) ^ (1 / ~n - 1)) / (~b * ~e - ~d * (~x) ^ q) ^ (1 / ~n + 1), ~x), ~x, ((~e * (~a + ~b * (~x) ^ ~n)) / (~c + ~d * (~x) ^ ~n)) ^ (1 / q))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (FractionQ(~p) && IntegerQ(1 / ~n))
    @apply_utils Antiderivative((~x) ^ ~(m') * ((~(e') * (~(a') + ~(b') * ~x)) / (~c + ~(d') * ~x)) ^ ~p, ~x) => With([q = Denominator(~p)], q * ~e * (~b * ~c - ~a * ~d) * Subst(Antiderivative(((~x) ^ (q * (~p + 1) - 1) * (-(~a) * ~e + ~c * (~x) ^ q) ^ ~m) / (~b * ~e - ~d * (~x) ^ q) ^ (~m + 2), ~x), ~x, ((~e * (~a + ~b * ~x)) / (~c + ~d * ~x)) ^ (1 / q))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (FractionQ(~p) && IntegerQ(~m))
    @apply_utils Antiderivative((~x) ^ ~(m') * ((~(e') * (~(a') + ~(b') * (~x) ^ ~(n'))) / (~c + ~(d') * (~x) ^ ~(n'))) ^ ~p, ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * ((~e * (~a + ~b * ~x)) / (~c + ~d * ~x)) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && IntegerQ(Simplify((~m + 1) / ~n))
    @apply_utils Antiderivative((~f * ~x) ^ ~m * ((~(e') * (~(a') + ~(b') * (~x) ^ ~(n'))) / (~c + ~(d') * (~x) ^ ~(n'))) ^ ~p, ~x) => Simp((~c * ~x) ^ ~m / (~x) ^ ~m) * Antiderivative((~x) ^ ~m * ((~e * (~a + ~b * (~x) ^ ~n)) / (~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && IntegerQ(Simplify((~m + 1) / ~n))
    @apply_utils Antiderivative((~u) ^ ~(r') * ((~(e') * (~(a') + ~(b') * (~x) ^ ~(n'))) / (~c + ~(d') * (~x) ^ ~(n'))) ^ ~p, ~x) => With([q = Denominator(~p)], ((q * ~e * (~b * ~c - ~a * ~d)) / ~n) * Subst(Antiderivative(SimplifyIntegrand((((~x) ^ (q * (~p + 1) - 1) * (-(~a) * ~e + ~c * (~x) ^ q) ^ (1 / ~n - 1)) / (~b * ~e - ~d * (~x) ^ q) ^ (1 / ~n + 1)) * ReplaceAll(~u, (~x->begin
                                                    #= none:1 =#
                                                    #= none:3 =#
                                                    (-(~a) * ~e + ~c * (~x) ^ q) ^ (1 / ~n) / (~b * ~e - ~d * (~x) ^ q) ^ (1 / ~n)
                                                end)) ^ ~r, ~x), ~x), ~x, ((~e * (~a + ~b * (~x) ^ ~n)) / (~c + ~d * (~x) ^ ~n)) ^ (1 / q))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolynomialQ(~u, ~x) && (FractionQ(~p) && (IntegerQ(1 / ~n) && IntegerQ(~r))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~u) ^ ~(r') * ((~(e') * (~(a') + ~(b') * (~x) ^ ~(n'))) / (~c + ~(d') * (~x) ^ ~(n'))) ^ ~p, ~x) => With([q = Denominator(~p)], ((q * ~e * (~b * ~c - ~a * ~d)) / ~n) * Subst(Antiderivative(SimplifyIntegrand((((~x) ^ (q * (~p + 1) - 1) * (-(~a) * ~e + ~c * (~x) ^ q) ^ ((~m + 1) / ~n - 1)) / (~b * ~e - ~d * (~x) ^ q) ^ ((~m + 1) / ~n + 1)) * ReplaceAll(~u, (~x->begin
                                                    #= none:1 =#
                                                    #= none:3 =#
                                                    (-(~a) * ~e + ~c * (~x) ^ q) ^ (1 / ~n) / (~b * ~e - ~d * (~x) ^ q) ^ (1 / ~n)
                                                end)) ^ ~r, ~x), ~x), ~x, ((~e * (~a + ~b * (~x) ^ ~n)) / (~c + ~d * (~x) ^ ~n)) ^ (1 / q))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolynomialQ(~u, ~x) && (FractionQ(~p) && (IntegerQ(1 / ~n) && IntegersQ(~m, ~r))))
    @apply_utils Antiderivative(~(u') * (~a + ~(b') / (~c + ~(d') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(~u * ((~b + ~a * ~c + ~a * ~d * (~x) ^ ~n) / (~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x)
    @apply_utils Antiderivative(~(u') * (~(e') * (~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(q') * (~c + ~(d') * (~x) ^ ~n) ^ ~(r')) ^ ~p, ~x) => Simp((~e * (~a + ~b * (~x) ^ ~n) ^ ~q * (~c + ~d * (~x) ^ ~n) ^ ~r) ^ ~p / ((~a + ~b * (~x) ^ ~n) ^ (~p * ~q) * (~c + ~d * (~x) ^ ~n) ^ (~p * ~r))) * Antiderivative(~u * (~a + ~b * (~x) ^ ~n) ^ (~p * ~q) * (~c + ~d * (~x) ^ ~n) ^ (~p * ~r), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p, ~q, ~r], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * (~(c') / ~x) ^ ~n) ^ ~p, ~x) => -(~c) * Subst(Antiderivative((~a + ~b * (~x) ^ ~n) ^ ~p / (~x) ^ 2, ~x), ~x, ~c / ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * (~(c') / ~x) ^ ~n) ^ ~p, ~x) => -((~c) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * (~x) ^ ~n) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, ~c / ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~p], ~x) && IntegerQ(~m)
    @apply_utils Antiderivative((~(d') * ~x) ^ ~m * (~(a') + ~(b') * (~(c') / ~x) ^ ~n) ^ ~p, ~x) => -(~c) * (~d * ~x) ^ ~m * (~c / ~x) ^ ~m * Subst(Antiderivative((~a + ~b * (~x) ^ ~n) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, ~c / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && Not(IntegerQ(~m))
    @apply_utils Antiderivative((~(a') + ~(b') * (~(d') / ~x) ^ ~n + ~(c') * (~(d') / ~x) ^ ~(n2')) ^ ~p, ~x) => -(~d) * Subst(Antiderivative((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p / (~x) ^ 2, ~x), ~x, ~d / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && EqQ(~n2, 2 * ~n)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~(d') / ~x) ^ ~n + ~(c') * (~(d') / ~x) ^ ~(n2')) ^ ~p, ~x) => -((~d) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, ~d / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && IntegerQ(~m))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * (~a + ~(b') * (~(d') / ~x) ^ ~n + ~(c') * (~(d') / ~x) ^ ~(n2')) ^ ~p, ~x) => -(~d) * (~e * ~x) ^ ~m * (~d / ~x) ^ ~m * Subst(Antiderivative((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, ~d / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && Not(IntegerQ(~m)))
    @apply_utils Antiderivative((~(a') + ~(b') * (~(d') / ~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => -(~d) * Subst(Antiderivative((~a + ~b * (~x) ^ ~n + (~c / (~d) ^ (2 * ~n)) * (~x) ^ (2 * ~n)) ^ ~p / (~x) ^ 2, ~x), ~x, ~d / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (EqQ(~n2, -2 * ~n) && IntegerQ(2 * ~n))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~a + ~(b') * (~(d') / ~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => -((~d) ^ (~m + 1)) * Subst(Antiderivative((~a + ~b * (~x) ^ ~n + (~c / (~d) ^ (2 * ~n)) * (~x) ^ (2 * ~n)) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, ~d / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (EqQ(~n2, -2 * ~n) && (IntegerQ(2 * ~n) && IntegerQ(~m)))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * (~a + ~(b') * (~(d') / ~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => -(~d) * (~e * ~x) ^ ~m * (~d / ~x) ^ ~m * Subst(Antiderivative((~a + ~b * (~x) ^ ~n + (~c / (~d) ^ (2 * ~n)) * (~x) ^ (2 * ~n)) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, ~d / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (EqQ(~n2, -2 * ~n) && (Not(IntegerQ(~m)) && IntegerQ(2 * ~n)))
    @apply_utils Antiderivative(~(u') * (~(e') * (~a + ~(b') * (~x) ^ ~(n')) ^ ~(r')) ^ ~p * (~(f') * (~c + ~(d') * (~x) ^ ~(n')) ^ ~s) ^ ~q, ~x) => (((~e * (~a + ~b * (~x) ^ ~n) ^ ~r) ^ ~p * (~f * (~c + ~d * (~x) ^ ~n) ^ ~s) ^ ~q) / ((~a + ~b * (~x) ^ ~n) ^ (~p * ~r) * (~c + ~d * (~x) ^ ~n) ^ (~q * ~s))) * Antiderivative(~u * (~a + ~b * (~x) ^ ~n) ^ (~p * ~r) * (~c + ~d * (~x) ^ ~n) ^ (~q * ~s), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n, ~p, ~q, ~r, ~s], ~x)
    @apply_utils Antiderivative((~u) ^ ~m, ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m, ~x) <-- FreeQ(~m, ~x) && (LinearQ(~u, ~x) && Not(LinearMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~v) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * ExpandToSum(~v, ~x) ^ ~n, ~x) <-- FreeQ([~m, ~n], ~x) && (LinearQ([~u, ~v], ~x) && Not(LinearMatchQ([~u, ~v], ~x)))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~v) ^ ~(n') * (~w) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * ExpandToSum(~v, ~x) ^ ~n * ExpandToSum(~w, ~x) ^ ~p, ~x) <-- FreeQ([~m, ~n, ~p], ~x) && (LinearQ([~u, ~v, ~w], ~x) && Not(LinearMatchQ([~u, ~v, ~w], ~x)))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~v) ^ ~(n') * (~w) ^ ~(p') * (~z) ^ ~(q'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * ExpandToSum(~v, ~x) ^ ~n * ExpandToSum(~w, ~x) ^ ~p * ExpandToSum(~z, ~x) ^ ~q, ~x) <-- FreeQ([~m, ~n, ~p, ~q], ~x) && (LinearQ([~u, ~v, ~w, ~z], ~x) && Not(LinearMatchQ([~u, ~v, ~w, ~z], ~x)))
    @apply_utils Antiderivative((~u) ^ ~p, ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ(~p, ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * (~u) ^ ~(p'), ~x) => Antiderivative((~c * ~x) ^ ~m * ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ([~c, ~m, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~u) ^ ~(p') * (~v) ^ ~(q'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~p * ExpandToSum(~v, ~x) ^ ~q, ~x) <-- FreeQ([~p, ~q], ~x) && (BinomialQ([~u, ~v], ~x) && (EqQ(BinomialDegree(~u, ~x) - BinomialDegree(~v, ~x), 0) && Not(BinomialMatchQ([~u, ~v], ~x))))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~u) ^ ~(p') * (~v) ^ ~(q'), ~x) => Antiderivative((~e * ~x) ^ ~m * ExpandToSum(~u, ~x) ^ ~p * ExpandToSum(~v, ~x) ^ ~q, ~x) <-- FreeQ([~e, ~m, ~p, ~q], ~x) && (BinomialQ([~u, ~v], ~x) && (EqQ(BinomialDegree(~u, ~x) - BinomialDegree(~v, ~x), 0) && Not(BinomialMatchQ([~u, ~v], ~x))))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~v) ^ ~(p') * (~w) ^ ~(q'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * ExpandToSum(~v, ~x) ^ ~p * ExpandToSum(~w, ~x) ^ ~q, ~x) <-- FreeQ([~m, ~p, ~q], ~x) && (BinomialQ([~u, ~v, ~w], ~x) && (EqQ(BinomialDegree(~u, ~x) - BinomialDegree(~v, ~x), 0) && (EqQ(BinomialDegree(~u, ~x) - BinomialDegree(~w, ~x), 0) && Not(BinomialMatchQ([~u, ~v, ~w], ~x)))))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~u) ^ ~(p') * (~v) ^ ~(q') * (~z) ^ ~(r'), ~x) => Antiderivative((~g * ~x) ^ ~m * ExpandToSum(~u, ~x) ^ ~p * ExpandToSum(~v, ~x) ^ ~q * ExpandToSum(~z, ~x) ^ ~r, ~x) <-- FreeQ([~g, ~m, ~p, ~q, ~r], ~x) && (BinomialQ([~u, ~v, ~z], ~x) && (EqQ(BinomialDegree(~u, ~x) - BinomialDegree(~v, ~x), 0) && (EqQ(BinomialDegree(~u, ~x) - BinomialDegree(~z, ~x), 0) && Not(BinomialMatchQ([~u, ~v, ~z], ~x)))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~u) ^ ~(p'), ~x) => Antiderivative((~c * ~x) ^ ~m * ~Pq * ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ([~c, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x))))
    @apply_utils Antiderivative((~u) ^ ~p, ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ(~p, ~x) && (GeneralizedBinomialQ(~u, ~x) && Not(GeneralizedBinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * (~u) ^ ~(p'), ~x) => Antiderivative((~c * ~x) ^ ~m * ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ([~c, ~m, ~p], ~x) && (GeneralizedBinomialQ(~u, ~x) && Not(GeneralizedBinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~u) ^ ~p, ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ(~p, ~x) && (QuadraticQ(~u, ~x) && Not(QuadraticMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~v) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * ExpandToSum(~v, ~x) ^ ~p, ~x) <-- FreeQ([~m, ~p], ~x) && (LinearQ(~u, ~x) && (QuadraticQ(~v, ~x) && Not(LinearMatchQ(~u, ~x) && QuadraticMatchQ(~v, ~x))))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~v) ^ ~(n') * (~w) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * ExpandToSum(~v, ~x) ^ ~n * ExpandToSum(~w, ~x) ^ ~p, ~x) <-- FreeQ([~m, ~n, ~p], ~x) && (LinearQ([~u, ~v], ~x) && (QuadraticQ(~w, ~x) && Not(LinearMatchQ([~u, ~v], ~x) && QuadraticMatchQ(~w, ~x))))
    @apply_utils Antiderivative((~u) ^ ~(p') * (~v) ^ ~(q'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~p * ExpandToSum(~v, ~x) ^ ~q, ~x) <-- FreeQ([~p, ~q], ~x) && (QuadraticQ([~u, ~v], ~x) && Not(QuadraticMatchQ([~u, ~v], ~x)))
    @apply_utils Antiderivative((~z) ^ ~(m') * (~u) ^ ~(p') * (~v) ^ ~(q'), ~x) => Antiderivative(ExpandToSum(~z, ~x) ^ ~m * ExpandToSum(~u, ~x) ^ ~p * ExpandToSum(~v, ~x) ^ ~q, ~x) <-- FreeQ([~m, ~p, ~q], ~x) && (LinearQ(~z, ~x) && (QuadraticQ([~u, ~v], ~x) && Not(LinearMatchQ(~z, ~x) && QuadraticMatchQ([~u, ~v], ~x))))
    @apply_utils Antiderivative(~Pq * (~u) ^ ~(p'), ~x) => Antiderivative(~Pq * ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ(~p, ~x) && (PolyQ(~Pq, ~x) && (QuadraticQ(~u, ~x) && Not(QuadraticMatchQ(~u, ~x))))
    @apply_utils Antiderivative((~u) ^ ~(m') * ~Pq * (~v) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * ~Pq * ExpandToSum(~v, ~x) ^ ~p, ~x) <-- FreeQ([~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (LinearQ(~u, ~x) && (QuadraticQ(~v, ~x) && Not(LinearMatchQ(~u, ~x) && QuadraticMatchQ(~v, ~x)))))
    @apply_utils Antiderivative((~u) ^ ~p, ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ(~p, ~x) && (TrinomialQ(~u, ~x) && Not(TrinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~u) ^ ~(p'), ~x) => Antiderivative((~d * ~x) ^ ~m * ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ([~d, ~m, ~p], ~x) && (TrinomialQ(~u, ~x) && Not(TrinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~u) ^ ~(q') * (~v) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~q * ExpandToSum(~v, ~x) ^ ~p, ~x) <-- FreeQ([~p, ~q], ~x) && (BinomialQ(~u, ~x) && (TrinomialQ(~v, ~x) && Not(BinomialMatchQ(~u, ~x) && TrinomialMatchQ(~v, ~x))))
    @apply_utils Antiderivative((~u) ^ ~(q') * (~v) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~q * ExpandToSum(~v, ~x) ^ ~p, ~x) <-- FreeQ([~p, ~q], ~x) && (BinomialQ(~u, ~x) && (BinomialQ(~v, ~x) && Not(BinomialMatchQ(~u, ~x) && BinomialMatchQ(~v, ~x))))
    @apply_utils Antiderivative((~(f') * ~x) ^ ~(m') * (~z) ^ ~(q') * (~u) ^ ~(p'), ~x) => Antiderivative((~f * ~x) ^ ~m * ExpandToSum(~z, ~x) ^ ~q * ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ([~f, ~m, ~p, ~q], ~x) && (BinomialQ(~z, ~x) && (TrinomialQ(~u, ~x) && Not(BinomialMatchQ(~z, ~x) && TrinomialMatchQ(~u, ~x))))
    @apply_utils Antiderivative((~(f') * ~x) ^ ~(m') * (~z) ^ ~(q') * (~u) ^ ~(p'), ~x) => Antiderivative((~f * ~x) ^ ~m * ExpandToSum(~z, ~x) ^ ~q * ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ([~f, ~m, ~p, ~q], ~x) && (BinomialQ(~z, ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~z, ~x) && BinomialMatchQ(~u, ~x))))
    @apply_utils Antiderivative(~Pq * (~u) ^ ~(p'), ~x) => Antiderivative(~Pq * ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ(~p, ~x) && (PolyQ(~Pq, ~x) && (TrinomialQ(~u, ~x) && Not(TrinomialMatchQ(~u, ~x))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~u) ^ ~(p'), ~x) => Antiderivative((~d * ~x) ^ ~m * ~Pq * ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ([~d, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (TrinomialQ(~u, ~x) && Not(TrinomialMatchQ(~u, ~x))))
    @apply_utils Antiderivative((~u) ^ ~p, ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ(~p, ~x) && (GeneralizedTrinomialQ(~u, ~x) && Not(GeneralizedTrinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~u) ^ ~(p'), ~x) => Antiderivative((~d * ~x) ^ ~m * ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ([~d, ~m, ~p], ~x) && (GeneralizedTrinomialQ(~u, ~x) && Not(GeneralizedTrinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative(~z * (~u) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~z, ~x) * ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ(~p, ~x) && (BinomialQ(~z, ~x) && (GeneralizedTrinomialQ(~u, ~x) && (EqQ(BinomialDegree(~z, ~x) - GeneralizedTrinomialDegree(~u, ~x), 0) && Not(BinomialMatchQ(~z, ~x) && GeneralizedTrinomialMatchQ(~u, ~x)))))
    @apply_utils Antiderivative((~(f') * ~x) ^ ~(m') * ~z * (~u) ^ ~(p'), ~x) => Antiderivative((~f * ~x) ^ ~m * ExpandToSum(~z, ~x) * ExpandToSum(~u, ~x) ^ ~p, ~x) <-- FreeQ([~f, ~m, ~p], ~x) && (BinomialQ(~z, ~x) && (GeneralizedTrinomialQ(~u, ~x) && (EqQ(BinomialDegree(~z, ~x) - GeneralizedTrinomialDegree(~u, ~x), 0) && Not(BinomialMatchQ(~z, ~x) && GeneralizedTrinomialMatchQ(~u, ~x)))))
end

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


text:
using Base.CoreLogging

import Base.CoreLogging: min_enabled_level, shouldlog, handle_message, catch_exceptions

export VisualDLLogger,
 ScalarLevel, HistogramLevel, EmbeddingLevel, TextLevel, AudioLevel, ImageLevel,
 as_mode, add_component, save,  with_logger,
 set_caption, start_sampling, finish_sampling,
 @log_scalar, @log_histogram, @log_image, @log_text

const ScalarLevel = LogLevel(880)
const HistogramLevel = LogLevel(881)
const EmbeddingLevel = LogLevel(882)
const TextLevel = LogLevel(883)
const AudioLevel = LogLevel(884)
const ImageLevel = LogLevel(885)

"""
`VisualDLLogger` is  a subtype of `AbstractLogger`. And all the necessary interfaces
are implemented. This means that you can use `Base.CoreLogging.disable_logging`
to control the log level and use the similar grammar like `@info`.

The provided log levels are:

```
const ScalarLevel = LogLevel(880)
const HistogramLevel = LogLevel(881)
const EmbeddingLevel = LogLevel(882)
const TextLevel = LogLevel(883)
const AudioLevel = LogLevel(884)
const ImageLevel = LogLevel(885)
```

which are between the `Info` level(`LogLevel(0)``) and `Warn` level(`LogLevel(1000)`).

There are two members in the `VisualDLLogger`, a `pylogger` which is a wrapper of
the python logger instance, and a `components` which is a `Dict` containing all the
components related to the `logger`.
"""
struct VisualDLLogger <: AbstractLogger
    pylogger::PyObject
    components::Dict{Symbol, PyObject}
end

"""
    VisualDLLogger(log_path, sync_cycle, mode="train")

Create a logger instance. `log_path` is used to specify where to store the log data.
`sync_cycle` specify how often should the system store data into the file system.
Typically adding a record requires 6 operations. System will save the data into
the file system once operations count reaches `sync_cycle`.
"""
function VisualDLLogger(log_path::AbstractString, sync_cycle::Int, mode::AbstractString="train")
    pylogger = LogWriter(log_path, sync_cycle)
    pylogger[:mode](mode)
    VisualDLLogger(pylogger, Dict{Symbol, PyObject}())
end

"""
    as_mode(logger, mode)

Clone a logger and reset the mode.
"""
function as_mode(logger::VisualDLLogger, mode::AbstractString)
    VisualDLLogger(logger.pylogger[:as_mode](mode), Dict{Symbol, PyObject}())
end

"""
    set_caption(logger, tag, caption)

Set the caption of a figure.
"""
function set_caption(logger::VisualDLLogger, tag::Symbol, caption::AbstractString)
    logger.components[tag][:set_caption](caption)
end

"""
    start_sampling(logger, tag)

This function is only valid for `@log_image` and `@log_audio`.
"""
function start_sampling(logger::VisualDLLogger, tag::Symbol)
    logger.components[tag][:start_sampling]()
end

"""
    finish_sampling(logger, tag)

This function is only valid for `@log_image` and `@log_audio`.
"""
function finish_sampling(logger::VisualDLLogger, tag::Symbol)
    logger.components[tag][:finish_sampling]()
end

"""
    save(logger)

Although the `logger` will automatically save data when the `sync_cycle`
is reached. It is better to force save in the end.
"""
function save(logger::VisualDLLogger)
    logger.pylogger[:save]()
end

"""
    add_component(logger, level, tag; kwargs...)

Althrough this package will automatically initial a component for you when you
add records to a component for the first time, it is strongly suggested to use
this function to initial a component manually.

For the detail of `kwargs`, please read the [python api](http://visualdl.paddlepaddle.org/docs/develop/visualdl/en/api/write_logs.html)
"""
function add_component(logger::VisualDLLogger, level::LogLevel, tag::Symbol; kwargs...)
    kwargs = Dict(kwargs)
    if level == ScalarLevel
        logger.components[tag] = logger.pylogger[:scalar](string(tag))
    elseif level == HistogramLevel
        logger.components[tag] = logger.pylogger[:histogram](string(tag), get(kwargs, :num_buckets, 10))
    elseif level == EmbeddingLevel
        logger.components[tag] = logger.pylogger[:embedding](string(tag))
    elseif level == TextLevel
        logger.components[tag] = logger.pylogger[:text](string(tag))
    elseif level == AudioLevel
        logger.components[tag] = logger.pylogger[:audio](
            string(tag),
            get(kwargs, :num_samples, 1),
            get(kwargs, :step_cycle, 1))
    elseif level == ImageLevel
        logger.components[tag] = logger.pylogger[:image](
            string(tag),
            get(kwargs, :num_samples, 1),
            get(kwargs, :step_cycle, 1))
    end
end

function add_record(logger::VisualDLLogger, level::LogLevel, tag::Symbol, args)
    if level == ScalarLevel
        logger.components[tag][:add_record](args...)
    elseif level == HistogramLevel
        logger.components[tag][:add_record](args...)
    elseif level == EmbeddingLevel
        logger.components[tag][:add_embeddings_with_word_dict](args...)
    elseif level == TextLevel
        logger.components[tag][:add_record](args...)
    elseif level == AudioLevel
        start_sampling(logger, tag)
        logger.components[tag][:add_sample](args...)
        finish_sampling(logger, tag)
    elseif level == ImageLevel
        start_sampling(logger, tag)
        if args isa Array{T, 3} where T <: Number
            logger.components[tag][:add_sample](size(args), collect(Iterators.flatten(args)))
        else
            logger.components[tag][:add_sample](args...)
        end
        finish_sampling(logger, tag)
    end
end

# interfaces
min_enabled_level(logger::VisualDLLogger) = ScalarLevel
shouldlog(logger::VisualDLLogger, level, _module, group, id) =  ScalarLevel ≤ level ≤ ImageLevel
catch_exceptions(logger::VisualDLLogger) = false

function handle_message(logger::VisualDLLogger, level, message, _module, group, id, file, line; kwargs...)
    for (key,val) in pairs(kwargs)
        haskey(logger.components, key) ? logger.components[key] : add_component(logger, level, key)
        add_record(logger, level, key, val)
    end
end

function preprocess(exs)
    exs_interp = Expr[]
    for ex in exs
        if ex isa Expr && ex.head === :(=)
            k,v = ex.args
            push!(exs_interp, Expr(:(=), k, esc(v)))
        else
            throw(ArgumentError("No component tag found in `$ex`, use `tagname=$ex` instead!"))
        end
    end
    exs_interp
end

macro log_scalar(exs...) :(@logmsg ScalarLevel "" $(preprocess(exs)...)) end
macro log_histogram(exs...) :(@logmsg HistogramLevel "" $(preprocess(exs)...)) end
macro log_embedding(exs...) :(@logmsg EmbeddingLevel "" $(preprocess(exs)...)) end
macro log_text(exs...) :(@logmsg TextLevel "" $(preprocess(exs)...)) end
macro log_audio(exs...) :(@logmsg AudioLevel "" $(preprocess(exs)...)) end
macro log_image(exs...) :(@logmsg ImageLevel "" $(preprocess(exs)...)) end

_log_comp_docs = """
    @log_scalar [tag_name=(args...) | ...]
    @log_histogram [tag_name=(args...) | ...]
    @log_embedding [tag_name=(args...) | ...]
    @log_text [tag_name=(args...) | ...]
    @log_audio [tag_name=(args...) | ...]
    @log_image [tag_name=(args...) | ...]

Add a record to a component which has the `tag_name`. If the component is not initialized yet,
a new component will be created automatically with default arguments. For the detail of `args`
of each component, please read the [doc](http://visualdl.paddlepaddle.org/docs/develop/visualdl/en/api/write_logs.html)
"""

# Logging macros share documentation
@eval @doc $_log_comp_docs :(@log_scalar)
@eval @doc $_log_comp_docs :(@log_histogram)
@eval @doc $_log_comp_docs :(@log_embedding)
@eval @doc $_log_comp_docs :(@log_text)
@eval @doc $_log_comp_docs :(@log_audio)
@eval @doc $_log_comp_docs :(@log_image)

ext: jl


text:
# This file investigates integration error with respect to free variables in 1D integration

using FractalTools
using Makie

# Construct interpolation data
f(x) = sin(2π * x) + 1
dx   = 0.1
xi   = 0.
xf   = 1.
x    = collect(xi : dx : xf)
y    = f.(x)
pts  = collect.(zip(x, y))

# Compute errors
fval = (cos(2π * xi) - cos(2π * xf)) / (2π) + (xf - xi)
freevars = 0.001 : 0.01 : 0.999
mse = map(freevars) do freevar
    ival = integrate(pts, Interp1D(freevar))
    abs(fval - ival)
end

# Plot mse
fig = Figure()
ax = fig[1, 1] = Axis(fig, xlabel="Free Variable", ylabel="Integration Error", title="1D Integration Error")
stem!(ax, freevars, mse, color=:black)
save(joinpath(@__DIR__, "integ1d_error.png"), fig)
display(fig)


text:
data = readdlm("RLBinding.csv",',') # julia uses "" to declare
# strings, '' for characters and # for comments

# one of the improvements over matlab - the use of squared brackets [] to
# index into arrays, to clearly distinguish it from function calls.
# Indices start at 1 in julia
t = data[:,1]; R = data[:,2]; L = data[:,3]; RL = data[:,4];

# similar to import in python - declare which package to use
using Winston # Winston is one of the three graphics packages available for julia at the moment

# FramedPlot can be used for plotting a single plot
p = FramedPlot("xlabel","time")
setattr(p,"ylabel","concentration") # attributes can be changed with setter functions
# similar to MATLABs `set` and `get`

# building a plot using Winston is like composing a collage:
# build one element after the other
# then add it all to the canvas (the FramedPlot)
r = Points(time,R,"type","diamond", "color","blue")
l = Points(time,L,"type","circle", "color","blue")
rl = Points(time,RL,"type","cross", "color","blue")

cr = Curve(time,R, "color","blue")
cl = Curve(time,L, "color","blue")
crl = Curve(time,RL, "color","blue")

setattr(r,"label","R")
setattr(l,"label","L")
setattr(rl,"label","RL")
lgnd = Legend(.8,.9,{r,l,rl})

# adding does not! need to be done at the end
add(p,r,l,rl,cr,cl,crl,lgnd)
# Winston/ijulia do not yet support the display of the generated plots,
# but it is possible to save the plot to a file
# eps, png and svg are supposed to be supported, but
# svg did not work on my machine out of the box
file(p,"julia_rlbinding.png")


text:
#Base.include(@__MODULE__,"../basis/BasisStructs.jl")

"""
  module Gradient
"""
module Gradient

include("GradHelpers.jl")

using JuliaChem.JCModules
using JuliaChem.JERI

using MPI
using JSON

function run(mol::Molecule, basis::Basis, rhf_energy;
  output="none")

  comm=MPI.COMM_WORLD

  if MPI.Comm_rank(comm) == 0 && output >= 2
      println("--------------------------------------------------------------------------------")
      println("                       ========================================                 ")
      println("                                RESTRICTED CLOSED-SHELL                         ")
      println("                                 HARTREE-FOCK GRADIENT                          ")
      println("                       ========================================                 ")
      println("")
  end

  #== initial setup ==#
  if true
    throw("JuliaChem.jl does not yet support RHF gradients!")
  else

  jeri_oei_grad_engine = JERI.OEIEngine(mol.mol_cxx,
    basis.basis_cxx, 1)
  jeri_prop_engine = JERI.PropEngine(mol.mol_cxx,
    basis.basis_cxx)

  W = rhf_energy["Energy-Weighted Density"]
  P = rhf_energy["Density"]

  #display(W); println()
  #display(P); println()

  #== compute nuclear repulsion gradient contribution==#
  nuclear_gradient = compute_nuc_grad(mol)
  println("NUCLEAR REPULSION: ")
  display(nuclear_gradient); println()

  #== compute overlap gradient contribution ==#
  overlap_gradient = compute_overlap_grad(mol, basis, W, jeri_oei_grad_engine)
  println("OVERLAP GRADIENT: ")
  display(overlap_gradient); println()

  #== compute kinetic gradient contribution ==#
  kinetic_gradient = compute_kinetic_grad(mol, basis, P, jeri_oei_grad_engine)
  println("KINETIC GRADIENT: ")
  display(kinetic_gradient); println()

  #== compute nuclear attraction gradient contribution ==#
  nuc_attr_gradient = compute_nuc_attr_grad(mol, basis, P, jeri_oei_grad_engine)
  println("NUC. ATTR. GRADIENT: ")
  display(nuc_attr_gradient); println()

  #== comput total gradient ==#
  total_gradient = nuclear_gradient .+ overlap_gradient .+ kinetic_gradient .+ nuc_attr_gradient
  #total_gradient = nuclear_gradient .+ overlap_gradient
  println("TOTAL GRADIENT: ")
  display(total_gradient); println()

  #== compute dipole moment ==#
  dipole = compute_dipole_moment(mol, basis, P, jeri_prop_engine)
  println("DIPOLE:")
  display(dipole); println()

  println("DIPOLE MOMENT:")
  display(sqrt(dipole[1]^2 + dipole[2]^2 + dipole[3]^2)); println()

  end

  if MPI.Comm_rank(comm) == 0 && output >= 2
    println("                       ========================================                 ")
    println("                              END RESTRICTED CLOSED-SHELL                       ")
    println("                                 HARTREE-FOCK GRADIENT                          ")
    println("                       ========================================                 ")
    println("--------------------------------------------------------------------------------")
  end

  #return nuclear_gradient
end
export run

end


text:
# some implementations for convenience
# maintained by Zach Sunberg

actions(mdp::MDP{S,Bool}) where S = (true, false)
actions(mdp::MDP{S,Bool}, s::S, t::Tuple{Bool,Bool}=(true,false)) where S = (true,false)
actions(mdp::POMDP{S,Bool,O}) where {S,O} = (true, false)
actions(mdp::POMDP{S,Bool,O}, s::S, t::Tuple{Bool,Bool}=(true,false)) where {S,O} = (true,false)
n_actions(mdp::MDP{S,Bool}) where S = 2
n_actions(mdp::POMDP{S,Bool,O}) where {S, O} = 2

rand(rng::AbstractRNG, t::Tuple{Bool, Bool}) = rand(rng, Bool)
rand(t::Tuple{Bool, Bool}) = rand(Bool)

support(s::AbstractVector) = s
support(s::Tuple) = s
support(r::AbstractRange) = r
support(g::Base.Generator) = g

states(mdp::MDP{Bool}) = (true, false)
states(mdp::POMDP{Bool}) = (true, false)
n_states(mdp::MDP{Bool}) = 2
n_states(mdp::POMDP{Bool}) = 2

observations(::POMDP{S,A,Bool}) where {S,A} = (true,false)
observations(::POMDP{S,A,Bool}, s::S) where {S,A} = (true,false)
n_observations(::POMDP{S,A,Bool}) where {S,A} = 2

stateindex(mdp::Union{MDP, POMDP}, s::Int) = s
actionindex(mdp::Union{MDP, POMDP}, a::Int) = a
obsindex(mdp::Union{MDP, POMDP}, o::Int) = o

stateindex(mdp::Union{MDP, POMDP}, s::Bool) = s+1
actionindex(mdp::Union{MDP, POMDP}, a::Bool) = a+1
obsindex(mdp::Union{MDP, POMDP}, o::Bool) = o+1


text:
"""
    struct GlobalMatrices
Stores global matrices and some other things.
Everything that is very expensive to compute which stays
constant throughout the simulation should be stored here.
    GlobalMatrices(basis::Basis, filter::Filter, dimensions)
Initialize GlobalMatrices for `basis`, `filter` and `dimensions`.
"""
 struct GlobalMatrices
     normalsigns::Dict{Face, Int64}
     normalidxs::Dict{Face, Int64}
     oppositefaces::Dict{Face, Face}

     function GlobalMatrices(dimensions)
        normalsigns = Dict(left => -1, right => 1, top => 1, bottom => -1)
        normalidxs = Dict(left => 1, right => 1, top => 2, bottom => 2)
        oppositefaces = Dict(left => right, right => left, top => bottom, bottom => top)

         new(normalsigns,
            normalidxs,
            oppositefaces)
     end
 end


text:
using Surrogates
using Distributions
using Test

#1D
lb = 0.0
ub = 5.0
n = 5
d = 1
Surrogates.sample(n,lb,ub,GridSample(0.1))
Surrogates.sample(n,lb,ub,UniformSample())
Surrogates.sample(n,lb,ub,SobolSample())
Surrogates.sample(n,lb,ub,LatinHypercubeSample())
Surrogates.sample(20,lb,ub,LowDiscrepancySample(10))
Surrogates.sample(5,d,Cauchy())
Surrogates.sample(5,d,Normal(0,4))

#ND
lb = [0.1,-0.5]
ub = [1.0,20.0]
n = 5
d = 2

#GridSample{T}
s = Surrogates.sample(n,lb,ub,GridSample([0.1,0.5]))
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

#UniformSample()
s = Surrogates.sample(n,lb,ub,UniformSample())
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

#SobolSample()
s = Surrogates.sample(n,lb,ub,SobolSample())
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

#LHS
s = Surrogates.sample(n,lb,ub,LatinHypercubeSample())
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

#LDS
s = Surrogates.sample(n,lb,ub,LowDiscrepancySample([10,3]))
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

#Distribution 1
s = Surrogates.sample(n,d,Cauchy())
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

#Distribution 2
s = Surrogates.sample(n,d,Normal(3,5))
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true


text:
@with_kw struct IZParameter
    a::SNNFloat = 0.01
    b::SNNFloat = 0.2
    c::SNNFloat = -65
    d::SNNFloat = 2

end

@with_kw mutable struct IZ
    param::IZParameter = IZParameter(a,b,c,d)
    N::SNNInt = 1
    v::Vector{SNNFloat} = fill(param.c, N)
    u::Vector{SNNFloat} = param.b * v
    fire::Vector{Bool} = zeros(Bool, N)
    I::Vector{SNNFloat} = zeros(N)
    records::Dict = Dict()
end

function integrate!(p::IZ, param::IZParameter, dt::SNNFloat)
    @unpack N, v, u, fire, I = p
    @unpack a, b, c, d = param
    @inbounds for i = 1:N
        #v[i] += 0.5f0dt * (0.04f0v[i]^2 + 5f0v[i] + 140f0 - u[i] + I[i])
        v[i] += 0.5f0dt * (0.04f0v[i]^2 + 5f0v[i] + 140f0 - u[i] + I[i])
        u[i] += dt * (a * (b * v[i] - u[i]))
    end

    @inbounds for i = 1:N
        fire[i] = v[i] > 30f0
        v[i] = ifelse(fire[i], c, v[i])
        u[i] += ifelse(fire[i], d, 0f0)
    end

end
#=
@inbounds for m = 1:N
    vT = v[m]+ 0.5f0dt * (k*(v[m] - vr)*(v[m] - vt)-u[m] + Iext[m])/C;
    v[m+1] = vT + 0.5f0dt  * (k*(v[m] - vr)*(v[m] - vt)-u[m] + Iext[m])/C;
    u[m+1] = u[m] + dt * a*(b*(v[m]-vr)-u[m]);
    u[m+1] = u[m] + dt * a*(b*(v[m+1]-vr)-u[m]);
    if v[m+1]>= vPeak:# % a spike is fired!
        v[m] = vPeak;# % padding the spike amplitude
        v[m+1] = c;# % membrane voltage reset
        u[m+1] = u[m+1] + d;# % recovery variable update
    end
end
=#

max_stars_repo_name: russelljjarvis/juliaSNN

max_issues_repo_name: russelljjarvis/juliaSNN

max_forks_repo_name: russelljjarvis/juliaSNN


text:
include("bpMCSetup.jl")

import MonteCarloMarkovKernels: estimateBM, kde
using StatsBase
import Statistics.var

using Plots
Plots.gr()

@load "data/bpMCSimple.jld" chainSimple sar simpleCovEstimate simpleElapsed
@load "data/bpMCRecycle.jld" chainRecycle rar recycleCovEstimate recycleElapsed
@load "data/bpMCBiased.jld" chainBiased bar biasedCovEstimate biasedElapsed

savefigures = false

vsSimple = (i->(x->x[i]).(chainSimple)).(1:d)
vsRecycle = (i->(x->x[i]).(chainRecycle)).(1:d)
vsBiased = (i->(x->x[i]).(chainBiased)).(1:d)

plot(kde(vsSimple[1], sar))
plot!(kde(vsRecycle[1], rar))
plot!(kde(vsBiased[1], bar))
savefigures && savefig("gk_kde1.png")

plot(kde(vsSimple[2], sar))
plot!(kde(vsRecycle[2], rar))
plot!(kde(vsBiased[2], bar))
savefigures && savefig("gk_kde2.png")

plot(kde(vsSimple[3], sar))
plot!(kde(vsRecycle[3], rar))
plot!(kde(vsBiased[3], bar))
savefigures && savefig("bp_kde3.png")

contour(kde(vsSimple[1], vsSimple[2], sar))
contour!(kde(vsRecycle[1], vsRecycle[2], rar))
savefigures && savefig("bp_kde12.png")

contour(kde(vsSimple[1], vsSimple[3], sar))
contour!(kde(vsRecycle[1], vsRecycle[3], rar))
savefigures && savefig("bp_kde13.png")
#
contour(kde(vsSimple[2], vsSimple[3], sar))
contour!(kde(vsRecycle[2], vsRecycle[3], rar))
savefigures && savefig("bp_kde23.png")

plot(autocor(vsSimple[1]))
plot!(autocor(vsRecycle[1]))
savefigures && savefig("pmmh_acf1.png")
plot(autocor(vsSimple[2]))
plot!(autocor(vsRecycle[2]))
savefigures && savefig("pmmh_acf2.png")
plot(autocor(vsSimple[3]))
plot!(autocor(vsRecycle[3]))
savefigures && savefig("pmmh_acf3.png")

avarsSimple = estimateBM.(vsSimple)
avarsRecycle = estimateBM.(vsRecycle)
avarsBiased = estimateBM.(vsBiased)

varsSimple = var.(vsSimple)
varsRecycle = var.(vsRecycle)
varsBiased = var.(vsBiased)

essesSimple = varsSimple./avarsSimple*chainLength
essesRecycle = varsRecycle./avarsRecycle*chainLength
essesBiased = varsBiased./avarsBiased*chainLength

estimatedSimpleElapsed = 4536000

essPerSecondSimple = essesSimple / estimatedSimpleElapsed
essPerSecondRecycle = essesRecycle / recycleElapsed
essPerSecondBiased = essesBiased / biasedElapsed

plot(kde(vsRecycle[1], rar), label="", xlabel="lambda", ylabel="density")
savefig("demo/bp_kde1.pdf")

plot(kde(vsRecycle[2], rar), label="", xlabel="k_on", ylabel="density")
savefig("demo/bp_kde2.pdf")

plot(kde(vsRecycle[3], rar), label="", xlabel="k_off", ylabel="density")
savefig("demo/bp_kde3.pdf")


text:
module AudioDisplay

export audioplayer

using WAV

function audioplayer(filepath)
  markup = """<audio controls="controls" {autoplay}>
              <source src="$filepath" />
              Your browser does not support the audio element.
              </audio>"""

  display(MIME("text/html") ,markup)
end

function audioplayer(s, fs)
    buf = IOBuffer()
    wavwrite(s, buf; Fs=fs)
    data = base64(bytestring(buf))
    markup = """<audio controls="controls" {autoplay}>
                <source src="data:audio/wav;base64,$data" type="audio/wav" />
                Your browser does not support the audio element.
                </audio>"""
    display(MIME("text/html") ,markup)
end

end

max_stars_repo_name: jfsantos/IJuliaNotebooks

max_issues_repo_name: jfsantos/IJuliaNotebooks

max_forks_repo_name: jfsantos/IJuliaNotebooks


text:
module ArraysTests

using Test

@testset "Interfaces" begin include("InterfaceTests.jl") end

@testset "BlockArraysCoo" begin include("BlockArraysCooTests.jl") end

@testset "VectorsOfBlockArrayCoo" begin include("VectorsOfBlockArrayCooTests.jl") end

@testset "CachedArrays" begin include("CachedArraysTests.jl") end

@testset "Kernels" begin include("KernelsTests.jl") end

@testset "Apply" begin include("ApplyTests.jl") end

@testset "CompressedArrays" begin include("CompressedArraysTests.jl") end

@testset "LocalToGlobalArrays" begin include("LocalToGlobalArraysTests.jl") end

@testset "LocalToGlobalPosNegArrays" begin include("LocalToGlobalPosNegArraysTests.jl") end

@testset "FilteredArraysTests" begin include("FilteredArraysTests.jl") end

@testset "Tables" begin include("TablesTests.jl") end

@testset "Reindex" begin include("ReindexTests.jl") end

@testset "IdentityVectors" begin include("IdentityVectorsTests.jl") end

@testset "SubVectors" begin include("SubVectorsTests.jl") end

@testset "ArrayPairs" begin include("ArrayPairsTests.jl") end

@testset "AppendedArrays" begin include("AppendedArraysTests.jl") end

@testset "AutodiffTests" begin include("AutodiffTests.jl") end

@testset "VectorWithEntryRemovedTests" begin include("VectorWithEntryRemovedTests.jl") end

@testset "VectorWithEntryInsertedTests" begin include("VectorWithEntryInsertedTests.jl") end

end # module


text:
using Gadfly

keys=["a.","b,","c/","d?","g;","h:","i'","l{","m|","p~","q!","r@","s#","t\$","u%","x*","y(","z_","A-","B=","C+"]

plot(y=collect(1:length(keys)), color=keys)


text:
include("air.jl")
UPPE.permittivity(w) = permittivity(w)
UPPE.permeability(w) = permeability(w)

grid = GridRT(
    GPU();
    ru = 1e-3,   # [m] unit of space in r direction
    rmax = 10,   # [ru] area in spatial domain
    Nr = 500,   # number of points in spatial domain
    rguard = 0,   # [ru] the width of the lossy slab at the end of r grid
    kguard = 90,   # [degrees] the cut-off angle for wave vectors
    tu = 1e-15,   # [s] unit of time
    tmin = -200,   # [tu] left boundary of t grid
    tmax = 200,   # [tu] right boundary of t grid
    Nt = 2^14,   # number of points on t grid
    tguard = 20,   # [tu] the width of the lossy slab at the ends of t grid
    wguard = 1e16,   # [1/s] the cut-off angular frequency
)

lam0 = 0.8e-6   # [m] central wavelength
a0 = 1e-3   # [m] initial beam radius
tau0 = 35e-15 / (2 * sqrt(log(2)))  # [s] initial pulse duration
I0 = 2e12 * 1e4   # [W/cm^2] initial intensity

w0 = 2 * pi * C0 / lam0   # central frequency
n0 = refractive_index(w0)
Iu = I0   # [W/m**2] unit of intensity
Eu = intensity2field(Iu, n0)

E = zeros((grid.Nr, grid.Nt))
for it=1:grid.Nt
for ir=1:grid.Nr
    E[ir,it] = sqrt(I0 / Iu) *
               exp(-0.5 * (grid.r[ir] * grid.ru)^2 / a0^2) *
               exp(-0.5 * (grid.t[it] * grid.tu)^2 / tau0^2) *
               cos(w0 * grid.t[it] * grid.tu)
end
end

field = Field(grid, w0, E; Eu, Iu)

ne = adapt_array(grid.arch, zeros(Float32, size(field.E)))
Pnl = zero(field.E)
J = zero(field.E)
response =  Response(ne, Pnl, J, nothing, nothing)

model = Model(
    grid, field, response;
    zu = 1,   # [m] unit of space in z direction
    z = 0,   # [m] initial propagation distance
    zmax = 7,   # [zu] propagation distance
    kparaxial = true,
    nonlinearity = false,
)

simulation = Simulation(
    model, n2, N0;
    prefix = "results/",
    dz0 = 1,   # [zu] initial z step
    dzout = model.zmax,   # [zu] z step for writing field into output file
)

run!(simulation)

# Compare to analytic solution:
E = collect(field.E)
(; ru, tu, Nr, Nt, r, t) = grid

zmax = model.zmax * model.zu
k0 = k_func(w0)
vg = group_velocity(w0)
zdiff = diffraction_length(w0, a0)
zdisp = dispersion_length(w0, tau0)
taud = tau0 * sqrt(1 - 1im * zmax / zdisp)
D = 1 + 1im * zmax / zdiff

Ea = zeros(ComplexF64, (Nr, Nt))
for it=1:Nt
for ir=1:Nr
    Ea[ir,it] = tau0 / taud / D *
                exp(-0.5 * (r[ir] * ru)^2 / a0^2 / D) *
                exp(-0.5 * (t[it] * tu)^2 / taud^2) *
                exp(1im * k0 * zmax - 1im * w0 * (t[it] * tu + zmax / vg))
end
end

@test isapprox(E, Ea; rtol=1e-2)


text:
using Distances, Multisets, StatsBase
export get_dist_dict, KL_error,  get_inner_len_dist, get_outer_len_dist

function get_dist_dict(output::SisMcmcOutput)
    d = Dict{Vector{Path}, Real}()
    N = length(output.sample)

    for p in output.sample
        myaddcounts!(d, p)
    end

    # normalise
    map!(x -> x/N, values(d))
    return d
end

function get_dist_dict(output::SimMcmcOutput)
    d = Dict{Multiset{Path}, Real}()
    N = length(output.sample)
    for p in output.sample
        myaddcounts!(d, Multiset(p))
    end
    # normalise
    map!(x -> x/N, values(d))
    return d
end

"""
`KL_error(output::SisMcmcOutput)`

Find the KL-divergence between the true distribution and the approximation obtained via
the MCMC sample `output`.

That is, if ``π(\\mathcal{S})`` is the true distribution and ``\\tilde{\\pi}(\\mathcal{S})`` is the approximation then this function will return

```math
\\sum_{\\mathcal{S}} \\tilde{\\pi}(\\mathcal{S}) \\log\\left( \\frac{\\tilde{\\pi}(\\mathcal{S})}{\\pi(\\mathcal{S})}\\right)
```

where we assume that ``0 \\cdot \\log(0) = 0``.
"""
function KL_error(output::SisMcmcOutput)
    d_approx = get_dist_dict(output)
    d_true = get_true_dist_dict(output.model)
    D = 0.0
    for key in keys(d_approx)
        D += d_approx[key] * ( log(d_approx[key]) - log(d_true[key]) )
    end
    return D
end

"""
`KL_error(output::SisMcmcOutput, true_dist::Dict)`

Takes the true distribution, as represented by a dictionary, as an additional parameter.
"""
function KL_error(output::SisMcmcOutput, true_dist::Dict)
    d_approx = get_dist_dict(output)
    D = 0.0
    for key in keys(d_approx)
        D += d_approx[key] * ( log(d_approx[key]) - log(true_dist[key]) )
    end
    return D
end

"""
`KL_error(output::SimMcmcOutput, true_dist::Dict)`

Takes the true distribution, as represented by a dictionary, as an additional parameter.
"""
function KL_error(output::SimMcmcOutput, true_dist::Dict)
    d_approx = get_dist_dict(output)
    D = 0.0
    for key in keys(d_approx)
        D += d_approx[key] * ( log(d_approx[key]) - log(true_dist[key]) )
    end
    return D
end

function get_outer_len_dist(mcmc_output::SisMcmcOutput, α::Float64)
    p = proportions(length.(mcmc_output.sample), 1:mcmc_output.model.K_outer)
    p .+= α
    p /= sum(p)
    return Categorical(p)
end

function get_inner_len_dist(mcmc_output::SisMcmcOutput, α::Float64)
    p = zeros(Int, mcmc_output.model.K_inner)
    for S in mcmc_output.sample
        addcounts!(p, length.(S), 1:mcmc_output.model.K_inner)
    end
    p /= sum(p)
    p .+= α
    p /= sum(p)
    return Categorical(p)
end


text:
include("branch_and_price.jl")

function solve_multi_depot(data, graph, approx, verbose)
    """Solve problem in the case of a single-depot dataset."""

    if approx
        println("Algorithm used : Branch-and-Price with an approximated solution")
    else
        println("Algorithm used : Branch-and-Price with an optimal solution")
    end
    println("Solving problem ...")
    (tours, solving_time) = branch_and_price(data, graph, approx, verbose)

    return(tours, solving_time)
end

max_stars_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_issues_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_forks_repo_name: TheoGuyard/Groomer-Optimization-Problem


text:
using WebIO
using Blink
using Observables
using Test

notinstalled = !AtomShell.isinstalled()

notinstalled && AtomShell.install()

"""
Execute function f() with a timeout of `timeout` seconds. Returns the
result of f() or `nothing` in the case of a timeout.
"""
function with_timeout(f::Function, timeout)
    c = Channel{Any}(1)
    @async begin
        put!(c, f())
    end
    @async begin
        sleep(timeout)
        put!(c, nothing)
    end
    take!(c)
end

if !(@isdefined open_window)
    include("./test-utils.jl")
end

# IMPORTANT: Cannot open Window's inside of @testsets.
# See https://github.com/JunoLab/Blink.jl/pull/193 for information.
w = open_window()
@testset "Blink mocks" begin
    body!(w, dom"div"("hello, blink"))
    sleep(5) # wait for it to render.

    content = Blink.@js(w, document.body.innerHTML)
    @test occursin("
hello, blink
", content)

    @testset "round-trip communication" begin
        scope = Scope()
        js_to_julia = Observable(scope, "inbox", "")
        output = Channel{String}(1)
        on(js_to_julia) do x
            # Put the result in a channel so we can watch for it
            put!(output, x)
        end

        julia_to_js = Observable(scope, "outbox", "")
        onjs(julia_to_js, js"""
        function (val) {
            _webIOScope.setObservableValue("inbox", val);
        }
        """)

        # Render the scope into the Blink window
        body!(w, scope)

        # Send a message and make sure we get it back from JS
        julia_to_js[] = "hello Blink"
        @test with_timeout(() -> take!(output), 5) == "hello Blink"
    end

    function scope_import(w::Window, url::AbstractString, use_iframe=false)
        scope = Scope(
            imports=[url]
        )
        js_to_julia = Observable{Any}(scope, "inbox", "")
        output = Channel{Any}(1)
        on(js_to_julia) do x
            # Put the result in a channel so we can watch for it
            put!(output, x)
        end
        onimport(
            scope,
            js"""
            function (mod) {
                _webIOScope.setObservableValue("inbox", mod.x);
            }
            """)
        if use_iframe
            body!(w, iframe(scope))
        else
            body!(w, scope)
        end

        return with_timeout(() -> take!(output), 5)
    end

    @testset "scope imports" begin
        for use_iframe in (false, true)
            @testset "local package, AssetRegistry" begin
                @test scope_import(w, joinpath(@__DIR__, "assets", "trivial_import.js"), use_iframe) == "ok"
            end

            @testset "global URL, no http:" begin
                # TODO: change this to a permanent URL because this CSAIL account
                # will eventually expire.
                @test scope_import(w, "//people.csail.mit.edu/rdeits/webio_tests/trivial_import.js", use_iframe) == "ok"
            end

            @testset "global URL, with http:" begin
                # TODO: change this to a permanent URL because this CSAIL account
                # will eventually expire.
                @test scope_import(w, "http://people.csail.mit.edu/rdeits/webio_tests/trivial_import.js", use_iframe) == "ok"
            end
        end
    end
end

window = open_window()
@testset "SVG/Namespaces" begin
    n = 10
    color = "yellow"
    h, w = 100, 100
    attributes = Dict(
        "fill" => color,
        "points" => join(["$(w/2*(1+sin(θ))),$(h/2*(1+cos(θ)))" for θ in 0:2π/n:2π], ' '),
    )
    my_svg = dom"svg:svg[width=$w, height=$h]"(
        dom"svg:polygon"(attributes=attributes),
        attributes=Dict("data-test-key" => "1234"),
    )
    body!(window, dom"div"(my_svg))
    sleep(5) # wait for it to render.
    child_count = @js window document.querySelector("svg[data-test-key=\"1234\"]").childElementCount
    @test child_count == 1
end

example_renderable_was_rendered = false

struct ExampleRenderableType
end

function WebIO.render(::ExampleRenderableType)
    global example_renderable_was_rendered
    example_renderable_was_rendered = true
    node(:div, "hello world")
end

@WebIO.register_renderable ExampleRenderableType

w = open_window()
@testset "register_renderable" begin
    body!(w, ExampleRenderableType())
    @test example_renderable_was_rendered
end

notinstalled && AtomShell.uninstall()


text:
"""
    module BookBayes

Library of functions given in Algorithms for Decision Making
"""
module BookBayes

using Reexport
using NamedTupleTools
using LightGraphs
using IterTools
using Distributions
using SpecialFunctions
using LinearAlgebra
using JuMP
using GLPK

export
        Variable,
        Factor,
        variablenames,
        namedtuple,
        assignments,
        normalize!,
        FactorTable,
        BayesianNetwork,
        probability

using Base.Iterators: product

"""
    const DAG = SimpleDiGraph

Data structure representing a directed acyclic graph
"""
const DAG = SimpleDiGraph

include("ch2.jl")
include("ch3.jl")
include("ch4.jl")
include("ch5.jl")
include("ch6.jl")
include("ch7_Exact_Solution_Methods.jl")
include("ch8_Approximate_Value_Functions.jl")
include("ch9_OnlinePlanning.jl")
include("ch10_policy_search.jl")
include("ch11_PolicyGradientEstimation.jl")
include("ch12_PolicyGradientOptimization.jl")

end # module


text:
"""
    get_min_type(n)

Return the minimum Integer type which can contain `n`
"""
function get_min_type(n::Integer)
    (n >= 64 ) && error("scale must be less than equal to 64")
    (n < 8)  && return UInt8
    (n < 16) && return UInt16
    (n < 32) && return UInt32
    (n < 64) && return UInt64
end

"""
    function kronecker_generator(scale, edgefactor; replicate, seed, A=0.57, B=0.19, C=0.19)

Generate a edge_list based on kronecker algorithm(http://graph500.org/?page_id=12#sec-3)
with kronecker parameters `A`, `B` and `C`. `scale` is the logarithm base two of the number
of vertices. `edgefactor` is the ratio of the graph’s edge count to its vertex count (i.e.,
half the average degree of a vertex in the graph)

###
References
- http://graph500.org/?page_id=12#alg:generator
"""

function kronecker_generator(
  scale::Integer,
  edgefactor::Integer;
  replicate::Bool=false,
  seed::AbstractArray=[],
  A::AbstractFloat=0.57,
  B::AbstractFloat=0.19,
  C::AbstractFloat=0.19
  )
    N  = 2^scale            # Set number of vertices
    M  = edgefactor * N     # Set number of edges

    # getting correct Inttype
    T = get_min_type(scale)

    # loop over each order of bit
    ab = A + B
    c_norm = C/(1 - (A + B))
    a_norm = A/(A + B)

    ij = ones(T, 2, M)    # Create index arrays
    T_one = T(1)
    T_two = T(2)

    # seeding the processors if replicate is true
    # if replicate
    #     for i in 1:nprocs()
    #         if replicate
    #             @spawnat i rng =  MersenneTwister(seed[i])
    #         else
    #             @spawnat i rng = MersenneTwister()
    #         end
    #     end
    # end

    temp =  @distributed (+) for ib in 1:T(scale)
        # Compare with probabilities and set bits of indices.
        if replicate
            rng = Random.seed!(seed[myid()])
        else
            rng = Random.seed!()
        end
        random_bits = falses(2, M)
        random_bits[1, :] = rand(rng, M) .> (ab)
        random_bits[2, :] = rand(rng, M) .> ( c_norm.*(random_bits[1, :]) + a_norm.*.!(random_bits[1, :]) )
        T_two^(ib-T_one).*random_bits
    end

    ij += temp

    # Permute vertex labels
    p = randperm(N)
    ij[1:2, :] = p[ij[1:2, :]]

    return ij
end


text:
# To pass incomplete specifications, we provide the "only_"
# family of functions. It allows the users to simply provide
# fg! for optimization or fj! for solving systems of equations
# and have NLSolversBase wrap them up properly so "F" can be
# calculated on its own and the save for partial derivatives.
# Note, that for TwiceDifferentiable we cannot provide con-
# structors if h == nothing, as that requires automatic dif-
# fferentiation of some sort.
struct InplaceObjective{DF, FDF, FGH}
    df::DF
    fdf::FDF
    fgh::FGH
end
struct NotInplaceObjective{DF, FDF, FGH}
    df::DF
    fdf::FDF
    fgh::FGH
end
# Mutating version
only_fg!(fg)   = InplaceObjective(nothing, fg,      nothing)
only_fgh!(fgh) = InplaceObjective(nothing, nothing, fgh)
only_fj!(fj)   = InplaceObjective(nothing, fj,      nothing)

# Non-mutating version
only_fg(fg)  = NotInplaceObjective(nothing, fg,      nothing)
only_fj(fj)  = NotInplaceObjective(nothing, fj,      nothing)

only_g_and_fg(g, fg) = NotInplaceObjective(g, fg, nothing)
only_j_and_fj(j, fj) = NotInplaceObjective(j, fj, nothing)

df(t::Union{InplaceObjective, NotInplaceObjective}) = t.df
fdf(t::Union{InplaceObjective, NotInplaceObjective}) = t.fdf

# Mutating version
make_f(t::InplaceObjective, x, F::Real) = x -> fdf(t)(F, nothing, x)
make_f(t::InplaceObjective, x, F) =  (F, x) -> fdf(t)(F, nothing, x)
make_df(t::InplaceObjective, x, F) = (DF, x) -> fdf(t)(nothing, DF, x)
make_fdf(t::InplaceObjective, x, F::Real) = (G, x) -> fdf(t)(F, G, x)
make_fdf(t::InplaceObjective, x, F) = fdf(t)

# Non-mutating version
# The contract with the user is that fdf returns (F, DF)
# and then we simply need to pick out the appropriate element
# of whatever fdf returns.
make_f(t::NotInplaceObjective, x, F::Real) = x -> fdf(t)(x)[1]
make_f(t::NotInplaceObjective, x, F) = (F, x) -> copyto!(F, fdf(t)(x)[1])
make_df(t::NotInplaceObjective{DF, TDF}, x, F) where {DF<:Nothing, TDF} = (DF, x) -> copyto!(DF, fdf(t)(x)[2])
make_df(t::NotInplaceObjective, x, F) = t.df
function make_fdf(t::NotInplaceObjective, x, F::Real)
    return function ffgg!(G, x)
        f, g = fdf(t)(x)
        copyto!(G, g)
        f
    end
end
function make_fdf(t::NotInplaceObjective, x, F)
    return function ffjj!(F, J, x)
        f, j = fdf(t)(x)
        copyto!(J, j)
        copyto!(F, f)
    end
end

# Constructors
function NonDifferentiable(t::Union{InplaceObjective, NotInplaceObjective}, x::AbstractArray, F::Real = real(zero(eltype(x))))
    f = make_f(t, x, F)
    NonDifferentiable(f, x, F)
end
# this would not be possible if we could mark f, g, ... as non-AbstractArrays
function NonDifferentiable(t::Union{InplaceObjective, NotInplaceObjective}, x::AbstractArray, F::AbstractArray)
    f = make_f(t, x, F)
    NonDifferentiable(f, x, F)
end

function TwiceDifferentiable(t::InplaceObjective{<: Nothing, <: Nothing, TH}, x::AbstractArray, F::Real = real(zero(eltype(x))), G::AbstractArray = similar(x), H = alloc_H(x)) where {TH}
    f   =     x  -> t.fgh(F, nothing, nothing, x)
    df  = (G, x) -> t.fgh(nothing, G, nothing, x)
    fdf = (G, x) -> t.fgh(F, G, nothing, x)
    h   = (H, x) -> t.fgh(F, nothing, H, x)
    TwiceDifferentiable(f, df, fdf, h, x, F, G, H)
end

function TwiceDifferentiable(t::InplaceObjective{<: Nothing, <: Nothing, TH}, x::AbstractVector, F::Real = real(zero(eltype(x))), G::AbstractVector = similar(x)) where {TH}

    H = alloc_H(x)
    f   =     x  -> t.fgh(F, nothing, nothing, x)
    df  = (G, x) -> t.fgh(nothing, G, nothing, x)
    fdf = (G, x) -> t.fgh(F, G, nothing, x)
    h   = (H, x) -> t.fgh(F, nothing, H, x)
    TwiceDifferentiable(f, df, fdf, h, x, F, G, H)
end


text:
"""
## Indexing
`pa` is an instance of `SeisForwExpt`.
* `pa[:data,i]` : get data for ith supersource (defaults to first), after `update!`
* `pa[:ageom]` : get `ageom` originally input while creating `pa`
* `pa[:srcwav]`: `srcwav` used to create `pa`
* `pa[:snaps,i]` : snapshots corresponding to `tsnaps` and ith source, after `update!`
"""
function Base.getindex(pa::PFdtd, s::Symbol, iss::Int=1)
	@assert iss≤length(pa.c.ageom[1])
	@assert s in [:snaps,:data, :ageom, :srcwav]
	if(s==:snaps)
		nzd,nxd=length.(pa.c.model.mgrid)
		snaps_all=zeros(nzd,nxd,length(pa.c.itsnaps))
		ip=findall(in.(iss,pa.sschunks))[1]
		issp=findall(x->x==iss,pa.sschunks[ip])[]
		@sync begin
			p=procs(pa.p)[ip]
			@sync remotecall_wait(p) do
				pap=localpart(pa.p)
				for i in eachindex(snaps_all)
					snaps_all[i]=pap.ss[issp].snaps[i]
				end
			end
		end
		return snaps_all
	elseif(s==:data)
		(iss==0) ?  (return pa.c.data[1]) :  (return pa.c.data[1][iss])
	elseif(s==:ageom)
		return pa.c.ageom
	elseif(s==:srcwav)
		return pa.c.srcwav
	end
end

function Base.show(io::Base.IO, pa::PFdtd)
	println(typeof(pa), "")
	println("pa[:snaps]")
	println("pa[:data] : modeled data after running `update!`")
	println("pa[:ageom]")
	println("pa[:srcwav]")
end

Base.show(io::Base.IO, pa::PFdtdp)=nothing
Base.show(io::Base.IO, pa::PFdtdc)=nothing
Base.show(io::Base.IO, pa::PFdtdss)=nothing


text:
using Base.Test

using MarketData
using TimeSeries

using MarketTechnicals

@testset "Volume" begin

@testset "obv" begin
    @test obv(ohlcv).values[1]     == 4783900 # TTR value is 4783900
    @test obv(ohlcv).values[12]    == 9390200 # TTR value is 2000-01-19   9390200
    @test obv(ohlcv).timestamp[12] == Date(2000,1,19)
end

@testset "vwap" begin
    # TTR value 2000-01-14 97.92154
    @test isapprox(vwap(ohlcv).values[1]   , 97.9215, atol=.01)
    # TTR value  21.44458
    @test isapprox(vwap(ohlcv).values[end] , 21.4446, atol=.01)
    @test vwap(ohlcv).timestamp[1] == Date(2000,1,14)
end

@testset "adl" begin
    """
    Quote from TTR

    > chaikinAD(x[,c('High', 'Low', 'Close')], x['Volume'])
        [1]     4288250.8     984498.1   -4126362.8  -10983262.8   -9113399.1
        [6]   -10015299.1  -11959861.3  -19338476.6  -16019132.6  -16929977.6
        [11]  -15868158.9  -14884408.6  -31233808.6  -33218729.3  -35989170.7
        [16]  -32269505.9  -34897972.6  -34897972.6  -37941049.9  -31849706.7
        [21]  -34405346.7  -35619584.9  -33370900.9  -31954278.7  -28196260.3
        [26]  -26405674.2  -28872966.5  -26697651.1  -28269906.1  -25051587.6
        [31]  -22431126.6  -23556859.9  -22318762.4  -24060310.9  -22592574.8
        [36]  -21272543.5  -21486881.4  -23543250.6  -22161637.3  -22562010.6
        [41]  -15427942.4  -17202733.1  -14476907.6  -15727198.9  -16807148.2
        [46]  -16131953.1  -15674416.0  -14855621.8  -16166047.8  -19805984.4
        [51]  -21014254.5  -18034460.1  -15310860.1  -16551649.3  -12990243.8
        [56]   -8069166.1  -11824243.5  -12798179.3  -13593182.8  -13919972.9
        [61]  -14450551.5  -17963286.2  -15275363.0  -15951899.0  -14184973.1
    """
    ta = adl(ohlcv)

    @test ta.meta           == ohlcv.meta
    @test ta.colnames       == ["adl"]
    @test ta.timestamp[1]   == ohlcv.timestamp[1]
    @test ta.timestamp[end] == ohlcv.timestamp[end]
    @test isapprox(ta.values[1], 4.2882507863e6, atol=.01)
    @test isapprox(ta.values[2], 984498.0822, atol=.01)
end

end  # @testset "Volume"


text:
FIXTURES = joinpath(@__DIR__, "fixtures")
SBML_FILE = joinpath(FIXTURES, "_model_v1.xml")

# test getmodel
@test_nowarn getmodel(SBML_FILE)
model = getmodel(SBML_FILE)

# test getparameters
parameters = getparameters(model)
trueparameters = Pair{ModelingToolkit.Num,Float64}[a0 => 1.0,
                                                   b0 => 0.0,
                                                   k1 => 0.8,
                                                   k2 => 0.6,
                                                   compartment => 1.0]
@test repr(parameters) == repr(trueparameters)

# test getinitialconditions
u0 = getinitialconditions(model)
true_u0 = Pair{ModelingToolkit.Num,Float64}[A => 1.0,
                                            B => 0.0]
@test repr(u0) == repr(true_u0)

# test getodes
eqs = getodes(model)
trueeqs = ModelingToolkit.Equation[
            D(B) ~ 1.0 * (compartment * k1 * A) - 1.0 * (compartment * k2 * B),
            D(A) ~ -1.0 * (compartment * k1 * A) + 1.0 * (compartment * k2 * B)
          ]
@test repr(eqs) == repr(trueeqs)
sys = ModelingToolkit.ODESystem(eqs)
@test repr(ModelingToolkit.get_iv(sys)) == "t"
@test repr(ModelingToolkit.get_states(sys)) == "Term{Real}[B(t), A(t)]"

# test sbml2odesystem
sys,u0,p = sbml2odesystem(SBML_FILE)
@test repr(ModelingToolkit.get_iv(sys)) == "t"
@test repr(ModelingToolkit.get_states(sys)) == "Term{Real}[B(t), A(t)]"

# test sbml2odeproblem
prob = sbml2odeproblem(SBML_FILE,(0.0,10.0))
@test prob.u0 == [0.0, 1.0]  # "Pair{ModelingToolkit.Num,Float64}[A => 1.0,B => 0.0]"
@test prob.tspan == (0.0, 10.0)
@test Set(prob.p) == Set([0.8, 0.6, 1.0])  # Todo: try removing the Set() when new ODEproblem version is released.
@test_nowarn OrdinaryDiffEq.solve(prob,OrdinaryDiffEq.Tsit5())
@test_nowarn sbml2odeproblem(SBML_FILE,(0.0,1.0),jac=false)

# test simulatesbml
sol1 = simulatesbml(SBML_FILE,(0.0,10.0))
@test sol1[end] ≈ [0.5714271572277619, 0.42857284277223784]
@test_nowarn simulatesbml(SBML_FILE,(0.0,1.0),saveat=0.5,jac=false,solver=OrdinaryDiffEq.Vern9())
@test_nowarn simulatesbml(joinpath(FIXTURES, "model_Vinod_FEBS2015.xml"),(0.0,1.0))


text:
Vk.@bitmask_flag MemoryAccess::UInt32 begin
    READ = 1
    WRITE = 2
    NO_ACCESS = 4
end

struct Attachment{IV<:ImageView}
    view::IV
    usage::MemoryAccess
end

@forward Attachment.view dims, format, samples, image_layout, image

function load_op(usage::MemoryAccess, clear::Bool)
    @match usage begin
        &READ => Vk.ATTACHMENT_LOAD_OP_LOAD
        &WRITE || &(READ | WRITE) => clear ? Vk.ATTACHMENT_LOAD_OP_CLEAR : Vk.ATTACHMENT_LOAD_OP_DONT_CARE
        _ => Vk.ATTACHMENT_LOAD_OP_DONT_CARE
    end
end

function store_op(usage::MemoryAccess)
    @match usage begin
        &WRITE || &(READ | WRITE) => Vk.ATTACHMENT_STORE_OP_STORE
        &READ => Vk.ATTACHMENT_STORE_OP_DONT_CARE
        _ => nothing
    end
end

function Vk.AttachmentDescription2(att::Attachment, clear::Bool, initial_layout::Vk.ImageLayout, final_layout::Vk.ImageLayout, aspect::Vk.ImageAspectFlag)
    _load_op, _store_op = if Vk.IMAGE_ASPECT_COLOR_BIT in aspect
        (load_op(att.usage, clear), store_op(att.usage))
    else
        (Vk.ATTACHMENT_LOAD_OP_DONT_CARE, Vk.ATTACHMENT_STORE_OP_DONT_CARE)
    end

    stencil_load_op, stencil_store_op = if Vk.IMAGE_ASPECT_STENCIL_BIT in aspect
        (load_op(att.usage, clear), store_op(att.usage))
    else
        (Vk.ATTACHMENT_LOAD_OP_DONT_CARE, Vk.ATTACHMENT_STORE_OP_DONT_CARE)
    end

    Vk.AttachmentDescription2(
        format(att),
        samples(att),
        _load_op,
        _store_op,
        stencil_load_op,
        stencil_store_op,
        initial_layout,
        final_layout,
    )
end

@auto_hash_equals struct TargetAttachments
    color::Vector{Symbol}
    depth::Vector{Symbol}
    stencil::Vector{Symbol}
end

TargetAttachments(color; depth = [], stencil = []) = TargetAttachments(color, depth, stencil)


text:
mutable struct PopulationData
	pops::Vector{String}
	npop::Int
	outgroup::Int
	f3::Matrix{Float64}
	cov::Array{Float64,4}
end

function PopulationData(
	f3file::String,
	covfile::String)

	f3dat  = CSV.read(f3file)
	covdat = CSV.read(covfile)

    # save populations
    pops = sort(vcat(unique(f3dat[:Outgroup]),unique(f3dat[:A])))
	npop = length(pops)
    popinds = Dict(zip(pops,1:npop))
    outgroup = findfirst(pops .== f3dat[:Outgroup][1])

    # get f3(Outgroup; A, B) stats
    f3 = zeros(npop, npop)
    for i in 1:DataFrames.nrow(f3dat)
        a = popinds[f3dat[i,:A]]
        b = popinds[f3dat[i,:B]]
        f3[a,b] = f3dat[i,:f3]
        f3[b,a] = f3dat[i,:f3]
    end

    # get covariance stats
    cov = zeros(npop, npop, npop, npop);
    for i in 1:DataFrames.nrow(covdat)
        a1 = min(popinds[covdat[i,:A1]],popinds[covdat[i,:B1]])
        b1 = max(popinds[covdat[i,:A1]],popinds[covdat[i,:B1]])
        a2 = min(popinds[covdat[i,:A2]],popinds[covdat[i,:B2]])
        b2 = max(popinds[covdat[i,:A2]],popinds[covdat[i,:B2]])
        cov[a1,b1,a2,b2] = covdat[i,:covariance]
        cov[b1,a1,b2,a2] = covdat[i,:covariance]
        cov[a2,b2,a1,b1] = covdat[i,:covariance]
        cov[b2,a2,b1,a1] = covdat[i,:covariance]
    end

    PopulationData(pops,npop,outgroup,f3,cov)
end

function Base.show(io::IO, pd::PopulationData; offset::String="")
    println(io, offset, string(typeof(pd)))
    println(io, offset, "    * npop: $(pd.npop)  (pops: $(pd.pops[1]), ..., $(pd.pops[end]))")
end


text:
function PhaseMinusCodeCycleSlipDetector(arg0::jdouble, arg1::jdouble, arg2::jint, arg3::jint)
    return PhaseMinusCodeCycleSlipDetector((jdouble, jdouble, jint, jint), arg0, arg1, arg2, arg3)
end


text:
kon = rand()
koff = rand()

c1 = chem2diff("""
               $(kon),$(koff): 2 A -> B
          """)
c2 = chem2diff("""
               $(kon),$(koff): 2A -> B
               """)

c3 = chem2diff("""
               $(kon),$(koff):A +A -> B
               """)

n = 2
z = 1.0*ones(n)
dz = zeros(n)

dz[1] = 2*(-kon*z[1]^2 + koff*z[2])
dz[2] = kon*z[1]^2 - koff*z[2]

#This works only because of the fixed usage of `y`,`dy`
y = ones(n)
dy = zeros(n)

eval(c1)
@test all(dy .≈ dz)

eval(c2)
@test all(dy .≈ dz)

eval(c3)
@test all(dy .≈ dz)


text:
using NLS_Solver
using Documenter

DocMeta.setdocmeta!(NLS_Solver, :DocTestSetup, :(using NLS_Solver); recursive=true)
doctest(NLS_Solver, fix=true)


text:
"""
Bridges `CP.SortPermutation` to MILP constraints by adding O(n²) binary
variables, with a transportation-like model.

The implemented model is highly similar to that of [`Sort2MILPBridge`](@ref).
"""
struct SortPermutation2MILPBridge{T} <: MOIBC.AbstractBridge
    vars_flow::Matrix{MOI.VariableIndex}
    vars_unicity::Matrix{MOI.VariableIndex}
    vars_unicity_bin::Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
    cons_transportation_x::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
    cons_transportation_y::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
    cons_unicity_x::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
    cons_unicity_y::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
    cons_flow_gt::Matrix{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}}
    cons_flow_lt::Matrix{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}
    cons_sort::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}}
    cons_perm::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
end

function MOIBC.bridge_constraint(
    ::Type{SortPermutation2MILPBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.SortPermutation,
) where {T}
    return MOIBC.bridge_constraint(
        SortPermutation2MILPBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{SortPermutation2MILPBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.SortPermutation,
) where {T}
    f_scalars = MOIU.scalarize(f)
    dim = s.dimension

    f_sorted = f_scalars[1:dim]
    f_to_sort = f_scalars[(dim + 1):(2 * dim)]
    f_perm = f_scalars[(2 * dim + 1):(3 * dim)]

    # For this formulation work, both lower and upper bounds are required on
    # the values of the array to sort.
    for i in 1:dim
        @assert CP.has_lower_bound(model, f_to_sort[i])
        @assert CP.has_upper_bound(model, f_to_sort[i])
    end

    # Create the new variables.
    vars_flow = Matrix{MOI.VariableIndex}(undef, dim, dim)
    vars_unicity = Matrix{MOI.VariableIndex}(undef, dim, dim)
    vars_unicity_bin = Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}(undef, dim, dim)
    for i in 1:dim
        vars_flow[i, :] = MOI.add_variables(model, dim)
        vars_unicity[i, :], vars_unicity_bin[i, :] = MOI.add_constrained_variables(model, [MOI.ZeroOne() for _ in 1:dim])
    end

    # Write the transportation equations.
    cons_transportation_x = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}[
        MOI.add_constraint(
            model,
            f_sorted[i] - sum(one(T) .* vars_flow[i, :]),
            MOI.EqualTo(zero(T))
        )
        for i in 1:dim
    ]
    cons_transportation_y = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}[
        MOI.add_constraint(
            model,
            f_to_sort[j] - sum(one(T) .* vars_flow[:, j]),
            MOI.EqualTo(zero(T))
        )
        for j in 1:dim
    ]

    # Unicity constraints, so that only one input and one output contribute to each flow.
    cons_unicity_x = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}[
        MOI.add_constraint(
            model,
            sum(one(T) .* vars_unicity[i, :]),
            MOI.EqualTo(one(T))
        )
        for i in 1:dim
    ]
    cons_unicity_y = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}[
        MOI.add_constraint(
            model,
            sum(one(T) .* vars_unicity[:, j]),
            MOI.EqualTo(one(T))
        )
        for j in 1:dim
    ]

    # Relate flows and unicity.
    cons_flow_gt = Matrix{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}}(undef, dim, dim)
    cons_flow_lt = Matrix{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}(undef, dim, dim)
    for i in 1:dim
        L = T(CP.get_lower_bound(model, f_to_sort[i]))
        U = T(CP.get_upper_bound(model, f_to_sort[i]))

        for j in 1:dim
            cons_flow_gt[i, j] = MOI.add_constraint(
                model,
                vars_flow[i, j] - U * vars_unicity[i, j],
                MOI.GreaterThan(zero(T))
            )
            cons_flow_lt[i, j] = MOI.add_constraint(
                model,
                vars_flow[i, j] - L * vars_unicity[i, j],
                MOI.LessThan(zero(T))
            )
        end
    end

    # Ensure that the end array is sorted.
    cons_sort = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}[
        MOI.add_constraint(
            model,
            f_sorted[i] - f_sorted[i + 1],
            MOI.GreaterThan(zero(T))
        )
        for i in 1:(dim - 1)
    ]

    # Constrain the variables of the permutation.
    cons_perm = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}[
        MOI.add_constraint(
            model,
            f_perm[i] - sum(T(j) * vars_unicity[i, j] for j in 1:dim),
            MOI.EqualTo(zero(T))
        )
        for i in 1:dim
    ]

    return SortPermutation2MILPBridge(
        vars_flow,
        vars_unicity,
        vars_unicity_bin,
        cons_transportation_x,
        cons_transportation_y,
        cons_unicity_x,
        cons_unicity_y,
        cons_flow_gt,
        cons_flow_lt,
        cons_sort,
        cons_perm,
    )
end

function MOI.supports_constraint(
    ::Type{SortPermutation2MILPBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.SortPermutation},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{SortPermutation2MILPBridge{T}}) where {T}
    return [(MOI.ZeroOne,)]
end

function MOIB.added_constraint_types(::Type{SortPermutation2MILPBridge{T}}) where {T}
    return [
        (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
    ]
end

function MOI.get(b::SortPermutation2MILPBridge, ::MOI.NumberOfVariables)
    return length(b.vars_flow) + length(b.vars_unicity)
end

function MOI.get(
    b::SortPermutation2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return length(b.vars_unicity)
end

function MOI.get(
    b::SortPermutation2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T},
    },
) where {T}
    return length(b.cons_flow_gt) + length(b.cons_sort)
end

function MOI.get(
    b::SortPermutation2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.LessThan{T},
    },
) where {T}
    return length(b.cons_flow_lt)
end

function MOI.get(
    b::SortPermutation2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return length(b.cons_transportation_x) + length(b.cons_transportation_y) + length(b.cons_unicity_x) + length(b.cons_unicity_y)
end

function MOI.get(
    b::SortPermutation2MILPBridge{T},
    ::MOI.ListOfVariableIndices,
) where {T}
    return vcat(vec(b.vars_flow), vec(b.vars_unicity))
end

function MOI.get(
    b::SortPermutation2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return copy(b.vars_unicity_bin)
end

function MOI.get(
    b::SortPermutation2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T},
    },
) where {T}
    return vcat(vec(b.cons_flow_gt), b.cons_sort)
end

function MOI.get(
    b::SortPermutation2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.LessThan{T},
    },
) where {T}
    return vec(b.cons_flow_lt)
end

function MOI.get(
    b::SortPermutation2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return vcat(b.cons_transportation_x, b.cons_transportation_y, b.cons_unicity_x, b.cons_unicity_y)
end


text:
function estimate(entity::FundamentalMatrix, method::DirectLinearTransform, 𝒟::Tuple{AbstractArray, Vararg{AbstractArray}})
    ℳ, ℳʹ =  𝒟
    N = length(ℳ)
    if (N != length(ℳʹ))
          throw(ArgumentError("There should be an equal number of points for each view."))
    end
    𝒪, 𝐓 = transform(HomogeneousCoordinates(),CanonicalToHartley(),ℳ)
    𝒪ʹ, 𝐓ʹ  = transform(HomogeneousCoordinates(),CanonicalToHartley(),ℳʹ)
    𝐀 = moments(FundamentalMatrix(), (𝒪, 𝒪ʹ))
    λ, f = smallest_eigenpair(Symmetric(𝐀))
    𝐅 = reshape(f,(3,3))
    𝐅 = enforce_ranktwo!(Array(𝐅))
    𝐅 = SMatrix{3,3,Float64,9}(𝐅 / norm(𝐅))
    # Transform estimate back to the original (unnormalised) coordinate system.
    𝐓ʹ'*𝐅*𝐓
end

# TODO fix numerical instability
# function estimate(entity::FundamentalMatrix, method::Taubin, matches...)
#     ℳ, ℳʹ = matches
#     N = length(ℳ)
#     if (N != length(ℳʹ))
#           throw(ArgumentError("There should be an equal number of points for each view."))
#     end
#     (ℳ,𝐓) = hartley_normalization(ℳ)
#     (ℳʹ,𝐓ʹ) = hartley_normalization(ℳʹ)
#     𝐀::Matrix{Float64} = moments(FundamentalMatrix(), ℳ, ℳʹ)
#     𝐁::Matrix{Float64} = mean_covariance(FundamentalMatrix(), ℳ, ℳʹ)
#     (λ::Float64, f::Vector{Float64}) = smallest_eigenpair(𝐀,𝐁)
#     𝐅::Matrix{Float64} = reshape(f,(3,3))
#     enforce_ranktwo!(𝐅)
#     # Transform estimate back to the original (unnormalised) coordinate system.
#     𝐅 = 𝐓ʹ'*𝐅*𝐓
# end

function estimate(entity::FundamentalMatrix, method::FundamentalNumericalScheme,  𝒞::Tuple{AbstractArray, Vararg{AbstractArray}}, 𝒟::Tuple{AbstractArray, Vararg{AbstractArray}})
    ℳ, ℳʹ = 𝒟
    Λ₀, Λ₀ʹ = 𝒞
    N = length(ℳ)
    if (N != length(ℳʹ))
          throw(ArgumentError("There should be an equal number of points for each view."))
    end
    if (N != length(Λ₀) || N != length(Λ₀ʹ) )
          throw(ArgumentError("There should be a covariance matrix for each point correspondence."))
    end
    # Map corresponding points to the normalized coordinate system.
    𝒪, 𝐓 = transform(HomogeneousCoordinates(),CanonicalToHartley(),ℳ)
    𝒪ʹ, 𝐓ʹ = transform(HomogeneousCoordinates(),CanonicalToHartley(),ℳʹ)
    # Map seed to the normalized coordinate system.
    𝛉 = (inv(𝐓') ⊗ inv(𝐓ʹ')) * method.𝛉₀
    # Map covariance matrices to the normalized coordinate system.
    Λ₁ = transform(CovarianceMatrices(), CanonicalToHartley(), Λ₀ , 𝐓)
    Λ₁ʹ = transform(CovarianceMatrices(), CanonicalToHartley(), Λ₀ʹ , 𝐓ʹ)
    for i = 1:method.max_iter
        𝐗 = X(AML(),FundamentalMatrix(), 𝛉, (Λ₁,Λ₁ʹ), (𝒪, 𝒪ʹ))
        λ, 𝛉⁺ = smallest_eigenpair(Symmetric(𝐗/N))
        𝛉 = reshape(𝛉⁺,length(𝛉⁺),1)
    end
    𝐅 = reshape(𝛉,(3,3))
    𝐅 = enforce_ranktwo!(Array(𝐅))
    𝐅 = SMatrix{3,3,Float64,9}(𝐅)
    # Transform estimate back to the original (unnormalised) coordinate system.
    𝐅 = 𝐓ʹ'*𝐅*𝐓
end

function estimate(entity::FundamentalMatrix, method::BundleAdjustment,  𝒟::Tuple{AbstractArray, Vararg{AbstractArray}})
    ℳ, ℳʹ = 𝒟
    N = length(ℳ)
    if (N != length(ℳʹ))
          throw(ArgumentError("There should be an equal number of points for each view."))
    end
    𝐅 = SMatrix{3,3,Float64,9}(reshape(method.𝛉₀,(3,3)))
    𝒳 = triangulate(DirectLinearTransform(),𝐅,(ℳ,ℳʹ))

    𝐏₁, 𝐏₂ = construct(ProjectionMatrix(),𝐅)

    # Construct a length-(12+3*N) vector consisting of the projection matrix associated
    # with the second view (the first twelve dimensions), as well as N three-dimensional points
    # (the remaining dimensions).
    𝛉 = pack(FundamentalMatrix(), 𝐏₂, 𝒳)

    index₁ = SVector(1,2)
    index₂ = SVector(3,4)
    pts = Matrix{Float64}(undef,4,N)
    for n = 1:N
        pts[index₁,n] = ℳ[n][index₁]
        pts[index₂,n] = ℳʹ[n][index₁]
    end
    fit = curve_fit(model_fundamental, jacobian_model,  𝐏₁, reshape(reinterpret(Float64,vec(pts)),(4*N,)) , 𝛉; show_trace = false)
    𝐏₂ = reshape(fit.param[1:12],(3,4))
    𝐅 = construct(FundamentalMatrix(), 𝐏₁, 𝐏₂)
    𝐅, fit
end

function model_fundamental(𝐏₁,𝛉)
    # Twelve parameters for the projection matrix, and 3 parameters per 3D point.
    N = Int((length(𝛉) - 12) / 3)
    index₁ = SVector(1,2)
    index₂ = SVector(3,4)
    reprojections = Matrix{Float64}(undef,4,N)
    𝛉v = @view 𝛉[1:12]
    𝐏₂ = SMatrix{3,4,Float64,12}(reshape(𝛉v,(3,4)))
    i = 13
    for n = 1:N
        # Extract 3D point and convert to homogeneous coordinates
        v = @view 𝛉[i:i+2]
        M = hom(SVector{3,Float64}(𝛉[i:i+2]))
        reprojections[index₁,n] = hom⁻¹(𝐏₁ * M)
        reprojections[index₂,n] = hom⁻¹(𝐏₂ * M)
        i = i + 3
    end
    reshape(reinterpret(Float64,vec(reprojections)),(4*N,))
end

function jacobian_model(𝐏₁,𝛉)
    # Twelve parameters for the projection matrix, and 3 parameters per 3D point.
    N = Int((length(𝛉) - 12) / 3)
    index₁ = SVector(1,2)
    index₂ = SVector(3,4)
    reprojections = Matrix{Float64}(undef,4,N)
    𝛉v = @view 𝛉[1:12]
    𝐏₂ = SMatrix{3,4,Float64,12}(reshape(𝛉v,(3,4)))
    𝐉 = zeros(4*N,12+3*N)
    # Create a view of the jacobian matrix 𝐉 and reshape it so that
    # it will be more convenient to index into the appropriate entires
    # whilst looping over all of the data points.
    𝐉v = reshape(reinterpret(Float64,𝐉), 4, N, 12+3*N)
    𝐀 = SMatrix{2,3,Float64,6}(1,0,0,1,0,0)
    𝐈₃ = SMatrix{3,3}(1.0I)
    i = 13
    for n = 1:N
        # Extract 3D point and convert to homogeneous coordinates.
        𝐌 = hom(SVector{3,Float64}(𝛉[i:i+2]))

        # Derivative of residual in first and second image w.r.t 3D point.
        ∂𝐫₁_d𝐌 = -𝐀 * ∂hom⁻¹(𝐏₁ * 𝐌) * 𝐏₁
        ∂𝐫₂_d𝐌 = -𝐀 * ∂hom⁻¹(𝐏₂ * 𝐌) * 𝐏₂

        # Derivative of residual in second image w.r.t projection martix
        # ∂𝐫₁_d𝐏₁ is the zero vector.
        ∂𝐫₂_d𝐏₂ = 𝐀 * ∂hom⁻¹(𝐏₂ * 𝐌) * (𝐌' ⊗ 𝐈₃)

        𝐉v[index₂,n,1:12] = ∂𝐫₂_d𝐏₂
        𝐉v[index₁,n,i:i+2] = ∂𝐫₁_d𝐌[:,1:3]
        𝐉v[index₂,n,i:i+2] = ∂𝐫₂_d𝐌[:,1:3]
        i = i + 3
    end
    𝐉
end

#𝛉 = reshape(𝛉₀,length(𝛉₀),1)
# function z(entity::FundamentalMatrix, 𝐦::Vector{Float64}, 𝐦ʹ::Vector{Float64})
# 𝐮 = 𝐦 ⊗ 𝐦ʹ
# 𝐮[1:end-1]
# end

function mean_covariance(entity::ProjectiveEntity, matches...)
    ℳ, ℳʹ = matches
    N = length(ℳ)
    if (N != length(ℳʹ))
          throw(ArgumentError("There should be an equal number of points for each view."))
    end
    𝐁 = fill(0.0,(9,9))
    𝚲 = eye(4)
    for correspondence in zip(ℳ, ℳʹ)
        m , mʹ = correspondence
        𝐦  = 𝑛(collect(Float64,m.coords))
        𝐦ʹ = 𝑛(collect(Float64,mʹ.coords))
        ∂ₓ𝐮 = ∂ₓu(entity, 𝐦 , 𝐦ʹ)
        𝐁 = 𝐁 + ∂ₓ𝐮 * 𝚲 * ∂ₓ𝐮'
    end
    𝐁/N
end

function enforce_ranktwo!(𝐅::AbstractArray)
    # Enforce the rank-2 constraint.
    U,S,V = svd(𝐅)
    S[end] = 0.0
    𝐅 = U*Matrix(Diagonal(S))*V'
end

# Construct a parameter vector consisting of a projection matrix and 3D points.
function pack(entity::FundamentalMatrix, 𝐏₂::AbstractArray, 𝒳::AbstractArray)
    N = length(𝒳)
    𝛉 = Vector{Float64}(undef,12+N*3)
    𝛉[1:12] = Array(𝐏₂[:])
    i = 13
    for n = 1:N
        𝛉[i:i+2] = 𝒳[n][1:3]
        i = i + 3
    end
    𝛉
end


text:
module MathTeXEngine
# Adapted from matplotlib mathtext engine
# https://github.com/matplotlib/matplotlib/blob/master/lib/matplotlib/_mathtext.py

using AbstractTrees
using Automa
using FreeTypeAbstraction
using LaTeXStrings

import Automa.RegExp: @re_str
import DataStructures: Stack
import FreeTypeAbstraction:
    ascender, boundingbox, descender, get_extent, hadvance, inkheight, inkwidth,
    leftinkbound, rightinkbound, topinkbound, bottominkbound
import GeometryBasics: Point2f
import REPL.REPLCompletions: latex_symbols
import RelocatableFolders

const re = Automa.RegExp

export TeXExpr, texparse
export TeXElement, TeXChar, VLine, HLine, generate_tex_elements
export get_font, get_fontpath

include("parser/texexpr.jl")
include("parser/commands_data.jl")
include("parser/commands_registration.jl")
include("parser/parser.jl")

include("engine/fonts.jl")
include("engine/layout_context.jl")
include("engine/texelements.jl")
include("engine/layout.jl")

end # module


text:
function LUDecomposition(arg0::RealMatrix)
    return LUDecomposition((RealMatrix,), arg0)
end

function LUDecomposition(arg0::RealMatrix, arg1::jdouble)
    return LUDecomposition((RealMatrix, jdouble), arg0, arg1)
end

function get_determinant(obj::LUDecomposition)
    return jcall(obj, "getDeterminant", jdouble, ())
end

function get_l(obj::LUDecomposition)
    return jcall(obj, "getL", RealMatrix, ())
end

function get_p(obj::LUDecomposition)
    return jcall(obj, "getP", RealMatrix, ())
end

function get_pivot(obj::LUDecomposition)
    return jcall(obj, "getPivot", Vector{jint}, ())
end

function get_solver(obj::LUDecomposition)
    return jcall(obj, "getSolver", DecompositionSolver, ())
end

function get_u(obj::LUDecomposition)
    return jcall(obj, "getU", RealMatrix, ())
end


text:
using Printf
import Flux: σ
using ModelingToolkit
using GalacticOptim
using Optim
using DiffEqFlux
using NeuralPDE
using Quadrature, Cubature, Cuba
using Plots

@parameters t,x
@variables c1(..),c2(..),c3(..),c4(..),c5(..),c6(..)
Dt = Differential(t)
Dx = Differential(x)

# Parameters
d_Fei = 1.
d_Cri = 1.
d_Nii = 1.
d_Sii = 1.
d_Fev = 1.
d_Crv = 1.
d_Niv = 1.
d_Siv = 1.
v = 1.
R1 = 0.
R2 = 0.
R3 = 0.
R4 = 0.
R5 = 0.
R6 = 0.

D = 0.1 # diffusion
α = 1.0
t_max = 2.0
x_min = 0.0
x_max = 1.0

f1 = α * c1(t, x) * (d_Fei * Dx(c3(x,t)) + d_Nii * Dx(c4(x,t)) + d_Cri * Dx(c5(x,t)) + d_Sii * Dx(c6(x,t)))+ Dx(c1(x,t)) * (d_Fei * c3(t, x) + d_Nii * c4(t, x) + d_Cri * c5(t, x)+ d_Sii * c6(t, x))
f2 = -α * c2(t, x) * (d_Fev * Dx(c3(x,t)) + d_Niv * Dx(c4(x,t)) + d_Crv * Dx(c5(x,t)) + d_Siv * Dx(c6(x,t)))+ Dx(c2(x,t)) * (d_Fev * c3(t, x) + d_Niv * c4(t, x) + d_Crv * c5(t, x)+ d_Siv * c6(t, x))
f3 = α * Dx(c3(x,t)) * (d_Fei * c1(t, x) + d_Fev * c2(t, x))  - c3(t, x) * (d_Fev * Dx(c2(x,t)) - d_Fei * Dx(c1(x,t)))
f4 = α * Dx(c4(x,t)) * (d_Nii * c1(t, x) + d_Niv * c2(t, x))  - c4(t, x) * (d_Niv * Dx(c2(x,t)) - d_Nii * Dx(c1(x,t)))
f5 = α * Dx(c5(x,t)) * (d_Cri * c1(t, x) + d_Crv * c2(t, x))  - c5(t, x) * (d_Crv * Dx(c2(x,t)) - d_Cri * Dx(c1(x,t)))
f6 = α * Dx(c6(x,t)) * (d_Sii * c1(t, x) + d_Siv * c2(t, x))  - c6(t, x) * (d_Siv * Dx(c2(x,t)) - d_Sii * Dx(c1(x,t)))

# Equations, initial and boundary conditions
eqs = [
Dt(c1(t, x)) ~  Dx(f1*c1(t,x)) + R1,
Dt(c2(t, x)) ~  Dx(f1*c2(t,x)) + R2,
Dt(c3(t, x)) ~  Dx(f1*c3(t,x)) + R3,
Dt(c4(t, x)) ~  Dx(f1*c4(t,x)) + R4,
Dt(c5(t, x)) ~  Dx(f1*c5(t,x)) + R5,
Dt(c6(t, x)) ~  Dx(f1*c6(t,x)) + R6 ]
bcs = [
        c1(0, x) ~  1.0,  c2(0, x) ~  1.0, c3(0, x) ~  0.68, c4(0, x) ~  0.2, c5(0, x) ~  0.1, c6(0, x) ~  0.02,
        c1(t, x_min) ~ c1(t, x_max),
        c2(t, x_min) ~ c2(t, x_max),
        Dx(c3(t, x_min)) ~ 0., Dx(c3(t, x_max)) ~ 0.,
        Dx(c4(t, x_min)) ~ 0., Dx(c4(t, x_max)) ~ 0.,
        Dx(c5(t, x_min)) ~ 0., Dx(c5(t, x_max)) ~ 0.,
        Dx(c6(t, x_min)) ~ 0., Dx(c6(t, x_max)) ~ 0.,

]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,t_max),
        x ∈ IntervalDomain(x_min,x_max)
]

# Discretization
nx = 32
dx = (x_max-x_min) / (nx - 1)
dt = 0.01

# Neural network
dim = length(domains)
output = length(eqs)
hidden = 32

# chain = FastChain( FastDense(dim, hidden, tanh),
#                     FastDense(hidden, hidden, tanh),
#                     FastDense(hidden, 1))

chain1 = FastChain(FastDense(dim,hidden,Flux.σ),FastDense(hidden,hidden,Flux.σ),FastDense(hidden,1))
chain2 = FastChain(FastDense(dim,hidden,Flux.σ),FastDense(hidden,hidden,Flux.σ),FastDense(hidden,1))
chain3 = FastChain(FastDense(dim,hidden,Flux.σ),FastDense(hidden,hidden,Flux.σ),FastDense(hidden,1))
chain4 = FastChain(FastDense(dim,hidden,Flux.σ),FastDense(hidden,hidden,Flux.σ),FastDense(hidden,1))
chain5 = FastChain(FastDense(dim,hidden,Flux.σ),FastDense(hidden,hidden,Flux.σ),FastDense(hidden,1))
chain6 = FastChain(FastDense(dim,hidden,Flux.σ),FastDense(hidden,hidden,Flux.σ),FastDense(hidden,1))

strategy = GridTraining(dx)

discretization = PhysicsInformedNN([chain1, chain2, chain3, chain4, chain5, chain6], strategy)

pde_system = PDESystem(eqs, bcs, domains, [t,x], [c1, c2, c3, c4, c5, c6])
prob = discretize(pde_system,discretization)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

res = GalacticOptim.solve(prob,Optim.BFGS(),cb=cb,maxiters=100)

max_stars_repo_name: drozzy/NeuralOceananigans

max_issues_repo_name: drozzy/NeuralOceananigans

max_forks_repo_name: drozzy/NeuralOceananigans


text:
:% s/DATE/\=strftime('%Y.%m.%d')/ge
:% s/YEAR/\=strftime('%Y')/ge
:% s/THEFUNCNAME/\=expand("%:t:r")/g
:% s/CLASS/\=expand("%:p:h:t:r")/g
:if search('<+CURSOR+>')
:  normal! "_da>
:endif
:finish
#
# CLASS/THEFUNCNAME
#
# Created: DATE
# Copyright: Steven E. Pav, YEAR
# Author: Steven E. Pav, steven@gilgamath.com

<+CURSOR+>

# doc
function funcname(blah::RealArray, p::Int)
  n = blah

  return ks
end

#for vim modeline: (do not edit)
# vim:ts=4:sw=4:tw=79:fdm=marker:fmr=FOLDUP,UNFOLD:cms=#%s:syn=julia:ft=julia:ai:si:et:cin:nu:fo=croql:cino=p0t0c5(0:

max_stars_repo_name: shabbychef/dotfiles

max_issues_repo_name: shabbychef/dotfiles

max_forks_repo_name: shabbychef/dotfiles


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCellDG
include("useful_routines.jl")
include("circular_bc_transformation_elasticity_solver.jl")

function compute_interface_L2_stress_error(polyorder, nelmts, width, penaltyfactor)
    K1, K2 = 247.0, 192.0
    mu1, mu2 = 126.0, 87.0
    lambda1 = lame_lambda(K1, mu1)
    lambda2 = lame_lambda(K2, mu2)
    stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)

    theta0 = -0.067
    stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)
    transfstress =
        CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

    meshwidth = [width, width]
    numqp = required_quadrature_order(polyorder) + 2

    dx = width / nelmts
    penalty = penaltyfactor / dx * 0.5 * (lambda1 + mu1 + lambda2 + mu2)

    basis = TensorProductBasis(2, polyorder)
    interfacecenter = [0.5, 0.5]
    interfaceradius = minimum(meshwidth) / 3.0
    outerradius = 2.0
    analyticalsolution = AnalyticalSolution(
        interfaceradius,
        outerradius,
        interfacecenter,
        lambda1,
        mu1,
        lambda2,
        mu2,
        theta0,
    )

    mesh, cellquads, facequads, interfacequads = construct_mesh_and_quadratures(
        meshwidth,
        nelmts,
        basis,
        interfacecenter,
        interfaceradius,
        numqp,
    )
    nodaldisplacement = nodal_displacement(
        mesh,
        basis,
        cellquads,
        facequads,
        interfacequads,
        stiffness,
        theta0,
        analyticalsolution,
        penalty,
    )
    parenterror, producterror = interface_stress_L2_error(
        nodaldisplacement,
        basis,
        interfacequads,
        stiffness,
        transfstress,
        theta0,
        mesh,
        x -> core_stress(analyticalsolution),
        x -> shell_stress(analyticalsolution, x),
    )

    parentstressnorm = integral_norm_on_interface(
        x -> core_stress(analyticalsolution),
        interfacequads,
        -1,
        mesh,
        4,
    )
    productstressnorm = integral_norm_on_interface(
        x -> shell_stress(analyticalsolution, x),
        interfacequads,
        +1,
        mesh,
        4,
    )

    return parenterror ./ parentstressnorm, producterror ./ productstressnorm
end

width = 1.0
polyorder = 2
penaltyfactor = 1e3
powers = [3,4,5]
nelmts = [2^p for p in powers]

dx = width ./ nelmts

stresserror = [compute_interface_L2_stress_error(polyorder, ne, width, penaltyfactor) for ne in nelmts]

parenterror = [st[1] for st in stresserror]
parenterror = Array(transpose(hcat(parenterror...)))

parentrates = mapslices(x->convergence_rate(dx,x),parenterror,dims=1)
parentrates = parentrates[:,[1,2,4]]

# @test all(parentrates .> 1.4)

producterror = [st[2] for st in stresserror]

producterror = Array(transpose(hcat(producterror...)))
productrates = mapslices(x->convergence_rate(dx,x),producterror,dims=1)

# @test all(productrates .> 1.4)