text:
int_rules_1_2_1_5 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.1.5*(a+b*x+c*x^2)^p*(d+e*x+f*x^2)^q =#
    @apply_utils Antiderivative((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~(p') * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~(q'), ~x) => (~c / ~f) ^ ~p * Antiderivative((~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~p + ~q), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~p, ~q], ~x) && (EqQ(~c * ~d - ~a * ~f, 0) && (EqQ(~b * ~d - ~a * ~e, 0) && ((IntegerQ(~p) || GtQ(~c / ~f, 0)) && (Not(IntegerQ(~q)) || LeafCount(~d + ~e * ~x + ~f * (~x) ^ 2) <= LeafCount(~a + ~b * ~x + ~c * (~x) ^ 2)))))
    @apply_utils Antiderivative((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~(q'), ~x) => (((~a) ^ IntPart(~p) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ FracPart(~p)) / ((~d) ^ IntPart(~p) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ FracPart(~p))) * Antiderivative((~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~p + ~q), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~p, ~q], ~x) && (EqQ(~c * ~d - ~a * ~f, 0) && (EqQ(~b * ~d - ~a * ~e, 0) && (Not(IntegerQ(~p)) && (Not(IntegerQ(~q)) && Not(GtQ(~c / ~f, 0))))))
    @apply_utils Antiderivative((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~(q'), ~x) => ((~a + ~b * ~x + ~c * (~x) ^ 2) ^ FracPart(~p) / ((4 * ~c) ^ IntPart(~p) * (~b + 2 * ~c * ~x) ^ (2 * FracPart(~p)))) * Antiderivative((~b + 2 * ~c * ~x) ^ (2 * ~p) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~p, ~q], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p)))
    @apply_utils Antiderivative((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(f') * (~x) ^ 2) ^ ~(q'), ~x) => ((~a + ~b * ~x + ~c * (~x) ^ 2) ^ FracPart(~p) / ((4 * ~c) ^ IntPart(~p) * (~b + 2 * ~c * ~x) ^ (2 * FracPart(~p)))) * Antiderivative((~b + 2 * ~c * ~x) ^ (2 * ~p) * (~d + ~f * (~x) ^ 2) ^ ~q, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~p, ~q], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p)))
    #= Antiderivative(((~a)+(~b')*(~x)+(~c')*(~x)^2)^(~p)*((~d')+(~e')*(~x)+(~f')*(~x)^2)^(~q),~x) :=   1/(2^(2*p+2*q+1)*c*(-c/(b^2-4*a*c))^p*(-f/(e^2-4*d*f))^q)* Subst(Antiderivative((1-x^2/(b^2-4*a*c))^p*(1+e*x^2/(b*(4*c*d-b*e)))^q,x),x, b+2*c*x) <-- FreeQ([a,b,c,d,e,f,p,q],x) && NeQ(b^2-4*a*c,0) &&  NeQ(e^2-4*d*f,0) && EqQ(c*e-b*f,0) && (IntegerQ(p) || GtQ(-c/(b^2-4*a*c),0)) && (IntegerQ(q) ||  GtQ(-f/(e^2-4*d*f),0)) =#
    #= Antiderivative(((~a)+(~b')*(~x)+(~c')*(~x)^2)^(~p)*((~d')+(~e')*(~x)+(~f')*(~x)^2)^(~q),~x) :=   (d+e*x+f*x^2)^q/(2^(2*p+2*q+1)*c*(-c/(b^2-4*a*c))^p*(-f*(d+e*x+f*x^2) /(e^2-4*d*f))^q)* Subst(Antiderivative((1-x^2/(b^2-4*a*c))^p*(1+e*x^2/(b*(4*c*d-b*e)))^q,x),x, b+2*c*x) <-- FreeQ([a,b,c,d,e,f,p,q],x) && NeQ(b^2-4*a*c,0) &&  NeQ(e^2-4*d*f,0) && EqQ(c*e-b*f,0) && (IntegerQ(p) || GtQ(-c/(b^2-4*a*c),0)) && Not(IntegerQ(q) ||  GtQ(-f/(e^2-4*d*f),0)) =#
    #= Antiderivative(((~a)+(~b')*(~x)+(~c')*(~x)^2)^(~p)*((~d')+(~e')*(~x)+(~f')*(~x)^2)^(~q),~x) :=   (a+b*x+c*x^2)^p*(d+e*x+f*x^2)^q/(2^(2*p+2*q+1)*c*(-c*(a+b*x+c*x^2)/( b^2-4*a*c))^p*(-f*(d+e*x+f*x^2)/(e^2-4*d*f))^q)* Subst(Antiderivative((1-x^2/(b^2-4*a*c))^p*(1+e*x^2/(b*(4*c*d-b*e)))^q,x),x, b+2*c*x) <-- FreeQ([a,b,c,d,e,f,p,q],x) && NeQ(b^2-4*a*c,0) &&  NeQ(e^2-4*d*f,0) && EqQ(c*e-b*f,0) =#
    @apply_utils Antiderivative((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q, ~x) => ((~b + 2 * ~c * ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q) / (((~b) ^ 2 - 4 * ~a * ~c) * (~p + 1)) - (1 / (((~b) ^ 2 - 4 * ~a * ~c) * (~p + 1))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~q - 1) * Simp(2 * ~c * ~d * (2 * ~p + 3) + ~b * ~e * ~q + (2 * ~b * ~f * ~q + 2 * ~c * ~e * (2 * ~p + ~q + 3)) * ~x + 2 * ~c * ~f * (2 * ~p + 2 * ~q + 3) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (LtQ(~p, -1) && (GtQ(~q, 0) && Not(IGtQ(~q, 0))))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(f') * (~x) ^ 2) ^ ~q, ~x) => ((~b + 2 * ~c * ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~f * (~x) ^ 2) ^ ~q) / (((~b) ^ 2 - 4 * ~a * ~c) * (~p + 1)) - (1 / (((~b) ^ 2 - 4 * ~a * ~c) * (~p + 1))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~f * (~x) ^ 2) ^ (~q - 1) * Simp(2 * ~c * ~d * (2 * ~p + 3) + (2 * ~b * ~f * ~q) * ~x + 2 * ~c * ~f * (2 * ~p + 2 * ~q + 3) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (LtQ(~p, -1) && (GtQ(~q, 0) && Not(IGtQ(~q, 0)))))
    @apply_utils Antiderivative((~(a') + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q, ~x) => ((2 * ~c * ~x) * (~a + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q) / ((-4 * ~a * ~c) * (~p + 1)) - (1 / ((-4 * ~a * ~c) * (~p + 1))) * Antiderivative((~a + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~q - 1) * Simp(2 * ~c * ~d * (2 * ~p + 3) + (2 * ~c * ~e * (2 * ~p + ~q + 3)) * ~x + 2 * ~c * ~f * (2 * ~p + 2 * ~q + 3) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f], ~x) && (NeQ((~e) ^ 2 - 4 * ~d * ~f) && (LtQ(~p, -1) && (GtQ(~q, 0) && Not(IGtQ(~q, 0)))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q, ~x) => (((2 * ~a * (~c) ^ 2 * ~e - (~b) ^ 2 * ~c * ~e) + (~b) ^ 3 * ~f + ~b * ~c * (~c * ~d - 3 * ~a * ~f) + ~c * ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (~b * ~e + 2 * ~a * ~f)) * ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~q + 1)) / (((~b) ^ 2 - 4 * ~a * ~c) * ((~c * ~d - ~a * ~f) ^ 2 - (~b * ~d - ~a * ~e) * (~c * ~e - ~b * ~f)) * (~p + 1)) - (1 / (((~b) ^ 2 - 4 * ~a * ~c) * ((~c * ~d - ~a * ~f) ^ 2 - (~b * ~d - ~a * ~e) * (~c * ~e - ~b * ~f)) * (~p + 1))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q * Simp(((2 * ~c * ((~c * ~d - ~a * ~f) ^ 2 - (~b * ~d - ~a * ~e) * (~c * ~e - ~b * ~f)) * (~p + 1) - ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (~b * ~e + 2 * ~a * ~f)) * (~a * ~f * (~p + 1) - ~c * ~d * (~p + 2))) - ~e * ((((~b) ^ 2 * ~c * ~e - 2 * ~a * (~c) ^ 2 * ~e) - (~b) ^ 3 * ~f) - ~b * ~c * (~c * ~d - 3 * ~a * ~f)) * (~p + ~q + 2)) + (2 * ~f * ((2 * ~a * (~c) ^ 2 * ~e - (~b) ^ 2 * ~c * ~e) + (~b) ^ 3 * ~f + ~b * ~c * (~c * ~d - 3 * ~a * ~f)) * (~p + ~q + 2) - ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (~b * ~e + 2 * ~a * ~f)) * (~b * ~f * (~p + 1) - ~c * ~e * (2 * ~p + ~q + 4))) * ~x + ~c * ~f * ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (~b * ~e + 2 * ~a * ~f)) * (2 * ~p + 2 * ~q + 5) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~q], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (LtQ(~p, -1) && (NeQ((~c * ~d - ~a * ~f) ^ 2 - (~b * ~d - ~a * ~e) * (~c * ~e - ~b * ~f), 0) && (Not(Not(IntegerQ(~p)) && ILtQ(~q, -1)) && Not(IGtQ(~q, 0)))))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(f') * (~x) ^ 2) ^ ~q, ~x) => (((~b) ^ 3 * ~f + ~b * ~c * (~c * ~d - 3 * ~a * ~f) + ~c * ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (2 * ~a * ~f)) * ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~f * (~x) ^ 2) ^ (~q + 1)) / (((~b) ^ 2 - 4 * ~a * ~c) * ((~b) ^ 2 * ~d * ~f + (~c * ~d - ~a * ~f) ^ 2) * (~p + 1)) - (1 / (((~b) ^ 2 - 4 * ~a * ~c) * ((~b) ^ 2 * ~d * ~f + (~c * ~d - ~a * ~f) ^ 2) * (~p + 1))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~f * (~x) ^ 2) ^ ~q * Simp((2 * ~c * ((~b) ^ 2 * ~d * ~f + (~c * ~d - ~a * ~f) ^ 2) * (~p + 1) - ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (2 * ~a * ~f)) * (~a * ~f * (~p + 1) - ~c * ~d * (~p + 2))) + (2 * ~f * ((~b) ^ 3 * ~f + ~b * ~c * (~c * ~d - 3 * ~a * ~f)) * (~p + ~q + 2) - ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (2 * ~a * ~f)) * (~b * ~f * (~p + 1))) * ~x + ~c * ~f * ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (2 * ~a * ~f)) * (2 * ~p + 2 * ~q + 5) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~q], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (LtQ(~p, -1) && (NeQ((~b) ^ 2 * ~d * ~f + (~c * ~d - ~a * ~f) ^ 2, 0) && (Not(Not(IntegerQ(~p)) && ILtQ(~q, -1)) && Not(IGtQ(~q, 0))))))
    @apply_utils Antiderivative((~(a') + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q, ~x) => ((2 * ~a * (~c) ^ 2 * ~e + ~c * (2 * (~c) ^ 2 * ~d - ~c * (2 * ~a * ~f)) * ~x) * (~a + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~q + 1)) / ((-4 * ~a * ~c) * (~a * ~c * (~e) ^ 2 + (~c * ~d - ~a * ~f) ^ 2) * (~p + 1)) - (1 / ((-4 * ~a * ~c) * (~a * ~c * (~e) ^ 2 + (~c * ~d - ~a * ~f) ^ 2) * (~p + 1))) * Antiderivative((~a + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q * Simp(((2 * ~c * ((~c * ~d - ~a * ~f) ^ 2 - (-(~a) * ~e) * (~c * ~e)) * (~p + 1) - (2 * (~c) ^ 2 * ~d - ~c * (2 * ~a * ~f)) * (~a * ~f * (~p + 1) - ~c * ~d * (~p + 2))) - ~e * (-2 * ~a * (~c) ^ 2 * ~e) * (~p + ~q + 2)) + (2 * ~f * (2 * ~a * (~c) ^ 2 * ~e) * (~p + ~q + 2) - (2 * (~c) ^ 2 * ~d - ~c * (2 * ~a * ~f)) * (-(~c) * ~e * (2 * ~p + ~q + 4))) * ~x + ~c * ~f * (2 * (~c) ^ 2 * ~d - ~c * (2 * ~a * ~f)) * (2 * ~p + 2 * ~q + 5) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~q], ~x) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (LtQ(~p, -1) && (NeQ(~a * ~c * (~e) ^ 2 + (~c * ~d - ~a * ~f) ^ 2, 0) && (Not(Not(IntegerQ(~p)) && ILtQ(~q, -1)) && Not(IGtQ(~q, 0))))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q, ~x) => (((~b * ~f * (3 * ~p + 2 * ~q) - ~c * ~e * (2 * ~p + ~q)) + 2 * ~c * ~f * (~p + ~q) * ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p - 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~q + 1)) / (2 * (~f) ^ 2 * (~p + ~q) * (2 * ~p + 2 * ~q + 1)) - (1 / (2 * (~f) ^ 2 * (~p + ~q) * (2 * ~p + 2 * ~q + 1))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p - 2) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q * Simp(((~b * ~d - ~a * ~e) * (~c * ~e - ~b * ~f) * (1 - ~p) * (2 * ~p + ~q) - (~p + ~q) * ((~b) ^ 2 * ~d * ~f * (1 - ~p) - ~a * (~f * (~b * ~e - 2 * ~a * ~f) * (2 * ~p + 2 * ~q + 1) + ~c * (2 * ~d * ~f - (~e) ^ 2 * (2 * ~p + ~q))))) + (2 * (~c * ~d - ~a * ~f) * (~c * ~e - ~b * ~f) * (1 - ~p) * (2 * ~p + ~q) - (~p + ~q) * (((~b) ^ 2 - 4 * ~a * ~c) * ~e * ~f * (1 - ~p) + ~b * (~c * ((~e) ^ 2 - 4 * ~d * ~f) * (2 * ~p + ~q) + ~f * ((2 * ~c * ~d - ~b * ~e) + 2 * ~a * ~f) * (2 * ~p + 2 * ~q + 1)))) * ~x + ((~c * ~e - ~b * ~f) ^ 2 * (1 - ~p) * ~p + ~c * (~p + ~q) * (~f * (~b * ~e - 2 * ~a * ~f) * ((4 * ~p + 2 * ~q) - 1) - ~c * (2 * ~d * ~f * (1 - 2 * ~p) + (~e) ^ 2 * ((3 * ~p + ~q) - 1)))) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~q], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (GtQ(~p, 1) && (NeQ(~p + ~q, 0) && (NeQ(2 * ~p + 2 * ~q + 1, 0) && (Not(IGtQ(~p, 0)) && Not(IGtQ(~q, 0))))))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(f') * (~x) ^ 2) ^ ~q, ~x) => ((~b * (3 * ~p + 2 * ~q) + 2 * ~c * (~p + ~q) * ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p - 1) * (~d + ~f * (~x) ^ 2) ^ (~q + 1)) / (2 * ~f * (~p + ~q) * (2 * ~p + 2 * ~q + 1)) - (1 / (2 * ~f * (~p + ~q) * (2 * ~p + 2 * ~q + 1))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p - 2) * (~d + ~f * (~x) ^ 2) ^ ~q * Simp((((~b) ^ 2 * ~d * (~p - 1) * (2 * ~p + ~q) - (~p + ~q) * ((~b) ^ 2 * ~d * (1 - ~p) - 2 * ~a * (~c * ~d - ~a * ~f * (2 * ~p + 2 * ~q + 1)))) - (2 * ~b * (~c * ~d - ~a * ~f) * (1 - ~p) * (2 * ~p + ~q) - 2 * (~p + ~q) * ~b * (2 * ~c * ~d * (2 * ~p + ~q) - (~c * ~d + ~a * ~f) * (2 * ~p + 2 * ~q + 1))) * ~x) + ((~b) ^ 2 * ~f * ~p * (1 - ~p) + 2 * ~c * (~p + ~q) * (~c * ~d * (2 * ~p - 1) - ~a * ~f * ((4 * ~p + 2 * ~q) - 1))) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~q], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (GtQ(~p, 1) && (NeQ(~p + ~q, 0) && (NeQ(2 * ~p + 2 * ~q + 1, 0) && (Not(IGtQ(~p, 0)) && Not(IGtQ(~q, 0)))))))
    @apply_utils Antiderivative((~(a') + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q, ~x) => (-(~c) * (~e * (2 * ~p + ~q) - 2 * ~f * (~p + ~q) * ~x) * (~a + ~c * (~x) ^ 2) ^ (~p - 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~q + 1)) / (2 * (~f) ^ 2 * (~p + ~q) * (2 * ~p + 2 * ~q + 1)) - (1 / (2 * (~f) ^ 2 * (~p + ~q) * (2 * ~p + 2 * ~q + 1))) * Antiderivative((~a + ~c * (~x) ^ 2) ^ (~p - 2) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q * Simp(-(~a) * ~c * (~e) ^ 2 * (1 - ~p) * (2 * ~p + ~q) + ~a * (~p + ~q) * (-2 * ~a * (~f) ^ 2 * (2 * ~p + 2 * ~q + 1) + ~c * (2 * ~d * ~f - (~e) ^ 2 * (2 * ~p + ~q))) + (2 * (~c * ~d - ~a * ~f) * (~c * ~e) * (1 - ~p) * (2 * ~p + ~q) + 4 * ~a * ~c * ~e * ~f * (1 - ~p) * (~p + ~q)) * ~x + (~p * (~c) ^ 2 * (~e) ^ 2 * (1 - ~p) - ~c * (~p + ~q) * (2 * ~a * (~f) ^ 2 * ((4 * ~p + 2 * ~q) - 1) + ~c * (2 * ~d * ~f * (1 - 2 * ~p) + (~e) ^ 2 * ((3 * ~p + ~q) - 1)))) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~q], ~x) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (GtQ(~p, 1) && (NeQ(~p + ~q, 0) && (NeQ(2 * ~p + 2 * ~q + 1, 0) && (Not(IGtQ(~p, 0)) && Not(IGtQ(~q, 0)))))))
    @apply_utils Antiderivative(1 / ((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([q = (((((~c) ^ 2 * (~d) ^ 2 - ~b * ~c * ~d * ~e) + ~a * ~c * (~e) ^ 2 + (~b) ^ 2 * ~d * ~f) - 2 * ~a * ~c * ~d * ~f) - ~a * ~b * ~e * ~f) + (~a) ^ 2 * (~f) ^ 2], (1 / q) * Antiderivative((((((~c) ^ 2 * ~d - ~b * ~c * ~e) + (~b) ^ 2 * ~f) - ~a * ~c * ~f) - ((~c) ^ 2 * ~e - ~b * ~c * ~f) * ~x) / (~a + ~b * ~x + ~c * (~x) ^ 2), ~x) + (1 / q) * Antiderivative((((~c * (~e) ^ 2 - ~c * ~d * ~f) - ~b * ~e * ~f) + ~a * (~f) ^ 2 + (~c * ~e * ~f - ~b * (~f) ^ 2) * ~x) / (~d + ~e * ~x + ~f * (~x) ^ 2), ~x) <-- NeQ(q, 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NeQ((~e) ^ 2 - 4 * ~d * ~f, 0))
    @apply_utils Antiderivative(1 / ((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * (~d + ~(f') * (~x) ^ 2)), ~x) => With([q = (((~c) ^ 2 * (~d) ^ 2 + (~b) ^ 2 * ~d * ~f) - 2 * ~a * ~c * ~d * ~f) + (~a) ^ 2 * (~f) ^ 2], (1 / q) * Antiderivative(((((~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~a * ~c * ~f) + ~b * ~c * ~f * ~x) / (~a + ~b * ~x + ~c * (~x) ^ 2), ~x) - (1 / q) * Antiderivative(((~c * ~d * ~f - ~a * (~f) ^ 2) + ~b * (~f) ^ 2 * ~x) / (~d + ~f * (~x) ^ 2), ~x) <-- NeQ(q, 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~f], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative(1 / ((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => -2 * ~e * Subst(Antiderivative(1 / (~e * (~b * ~e - 4 * ~a * ~f) - (~b * ~d - ~a * ~e) * (~x) ^ 2), ~x), ~x, (~e + 2 * ~f * ~x) / sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && EqQ(~c * ~e - ~b * ~f, 0)))
    #= Antiderivative(1/(((~a)+(~b')*(~x)+(~c')*(~x)^2)*sqrt((~d')+(~e')*(~x)+(~f')*(~x)^2)),~x) :=   With([k=Rt((a/c-d/f)^2+(b/c-e/f)*(b*d/(c*f)-a*e/(c*f)),2)], -2*(c*d-a*f+c*f*k+(c*e-b*f)*x)*sqrt((d+e*x+f*x^2)*((c*f*k)/(c*d-a*f+ c*f*k+(c*e-b*f)*x))^2)/(c*sqrt(d+e*x+f*x^2))* Subst(Antiderivative((1-x)/( (b*d-a*e-b*f*k-(c*d-a*f-c*f*k)^2/(c*e-b*f)+(b*d-a*e+b*f*k-(a*f- c*d-c*f*k)^2/(c*e-b*f))*x^2)* sqrt(-f*((b*d-a*e-c*e*k)/(c*e-b*f)-(c*d-a*f-c*f*k)^2/(c*e-b*f)^ 2)-f*((b*d-a*e+c*e*k)/(c*e-b*f)-(a*f-c*d-c*f*k)^2/(c*e-b*f)^2)*x^2)), x),x, (c*d-a*f-c*f*k+(c*e-b*f)*x)/(c*d-a*f+c*f*k+(c*e-b*f)*x))) <--  FreeQ([a,b,c,d,e,f],x) && RationalQ(a,b,c,d,e,f) && NeQ(b^2-4*a*c,0) &&  NeQ(e^2-4*d*f,0) && NeQ(c*e-b*f,0) && LtQ(b^2-4*a*c,0) =#
    #= Antiderivative(1/(((~a)+(~c')*(~x)^2)*sqrt((~d')+(~e')*(~x)+(~f')*(~x)^2)),~x) :=  With([k=Rt((a/c-d/f)^2+a*e^2/(c*f^2),2)], -2*(c*d-a*f+c*f*k+c*e*x)*sqrt((d+e*x+f*x^2)*((c*f*k)/(c*d-a*f+c*f*k+ c*e*x))^2)/(c*sqrt(d+e*x+f*x^2))* Subst(Antiderivative((1-x)/( (-a*e-(c*d-a*f-c*f*k)^2/(c*e)+(-a*e-(a*f-c*d-c*f*k)^2/(c*e))*x^ 2)* sqrt(-f*((-a*e-c*e*k)/(c*e)-(c*d-a*f-c*f*k)^2/(c*e)^2)-f*((-a*e+ c*e*k)/(c*e)-(a*f-c*d-c*f*k)^2/(c*e)^2)*x^2)),x),x, (c*d-a*f-c*f*k+(c*e)*x)/(c*d-a*f+c*f*k+(c*e)*x))) <--  FreeQ([a,c,d,e,f],x) && RationalQ(a,c,d,e,f) && NeQ(e^2-4*d*f,0) &&  LtQ(-a*c,0) =#
    #= Antiderivative(1/(((~a)+(~b')*(~x)+(~c')*(~x)^2)*sqrt((~d')+(~f')*(~x)^2)),~x) :=  With([k=Rt((a/c-d/f)^2+b^2*d/(c^2*f),2)], -2*(c*d-a*f+c*f*k-b*f*x)*sqrt((d+f*x^2)*((c*f*k)/(c*d-a*f+c*f*k-b*f* x))^2)/(c*sqrt(d+f*x^2))* Subst(Antiderivative((1-x)/( (b*d-b*f*k+(c*d-a*f-c*f*k)^2/(b*f)+(b*d+b*f*k+(a*f-c*d-c*f*k)^2/ (b*f))*x^2)* sqrt(-f*(-d/f-(c*d-a*f-c*f*k)^2/(b*f)^2)-f*(-d/f-(a*f-c*d-c*f*k) ^2/(b*f)^2)*x^2)),x),x, (c*d-a*f-c*f*k+(-b*f)*x)/(c*d-a*f+c*f*k+(-b*f)*x))) <--  FreeQ([a,b,c,d,f],x) && RationalQ(a,b,c,d,f) && NeQ(b^2-4*a*c,0) &&  LtQ(b^2-4*a*c,0) =#
    @apply_utils Antiderivative(1 / ((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c) / q) * Antiderivative(1 / (((~b - q) + 2 * ~c * ~x) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) - ((2 * ~c) / q) * Antiderivative(1 / ((~b + q + 2 * ~c * ~x) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (NeQ(~c * ~e - ~b * ~f, 0) && PosQ((~b) ^ 2 - 4 * ~a * ~c))))
    @apply_utils Antiderivative(1 / ((~a + ~(c') * (~x) ^ 2) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => (1 / 2) * Antiderivative(1 / ((~a - Rt(-(~a) * ~c, 2) * ~x) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) + (1 / 2) * Antiderivative(1 / ((~a + Rt(-(~a) * ~c, 2) * ~x) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f], ~x) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && PosQ(-(~a) * ~c))
    @apply_utils Antiderivative(1 / ((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~d + ~(f') * (~x) ^ 2)), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c) / q) * Antiderivative(1 / (((~b - q) + 2 * ~c * ~x) * sqrt(~d + ~f * (~x) ^ 2)), ~x) - ((2 * ~c) / q) * Antiderivative(1 / ((~b + q + 2 * ~c * ~x) * sqrt(~d + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~f], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && PosQ((~b) ^ 2 - 4 * ~a * ~c))
    @apply_utils Antiderivative(1 / ((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([q = Rt((~c * ~d - ~a * ~f) ^ 2 - (~b * ~d - ~a * ~e) * (~c * ~e - ~b * ~f), 2)], (1 / (2q)) * Antiderivative(((~c * ~d - ~a * ~f) + q + (~c * ~e - ~b * ~f) * ~x) / ((~a + ~b * ~x + ~c * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) - (1 / (2q)) * Antiderivative((((~c * ~d - ~a * ~f) - q) + (~c * ~e - ~b * ~f) * ~x) / ((~a + ~b * ~x + ~c * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (NeQ(~c * ~e - ~b * ~f, 0) && NegQ((~b) ^ 2 - 4 * ~a * ~c))))
    @apply_utils Antiderivative(1 / ((~(a') + ~(c') * (~x) ^ 2) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([q = Rt((~c * ~d - ~a * ~f) ^ 2 + ~a * ~c * (~e) ^ 2, 2)], (1 / (2q)) * Antiderivative(((~c * ~d - ~a * ~f) + q + ~c * ~e * ~x) / ((~a + ~c * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) - (1 / (2q)) * Antiderivative((((~c * ~d - ~a * ~f) - q) + ~c * ~e * ~x) / ((~a + ~c * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~c, ~d, ~e, ~f], ~x) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && NegQ(-(~a) * ~c))
    @apply_utils Antiderivative(1 / ((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~(d') + ~(f') * (~x) ^ 2)), ~x) => With([q = Rt((~c * ~d - ~a * ~f) ^ 2 + (~b) ^ 2 * ~d * ~f, 2)], (1 / (2q)) * Antiderivative(((~c * ~d - ~a * ~f) + q + (-(~b) * ~f) * ~x) / ((~a + ~b * ~x + ~c * (~x) ^ 2) * sqrt(~d + ~f * (~x) ^ 2)), ~x) - (1 / (2q)) * Antiderivative((((~c * ~d - ~a * ~f) - q) + (-(~b) * ~f) * ~x) / ((~a + ~b * ~x + ~c * (~x) ^ 2) * sqrt(~d + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~f], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NegQ((~b) ^ 2 - 4 * ~a * ~c))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * ~x + ~(c') * (~x) ^ 2) / (~d + ~(e') * ~x + ~(f') * (~x) ^ 2), ~x) => (~c / ~f) * Antiderivative(1 / sqrt(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) - (1 / ~f) * Antiderivative(((~c * ~d - ~a * ~f) + (~c * ~e - ~b * ~f) * ~x) / (sqrt(~a + ~b * ~x + ~c * (~x) ^ 2) * (~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NeQ((~e) ^ 2 - 4 * ~d * ~f, 0))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * ~x + ~(c') * (~x) ^ 2) / (~d + ~(f') * (~x) ^ 2), ~x) => (~c / ~f) * Antiderivative(1 / sqrt(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) - (1 / ~f) * Antiderivative(((~c * ~d - ~a * ~f) - ~b * ~f * ~x) / (sqrt(~a + ~b * ~x + ~c * (~x) ^ 2) * (~d + ~f * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative(sqrt(~a + ~(c') * (~x) ^ 2) / (~d + ~(e') * ~x + ~(f') * (~x) ^ 2), ~x) => (~c / ~f) * Antiderivative(1 / sqrt(~a + ~c * (~x) ^ 2), ~x) - (1 / ~f) * Antiderivative(((~c * ~d - ~a * ~f) + ~c * ~e * ~x) / (sqrt(~a + ~c * (~x) ^ 2) * (~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f], ~x) && NeQ((~e) ^ 2 - 4 * ~d * ~f, 0)
    @apply_utils Antiderivative(1 / (sqrt(~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~d + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([r = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((sqrt(~b + r + 2 * ~c * ~x) * sqrt(2 * ~a + (~b + r) * ~x)) / sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)) * Antiderivative(1 / (sqrt(~b + r + 2 * ~c * ~x) * sqrt(2 * ~a + (~b + r) * ~x) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NeQ((~e) ^ 2 - 4 * ~d * ~f, 0))
    @apply_utils Antiderivative(1 / (sqrt(~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~d + ~(f') * (~x) ^ 2)), ~x) => With([r = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((sqrt(~b + r + 2 * ~c * ~x) * sqrt(2 * ~a + (~b + r) * ~x)) / sqrt(~a + ~b * ~x + ~c * (~x) ^ 2)) * Antiderivative(1 / (sqrt(~b + r + 2 * ~c * ~x) * sqrt(2 * ~a + (~b + r) * ~x) * sqrt(~d + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~f], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q, ~x) => Unintegrable((~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~p, ~q], ~x) && (Not(IGtQ(~p, 0)) && Not(IGtQ(~q, 0)))
    @apply_utils Antiderivative((~a + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q, ~x) => Unintegrable((~a + ~c * (~x) ^ 2) ^ ~p * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~p, ~q], ~x) && (Not(IGtQ(~p, 0)) && Not(IGtQ(~q, 0)))
    @apply_utils Antiderivative((~(a') + ~(b') * ~u + ~(c') * (~u) ^ 2) ^ ~(p') * (~(d') + ~(e') * ~u + ~(f') * (~u) ^ 2) ^ ~(q'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~p, ~q], ~x) && (LinearQ(~u, ~x) && NeQ(~u, ~x))
    @apply_utils Antiderivative((~(a') + ~(c') * (~u) ^ 2) ^ ~(p') * (~(d') + ~(e') * ~u + ~(f') * (~u) ^ 2) ^ ~(q'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~a + ~c * (~x) ^ 2) ^ ~p * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q, ~x), ~x, ~u) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~p, ~q], ~x) && (LinearQ(~u, ~x) && NeQ(~u, ~x))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.1 Quadratic/1.2.1.5 (a+b x+c x^2)^p (d+e x+f x^2)^q.jl


text:
"""
prime_check(number)

Checks to see if a number is a prime or not

A number is prime if it has exactly two factors: 1 and itself.

# Example

```julia
prime_check(2) # returns true
prime_check(3) # returns true
prime_check(5) # returns true
prime_check(7) # returns true
prime_check(11) # returns true
prime_check(13) # returns true
prime_check(17) # returns true
prime_check(19) # returns true
prime_check(23) # returns true
prime_check(29) # returns true
prime_check(30) # returns false
```

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee) and [Rratic](https://github.com/Rratic)
"""
function prime_check(number::N) where {N<:Integer}
    if 1 < number < 4
        # 2, 3 are primes
        return true
    elseif number < 2 || number % 2 == 0
        # Negative ,0,1, and even numbers except 2 are not primes
        return false
    end
    for i in 3:2:(isqrt(number)+1)
        if number % i == 0
            return false
        end
    end
    return true
end

ext: jl

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
include("SppGenerator.jl")
#
# function profFunction(var::Int64, c::Vector{Int64}, M::Array{Int64, 2})
#     return c[var]/sum(M[:,var])
# end

# function findVariable(c::Vector{Int64}, M::Array{Int64, 2}, nbVar::Int64)
#     varMax = 1
#     maxProf = profFunction(varMax, c, M)
#     for var = 2:nbVar
#         prof = profFunction(var, c, M)
#         if prof > maxProf
#             varMax = var
#             maxProf = prof
#         end
#     end
#     varMax = maxProf == 0 ? -1 : varMax
#     return varMax
# end

# function max(dict::Dict{Any, Int64})
#     valMax = -Inf
#     varMax = -1
#     for pair in dict
#         if pair.second > valMax
#             varMax = pair.first
#             valMax = pair.second
#         end
#     end
#     return (indMax, valMax)
# end

function profFunction(var::VI, c::Dict{VI, Int64}, occur::Dict{VI, Int64})
    return (haskey(c, var) && haskey(occur, var) ? c[var]/occur[var] : -Inf)
end

function findVariable(c::Dict{VI, Int64}, occurVar::Dict{VI, Int64}, variables::Vector{VI})
    varMax = VI(-1)
    valMax = typemin(Float64)
    for var in variables

        val = profFunction(var, c, occurVar)
        if val > valMax
            varMax = var
            valMax = val
        end
    end
    return varMax
end

function z(x::Dict{VI, Int64}, c::Dict{VI, Int64})
    val::Int64 = 0
    for (x_i, b) in x
        val += b*c[x_i]
    end
    return val
end

function heuristic(c::Dict{VI, Int64}, M::Dict{CI, Vector{VI}}, listOccur::Dict{VI, Dict{VI, Int64}}, listLiaison::Dict{VI, Vector{VI}}, occurVar::Dict{VI, Int64})
    nbVar = length(c)
    nbConst = length(M)
    x0::Dict{VI, Int64} = Dict(var => 0 for var=VI.(1:nbVar))
    # relationList::Vector{Vector{Int64}} = Vector{Vector{Int64}}(undef, nbrVar)
    variables = VI.(1:nbVar)
    while variables != []

        x_i = findVariable(c, occurVar, variables)
        filterList::Vector{VI} = haskey(listLiaison, x_i) ? listLiaison[x_i] : [x_i]
        filter!(x -> !(x in filterList), variables)
        x0[x_i] = 1
        # Reduire OccurVar

        for variable in variables
            if haskey(listOccur, variable)
                for pair in listOccur[variable]
                    # println(occurVar[variable])
                    # println(listOccur[variable])
                    # println(listOccur[variable][supprVar])
                    if pair.first in filterList
                        occurVar[variable] -= pair.second
                    end
                end
            end
        end

    end
    x = zeros(nbVar)
    for pair in x0
        x[pair.first.value] = pair.second
    end
    return x, z(x0, c)
end

function heuristic(nbVar::Int64, nbConst::Int64)
    c, M, listOccur, listLiaison, occurVar = generate(nbVar, nbConst)
    return heuristic(c, M, listOccur, listLiaison, occurVar)
end

function heuristic(fname::String)
    c, M, listOccur, listLiaison, occurVar = loadSPP(fname)
    return heuristic(c, M, listOccur, listLiaison, occurVar)
end

    # Faux ! Ecrit ton heuristique sur papier avant de la coder !
function main(fname::String)
    heuristic("Data/test.txt")
    return @time heuristic(fname)
end


text:
module FISH_MDS

    using Ipopt
    using ArgParse
    using Grid

    # exports...

    export run_mds
    export mds_main

    # things to run w/o command line?
    export load_data
    export new_mds
    export make_ipopt_problem
    export solveProblem
    export remove_infs
    export get_removed_indices
    export output_txt
    export interp_3d

    include("graph.jl")
    include("mds_metric.jl")
    include("mds_approximate_fish.jl")
    include("output.jl")
    include("interpolate.jl")
    include("mds_main.jl")

end # module


text:
import gmsh

gmsh.initialize(ARGS)
gmsh.option.setNumber("General.Terminal", 1)

gmsh.model.add("test");

# add discrete surface with tag 1
gmsh.model.addDiscreteEntity(2, 1)

# add 4 mesh nodes
gmsh.model.mesh.addNodes(2, 1,
                         [1, 2, 3, 4], # node tags: 1, 2, 3, and 4
                         [0., 0., 0., # coordinates of node 1
                          1., 0., 0., # coordinates of node 2
                          1., 1., 0., # ...
                          0., 1., 0.])

# add 2 triangles
gmsh.model.mesh.addElements(2, 1,
                            [2], # single type : 3-node triangle
                            [[1, 2]], # triangle tags: 1 and 2
                            [[1, 2, 3, # triangle 1: nodes 1, 2, 3
                              1, 3, 4]]) # triangle 2: nodes 1, 3, 4

# export the mesh ; use explore.py to read and examine the mesh
gmsh.write("test.msh")

gmsh.finalize()

max_stars_repo_name: Christophe-Foyer/PyFEA

max_issues_repo_name: Christophe-Foyer/pyFEA

max_forks_repo_name: Christophe-Foyer/pyFEA


text:
function velocityphi(r,x,b,mode,d)

#   evaluates the velocity from the solution
#
#   input
#   r        - radial grid (excludes boundaries)
#   x        - cos(theta) grid (excludes poles and equator)
#   b        - radial magnetic field (mT)
#   eta      - dimensionless diffusivity
#   mode     - solution vector (eigenfunction)
#   d        - eigenvalue  (frequency = d / 1im)

#   some constants
    omega = 0.7292e-4;       # rotation rate (1/s)
    sectoyear = 3.155e7;     # seconds in a year
    mu = 4 * pi * 1.0e-7;    # permeability
    eta = 0.8;            # magnetic diffusivity  (m^2/s)
    rho = 1.0e4;             # density
    L = 3.480e6;             # core radius (length scale)
    bscale = sqrt(rho * mu) * omega * L;

#   dimensionless quantities
    B = b * 1.0e-3 / bscale;
    Eta = eta / (omega * L^2);

#  allocate memory for velocity field
   nr = length(r);
   nx = length(x);
   ndim = nr * nx;
   dvdr = Array{ComplexF64}(undef,ndim);
   vphi = Array{ComplexF64}(undef,ndim);

#  radial magnetic field (dimensionless)
   br = magnetic_field(B,x);

#  retrieve required matrices
   dsp = diffusion(r,x,eta);

#  rate of change of bphi
   dbdt = d * mode[1:ndim];

#  account for diffusion
   dbdr2 = dsp[ndim+1:2*ndim,ndim+1:2*ndim] * mode[1:ndim];

#  combine to evaluate velocity gradient
   f = dbdt - dbdr2;
   for l = 1 : nx
       for k = 1 : nr
           ipos = (l-1) * nr + k;
           dvdr[ipos] = f[ipos] / br[l];
       end
   end

#  integrate for vphi
   dr = r[2] - r[1];
   for l = 1 : nx
       for k = 2 : nr
           ipos = (l-1) * nr + k;
           vphi[ipos] = vphi[ipos-1]+ 0.5*(dvdr[ipos] + dvdr[ipos-1])*dr;
       end
   end

#  extract CMB value
   vphicmb = Array{ComplexF64}(undef,nx)
   bphicmb = Array{ComplexF64}(undef,nx);

   for l = 1 : nx
       k = nr;
       ipos = (l-1)*nr + k;
#       println(" l ", l, " ipos ",ipos," vphi ",vphi[ipos]," ",vphi[ipos-1])
       vphicmb[l] = vphi[ipos];
       bphicmb[l] = mode[ipos];
#       println(" l ", vphicmb[l])
   end

   return vphicmb,bphicmb

end

max_stars_repo_name: bbuffett/MacWaves

max_issues_repo_name: bbuffett/MacWaves

max_forks_repo_name: bbuffett/MacWaves


text:
# coding: utf-8

# In[1]:

using Random
using Plots

# In[2]:

#bogo, what the heck does it even mean
#i d call it bullocks sort
#it is not really sorting
#it keeps shuffling the list until the list is sorted
#kinda feel like monte carlo simulation
#the time complexity is (n+1)!
#for a large list, we may never wait for the day it gets sorted

# In[3]:

function bogo_sort(arr)

    #deep copy is crucial
    arr_copy=deepcopy(arr)
    sort!(arr_copy)
    counter=0

    while arr_copy!=arr

        shuffle!(arr)
        counter+=1

    end

    return counter

end

# In[4]:

test_arr=rand(0:100,5)

#dont forget to deep copy
test_stats=[bogo_sort(deepcopy(test_arr)) for i in 1:100000]

# In[5]:

#visualize the ridiculous result of bogo sort by histogram
#it can go as far as 1400 times of shuffle to get it sorted from my test
#this is merely a five-element list

# In[6]:

#just out of curiosity
#could this be a power law distribution
histogram(test_stats,bins=10000,
    normed=:true,fill=(:blue, true),
    la=0,legend=:none,
    xlabel="iterations to reach sorted state",
    ylabel="frequency in 100k simulations")

max_stars_repo_path: bogo sort.jl

max_stars_repo_name: je-suis-tm/search-and-sort

max_issues_repo_path: bogo sort.jl

max_issues_repo_name: HafizSabia/search-and-sort

max_forks_repo_path: bogo sort.jl

max_forks_repo_name: HafizSabia/search-and-sort


text:
module Fetch
    using Random: randstring
    using HTTP
    using JSON3
    using StructTypes
    using Base64

    include("gdrive.jl")
    include("kaggle.jl")
end


text:
is_root(node) = isnull(node.parent)

is_leaf(node) = length(node.children) == 0

all_children_explored(node::Node{ACTIONC}) where ACTIONC =
    length(node.children)==ACTIONC

all_children_done(tree, node) =
    all_children_explored(node) && all(tree.dones[c.id] for c in values(node.children))

function get_rank(node)
    rank=1
    while true
        isnull(node.parent) && break
        node = get(node.parent)
        rank += 1
    end
    rank
end

get_best_leaf(tree, gamma) = tree.nodes[argmax(tree.accumulated_rewards + gamma*[maximum(qs) for qs in tree.children_qs])]

get_tree_rank(tree, gamma) = get_rank(get_best_leaf(tree, gamma))

get_accumulated_reward(tree) = maximum(tree.accumulated_rewards)

function calc_visitedc(tree, maxval=-1)
    # calculate number of times each node was visited
    # iterate in reverse, so parents have all children
    # already updated when we get to them
    visitedc = zeros(Int32, length(tree.nodes))
    for nix in length(tree.nodes):-1:1
        node = tree.nodes[nix]
        for ch in values(node.children)
            visitedc[nix] += visitedc[ch.id]
            if maxval > 0 && visitedc[nix] > maxval
                visitedc[nix] = maxval
            end
        end
        visitedc[nix] += 1
    end
    visitedc
end

function calc_action_visitedc(tree, visited_threshold=-1, nonexplored_value=0)
    visitedc = calc_visitedc(tree, visited_threshold)
    avc = Vector{Int32}[]
    for node in tree.nodes
        _avc = Int32(nonexplored_value)*ones(Int32, length(tree.children_qs[1]))
        for (ix,ch) in node.children
            _avc[ix] = visitedc[ch.id]
        end
        push!(avc, _avc)
    end
    avc
end

function backprop_weighted_q!(tree, gamma, visited_threshold=-1)
    avisitedc = calc_action_visitedc(tree, visited_threshold)
    # calculate weighted Qs, root is done separately
    # iterae in reverse, so parents have all children
    # already updated when we get to them
    for nix in length(tree.nodes):-1:2
        node = tree.nodes[nix]
        tree.dones[node.id] = all_children_done(tree, node)
        # update Q at parent
        parent = get(node.parent)
        node_reward = tree.accumulated_rewards[node.id] - tree.accumulated_rewards[parent.id]
        mn = if is_leaf(node)
            mean(tree.children_qs[nix])
        else
            sum(avisitedc[nix] .* tree.children_qs[nix]) ./ sum(avisitedc[nix])
        end
        tree.children_qs[parent.id][node.action_ix] = node_reward + gamma*mn
    end
    # update root
    tree.dones[1] = all_children_done(tree, tree.nodes[1])
end

function backprop_max_q!(tree, gamma)
    # backprop everybody in reverse, so parents have all children
    # already updated when we get to them
    for nix in length(tree.nodes):-1:2
        node = tree.nodes[nix]
        parent = get(node.parent)
        tree.dones[node.id] = all_children_done(tree, node)
        # update Q at parent
        mx::Float32 = maximum(tree.children_qs[node.id])
        node_reward = tree.accumulated_rewards[node.id] - tree.accumulated_rewards[parent.id]
        tree.children_qs[parent.id][node.action_ix] = node_reward + gamma*mx
    end
    # update root
    tree.dones[1] = all_children_done(tree, tree.nodes[1])
end

function expand!(heap, tree, parent_node::Node{ACTIONC}, w, action_ix, env, gamma) where ACTIONC
    actual_action = action_ix_to_action(env, action_ix)
    new_sim_state, reward, done = sim!(env, tree.states[parent_node.id], actual_action)
    new_sensors = get_sensors(env, new_sim_state)
    children_qs = if done
        zeros(Float32, ACTIONC)
    else
        network_predict(env, w,new_sensors)
    end
    id::Int32 = length(tree.states) + 1
    new_node = Node(Int32(action_ix), Nullable(parent_node), id)
    parent_node.children[action_ix] = new_node
    accumulated_reward::Float32 = tree.accumulated_rewards[parent_node.id] + Float32(reward)
    push!(tree.nodes, new_node)
    push!(tree.sensors, new_sensors)
    push!(tree.children_qs, children_qs)
    push!(tree.states, new_sim_state)
    push!(tree.accumulated_rewards, accumulated_reward)
    push!(tree.dones, done)
    if !done
        for (aix, q) in enumerate(children_qs)
            score::Float32 = accumulated_reward + Float32(gamma)*q
            push!(heap, Int32(aix), new_node.id, score)
        end
    end
end

function build_tree(w, env, root_state, stepc, epsilon, gamma)
    # build the root
    root_sensors = get_sensors(env, root_state)
    root_children_qs = network_predict(env, w, root_sensors)
    ACTIONC = length(root_children_qs)
    null_parent = Nullable{Node{ACTIONC}}()
    root_action_ix::Int32 = -1 # nonsensical action
    root_id::Int32 = 1
    root_accumulated_reward::Float32 = 0.0f0
    root_done::Bool = false
    root = Node(root_action_ix, null_parent, root_id)

    # build the tree
    nodes = Node{ACTIONC}[]
    sensors = typeof(root_sensors)[]
    children_qs = Vector{Float32}[]
    states = typeof(root_state)[]
    accumulated_rewards = Float32[]
    dones = Bool[]
    heap = Heap()
    tree = Tree(env, heap, nodes, sensors, children_qs, states, accumulated_rewards, dones)

    # update the tree with root
    push!(tree.nodes, root)
    push!(tree.sensors, root_sensors)
    push!(tree.children_qs, root_children_qs)
    push!(tree.states, root_state)
    push!(tree.accumulated_rewards, root_accumulated_reward)
    push!(tree.dones, root_done)
    for (aix, q) in enumerate(tree.children_qs[1])
        score::Float32 = Float32(gamma)*q
        push!(tree.heap, Int32(aix), root.id, score)
    end

    # add new nodes in a loop
    for i in 1:stepc
        if i % 200 == 0
            print(".")
            flush(stdout)
        end

        # choose parent and action
        tree.dones[1] && return tree # root is done
        length(heap) == 0 && return tree
        action_ix, parent_id = if rand()>epsilon
            pop_max!(heap)
        else
            pop_rand!(heap)
        end
        @assert !tree.dones[parent_id]
        parent = tree.nodes[parent_id]
        @assert !haskey(parent.children, action_ix)

        # add a new node to the tree
        expand!(heap, tree, parent, w, action_ix, env, gamma)
    end
    println()
    return tree
end


text:
export PetscVec, VecSetType, VecSetValues, VecAssemblyBegin, VecAssemblyEnd, VecSetSizes, VecGetSize, VecGetLocalSize, VecNorm, VecGetValues, VecGetOwnershipRange, VecGetArray, VecRestoreArray, VecGetArrayRead, VecRestoreArrayRead, VecSet, VecSqrtAbs, VecLog, VecExp, VecAbs, VecMax, VecMin, VecCopy, VecDuplicate, VecAXPY, VecAXPBY, VecAYPX, VecWAXPY, VecMAXPY, VecAXPBYPCZ, VecScale, VecDot, VecTDot, VecSum, VecSwap, VecReciprocal, VecShift, VecPointwiseMult, VecPointwiseDivide

export getLocalIndices

"""
Petsc Vector type.

Not a subtype of `AbstractArray` because Petsc vectors do not conform to that
API.
"""
type PetscVec
  pobj::Ptr{Void}

  """
    Constructor

    **Inputs**

     * comm: MPI Communicator
  """
  function PetscVec(comm::MPI_Comm)
#    comm = PETSC_COMM_SELF();
    vec = Array{Ptr{Void}}(1)
    err = ccall(( :VecCreate, libpetsclocation ), PetscErrorCode,(comm_type,Ptr{Void}),comm,vec);
    vec = new(vec[1])
#    finalizer(vec,PetscDestroy)
    # does not seem to be called immediately when vec is no longer visible, is it called later during garbage collection? - yes
    return vec

  end

  """
    Constructor from existing pointer to a Petsc Vec

    **Inputs**

     pobj: a Ptr{Void} to an existing Petsc object
  """
  function PetscVec(pobj::Ptr{Void})  # default constructor
    return new(pobj)
  end

end

"""
Union{AbstractVector, PetscVec}
"""
const AllVectors = Union{AbstractVector, PetscVec}

"""
  Free a Petsc vec.  Safe to call multiple times
"""
function PetscDestroy(vec::PetscVec)
  if (vec.pobj != 0)
    err = ccall(( :VecDestroy, libpetsclocation), PetscErrorCode, (Ptr{Ptr{Void}},), &vec.pobj);
  end
  vec.pobj = 0  # unnecessary? vec no longer has any references to it
#    println("VecDestroy called")
end

"""
  VecSetType
"""
function VecSetType(vec::PetscVec,name)
  err = ccall((:VecSetType,  libpetsclocation), PetscErrorCode, (Ptr{Void}, Cstring), vec.pobj,name);
end

  #TODO: make this better: use VecFromArray
function PetscVec(array::Array{PetscScalar})
  vec = PetscVec()
  err = ccall(( :VecSetType,  libpetsclocation), PetscErrorCode,(Ptr{Void},Cstring), vec.pobj,"seq");
  err = ccall( (:VecSetSizes,  libpetsclocation), PetscErrorCode, (Ptr{Void}, PetscInt, PetscInt), vec.pobj,length(array),length(array));
  # want a PetscInt array so build it ourselves
  idx = Array{PetscInt}(length(array));
  for i=1:length(array);  idx[i] = i-1;  end
  err = ccall( ( :VecSetValues,  libpetsclocation), PetscErrorCode,(Ptr{Void},PetscInt, Ptr{PetscInt},Ptr{PetscScalar},Int32), vec.pobj,length(idx),idx,array,INSERT_VALUES);
  err = ccall( ( :VecAssemblyBegin,  libpetsclocation), PetscErrorCode,(Ptr{Void},), vec.pobj);
  err = ccall( ( :VecAssemblyEnd,  libpetsclocation), PetscErrorCode, (Ptr{Void},), vec.pobj);
  return vec
end

"""
   VecSetValues
"""
function VecSetValues(vec::PetscVec,idx::Array{PetscInt},array::Array{PetscScalar},flag::Integer)

  err = ccall( ( :VecSetValues,  libpetsclocation), PetscErrorCode,
               (Ptr{Void},PetscInt ,Ptr{PetscInt},Ptr{PetscScalar},Int32),
               vec.pobj,length(idx), idx,array,flag);
  return err
end

"""
  VecSetValues method that implicitly uses `INSERT_VALUES`
"""
function VecSetValues(vec::PetscVec,idx::Array{PetscInt},array::Array{PetscScalar})
  VecSetValues(vec,idx,array,INSERT_VALUES)
end

#=
function VecSetValues(vec::PetscVec,array::Array{PetscScalar})
  idx = Array{PetscInt}(ength(array))
  for i=1:length(array);  idx[i] = i-1;  end
  VecSetValues(vec,idx,array,INSERT_VALUES)
end
=#

"""
  VecAssemblyBegin
"""
function VecAssemblyBegin(obj::PetscVec)
  err = ccall( ( :VecAssemblyBegin,  libpetsclocation), PetscErrorCode, (Ptr{Void},), obj.pobj);
end

"""
  VecAssemblyEnd
"""
function VecAssemblyEnd(obj::PetscVec)
  err = ccall( ( :VecAssemblyEnd,  libpetsclocation), PetscErrorCode,(Ptr{Void},), obj.pobj);
end

"""
  Convenience function for calling VecAssemblyBegin, and VecAssemblyEnd, in
  one go.
"""
function VecAssemble(obj::PetscVec)
  VecAssemblyBegin(obj)
  VecAssemblyEnd(obj)
end

"""
  VecSetSizes
"""
function VecSetSizes(vec::PetscVec,n::Integer, N::Integer)
  err = ccall( ( :VecSetSizes,  libpetsclocation), PetscErrorCode, (Ptr{Void}, PetscInt, PetscInt), vec.pobj,n,N);
end

"""
  PetscView for Petsc vector
"""
function PetscView(obj::PetscVec,viewer)
  err = ccall( ( :VecView,  libpetsclocation), PetscErrorCode, (Ptr{Void},Int64),obj.pobj,0);
end

"""
  VeGetSize

  **Inputs**

   * vec: a Petsc vector

  **Outputs**

   * the size
"""
function VecGetSize(obj::PetscVec)
  n = Array{PetscInt}(1)
  err = ccall( ( :VecGetSize,  libpetsclocation), PetscErrorCode, (Ptr{Void},Ptr{PetscInt}), obj.pobj,n);
  return n[1]
end

"""
  VecGetLocalSize

  **Inputs**

   * vec: a Petsc vector

  **Outputs**

   * the size

"""
function VecGetLocalSize(arg1::PetscVec)
  arg2 = Ref{PetscInt}()
  ccall((:VecGetLocalSize,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscInt}),arg1.pobj,arg2)
  return arg2[]
end

"""
  VecNorm

  **Inputs**

   * obj: Petsc vector
   * normtype: the Petsc enum for the norm type

  **Output**

   * the norm value
"""
function VecNorm(obj::PetscVec,normtype::Integer)
  n = Array{PetscReal}()
  err = ccall( ( :VecNorm,  libpetsclocation), PetscScalar, (Ptr{Void},Int32,Ptr{PetscReal}), obj.pobj,normtype, n);
  return n[1]
end

#=
function VecNorm(obj::PetscVec)
  return VecNorm(obj,PETSC_NORM_2)
end
=#
"""
  VecGetValues
"""
function VecGetValues(vec::PetscVec, ni::Integer, ix::AbstractArray{PetscInt,1}, y::AbstractArray{PetscScalar,1})

     # need indices to be PetscInt
#     ix_local = Array{PetscInt}(ni)
#     for i=1:ni
#       ix_local[i] = ix[i]
#     end

    err = ccall((:VecGetValues,petsc),PetscErrorCode,(Ptr{Void},PetscInt,Ptr{PetscInt},Ptr{PetscScalar}),vec.pobj, ni, ix, y)

    return nothing
end

"""
  VecGetValues with length of idx inferred from idx

  **Inputs**

   * vec: the Petsc vector
   * idx: array of PetscInt indices
   * y: array of PetscScalar values
"""
function VecGetValues(vec::PetscVec,idx::AbstractArray{PetscInt,1}, y::AbstractArray{PetscScalar, 1})

  VecGetValues(vec, length(idx), idx, y)
end

"""
  VecGetOwnershipRange

  **Inputs**

   * vec: Petsc vector

  **Outputs**

   * low: lowest index (zero-based) that is owned
   * high: highest index + 1 that is owned
"""
function VecGetOwnershipRange(vec::PetscVec)
    low = Array{PetscInt}(1)
    high = Array{PetscInt}(1)

    ccall((:VecGetOwnershipRange,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscInt},Ptr{PetscInt}),vec.pobj, low, high)

  return low[1], high[1]
end

# new functions
"""
  VecGetArray.  Users must call `VecRestoreArray` when finished.

  **Inputs**

   * vec: the Petsc vector

  **Outputs**

   * arr: a Julia Array{PetscScalar, 1}
"""
function VecGetArray(vec::PetscVec)
# gets a pointer to the data underlying a Petsc vec, turns it into a Julia
# array
# ptr_arr must be passed into PetscVecRestoreArray
    ptr_arr = Array{Ptr{PetscScalar}}(1)
    ccall((:VecGetArray,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Ptr{PetscScalar}}),vec.pobj, ptr_arr)

    first, last = VecGetOwnershipRange(vec)
    len = last - first
    arr = unsafe_wrap(Array, ptr_arr[1], len)
    return arr
end

"""
  VecRestoreArray.  Users must not access the array after calling this function.

  **Inputs**

   * vec: the PetscVector passed into `VecGetArray`
   * arr: the array returned by `VecGetArray
"""
function VecRestoreArray(vec::PetscVec, arr::Array{PetscScalar, 1})
    ptr_arr = Ref{Ptr{PetscScalar}}(pointer(arr))
    ccall((:VecRestoreArray,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Ptr{PetscScalar}}),vec.pobj, ptr_arr)
end

"""
  Similar to `VecGetArray`, but produces the array must not be written to.
"""
function VecGetArrayRead(vec::PetscVec)
    ptr_arr = Array{Ptr{PetscScalar}}(1)
    ccall((:VecGetArrayRead,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Ptr{PetscScalar}}),vec.pobj, ptr_arr)

    first, last = VecGetOwnershipRange(vec)
    len = last - first
    arr = unsafe_wrap(Array, ptr_arr[1], len)
    return arr
end

"""
  Similar to `VecRestoreArray`, but corresponds to `VecGetArrayRead`
"""
function VecRestoreArrayRead(vec::PetscVec, arr::Array{PetscScalar, 1})

    ptr_arr = Ref{Ptr{PetscScalar}}(pointer(arr))
    ccall((:VecRestoreArrayRead,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Ptr{PetscScalar}}),vec.pobj, ptr_arr)
#    ccall((:VecRestoreArrayRead,petsc),PetscErrorCode,(Vec,Ptr{Ptr{PetscScalar}}),arg1,arg2)
end

"""
  VecSet
"""
function VecSet(vec::PetscVec, val::PetscScalar)
    ccall((:VecSet,petsc),PetscErrorCode,(Ptr{Void},PetscScalar), vec.pobj, val)
end

"""
  VecSqrtAbs
"""
function VecSqrtAbs(vec::PetscVec)
    ccall((:VecSqrtAbs,petsc),PetscErrorCode,(Ptr{Void},), vec.pobj)
end

"""
  VecLog
"""
function VecLog(vec::PetscVec)
    ccall((:VecLog,petsc),PetscErrorCode,(Ptr{Void},),vec.pobj)
end

"""
  VecExp
"""
function VecExp(vec::PetscVec)
    ccall((:VecExp,petsc),PetscErrorCode,(Ptr{Void},),vec.pobj)
end

function VecAbs(vec::PetscVec)
    ccall((:VecAbs,petsc),PetscErrorCode,(Ptr{Void},),vec.pobj)
end

"""
  VecMax

  **Inputs**

   * vec: PetscVec

  **Output**

   * r: the maximum value
   * idx: the (zero-based) index of the maximum value
"""
function VecMax(vec::PetscVec)
    r = Array{PetscReal}(1) # max value
    idx = Array{PetscInt}(1)  # index of max value
    ccall((:VecMax,petsc),PetscErrorCode,(Ptr{Void}, Ptr{PetscInt},Ptr{PetscReal}), vec.pobj, idx, r)

    return r[1], idx[1]
end

"""
  VecMin. Same interface as `VecMax`
"""
function VecMin(vec::PetscVec)
    r = Array{PetscReal}(1) # min value
    idx = Array{PetscInt}(1)  # index of min value

    ccall((:VecMin,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscInt},Ptr{PetscReal}), vec.pobj, idx, r)

    return r[1], idx[1]
end

"""
  VecReciprocal
"""
function VecReciprocal(vec::PetscVec)
    ccall((:VecReciprocal,petsc),PetscErrorCode,(Ptr{Void},),vec.pobj)
end

"""
  VecShift
"""
function VecShift(vec::PetscVec, a::PetscScalar)
    ccall((:VecShift,petsc),PetscErrorCode,(Ptr{Void},PetscScalar), vec.pobj, a)
end

"""
  VecPointwiseMult
"""
function VecPointwiseMult(w::PetscVec, x::PetscVec,y::PetscVec)
    ccall((:VecPointwiseMult,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Void}, Ptr{Void}), w.pobj, x.pobj, y.pobj)
end

"""
  VecPointwiseDivide
"""
function VecPointwiseDivide(w::PetscVec, x::PetscVec, y::PetscVec)
    ccall((:VecPointwiseDivide,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Void}, Ptr{Void}), w.pobj, x.pobj, y.pobj)
end

"""
  VecCopy
"""
function VecCopy(vec::PetscVec , vec2::PetscVec)
    ccall((:VecCopy,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Void}), vec.pobj, vec2.pobj)
end

"""
  VecDulicate
"""
function VecDuplicate( vec::PetscVec)
    ptr_arr = Array{Ptr{Void}}(1)

    ccall((:VecDuplicate,petsc),PetscErrorCode,( Ptr{Void}, Ptr{Ptr{Void}}), vec.pobj, ptr_arr)

    return PetscVec(ptr_arr[1])
end

# Some vector linear algebra
"""
  VecAXPY
"""
function VecAXPY( vec1::PetscVec, a::PetscScalar, vec2::PetscVec)
    ccall((:VecAXPY,petsc),PetscErrorCode,(Ptr{Void},PetscScalar,Ptr{Void}), vec1.pobj, a, vec2.pobj)
end

"""
  VecAXPBY
"""
function VecAXPBY(vec1::PetscVec, a::PetscScalar, b::PetscScalar, vec2::PetscVec)
    ccall((:VecAXPBY,petsc),PetscErrorCode,(Ptr{Void},PetscScalar,PetscScalar,Ptr{Void}),vec1.pobj, a, b, vec2.pobj)
end

"""
  VecMAXPY
"""
function VecMAXPY(vec1::PetscVec, n::Integer, a::AbstractArray{PetscScalar, 1}, x::AbstractArray{Ptr{Void}, 1})
# the vector x must contains the pointers from the PetscVec objects, not the PetscVec objects themselves

    ccall((:VecMAXPY,petsc),PetscErrorCode,(Ptr{Void},PetscInt,Ptr{PetscScalar},Ptr{Ptr{Void}}),vec1.pobj, n, a, x)
end

"""
  VecAYPX
"""
function VecAYPX(vec1::PetscVec, a::PetscScalar, vec2::PetscVec)
    ccall((:VecAYPX,petsc),PetscErrorCode,(Ptr{Void},PetscScalar,Ptr{Void}),vec1.pobj, a, vec2.pobj)
end

"""
  VecWAXPY
"""
function VecWAXPY(w::PetscVec, a::PetscScalar, x::PetscVec, y::PetscVec)
    ccall((:VecWAXPY,petsc),PetscErrorCode,(Ptr{Void},PetscScalar,Ptr{Void},Ptr{Void}), w.pobj, a, x.pobj, y.pobj)
end

"""
  VecAXPBYPCZ
"""
function VecAXPBYPCZ(z::PetscVec, alpha::PetscScalar, beta::PetscScalar, gamma::PetscScalar, x::PetscVec, y::PetscVec)
    ccall((:VecAXPBYPCZ,petsc),PetscErrorCode,(Ptr{Void},PetscScalar,PetscScalar,PetscScalar,Ptr{Void},Ptr{Void}), z.pobj, alpha,beta, gamma, x.pobj, y.pobj)
end

"""
  VecScale
"""
function VecScale(vec::PetscVec, a::PetscScalar)
    ccall((:VecScale,petsc),PetscErrorCode,(Ptr{Void},PetscScalar), vec.pobj, a)
end

"""
  VecDot
"""
function VecDot(x::PetscVec, y::PetscVec)
    r = Array{PetscScalar}(1)
    err = ccall((:VecDot,petsc),PetscErrorCode,( Ptr{Void}, Ptr{Void}, Ptr{PetscScalar}), x.pobj, y.pobj, r)

    println("PetscVecDot error code = ", err)
    return r[1]
end

"""
  VecTDot
"""
function VecTDot(x::PetscVec, y::PetscVec)
    r = Array{PetscScalar}(1)
    ccall((:VecTDot,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Void},Ptr{PetscScalar}), x.pobj, y.pobj, r)

    return r[1]
end

"""
  VecSum
"""
function VecSum(vec::PetscVec)
    r = Array{PetscScalar}(1)
    ccall((:VecSum,petsc),PetscErrorCode,(Ptr{Void},Ptr{PetscScalar}),vec.pobj, r)
    return r[1]
end

"""
  VecSwap
"""
function VecSwap(x::PetscVec, y::PetscVec)
    ccall((:VecSwap,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Void}), x.pobj, y.pobj)

end

include("vec_interface.jl")


text:
# generated
type = "champion"
format = "standAloneComplex"
version = v"11.17.1"
for locale in ("en_US", "ko_KR")
    gendir = normpath(@__DIR__, "11.17.1", "generated", locale)
    include(normpath(gendir, "module.jl"))
end


text:
using LinearAlgebra

function _projecthermitian!(a::StridedMatrix)
    L = size(a,1)
    for J = 1:16:L
        for I = 1:16:L
            if I == J
                for j = J:min(J+15,L)
                    for i = J:j-1
                        x = (a[i,j] + conj(a[j,i]))/2
                        a[i,j] = x
                        a[j,i] = conj(x)
                    end
                    a[j,j] = (a[j,j] + conj(a[j,j]))/2
                end
            else
                for j = J:min(J+15,L)
                    for i = I:min(I+15,L)
                        x = (a[i,j] + conj(a[j,i]))/2
                        a[i,j] = x
                        a[j,i] = conj(x)
                    end
                end
            end
        end
    end
    return a
end
function _projecthermitian!(a::AbstractMatrix)
    a .= (a .+ a')./2
    return a
end

function _projectantihermitian!(a::StridedMatrix)
    L = size(a,1)
    for J = 1:16:L
        for I = 1:16:L
            if I == J
                for j = J:min(J+15,L)
                    for i = J:j-1
                        x = (a[i,j] - conj(a[j,i]))/2
                        a[i,j] = x
                        a[j,i] = -conj(x)
                    end
                    a[j,j] = (a[j,j] - conj(a[j,j]))/2
                end
            else
                for j = J:min(J+15,L)
                    for i = I:min(I+15,L)
                        x = (a[i,j] - conj(a[j,i]))/2
                        a[i,j] = x
                        a[j,i] = -conj(x)
                    end
                end
            end
        end
    end
    return a
end
function _projectantihermitian!(a::AbstractMatrix)
    a .= (a .- a')./2
    return a
end

function _addone!(a::AbstractMatrix)
    view(a, diagind(a)) .= view(a, diagind(a)) .+ 1
    return a
end
function _subtractone!(a::AbstractMatrix)
    view(a, diagind(a)) .= view(a, diagind(a)) .- 1
    return a
end
function _one!(A::AbstractMatrix)
    m, n = size(A)
    @inbounds for j = 1:n
        for i = 1:m
            A[i,j] = i==j
        end
    end
    return A
end

function _polarsdd!(A::StridedMatrix)
    U, S, V = svd!(A; alg = LinearAlgebra.DivideAndConquer())
    return mul!(A, U, V')
end
function _polarsvd!(A::StridedMatrix)
    U, S, V = svd!(A; alg = LinearAlgebra.QRIteration())
    return mul!(A, U, V')
end
function _polarnewton!(A::StridedMatrix; tol = 10*eleps(A), maxiter = 5)
    m, n = size(A)
    @assert m >= n
    A2 = copy(A)
    Q, R = qr!(A2)
    Ri = ldiv!(UpperTriangular(R)', _one!(similar(R)))
    R, Ri = _avgdiff!(R, Ri)
    i = 1
    R2 = view(A, 1:n, 1:n)
    fill!(view(A, n+1:m, 1:n), zero(eltype(A)))
    copyto!(R2, R)
    while maximum(abs, Ri) > tol
        if i == maxiter # if not converged by now, fall back to sdd
            _polarsdd!(Ri)
            break
        end
        Ri = ldiv!(lu!(R2)', _one!(Ri))
        R, Ri = _avgdiff!(R, Ri)
        copyto!(R2, R)
        i += 1
    end
    return lmul!(Q, A)
end
# in place computation of the average and difference of two arrays
function _avgdiff!(A::AbstractArray, B::AbstractArray)
    axes(A) == axes(B) || throw(DimensionMismatch())
    @simd for I in eachindex(A, B)
        @inbounds begin
            a = A[I]
            b = B[I]
            A[I] = (a+b)/2
            B[I] = b-a
        end
    end
    return A, B
end

# _stiefelexp(W, A, Z, α)
# given an isometry W, and a Stiefel tangent vector Δ = W*A + Z, compute the building blocks
# W′, Q′, R′ of the geodesic with respect to the canonical metric in the direction of α*Δ.
# Here, W′is the new isometry, and the local tangent vector is given by
# Δ′ = W′ * A + Z′ with Z′ = Q′*R′
# Here, Q′ is a set of orthogonal columns to the colums in W′.
function _stiefelexp(W::StridedMatrix, A::StridedMatrix, Z::StridedMatrix, α)
    n, p = size(W)
    r = min(2*p, n)
    QQ, _ = qr!([W Z])
    Q = similar(W, n, r-p)
    @inbounds for j = Base.OneTo(r-p)
        for i = Base.OneTo(n)
            Q[i,j] = (i == p+j)
        end
    end
    Q = lmul!(QQ, Q)
    R = Q'*Z
    A2 = similar(A, min(2*p, n), min(2*p, n))
    A2[1:p, 1:p] .= α .* A
    A2[p+1:end, 1:p] .= α .* R
    A2[1:p, p+1:end] .= (-α) .* (R')
    A2[p+1:end, p+1:end] .= 0
    U = [W Q]*exp(A2)
    U = _polarnewton!(U)
    W′ = U[:,1:p]
    Q′ = U[:,p+1:end]
    R′ = R
    return W′, Q, Q′, R′
end

function _stiefellog(Wold::StridedMatrix, Wnew::StridedMatrix;
                        tol = 10*eleps(Wold), maxiter = 100)
    n, p = size(Wold)
    r = min(2*p, n)
    P = Wold'*Wnew
    dW = Wnew - Wold*P
    QQ, _ = qr!([Wold dW])
    Q = similar(Wold, n, r-p)
    @inbounds for j = Base.OneTo(r-p)
        for i = Base.OneTo(n)
            Q[i,j] = (i == p+j)
        end
    end
    Q = lmul!(QQ, Q)
    R = Q'*dW
    Wext = [Wold Q];
    F = qr!([P; R])
    U = lmul!(F.Q, _one!(similar(P, r, r)))
    U[1:p, 1:p] .= P
    U[p+1:r, 1:p] .= R
    X = view(U, 1:p, p+1:r)
    Y = view(U, p+1:r, p+1:r)
    if p < n
        YSVD = svd!(Y)
        mul!(X, X*(YSVD.V), (YSVD.U)')
        UsqrtS = YSVD.U
        @inbounds for j = 1:size(UsqrtS, 2)
            s = sqrt(YSVD.S[j])
            @simd for i = 1:size(UsqrtS, 1)
                UsqrtS[i,j] *= s
            end
        end
        mul!(Y, UsqrtS, UsqrtS')
    end
    logU = _projectantihermitian!(log(U))
    if eltype(U) <: Real
        @assert mapreduce(abs ∘ imag, max, logU; init = abs(zero(eltype(logU))) ) <= tol
        K = real(logU)
    else
        K = logU
    end
    C = view(K, p+1:r, p+1:r)
    i = 1
    τ = mapreduce(abs, max, C; init = abs(zero(eltype(C))))
    while τ > tol
        if i > maxiter
            @warn "Stiefel logarithm: not converged in $maxiter iterations, τ = $τ"
            break
        end
        eC = exp(rmul!(C,-1))
        X .= X*eC
        Y .= Y*eC
        logU = _projectantihermitian!(log(U))
        if eltype(U) <: Real
            @assert mapreduce(abs ∘ imag, max, logU; init = abs(zero(eltype(logU))) ) <= tol
            K .= real.(logU)
        else
            K .= logU
        end
        τ = maximum(abs, C)
        i += 1
    end
    return K[1:p,1:p], Q, K[p+1:r, 1:p]
end


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

"""
Optional private token for additional security. Use of token prompts seed rehashing to
obstruct client-side interpretability.
"""
struct Token
    data::Any
end
Base.show(io::IO, ::Token) = print(io, "Token()")

"""
Salt random seed with private token for additional security.
"""
reseed(seed::UInt32, token::Token) = UInt32(hash((seed, token.data)) % typemax(UInt32))

"""
    ASTServer

Simulation-side server. Interacts with client via TCP.
"""
Base.@kwdef mutable struct ASTServer
    ip::IPAddr                              = ip"::1"  # address(es) to listen on
    port::Int64                             = 1812     # port to listen on
    serv::Union{Sockets.TCPServer, Nothing} = nothing
    mdp::ASTMDP
    token::Union{Token, Nothing}            = nothing
    tunnel::Bool                            = false
    verbose::Bool                           = false
    presample::Bool                         = true
end

"""
    ASTServer(mdp::ASTMDP; kwargs...)

Constructor for server.
"""
function ASTServer(mdp::ASTMDP; kwargs...)
    ASTServer(; mdp=mdp, kwargs...)
end

"""
    set_token(server::ASTServer, token)

Set private token to encrypt seeds. Recommended to use `set_password` instead, to avoid
leaking token through side effects.
"""
function set_token(server::ASTServer, token)
    if action_type(server.mdp) == SeedAction
        server.token = Token(token)
        @info "Private token set."
    else
        @info "Private token can only be set for seed-action (blackbox) simulators."
    end
    return
end

"""
    set_password(server::ASTServer)

Set private token to string input (password). Avoids explicit representation of password in
storage or side effects.
"""
function set_password(server::ASTServer)
    io = Base.getpass("Enter password")
    Base.shred!(io) do io
        set_token(server, read(io, String))
    end
    return
end

"""
Respond to echo request from ASTClient.
"""
_echo(::ASTMDP, x) = x

"""
Respond to information request from ASTClient.
"""
function _info(mdp::ASTMDP{S, A}) where {S, A}
    format = t -> Symbol(last(split("$t", ".")))
    Dict(
        :S => format(S),
        :A => format(A),
        :episodic => mdp.episodic,
        :num_steps => mdp.num_steps
    )
end

"""
Process simulation request from client and construct response.
"""
function respond(server::ASTServer, request::Dict)
    # interpret request
    sym = FMAP[request[:f]]
    f = getproperty(Interface, sym)
    args = request[:a]
    kwargs = request[:k]

    # reseeding
    if server.token !== nothing && sym == :act! && action_type(server.mdp) == SeedAction
        seed = reseed(args[1], server.token)
        args = (seed,)
    end

    if haskey(kwargs, :batch) && kwargs.batch
        r = sum(f.(Ref(server.mdp), args...)) # batch processing
    else
        r = f(server.mdp, args...)
    end

    # server-side sampling
    if server.presample && sym == :actions && action_type(server.mdp) == SampleAction
        r = rand(r; flat=true)
    end

    return Dict(:r => r)
end

"""
Listen for incoming requests and execute function calls on MDP. Can handle multiple client
connections asynchronously, but server currently holds a single MDP instance.
"""
function run(server::ASTServer)
    @async while true
        # listen for incoming connections
        conn = accept(server.serv)
        Sockets.nagle(conn, false)
        @info "Connected to AST client." conn
        @async while true
            request = BSON.load(conn)
            server.verbose && @info "Received request from client:" request
            response = try
                respond(server, request)
            catch e
                bson(conn, Dict(:e => true))
                throw(e)
            end
            bson(conn, response)
            server.verbose && @info "Sent response to client." response
        end
    end
end

"""
    connect!(server::ASTServer; remote::String="", remote_port::Int64=1812)

Connect server, optionally through SSH tunnel. Optional keyword argument `remote` should be
of the form `user@machine`.
"""
function connect!(server::ASTServer; remote::String="", remote_port::Int64=1812)
    disconnect!(server)
    !isempty(remote) && open_tunnel(server, remote, remote_port)
    server.serv = listen(server.ip, server.port)
    run(server)
    return
end

"""
    disconnect!(server::ASTServer)

Disconnect server.
"""
function disconnect!(server::ASTServer)
    if server.serv !== nothing
        close(server.serv)
        server.serv = nothing
    end
    server.tunnel && close_tunnel()
    server.tunnel = false
    return
end


text:
function Invariants.errormessage(io::IO, invs::AllInvariant, msgs)
    __combinator_errormessage(io, invs, msgs, map(__getmarker, msgs),
        faint("All invariants listed below should be satisfied:\n\n "))
end

function InvariantsCore.errormessage(io::IO, invs::AnyInvariant, msgs)
    __combinator_errormessage(io, invs, msgs, map(m -> isnothing(m) ? PASS : FAIL, msgs),
        faint("At least one of the invariants listed below should be satisfied:\n\n"))
end

const PASS = "\e[32m✔\e[0m\e[2m"
const FAIL = "\e[31m⨯\e[0m"
const UNKNOWN = "\e[33m?\e[0m"

function __combinator_errormessage(io::IO, invs, msgs, markers, msg)
    showdescription(io, invs)
    print(io, msg)
    for (inv, message, marker) in zip(invs.invariants, msgs, markers)
        __errormessage_child(io, inv; marker, message)
    end
end

function __errormessage_child(
        io,
        inv;
        marker = "o",
        message=nothing,
        indent = "    "
        )
    print(io, marker, " ")
    showtitle(io, inv)
    print(io, "\e[0m")
    println(io)

    if !(isnothing(message) || ismissing(message))
        ioi = WrapIO(io; indent)
        errormessage(ioi, inv, message)
    end
end

__getmarker(::Nothing) = PASS
__getmarker(::Missing) = UNKNOWN
__getmarker(_) = FAIL


text:
"""
    conductance(h::Hypergraph, subset::Set{Int})::Float64

Calculate unweighted hypergraph conductance of `subset`.
note: ∅ ⊊ `subset` ⊊ `1:nhv(h)`

For more information see `1. Introduction` at:
Generalizing the Hypergraph Laplacian via a Diffusion Processwith Mediators,
auhtors: T-H. Hubert Chan, Xhibin Liang.
"""
function conductance(h::Hypergraph, subset::Set{Int})::Float64
  if isempty(subset) error("`subset` is not allowed empty.") end
  if subset == Set(1:nhv(h)) error("`subset` is not allowed true subset of `h`.") end
  subset2 = setdiff(Set(1:nhv(h)), subset)
  wₛ = sum([length(gethyperedges(h, node)) for node in subset])
  w∂s = 0
  for he in 1:nhe(h)
    he_vertices = keys(getvertices(h, he))
    if isempty(intersect(he_vertices, subset)) continue end
    if isempty(intersect(he_vertices, subset2)) continue end
    w∂s += length(getvertices(h, he))
  end
  return w∂s / wₛ
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, PosDefManifoldML

makedocs(
   sitename="PosDefManifoldML",
   authors="Saloni Jain, Marco Congedo",
   modules=[PosDefManifoldML],
   pages =
   [
      "index.md",
      "Tutorial" => "tutorial.md",
      "Main Module" => "MainModule.md",
      "Tools" => "tools.md",

      "ML Models: PD Manifold" => Any[
                  "Minimum Distance to Mean" => "mdm.md",
      ],
      "ML Models: PD Tangent Space" => Any[
                  "Elastic-Net Logistic Regression" => "enlr.md",
      ],
      "fit, predict, cv" => "cv.md",

   ]
)

deploydocs(
   # root
   # target = "build", # add this folder to .gitignore!
   repo = "github.com/Marco-Congedo/PosDefManifoldML.jl.git",
   # branch = "gh-pages",
   # osname = "linux",
   # deps = Deps.pip("pygments", "mkdocs"),
   # devbranch = "dev",
   # devurl = "dev",
   # versions = ["stable" => "v^", "v#.#", devurl => devurl],
)


text:
push!(LOAD_PATH, "./")
using Rando
#=
feed a Matrix with random data
=#
function RandomFeedMatrix(M,R::RandoData)
    n,m=size(M)
    @inbounds for i = 1:n
        @inbounds for j = 1:m
            M[i,j]=Rando.fv!(R)
        end
    end
end

#= factor a matrix od size (n,n+1)
   size is not checked
=#
function factorMatrix!(M::Array{Float64, 2})
    n, m = size(M)
    for line = 1:n
        # find pivot
         @inbounds(cmax= line-1+indmax(abs.(M[line:n,line])))
        # exchange rows if necessary
        if cmax != line
            @simd for j = line:m
                @inbounds M[line, j], M[cmax, j] = M[cmax, j], M[line, j]
            end
        end

        pivot = M[line, line]
        for k = line+1:m
            @inbounds tmp = M[line, k] / pivot
            @simd for j = line+1:n
                @inbounds M[j, k] -= M[j, line] * tmp
            end
        end

        M[line+1:n, line] .= 0.
    end

end

function doall(n::Int64)
    Ro = RandoData()
    M = Array{Float64}(n, n+1)
    RandomFeedMatrix(M, Ro)

    t1 = time_ns()
    factorMatrix!(M)
    t2 = time_ns() - t1

end

# main starts here.

MatrixSize = 2048
liminf = 2
D = Dict{Integer, Real}()

# be sure to run once before actually running the benchmark!
T = doall(MatrixSize)
#

println("start")

#Profile.clear()

while MatrixSize > liminf
    T = doall(MatrixSize)
    println(MatrixSize)
    D[MatrixSize] = 1.e-9 * T
    MatrixSize = convert(Int64, MatrixSize/2)
end

# Results:
S = sort(collect(zip(keys(D), values(D))), rev=true)

open("RunningOn"*gethostname(),"w") do f
    for k in S
        write(f,string(k[1])," ",string(k[2]),"\n")
    end
end

#Profile.print()

max_stars_repo_name: wolfv/BenchmarksPythonJuliaAndCo

max_issues_repo_name: wolfv/BenchmarksPythonJuliaAndCo

max_forks_repo_name: wolfv/BenchmarksPythonJuliaAndCo


text:
# Dense GPU Array
abstract type GPUArray{T, N} <: DenseArray{T, N} end

# Sampler type that acts like a texture/image and allows interpolated access
abstract type Sampler{T, N} <: DenseArray{T, N} end

const GPUVector{T} = GPUArray{T, 1}
const GPUMatrix{T} = GPUArray{T, 2}
const GPUVecOrMat{T} = Union{GPUArray{T, 1}, GPUArray{T, 2}}

# GPU Local Memory
struct LocalMemory{T} <: GPUArray{T, 1}
    size::Int
    LocalMemory{T}(x::Integer) where T = new{T}(x)
end

# input/output

## serialization

import Serialization: AbstractSerializer, serialize, deserialize, serialize_type

function serialize(s::AbstractSerializer, t::T) where T <: GPUArray
    serialize_type(s, T)
    serialize(s, Array(t))
end
function deserialize(s::AbstractSerializer, ::Type{T}) where T <: GPUArray
    A = deserialize(s)
    T(A)
end

function to_cartesian(A, indices::Tuple)
    start = CartesianIndex(ntuple(length(indices)) do i
        val = indices[i]
        isa(val, Integer) && return val
        isa(val, UnitRange) && return first(val)
        isa(val, Colon) && return 1
        isa(val, Base.Slice{Base.OneTo{Int}}) && return 1
        error("GPU indexing only defined for integers or unit ranges. Found: $val")
    end)
    stop = CartesianIndex(ntuple(length(indices)) do i
        val = indices[i]
        isa(val, Integer) && return val
        isa(val, UnitRange) && return last(val)
        isa(val, Colon) && return size(A, i)
        isa(val, Base.Slice{Base.OneTo{Int}}) && return size(A, i)
        error("GPU indexing only defined for integers or unit ranges. Found: $val")
    end)
    CartesianIndices(start, stop)
end

## convert to CPU (keeping wrapper type)

Adapt.adapt_storage(::Type{<:Array}, xs::AbstractArray) = convert(Array, xs)
convert_to_cpu(xs) = adapt(Array, xs)

## showing

for (W, ctor) in (:AT => (A,mut)->mut(A), Adapt.wrappers...)
    @eval begin
        # display
        Base.print_array(io::IO, X::$W where {AT <: GPUArray}) =
            Base.print_array(io, $ctor(X, convert_to_cpu))

        # show
        Base._show_nonempty(io::IO, X::$W where {AT <: GPUArray}, prefix::String) =
            Base._show_nonempty(io, $ctor(X, convert_to_cpu), prefix)
        Base._show_empty(io::IO, X::$W where {AT <: GPUArray}) =
            Base._show_empty(io, $ctor(X, convert_to_cpu))
        Base.show_vector(io::IO, v::$W where {AT <: GPUArray}, args...) =
            Base.show_vector(io, $ctor(v, convert_to_cpu), args...)
    end
end

## collect to CPU (discarding wrapper type)

collect_to_cpu(xs::AbstractArray) = collect(convert_to_cpu(xs))

for (W, ctor) in (:AT => (A,mut)->mut(A), Adapt.wrappers...)
    @eval begin
        Base.collect(X::$W where {AT <: GPUArray}) = collect_to_cpu(X)
    end
end

# memory copying

## basic linear copies of identically-typed memory

materialize(x::AbstractArray) = Array(x)
materialize(x::GPUArray) = x

for (D, S) in ((GPUArray, AbstractArray), (Array, GPUArray), (GPUArray, GPUArray))
    @eval begin
        function Base.copyto!(dest::$D{T, N}, rdest::NTuple{N, UnitRange},
                              src::$S{T, N}, ssrc::NTuple{N, UnitRange}) where {T, N}
            drange = CartesianIndices(rdest)
            srange = CartesianIndices(ssrc)
            copyto!(dest, drange, src, srange)
        end

        function Base.copyto!(dest::$D{T}, d_range::CartesianIndices{1},
                              src::$S{T}, s_range::CartesianIndices{1}) where T
            amount = length(d_range)
            if length(s_range) != amount
                throw(ArgumentError("Copy range needs same length. Found: dest: $amount, src: $(length(s_range))"))
            end
            amount == 0 && return dest
            d_offset = first(d_range)[1]
            s_offset = first(s_range)[1]
            copyto!(dest, d_offset, materialize(src), s_offset, amount)
        end

        function Base.copyto!(dest::$D{T, N}, src::$S{T, N}) where {T, N}
            len = length(src)
            len == 0 && return dest
            copyto!(dest, 1, materialize(src), 1, len)
        end
    end
end

## generalized blocks of heterogeneous memory

Base.copyto!(dest::GPUArray, src::GPUArray) =
    copyto!(dest, CartesianIndices(dest), src, CartesianIndices(src))

function copy_kernel!(state, dest, dest_offsets, src, src_offsets, shape, shape_dest, shape_source, length)
    i = linear_index(state)
    if i <= length
        # TODO can this be done faster and smarter?
        idx = gpu_ind2sub(shape, i)
        dest_idx = gpu_sub2ind(shape_dest, idx .+ dest_offsets)
        src_idx = gpu_sub2ind(shape_source, idx .+ src_offsets)
        @inbounds dest[dest_idx] = src[src_idx]
    end
    return
end

function Base.copyto!(dest::GPUArray{T, N}, destcrange::CartesianIndices{N},
                      src::GPUArray{U, N}, srccrange::CartesianIndices{N}) where {T, U, N}
    shape = size(destcrange)
    if shape != size(srccrange)
        throw(DimensionMismatch("Ranges don't match their size. Found: $shape, $(size(srccrange))"))
    end
    len = length(destcrange)

    dest_offsets = first.(destcrange.indices) .- 1
    src_offsets = first.(srccrange.indices) .- 1
    gpu_call(copy_kernel!, dest,
             (dest, dest_offsets, src, src_offsets, shape, size(dest), size(src), len),
             len)
    dest
end

function Base.copyto!(dest::GPUArray{T, N}, destcrange::CartesianIndices{N},
                      src::AbstractArray{T, N}, srccrange::CartesianIndices{N}) where {T, N}
    # Is this efficient? Maybe!
    # TODO: compare to a pure intrinsic copyto implementation!
    # this would mean looping over linear sections of memory and
    # use copyto!(dest, offset::Integer, buffer(src), offset::Integer, amout::Integer)
    src_gpu = typeof(dest)(map(idx-> src[idx], srccrange))
    nrange = CartesianIndices(size(src_gpu))
    copyto!(dest, destcrange, src_gpu, nrange)
    dest
end

function Base.copyto!(dest::AbstractArray{T, N}, destcrange::CartesianIndices{N},
                      src::GPUArray{T, N}, srccrange::CartesianIndices{N}) where {T, N}
    # Is this efficient? Maybe!
    dest_gpu = similar(src, size(destcrange))
    nrange = CartesianIndices(size(dest_gpu))
    copyto!(dest_gpu, nrange, src, srccrange)
    copyto!(dest, destcrange, Array(dest_gpu), nrange)
    dest
end

## other

Base.copy(x::GPUArray) = identity.(x)

Base.deepcopy(x::GPUArray) = copy(x)

# reinterpret

#=
copied from julia base/array.jl
Copyright (c) 2009-2016: Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and other contributors:

https://github.com/JuliaLang/julia/contributors

Permission is hereby granted, free of charge, to any person obtaining a copie of this
software and associated documentation files (the "Software"), to deal in the Software
without restriction, including without limitation the rights to use, copy, modify, merge,
publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons
to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or
substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE
FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.
=#

import Base.reinterpret

"""
Unsafe reinterpret for backends to overload.
This makes it easier to do checks just on the high level.
"""
function unsafe_reinterpret end

function reinterpret(::Type{T}, a::GPUArray{S,1}) where T where S
    nel = (length(a)*sizeof(S)) ÷ sizeof(T)
    # TODO: maybe check that remainder is zero?
    return reinterpret(T, a, (nel,))
end

function reinterpret(::Type{T}, a::GPUArray{S}) where T where S
    if sizeof(S) != sizeof(T)
        throw(ArgumentError("result shape not specified"))
    end
    reinterpret(T, a, size(a))
end

function reinterpret(::Type{T}, a::GPUArray{S}, dims::NTuple{N, Integer}) where T where S where N
    if !isbitstype(T)
        throw(ArgumentError("cannot reinterpret Array{$(S)} to ::Type{Array{$(T)}}, type $(T) is not a bits type"))
    end
    if !isbitstype(S)
        throw(ArgumentError("cannot reinterpret Array{$(S)} to ::Type{Array{$(T)}}, type $(S) is not a bits type"))
    end
    nel = div(length(a)*sizeof(S),sizeof(T))
    if prod(dims) != nel
        throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $(nel)"))
    end
    unsafe_reinterpret(T, a, dims)
end

function Base._reshape(A::GPUArray{T}, dims::Dims) where T
    n = length(A)
    prod(dims) == n || throw(DimensionMismatch("parent has $n elements, which is incompatible with size $dims"))
    return unsafe_reinterpret(T, A, dims)
end
#ambig
function Base._reshape(A::GPUArray{T, 1}, dims::Tuple{Integer}) where T
    n = Base._length(A)
    prod(dims) == n || throw(DimensionMismatch("parent has $n elements, which is incompatible with size $dims"))
    return unsafe_reinterpret(T, A, dims)
end

# filtering

# TODO: filter!

# revert of JuliaLang/julia#31929
Base.filter(f, As::GPUArray) = As[map(f, As)::GPUArray{Bool}]