text:
int_rules_4_7_4 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.7.4*(c*trig)^m*(d*trig)^n =#
    @apply_utils Antiderivative(sin(~(a') + ~(b') * ~x) * sin(~(c') + ~(d') * ~x), ~x) => sin((~a - ~c) + (~b - ~d) * ~x) / (2 * (~b - ~d)) - sin(~a + ~c + (~b + ~d) * ~x) / (2 * (~b + ~d)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~b) ^ 2 - (~d) ^ 2, 0)
    @apply_utils Antiderivative(cos(~(a') + ~(b') * ~x) * cos(~(c') + ~(d') * ~x), ~x) => sin((~a - ~c) + (~b - ~d) * ~x) / (2 * (~b - ~d)) + sin(~a + ~c + (~b + ~d) * ~x) / (2 * (~b + ~d)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~b) ^ 2 - (~d) ^ 2, 0)
    @apply_utils Antiderivative(sin(~(a') + ~(b') * ~x) * cos(~(c') + ~(d') * ~x), ~x) => -(cos((~a - ~c) + (~b - ~d) * ~x)) / (2 * (~b - ~d)) - cos(~a + ~c + (~b + ~d) * ~x) / (2 * (~b + ~d)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~b) ^ 2 - (~d) ^ 2, 0)
    @apply_utils Antiderivative(cos(~(a') + ~(b') * ~x) ^ 2 * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (1 / 2) * Antiderivative((~g * sin(~c + ~d * ~x)) ^ ~p, ~x) + (1 / 2) * Antiderivative(cos(~c + ~d * ~x) * (~g * sin(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && IGtQ(~p / 2, 0)))
    @apply_utils Antiderivative(sin(~(a') + ~(b') * ~x) ^ 2 * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (1 / 2) * Antiderivative((~g * sin(~c + ~d * ~x)) ^ ~p, ~x) - (1 / 2) * Antiderivative(cos(~c + ~d * ~x) * (~g * sin(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && IGtQ(~p / 2, 0)))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~(m') * sin(~(c') + ~(d') * ~x) ^ ~(p'), ~x) => (2 ^ ~p / (~e) ^ ~p) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m + ~p) * sin(~a + ~b * ~x) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && IntegerQ(~p)))
    @apply_utils Antiderivative((~(f') * sin(~(a') + ~(b') * ~x)) ^ ~(n') * sin(~(c') + ~(d') * ~x) ^ ~(p'), ~x) => (2 ^ ~p / (~f) ^ ~p) * Antiderivative(cos(~a + ~b * ~x) ^ ~p * (~f * sin(~a + ~b * ~x)) ^ (~n + ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~n], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && IntegerQ(~p)))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~m * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => ((~e) ^ 2 * (~e * cos(~a + ~b * ~x)) ^ (~m - 2) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~p + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~g, ~m, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && EqQ((~m + ~p) - 1, 0))))
    @apply_utils Antiderivative((~(e') * sin(~(a') + ~(b') * ~x)) ^ ~m * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-((~e) ^ 2) * (~e * sin(~a + ~b * ~x)) ^ (~m - 2) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~p + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~g, ~m, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && EqQ((~m + ~p) - 1, 0))))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~(m') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-((~e * cos(~a + ~b * ~x)) ^ ~m) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (~b * ~g * ~m) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~g, ~m, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && EqQ(~m + 2 * ~p + 2, 0))))
    @apply_utils Antiderivative((~(e') * sin(~(a') + ~(b') * ~x)) ^ ~(m') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => ((~e * sin(~a + ~b * ~x)) ^ ~m * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (~b * ~g * ~m) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~g, ~m, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && EqQ(~m + 2 * ~p + 2, 0))))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~m * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => ((~e) ^ 2 * (~e * cos(~a + ~b * ~x)) ^ (~m - 2) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~p + 1)) + (((~e) ^ 4 * ((~m + ~p) - 1)) / (4 * (~g) ^ 2 * (~p + 1))) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m - 4) * (~g * sin(~c + ~d * ~x)) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~m, 2) && (LtQ(~p, -1) && ((GtQ(~m, 3) || EqQ(~p, -3 / 2)) && IntegersQ(2 * ~m, 2 * ~p)))))))
    @apply_utils Antiderivative((~(e') * sin(~(a') + ~(b') * ~x)) ^ ~m * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-((~e) ^ 2) * (~e * sin(~a + ~b * ~x)) ^ (~m - 2) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~p + 1)) + (((~e) ^ 4 * ((~m + ~p) - 1)) / (4 * (~g) ^ 2 * (~p + 1))) * Antiderivative((~e * sin(~a + ~b * ~x)) ^ (~m - 4) * (~g * sin(~c + ~d * ~x)) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~m, 2) && (LtQ(~p, -1) && ((GtQ(~m, 3) || EqQ(~p, -3 / 2)) && IntegersQ(2 * ~m, 2 * ~p)))))))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~m * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => ((~e * cos(~a + ~b * ~x)) ^ ~m * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~p + 1)) + (((~e) ^ 2 * (~m + 2 * ~p + 2)) / (4 * (~g) ^ 2 * (~p + 1))) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m - 2) * (~g * sin(~c + ~d * ~x)) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~m, 1) && (LtQ(~p, -1) && (NeQ(~m + 2 * ~p + 2, 0) && ((LtQ(~p, -2) || EqQ(~m, 2)) && IntegersQ(2 * ~m, 2 * ~p))))))))
    @apply_utils Antiderivative((~(e') * sin(~(a') + ~(b') * ~x)) ^ ~m * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-((~e * sin(~a + ~b * ~x)) ^ ~m) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~p + 1)) + (((~e) ^ 2 * (~m + 2 * ~p + 2)) / (4 * (~g) ^ 2 * (~p + 1))) * Antiderivative((~e * sin(~a + ~b * ~x)) ^ (~m - 2) * (~g * sin(~c + ~d * ~x)) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~m, 1) && (LtQ(~p, -1) && (NeQ(~m + 2 * ~p + 2, 0) && ((LtQ(~p, -2) || EqQ(~m, 2)) && IntegersQ(2 * ~m, 2 * ~p))))))))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~m * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => ((~e) ^ 2 * (~e * cos(~a + ~b * ~x)) ^ (~m - 2) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~m + 2 * ~p)) + (((~e) ^ 2 * ((~m + ~p) - 1)) / (~m + 2 * ~p)) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m - 2) * (~g * sin(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~g, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~m, 1) && (NeQ(~m + 2 * ~p, 0) && IntegersQ(2 * ~m, 2 * ~p))))))
    @apply_utils Antiderivative((~(e') * sin(~(a') + ~(b') * ~x)) ^ ~m * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-((~e) ^ 2) * (~e * sin(~a + ~b * ~x)) ^ (~m - 2) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~m + 2 * ~p)) + (((~e) ^ 2 * ((~m + ~p) - 1)) / (~m + 2 * ~p)) * Antiderivative((~e * sin(~a + ~b * ~x)) ^ (~m - 2) * (~g * sin(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~g, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~m, 1) && (NeQ(~m + 2 * ~p, 0) && IntegersQ(2 * ~m, 2 * ~p))))))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~m * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-((~e * cos(~a + ~b * ~x)) ^ ~m) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~m + ~p + 1)) + ((~m + 2 * ~p + 2) / ((~e) ^ 2 * (~m + ~p + 1))) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m + 2) * (~g * sin(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~g, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (LtQ(~m, -1) && (NeQ(~m + 2 * ~p + 2, 0) && (NeQ(~m + ~p + 1, 0) && IntegersQ(2 * ~m, 2 * ~p)))))))
    @apply_utils Antiderivative((~(e') * sin(~(a') + ~(b') * ~x)) ^ ~m * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => ((~e * sin(~a + ~b * ~x)) ^ ~m * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~m + ~p + 1)) + ((~m + 2 * ~p + 2) / ((~e) ^ 2 * (~m + ~p + 1))) * Antiderivative((~e * sin(~a + ~b * ~x)) ^ (~m + 2) * (~g * sin(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~g, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (LtQ(~m, -1) && (NeQ(~m + 2 * ~p + 2, 0) && (NeQ(~m + ~p + 1, 0) && IntegersQ(2 * ~m, 2 * ~p)))))))
    @apply_utils Antiderivative(cos(~(a') + ~(b') * ~x) * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (2 * sin(~a + ~b * ~x) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~d * (2 * ~p + 1)) + ((2 * ~p * ~g) / (2 * ~p + 1)) * Antiderivative(sin(~a + ~b * ~x) * (~g * sin(~c + ~d * ~x)) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~p, 0) && IntegerQ(2 * ~p)))))
    @apply_utils Antiderivative(sin(~(a') + ~(b') * ~x) * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-2 * cos(~a + ~b * ~x) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~d * (2 * ~p + 1)) + ((2 * ~p * ~g) / (2 * ~p + 1)) * Antiderivative(cos(~a + ~b * ~x) * (~g * sin(~c + ~d * ~x)) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~p, 0) && IntegerQ(2 * ~p)))))
    @apply_utils Antiderivative(cos(~(a') + ~(b') * ~x) * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (cos(~a + ~b * ~x) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~p + 1)) + ((2 * ~p + 3) / (2 * ~g * (~p + 1))) * Antiderivative(sin(~a + ~b * ~x) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (LtQ(~p, -1) && IntegerQ(2 * ~p)))))
    @apply_utils Antiderivative(sin(~(a') + ~(b') * ~x) * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-(sin(~a + ~b * ~x)) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~p + 1)) + ((2 * ~p + 3) / (2 * ~g * (~p + 1))) * Antiderivative(cos(~a + ~b * ~x) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (LtQ(~p, -1) && IntegerQ(2 * ~p)))))
    @apply_utils Antiderivative(cos(~(a') + ~(b') * ~x) / sqrt(sin(~(c') + ~(d') * ~x)), ~x) => -(asin(cos(~a + ~b * ~x) - sin(~a + ~b * ~x))) / ~d + log(cos(~a + ~b * ~x) + sin(~a + ~b * ~x) + sqrt(sin(~c + ~d * ~x))) / ~d <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && EqQ(~d / ~b, 2))
    @apply_utils Antiderivative(sin(~(a') + ~(b') * ~x) / sqrt(sin(~(c') + ~(d') * ~x)), ~x) => -(asin(cos(~a + ~b * ~x) - sin(~a + ~b * ~x))) / ~d - log(cos(~a + ~b * ~x) + sin(~a + ~b * ~x) + sqrt(sin(~c + ~d * ~x))) / ~d <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && EqQ(~d / ~b, 2))
    @apply_utils Antiderivative((~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p / cos(~(a') + ~(b') * ~x), ~x) => 2 * ~g * Antiderivative(sin(~a + ~b * ~x) * (~g * sin(~c + ~d * ~x)) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~g, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && IntegerQ(2 * ~p))))
    @apply_utils Antiderivative((~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p / sin(~(a') + ~(b') * ~x), ~x) => 2 * ~g * Antiderivative(cos(~a + ~b * ~x) * (~g * sin(~c + ~d * ~x)) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~g, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && IntegerQ(2 * ~p))))
    #= Antiderivative(((~e')*cos((~a')+(~b')*(~x)))^(~m)*((~g')*sin((~c')+(~d')*(~x)))^(~p),~x) :=   -(e*cos(a+b*x))^(m+1)*sin(a+b*x)*(g*sin(c+d*x))^p/(b*e*(m+p+1)*(sin( a+b*x)^2)^((p+1)/2))* HypergeometricFunctions._₂F₁(-(p-1)/2,(m+p+1)/2,(m+p+3)/2,cos(a+b*x)^2) <--  FreeQ([a,b,c,d,e,g,m,p],x) && EqQ(b*c-a*d,0) && EqQ(d/b,2) &&  Not(IntegerQ(p)) && Not(IntegerQ(m+p)) =#
    #= Antiderivative(((~f')*sin((~a')+(~b')*(~x)))^(~n')*((~g')*sin((~c')+(~d')*(~x)))^(~p),~x) :=   -cos(a+b*x)*(f*sin(a+b*x))^(n+1)*(g*sin(c+d*x))^p/(b*f*(p+1)*(sin(a+ b*x)^2)^((n+p+1)/2))* HypergeometricFunctions._₂F₁(-(n+p-1)/2,(p+1)/2,(p+3)/2,cos(a+b*x)^2) <--  FreeQ([a,b,c,d,f,g,n,p],x) && EqQ(b*c-a*d,0) && EqQ(d/b,2) &&  Not(IntegerQ(p)) && Not(IntegerQ(n+p)) =#
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~(m') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => ((~g * sin(~c + ~d * ~x)) ^ ~p / ((~e * cos(~a + ~b * ~x)) ^ ~p * sin(~a + ~b * ~x) ^ ~p)) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m + ~p) * sin(~a + ~b * ~x) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~g, ~m, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~(f') * sin(~(a') + ~(b') * ~x)) ^ ~(n') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => ((~g * sin(~c + ~d * ~x)) ^ ~p / (cos(~a + ~b * ~x) ^ ~p * (~f * sin(~a + ~b * ~x)) ^ ~p)) * Antiderivative(cos(~a + ~b * ~x) ^ ~p * (~f * sin(~a + ~b * ~x)) ^ (~n + ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~g, ~n, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative(cos(~(a') + ~(b') * ~x) ^ 2 * sin(~(a') + ~(b') * ~x) ^ 2 * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (1 / 4) * Antiderivative((~g * sin(~c + ~d * ~x)) ^ ~p, ~x) - (1 / 4) * Antiderivative(cos(~c + ~d * ~x) ^ 2 * (~g * sin(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && IGtQ(~p / 2, 0)))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~(m') * (~(f') * sin(~(a') + ~(b') * ~x)) ^ ~(n') * sin(~(c') + ~(d') * ~x) ^ ~(p'), ~x) => (2 ^ ~p / ((~e) ^ ~p * (~f) ^ ~p)) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m + ~p) * (~f * sin(~a + ~b * ~x)) ^ (~n + ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && IntegerQ(~p)))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~(m') * (~(f') * sin(~(a') + ~(b') * ~x)) ^ ~(n') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (~e * (~e * cos(~a + ~b * ~x)) ^ (~m - 1) * (~f * sin(~a + ~b * ~x)) ^ (~n + 1) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~b * ~f * (~n + ~p + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && EqQ(~m + ~p + 1, 0))))
    @apply_utils Antiderivative((~(e') * sin(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * cos(~(a') + ~(b') * ~x)) ^ ~n * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-(~e) * (~e * sin(~a + ~b * ~x)) ^ (~m - 1) * (~f * cos(~a + ~b * ~x)) ^ (~n + 1) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~b * ~f * (~n + ~p + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && EqQ(~m + ~p + 1, 0))))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~(m') * (~(f') * sin(~(a') + ~(b') * ~x)) ^ ~(n') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-((~e * cos(~a + ~b * ~x)) ^ (~m + 1)) * (~f * sin(~a + ~b * ~x)) ^ (~n + 1) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~b * ~e * ~f * (~m + ~p + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (EqQ(~m + ~n + 2 * ~p + 2, 0) && NeQ(~m + ~p + 1, 0)))))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * sin(~(a') + ~(b') * ~x)) ^ ~n * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => ((~e) ^ 2 * (~e * cos(~a + ~b * ~x)) ^ (~m - 2) * (~f * sin(~a + ~b * ~x)) ^ ~n * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~n + ~p + 1)) + (((~e) ^ 4 * ((~m + ~p) - 1)) / (4 * (~g) ^ 2 * (~n + ~p + 1))) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m - 4) * (~f * sin(~a + ~b * ~x)) ^ ~n * (~g * sin(~c + ~d * ~x)) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~m, 3) && (LtQ(~p, -1) && (NeQ(~n + ~p + 1, 0) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p)))))))
    @apply_utils Antiderivative((~(e') * sin(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * cos(~(a') + ~(b') * ~x)) ^ ~n * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-((~e) ^ 2) * (~e * sin(~a + ~b * ~x)) ^ (~m - 2) * (~f * cos(~a + ~b * ~x)) ^ ~n * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~n + ~p + 1)) + (((~e) ^ 4 * ((~m + ~p) - 1)) / (4 * (~g) ^ 2 * (~n + ~p + 1))) * Antiderivative((~e * sin(~a + ~b * ~x)) ^ (~m - 4) * (~f * cos(~a + ~b * ~x)) ^ ~n * (~g * sin(~c + ~d * ~x)) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~m, 3) && (LtQ(~p, -1) && (NeQ(~n + ~p + 1, 0) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p)))))))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * sin(~(a') + ~(b') * ~x)) ^ ~(n') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => ((~e * cos(~a + ~b * ~x)) ^ ~m * (~f * sin(~a + ~b * ~x)) ^ ~n * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~n + ~p + 1)) + (((~e) ^ 2 * (~m + ~n + 2 * ~p + 2)) / (4 * (~g) ^ 2 * (~n + ~p + 1))) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m - 2) * (~f * sin(~a + ~b * ~x)) ^ ~n * (~g * sin(~c + ~d * ~x)) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~m, 1) && (LtQ(~p, -1) && (NeQ(~m + ~n + 2 * ~p + 2, 0) && (NeQ(~n + ~p + 1, 0) && (IntegersQ(2 * ~m, 2 * ~n, 2 * ~p) && (LtQ(~p, -2) || (EqQ(~m, 2) || EqQ(~m, 3)))))))))))
    @apply_utils Antiderivative((~(e') * sin(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * cos(~(a') + ~(b') * ~x)) ^ ~(n') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-((~e * sin(~a + ~b * ~x)) ^ ~m) * (~f * cos(~a + ~b * ~x)) ^ ~n * (~g * sin(~c + ~d * ~x)) ^ (~p + 1)) / (2 * ~b * ~g * (~n + ~p + 1)) + (((~e) ^ 2 * (~m + ~n + 2 * ~p + 2)) / (4 * (~g) ^ 2 * (~n + ~p + 1))) * Antiderivative((~e * sin(~a + ~b * ~x)) ^ (~m - 2) * (~f * cos(~a + ~b * ~x)) ^ ~n * (~g * sin(~c + ~d * ~x)) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~m, 1) && (LtQ(~p, -1) && (NeQ(~m + ~n + 2 * ~p + 2, 0) && (NeQ(~n + ~p + 1, 0) && (IntegersQ(2 * ~m, 2 * ~n, 2 * ~p) && (LtQ(~p, -2) || (EqQ(~m, 2) || EqQ(~m, 3)))))))))))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * sin(~(a') + ~(b') * ~x)) ^ ~n * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (~e * (~e * cos(~a + ~b * ~x)) ^ (~m - 1) * (~f * sin(~a + ~b * ~x)) ^ (~n + 1) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~b * ~f * (~n + ~p + 1)) + (((~e) ^ 2 * ((~m + ~p) - 1)) / ((~f) ^ 2 * (~n + ~p + 1))) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m - 2) * (~f * sin(~a + ~b * ~x)) ^ (~n + 2) * (~g * sin(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~m, 1) && (LtQ(~n, -1) && (NeQ(~n + ~p + 1, 0) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p)))))))
    @apply_utils Antiderivative((~(e') * sin(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * cos(~(a') + ~(b') * ~x)) ^ ~n * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-(~e) * (~e * sin(~a + ~b * ~x)) ^ (~m - 1) * (~f * cos(~a + ~b * ~x)) ^ (~n + 1) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~b * ~f * (~n + ~p + 1)) + (((~e) ^ 2 * ((~m + ~p) - 1)) / ((~f) ^ 2 * (~n + ~p + 1))) * Antiderivative((~e * sin(~a + ~b * ~x)) ^ (~m - 2) * (~f * cos(~a + ~b * ~x)) ^ (~n + 2) * (~g * sin(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~m, 1) && (LtQ(~n, -1) && (NeQ(~n + ~p + 1, 0) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p)))))))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * sin(~(a') + ~(b') * ~x)) ^ ~(n') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (~e * (~e * cos(~a + ~b * ~x)) ^ (~m - 1) * (~f * sin(~a + ~b * ~x)) ^ (~n + 1) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~b * ~f * (~m + ~n + 2 * ~p)) + (((~e) ^ 2 * ((~m + ~p) - 1)) / (~m + ~n + 2 * ~p)) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m - 2) * (~f * sin(~a + ~b * ~x)) ^ ~n * (~g * sin(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~m, 1) && (NeQ(~m + ~n + 2 * ~p, 0) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p))))))
    @apply_utils Antiderivative((~(e') * sin(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * cos(~(a') + ~(b') * ~x)) ^ ~(n') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-(~e) * (~e * sin(~a + ~b * ~x)) ^ (~m - 1) * (~f * cos(~a + ~b * ~x)) ^ (~n + 1) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~b * ~f * (~m + ~n + 2 * ~p)) + (((~e) ^ 2 * ((~m + ~p) - 1)) / (~m + ~n + 2 * ~p)) * Antiderivative((~e * sin(~a + ~b * ~x)) ^ (~m - 2) * (~f * cos(~a + ~b * ~x)) ^ ~n * (~g * sin(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (GtQ(~m, 1) && (NeQ(~m + ~n + 2 * ~p, 0) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p))))))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * sin(~(a') + ~(b') * ~x)) ^ ~(n') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-(~f) * (~e * cos(~a + ~b * ~x)) ^ (~m + 1) * (~f * sin(~a + ~b * ~x)) ^ (~n - 1) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~b * ~e * (~m + ~n + 2 * ~p)) + ((2 * ~f * ~g * ((~n + ~p) - 1)) / (~e * (~m + ~n + 2 * ~p))) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m + 1) * (~f * sin(~a + ~b * ~x)) ^ (~n - 1) * (~g * sin(~c + ~d * ~x)) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (LtQ(~m, -1) && (GtQ(~n, 0) && (GtQ(~p, 0) && (NeQ(~m + ~n + 2 * ~p, 0) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p))))))))
    @apply_utils Antiderivative((~(e') * sin(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * cos(~(a') + ~(b') * ~x)) ^ ~(n') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (~f * (~e * sin(~a + ~b * ~x)) ^ (~m + 1) * (~f * cos(~a + ~b * ~x)) ^ (~n - 1) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~b * ~e * (~m + ~n + 2 * ~p)) + ((2 * ~f * ~g * ((~n + ~p) - 1)) / (~e * (~m + ~n + 2 * ~p))) * Antiderivative((~e * sin(~a + ~b * ~x)) ^ (~m + 1) * (~f * cos(~a + ~b * ~x)) ^ (~n - 1) * (~g * sin(~c + ~d * ~x)) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (LtQ(~m, -1) && (GtQ(~n, 0) && (GtQ(~p, 0) && (NeQ(~m + ~n + 2 * ~p, 0) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p))))))))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * sin(~(a') + ~(b') * ~x)) ^ ~(n') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-((~e * cos(~a + ~b * ~x)) ^ (~m + 1)) * (~f * sin(~a + ~b * ~x)) ^ (~n + 1) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~b * ~e * ~f * (~m + ~p + 1)) + ((~f * (~m + ~n + 2 * ~p + 2)) / (2 * ~e * ~g * (~m + ~p + 1))) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m + 1) * (~f * sin(~a + ~b * ~x)) ^ (~n - 1) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (LtQ(~m, -1) && (GtQ(~n, 0) && (LtQ(~p, -1) && (NeQ(~m + ~n + 2 * ~p + 2, 0) && (NeQ(~m + ~p + 1, 0) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p)))))))))
    @apply_utils Antiderivative((~(e') * sin(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * cos(~(a') + ~(b') * ~x)) ^ ~(n') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => ((~e * sin(~a + ~b * ~x)) ^ (~m + 1) * (~f * cos(~a + ~b * ~x)) ^ (~n + 1) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~b * ~e * ~f * (~m + ~p + 1)) + ((~f * (~m + ~n + 2 * ~p + 2)) / (2 * ~e * ~g * (~m + ~p + 1))) * Antiderivative((~e * sin(~a + ~b * ~x)) ^ (~m + 1) * (~f * cos(~a + ~b * ~x)) ^ (~n - 1) * (~g * sin(~c + ~d * ~x)) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (LtQ(~m, -1) && (GtQ(~n, 0) && (LtQ(~p, -1) && (NeQ(~m + ~n + 2 * ~p + 2, 0) && (NeQ(~m + ~p + 1, 0) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p)))))))))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * sin(~(a') + ~(b') * ~x)) ^ ~(n') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => (-((~e * cos(~a + ~b * ~x)) ^ (~m + 1)) * (~f * sin(~a + ~b * ~x)) ^ (~n + 1) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~b * ~e * ~f * (~m + ~p + 1)) + ((~m + ~n + 2 * ~p + 2) / ((~e) ^ 2 * (~m + ~p + 1))) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m + 2) * (~f * sin(~a + ~b * ~x)) ^ ~n * (~g * sin(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (LtQ(~m, -1) && (NeQ(~m + ~n + 2 * ~p + 2, 0) && (NeQ(~m + ~p + 1, 0) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p)))))))
    @apply_utils Antiderivative((~(e') * sin(~(a') + ~(b') * ~x)) ^ ~m * (~(f') * cos(~(a') + ~(b') * ~x)) ^ ~(n') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => ((~e * sin(~a + ~b * ~x)) ^ (~m + 1) * (~f * cos(~a + ~b * ~x)) ^ (~n + 1) * (~g * sin(~c + ~d * ~x)) ^ ~p) / (~b * ~e * ~f * (~m + ~p + 1)) + ((~m + ~n + 2 * ~p + 2) / ((~e) ^ 2 * (~m + ~p + 1))) * Antiderivative((~e * sin(~a + ~b * ~x)) ^ (~m + 2) * (~f * cos(~a + ~b * ~x)) ^ ~n * (~g * sin(~c + ~d * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && (Not(IntegerQ(~p)) && (LtQ(~m, -1) && (NeQ(~m + ~n + 2 * ~p + 2, 0) && (NeQ(~m + ~p + 1, 0) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p)))))))
    #= Antiderivative(((~e')*cos((~a')+(~b')*(~x)))^(~m)*((~f')*sin((~a')+(~b')*(~x)))^(~n')*((~g')*sin((~c') +(~d')*(~x)))^(~p),~x) :=  -(e*cos(a+b*x))^(m+1)*(f*sin(a+b*x))^(n+1)*(g*sin(c+d*x))^p/(b*e*f*(m+ p+1)*(sin(a+b*x)^2)^((n+p+1)/2))* HypergeometricFunctions._₂F₁(-(n+p-1)/2,(m+p+1)/2,(m+p+3)/2,cos(a+b*x)^2) <--  FreeQ([a,b,c,d,e,f,g,m,n,p],x) && EqQ(b*c-a*d,0) && EqQ(d/b,2) &&  Not(IntegerQ(p)) && Not(IntegerQ(m+p)) && Not(IntegerQ(n+p)) =#
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~(m') * (~(f') * sin(~(a') + ~(b') * ~x)) ^ ~(n') * (~(g') * sin(~(c') + ~(d') * ~x)) ^ ~p, ~x) => ((~g * sin(~c + ~d * ~x)) ^ ~p / ((~e * cos(~a + ~b * ~x)) ^ ~p * (~f * sin(~a + ~b * ~x)) ^ ~p)) * Antiderivative((~e * cos(~a + ~b * ~x)) ^ (~m + ~p) * (~f * sin(~a + ~b * ~x)) ^ (~n + ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && (EqQ(~d / ~b, 2) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~(e') * cos(~(a') + ~(b') * ~x)) ^ ~(m') * sin(~(c') + ~(d') * ~x), ~x) => (-((~m + 2)) * (~e * cos(~a + ~b * ~x)) ^ (~m + 1) * cos((~m + 1) * (~a + ~b * ~x))) / (~d * ~e * (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (EqQ(~b * ~c - ~a * ~d, 0) && EqQ(~d / ~b, Abs(~m + 2)))
end

max_stars_repo_path: src/rules/4 Trig functions/4.7 Miscellaneous/4.7.4 (c trig)^m (d trig)^n.jl


text:
"""
    append_to_file(filename, str)

Opens filename, appends str to it, and closes filename.

If filename is not a string but us type Base.PipeEndpoint (an IO stream)
then simply prints to it, without trying to open or close
"""
function append_to_file(filename, str)
    if typeof(filename)<:IO
        @printf(filename, "%s", str)
    else
        fstr = open(filename, "a")
        @printf(fstr, "%s", str)
        close(fstr)
    end
end

"""
evaluated_expression = replacer(P::String, mypars)

Given a string representing an expression to be evaluated, and a dictionary of
keys that are strings representing variables with values that are numbers,
parses the string into an expression, and substitutes any variables matching
the keys in mypars with the corresponding numeric values; finally, evaluates the
expression and returns the result.

# PARAMETERS:

- P::String   The expression to be evaluated, for example "t1*10 + sqrt(t2)"

- mypars      A dictionary mapping variable names to values, for example
                Dict("t1"=>5, "t2"=>100)

# RETURNS:

The result of evaluating the corresponding expression. Any variables that cannot
be instantiated into values will result in an Undefvar error

# EXAMPLE:

```jldoctest
julia>  replacer("t1*10+sqrt(t2)", Dict("t"=>-3, "t1"=>5, "t2"=>100))

15

```

"""
function replacer(P::String, mypars)
    return replacer(parse(P), mypars)
end

"""
evaluated_expression = replacer(P::Expr, mypars)

Given an expression to be evaluated, and a dictionary of
keys that are strings representing variables with values that are numbers,
substitutes any variables matching
the keys in mypars with the corresponding numeric values; finally, evaluates the
expression and returns the result.

# PARAMETERS:

- P::String   The expression to be evaluated, for example parse("t1*10 + sqrt(t2))"

- mypars      A dictionary mapping variable names to values, for example
                Dict("t1"=>5, "t2"=>100)

# RETURNS:

The result of evaluating the corresponding expression. Any variables that cannot
be instantiated into values will result in an Undefvar error

# EXAMPLE:

```jldoctest
julia>  replacer(parse("t1*10+sqrt(t2))", Dict("t"=>-3, "t1"=>5, "t2"=>100))

15

```

"""
function replacer(P, mypars)   # run through an expression tree, replacing known symbols with their values, then evaluate
    mypars = symbol_key_ize(mypars)
    ks = collect(keys(mypars))

    if typeof(P)<:Symbol
        idx = find(ks .== P)
        if length(idx)>0
            P = mypars[ks[idx[1]]]
        end
        return P
    end
    for i=1:length(P.args)
        if typeof(P.args[i])<:Expr || typeof(P.args[i])<:Symbol
            P.args[i] = replacer(P.args[i], mypars)
        end
    end
    # @printf("P = \n"); print(P)
    return eval(P)
end

"""
function print_vector(vec)

Takes a vector and uses @printf to put it on the screen with [%.3f, %.3f] format.

If passed a symbol (which must evaluate to a vector), then prints the string for that symbol,
an equals sign, the vector, and ends by adding a carriage return \n.
"""
function print_vector(vec)
    print_vector(STDOUT, vec)
end

"""
function print_vector(stream::IO, vec)

Takes a vector and uses @sprintf to put it on stream IO with [%.3f, %.3f] format.

If passed a symbol (which must evaluate to a vector), then prints the string for that symbol,
an equals sign, the vector, and ends by adding a carriage return \n.
"""
function print_vector(stream::IO, vec)

    if typeof(vec)==Symbol
        mystr = string(vec)
        @printf(stream, "%s = ", mystr);
        print_vector(stream, eval(vec))
        @printf(stream, "\n");
        return
    end

    @printf stream "["
    for p in [1:length(vec);]
        @printf(stream, "%.3f", vec[p])
        if p < length(vec) @printf(stream, ", "); end
    end
    @printf(stream, "]")
end

"""
function print_vector_g(vec)

Takes a vector and uses @printf to put it on the screen with [%g, %g] format.

If passed a symbol (which must evaluate to a vector), then prints the string for that symbol,
an equals sign, the vector, and ends by adding a carriage return \n.
"""
function print_vector_g(vec)
    print_vector_g(STDOUT, vec)
end

"""
function print_vector_g(stream::IO, vec)

Takes a vector and uses @printf to put it on stream with [%g, %g] format.

If passed a symbol (which must evaluate to a vector), then prints the string for that symbol,
an equals sign, the vector, and ends by adding a carriage return \n.
"""
function print_vector_g(stream::IO, vec)

    if typeof(vec)==Symbol
        mystr = string(vec)
        @printf(stream, "%s = ", mystr);
        print_vector_g(stream, eval(vec))
        @printf(stream, "\n");
        return
    end

    @printf stream "["
    for p in [1:length(vec);]
        @printf(stream, "%g", vec[p])
        if p < length(vec) @printf(stream, ", "); end
    end
    @printf(stream, "]")
end

"""
y = two_level_copy(x)

Like copy(x), but can go down a level. Can handle both Arrays and Dicts, otherwise gets confused.

EXAMPLE:

p = [1, 2, 3]
z = [p, 4]
y = Dict(:my=>p)

c = copy(z)
d = copy(y)

alpha = two_level_copy(z)
beta  = two_level_copy(y)
p[1]=1000

print("The inner levels of c and d are affected by the change to p:\n")
print(c); print("\n")
print(d); print("\n")
print("But the inner levels of alpha and beta are not:\n")
print(alpha); print("\n")
print(beta); print("\n")

"""
function two_level_copy(x)
    if typeof(x)<:Array
        y = copy(x)
        for i=1:length(x)
            if typeof(x[i])<:Tuple; y[i]=x[i];
            else y[i] = copy(x[i]) end;
        end
    elseif typeof(x)<:Dict
        y = copy(x)
        allkeys = keys(x)
        for i in allkeys
            if typeof(x[i])<:Tuple; y[i]=x[i];
            else y[i] = copy(x[i]) end;
        end
    else
        error(@sprintf("two_level_copy: Don't know how to handle type %s\n", typeof(x)))
    end
    return y
end

"""
new_fname = next_file(fbasename, ndigits)

Returns a numbered and presumably unused filename starting with the string fbasename, followed by an integer
digit. The returned integer will be one higher than the number of existing filenames starting with fbasename,
and will be written with ndigits numbers, using leading zeros if necessary.

# EXAMPLE:

If there are already 8 files starting with "Mydir/model" then

> next_file("Mydir/model_", 4)

"Mydir/model_0009"
"""
function next_file(fbasename, ndigits)
    mydir  = dirname(fbasename)
    myfile = basename(fbasename)
    if length(mydir)>0
        fnames = readdir(mydir)
    else
        fnames = readdir()
    end
    matched_filenames = Array{Bool}(length(fnames))
    for i=1:length(fnames)
        matched_filenames[i] = ismatch(Regex(@sprintf("^%s", myfile)), fnames[i])
    end

    mynum = length(find(matched_filenames))+1
    myname = @sprintf("%d", mynum)
    while length(myname)<ndigits
        myname = "0" * myname
    end

    if length(mydir)>0
        return mydir * "/" * myfile * myname
    else
        return myfile * myname
    end
end

"""
    fstring = num2fixed_string(n, ndigits)

Returns a string version of a positive integer, with
however many leading zeros are necessary to have
ndigits characters.

"""
function num2fixed_string(n, ndigits)
    if ndigits<=0
        error("ndigits must be bigger than zero")
    end

    if n<0
        error("n must be positive")
    end

    myname = @sprintf("%d", n)
    while length(myname)<ndigits
        myname = "0"*myname
    end

    return myname
end

"""
ad = ascii_key_ize(d)

Given a dictionary that has keys that can be converted to strings, returns a copy with all
keys converted to strings
"""
function ascii_key_ize(d)
    ad = Dict()
    for k in keys(d)
        get!(ad, string(k), d[k])
    end
    return ad
end

"""
sd = symbol_key_ize(d)

Given a dictionary that has keys that can be converted to Symbols, returns a copy with all
keys converted to Symbols
"""
function symbol_key_ize(d)
    sd = Dict()
    for k in keys(d)
        get!(sd, Symbol(k), d[k])
    end
    return sd
end

"""
vks = vectorize_dict(dictionary, ks)

Given a dictionary (in which) all keys are either strings or Symbols, and all values are Float64s),
and an array ks of keys into that dictionary, returns a Float64 array the same size as ks containing
the values. Each key is checked as either itself or the string version of itself or the Symbol version
of itself.

Thus the following all return the same

a = Dict(:this=>33.4, "that"=>28.7)

vectorize_dict(a, ["this", "that"])
vectorize_dict(a, [:this, "that"])
vectorize_dict(a, [:this, :that])
"""
function vectorize_dict(dictionary, ks)
    output = Array{Float64}(size(ks))
    for i=1:length(ks)
        if haskey(dictionary, ks[i])
            output[i] = dictionary[ks[i]]
        elseif typeof(ks[i])<:Symbol && haskey(dictionary, string(ks[i]))
            output[i] = dictionary[string(ks[i])]
        elseif typeof(ks[i])<:String && haskey(dictionary, Symbol(ks[i]))
            output[i] = dictionary[Symbol(ks[i])]
        else
            print("Troublesome key: "); print(ks[i]); print("\n")
            error("Found neither key nor string(key) nor Symbol(key) in the dictionary")
        end
    end
    return output
end

ext: jl


text:
_is_incremental(::Type{GLPK.Optimizer}) = true

function set_default_config!(ext::ExtensionType, d::GlobalOptimizer, m::Incremental{GLPK.Optimizer}, local_solver::Bool)
    return
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# L2,1 norm/Sum of norms of columns or rows (times a constant)

export NormL21

"""
**Sum-of-``L_2`` norms**

    NormL21(λ=1.0, dim=1)

Returns the function
```math
f(X) = λ⋅∑_i\\|x_i\\|
```
for a nonnegative `λ`, where ``x_i`` is the ``i``-th column of ``X`` if `dim == 1`, and the ``i``-th row of ``X`` if `dim == 2`.
In words, it is the sum of the Euclidean norms of the columns or rows.
"""
struct NormL21{R <: Real, I <: Integer} <: ProximableFunction
    lambda::R
    dim::I
    function NormL21{R,I}(lambda::R, dim::I) where {R <: Real, I <: Integer}
        if lambda < 0
            error("parameter λ must be nonnegative")
        else
            new(lambda, dim)
        end
    end
end

is_convex(f::NormL21) = true

NormL21(lambda::R=1.0, dim::I=1) where {R <: Real, I <: Integer} = NormL21{R, I}(lambda, dim)

function (f::NormL21)(X::AbstractArray{T,2}) where T <: RealOrComplex
    nslice = 0.0
    n21X = 0.0
    if f.dim == 1
        for j = 1:size(X,2)
            nslice = 0.0
            for i = 1:size(X,1)
                nslice += abs(X[i,j])^2
            end
            n21X += sqrt(nslice)
        end
    elseif f.dim == 2
        for i = 1:size(X,1)
            nslice = 0.0
            for j = 1:size(X,2)
                nslice += abs(X[i,j])^2
            end
            n21X += sqrt(nslice)
        end
    end
    return f.lambda*n21X
end

function prox!(Y::AbstractArray{T,2}, f::NormL21, X::AbstractArray{T,2}, gamma::Real=1.0) where T <: RealOrComplex
    gl = gamma*f.lambda
    nslice = zero(Float64)
    n21X = zero(Float64)
    if f.dim == 1
        for j = 1:size(X,2)
            nslice = 0.0
            for i = 1:size(X,1)
                nslice += abs(X[i,j])^2
            end
            nslice = sqrt(nslice)
            scal = 1-gl/nslice
            scal = scal <= 0.0 ? 0.0 : scal
            for i = 1:size(X,1)
                Y[i,j] = scal*X[i,j]
            end
            n21X += scal*nslice
        end
    elseif f.dim == 2
        for i = 1:size(X,1)
            nslice = 0.0
            for j = 1:size(X,2)
                nslice += abs(X[i,j])^2
            end
            nslice = sqrt(nslice)
            scal = 1-gl/nslice
            scal = scal <= 0.0 ? 0.0 : scal
            for j = 1:size(X,2)
                Y[i,j] = scal*X[i,j]
            end
            n21X += scal*nslice
        end
    end
    return f.lambda*n21X
end

fun_name(f::NormL21) = "sum of Euclidean norms"
fun_dom(f::NormL21) = "AbstractArray{Real,2}, AbstractArray{Complex,2}"
fun_expr(f::NormL21) = "x ↦ λsum(||x_i||)"
fun_params(f::NormL21) = "λ = $(f.lambda), dim = $(f.dim)"

function prox_naive(f::NormL21, X::AbstractArray{T,2}, gamma::Real=1.0) where T <: RealOrComplex
    Y = max.(0, 1.0 .- f.lambda*gamma./sqrt.(sum(abs.(X).^2, dims=f.dim))).*X
    return Y, f.lambda*sum(sqrt.(sum(abs.(Y).^2, dims=f.dim)))
end


text:
module Wrap

using MacroTools: isexpr, prewalk

using ..Inner, ..IRTools
import ..Inner: IR, Variable, Statement, Branch, BasicBlock, Meta, block!,
  unreachable, varmap, argument!, branch!, return!
import Core: CodeInfo, GotoNode, SSAValue
import Core.Compiler: IRCode, CFG, GotoIfNot, ReturnNode, StmtRange

@static if VERSION > v"1.6-"
  import Core.Compiler: InstructionStream
end

unvars(ex) = prewalk(x -> x isa Variable ? SSAValue(x.id) : x, ex)

function IRCode(ir::IR)
  defs = Dict()
  stmts, types, lines = [], [], Int32[]
  index = Int[]
  for b in IRTools.blocks(ir)
    if b.id == 1
      for (i, arg) in enumerate(IRTools.arguments(b))
        defs[arg] = Core.SlotNumber(i)
      end
    else
      @assert isempty(BasicBlock(b).args)
    end
    for (v, st) in b
      defs[v] = Variable(length(stmts)+1)
      ex = varmap(x -> get(defs, x, x), st.expr) |> unvars
      push!(stmts, ex)
      push!(types, st.type)
      push!(lines, st.line)
    end
    for br in BasicBlock(b).branches
      if IRTools.isreturn(br)
        x = get(defs, br.args[1], br.args[1]) |> unvars
        push!(stmts, ReturnNode(x))
      elseif br == unreachable
        push!(stmts, ReturnNode())
      elseif br.condition == nothing
        push!(stmts, GotoNode(br.block))
      else
        cond = get(defs, br.condition, br.condition) |> unvars
        push!(stmts, GotoIfNot(cond, br.block))
      end
      push!(types, Any); push!(lines, 0)
    end
    push!(index, length(stmts)+1)
  end
  ranges = StmtRange.([1, index[1:end-1]...], index.-1)
  succs = map.(x -> x.id, IRTools.successors.(IRTools.blocks(ir)))
  preds = map.(x -> x.id, IRTools.predecessors.(IRTools.blocks(ir)))
  bs = Core.Compiler.BasicBlock.(ranges, preds, succs)
  cfg = CFG(bs, index)
  flags = [0x00 for _ in stmts]
  sps = VERSION > v"1.2-" ? [] : Core.svec()

  @static if VERSION > v"1.6-"
    stmtinfo = Any[nothing for _ in 1:length(length(stmts))]
    stmts = InstructionStream(stmts, types, stmtinfo, lines, flags)
    meta = @static VERSION < v"1.9.0-DEV.472" ? [] : Expr[]
    IRCode(stmts, cfg, ir.lines, ir.blocks[1].argtypes, meta, sps)
  else
    IRCode(stmts, types, lines, flags, cfg, ir.lines, ir.blocks[1].argtypes, [], sps)
  end
end

if VERSION >= v"1.6.0-DEV.272"
    isgotoifnot(@nospecialize(ex)) = (@assert !isexpr(ex, :gotoifnot); ex isa Core.GotoIfNot)
    destruct_gotoifnot(@nospecialize(ex)) = Any[ex.cond, ex.dest]

    isreturn(@nospecialize(ex)) = (@assert !isexpr(ex, :return); ex isa ReturnNode)
    retval(@nospecialize(ex)) = ex.val
else
    isgotoifnot(@nospecialize(ex)) = isexpr(ex, :gotoifnot)
    destruct_gotoifnot(@nospecialize(ex)) = ex.args

    isreturn(@nospecialize(ex)) = isexpr(ex, :return)
    retval(@nospecialize(ex)) = ex.args[1]
end

function blockstarts(ci::CodeInfo)
  bs = Int[]
  terminator = false
  for i = 1:length(ci.code)
    ex = ci.code[i]
    if isgotoifnot(ex)
      _, dest = destruct_gotoifnot(ex)
      push!(bs, dest)
      terminator = true
    elseif ex isa GotoNode || isreturn(ex)
      ex isa GotoNode && push!(bs, ex.label)
      i < length(ci.code) && push!(bs, i+1)
      terminator = false
    elseif terminator
      push!(bs, i)
      terminator = false
    end
  end
  return sort(unique(bs))
end

# TODO more informative names, while avoiding clashes.
slotname(ci, s) = Symbol(:_, s.id)

function IR(ci::CodeInfo, nargs::Integer; meta = nothing)
  bs = blockstarts(ci)
  ir = IR(Core.LineInfoNode[ci.linetable...], meta = meta)
  _rename = Dict()
  rename(ex) = prewalk(ex) do x
    haskey(_rename, x) && return _rename[x]
    x isa Core.SlotNumber && return Inner.Slot(slotname(ci, x))
    return x
  end
  for i = 1:nargs
    _rename[Core.SlotNumber(i)] = argument!(ir)
  end
  for i = 1:length(ci.code)
    i in bs && block!(ir)
    ex = ci.code[i]
    if ex isa Core.NewvarNode
      continue
    elseif isexpr(ex, :enter)
      _rename[Core.SSAValue(i)] = push!(ir, Expr(:enter, findfirst(==(ex.args[1]), bs)+1))
    elseif ex isa GotoNode
      branch!(ir, findfirst(==(ex.label), bs)+1)
    elseif isgotoifnot(ex)
      cond, dest = destruct_gotoifnot(ex)
      branch!(ir, findfirst(==(dest), bs)+1,
              unless = rename(cond))
    elseif isreturn(ex)
      return!(ir, rename(retval(ex)))
    else
      _rename[Core.SSAValue(i)] = push!(ir, IRTools.stmt(rename(ex), line = ci.codelocs[i]))
    end
  end
  return ir
end

function IR(meta::Meta; slots = false, prune = true)
  ir = IR(meta.code, meta.nargs, meta = meta)
  slots && return ir
  ir = IRTools.ssa!(ir)
  if prune
    ir = ir |> IRTools.prune! |> IRTools.renumber
  end
  return ir
end

function IR(Ts::Type...; slots = false, prune = true)
  m = IRTools.meta(Tuple{Ts...})
  m == nothing && return
  IR(m, slots = slots, prune = prune)
end

end


text:
abstract type Results end

get_results(result::Results) = nothing
struct OperationsProblemResults <: Results
    variables::Dict{Symbol, DataFrames.DataFrame}
    total_cost::Dict
    optimizer_log::Dict
    time_stamp::DataFrames.DataFrame
end

get_duals(result::OperationsProblemResults) = nothing

function make_results(variables::Dict,
                      total_cost::Dict,
                      optimizer_log::Dict,
                      time_stamp::DataFrames.DataFrame)
    return OperationsProblemResults(variables, total_cost, optimizer_log, time_stamp)
end

function make_results(variables::Dict,
                      total_cost::Dict,
                      optimizer_log::Dict,
                      time_stamp::DataFrames.DataFrame,
                      duals::Dict)
    return AggregatedResults(variables, total_cost, optimizer_log, time_stamp, duals)
end
function get_variable(res_model::OperationsProblemResults, key::Symbol)
        try
            !isnothing(results.variables)
            return get(results.variables, key, nothing)
        catch
            error("No variable with key $(key) has been found.")
        end
end

function get_optimizer_log(results::OperationsProblemResults)
    return results.optimizer_log
end

function get_time_stamps(results::OperationsProblemResults, key::Symbol)
    return results.time_stamp
end

"""
    results = load_operation_results(path)

This function can be used to load results from a folder
of results from a single-step problem, or for a single foulder
within a simulation.

# Arguments
-`path::AbstractString = folder path`
-`directory::AbstractString = "2019-10-03T09-18-00"`: the foulder name that contains
feather files of the results.

# Example
```julia
results = load_operation_results("/Users/test/2019-10-03T09-18-00")
```
"""
function load_operation_results(folder_path::AbstractString)

    if isfile(folder_path)
        @error("not a folder path")
    end
    files_in_folder = collect(readdir(folder_path))
    variable_list = setdiff(files_in_folder, ["time_stamp.feather", "optimizer_log.json"])
    variables = Dict{Symbol, DataFrames.DataFrame}()

    for name in 1:length(variable_list)
        variable_name = splitext(variable_list[name])[1]
        file_path = joinpath(folder_path,variable_list[name])
        variables[Symbol(variable_name)] = Feather.read(file_path) #change key to variable
    end
    optimizer = JSON.parse(open(joinpath(folder_path, "optimizer_log.json")))
    time_stamp = Feather.read(joinpath(folder_path,"time_stamp.feather"))
    time_stamp = shorten_time_stamp(time_stamp)
    obj_value = Dict{Symbol, Any}(:OBJECTIVE_FUNCTION => optimizer["obj_value"])
    results = make_results(variables, obj_value, optimizer, time_stamp)
    return results
end

function shorten_time_stamp(time::DataFrames.DataFrame)
    time = time[1:(size(time,1)-1),:]
    return time
end


text:
export Halley

function Halley(nlp_stop,
				   PData :: PDataFact,
				   dₙ :: Vector,
				   gt :: Vector)
	# println("on est dans Halley")
	nlp_at_x = nlp_stop.current_state
	x = nlp_at_x.x
	# println("on a nlp_at_x et x")

	n = length(x)
	# @show n

	e₁ = zeros(n)
	# @show e₁
    e₁[1] = 1.0
	# @show e₁
    ∇³fdₙ = ARCTR.∇f³xuv(nlp_stop.pb, x, dₙ, e₁)
	# @show ∇³fdₙ
	# println("avant de rentrer dans le for")
	# @show n
    for i = 2:n
		# @show i
        eᵢ = zeros(n)
        eᵢ[i] = 1.0
        ∇³fdₙ = hcat(∇³fdₙ, ∇f³xuv(nlp_stop.pb, x, dₙ, eᵢ))
    end

	# @show PData.λ
	# @show PData.λ .== 0.0
	if PData.λ .== 0.0
		hess_plus_tenseur = (nlp_at_x.Hx .+ 0.5 .* ∇³fdₙ)
	else
		hess_plus_tenseur = (nlp_at_x.Hx .+ 0.5 .* ∇³fdₙ) .+ (PData.λ .* Matrix(1.0I, n, n))
	end
	# hess_plus_tenseur = (nlp_at_x.Hx .+ 0.5 .* ∇³fdₙ)
	# @show isposdef(hess_plus_tenseur)
	if isposdef(hess_plus_tenseur)
        (L, D, pp, ρ, ncomp) = ldlt_symm(hess_plus_tenseur, 'r')
   	else
		# println("on est ici")
	   	return dₙ
   	end
	# catch
 	# 	println("*******   Problem in LDLt Halley")
    #    	#res = PDataLDLt()
    #    	#res.OK = false
    #    	#return res
    #    	d = NaN * zeros(n)
    #    	return d, xdemi
   	# end

	if true in isnan.(D)
   	   println("*******   Problem in D from LDLt: NaN")
	   println(" cond (H) = $(cond(H))")
	   #res = PDataLDLt()
	   #res.OK = false
	   #return res
	   d = NaN * zeros(n)
	   return d, xdemi
   	end

	D = Hermitian(D)
   	X = eigen(D)
  	Δ = X.values
  	Q =  X.vectors

	ϵ2 =  sqrt(eps(eltype(gt)))
   	Γ = max.(abs.(Δ), ϵ2)

	d̃ = L \ nlp_at_x.gx[pp]
    d̂ = L' \ (Q * (Q' * d̃ ./ Γ))
    dH = -d̂[invperm(pp)]

	return dH
end


text:
"""
    ScaleTransform(l::Real)

Transformation that multiplies the input elementwise with `l`.

# Examples

```jldoctest
julia> l = rand(); t = ScaleTransform(l); X = rand(100, 10);

julia> map(t, ColVecs(X)) == ColVecs(l .* X)
true
```
"""
struct ScaleTransform{T<:Real} <: Transform
    s::Vector{T}
end

function ScaleTransform(s::T=1.0) where {T<:Real}
    return ScaleTransform{T}([s])
end

@functor ScaleTransform

set!(t::ScaleTransform, ρ::Real) = t.s .= [ρ]

(t::ScaleTransform)(x) = only(t.s) * x

_map(t::ScaleTransform, x::AbstractVector{<:Real}) = only(t.s) .* x
_map(t::ScaleTransform, x::ColVecs) = ColVecs(only(t.s) .* x.X)
_map(t::ScaleTransform, x::RowVecs) = RowVecs(only(t.s) .* x.X)

Base.isequal(t::ScaleTransform, t2::ScaleTransform) = isequal(only(t.s), only(t2.s))

Base.show(io::IO, t::ScaleTransform) = print(io, "Scale Transform (s = ", only(t.s), ")")

# Helpers

"""
    median_heuristic_transform(distance, x::AbstractVector)

Create a [`ScaleTransform`](@ref) that divides the input elementwise by the median
`distance` of the data points in `x`.

The `distance` has to support pairwise evaluation with `KernelFunctions.pairwise`. All
`PreMetric`s of the package [Distances.jl](https://github.com/JuliaStats/Distances.jl) such
as `Euclidean` satisfy this requirement automatically.

# Examples

```jldoctest
julia> using Distances, Statistics

julia> x = ColVecs(rand(100, 10));

julia> t = median_heuristic_transform(Euclidean(), x);

julia> y = map(t, x);

julia> median(euclidean(y[i], y[j]) for i in 1:10, j in 1:10 if i != j) ≈ 1
true
```
"""
function median_heuristic_transform(f, x::AbstractVector)
    # Compute pairwise distances between **different** elements
    n = length(x)
    distances = vec(pairwise(f, x))
    deleteat!(distances, 1:(n + 1):(n^2))

    return ScaleTransform(inv(median!(distances)))
end


text:
### RLEVectorsType

"""
# RLEVectors
`RLEVectors` is an alternate implementation of the Rle type from Bioconductor's
IRanges package by H. Pages, P. Aboyoun and M. Lawrence. RLEVectors represent a
vector with repeated values as the ordered set of values and repeat extents. In
the field of genomics, data of various types measured across the ~3 billion
letters in the human genome can often be represented in a few thousand runs.
It is useful to know the bounds of genome regions covered by these runs, the
values associated with these runs, and to be able to perform various
mathematical operations on these values.

`RLEVectors` can be created from a single vector or a vector of values and a
vector of run ends. In either case runs of values or zero length runs will
be compressed out. RLEVectors can be expanded to a full vector with `collect`.

## Aliases
Several aliases are defined for specific types of RLEVector (or collections thereof).

    FloatRle              RLEVector{Float64,UInt32}
    IntegerRle            RLEVector{Int64,UInt32}
    BoolRle               RLEVector{Bool,UInt32}
    StringRle             RLEVector{String,UInt32}
    RLEVectorList{T1,T2}  Vector{ RLEVector{T1,T2} }

## Constructors
`RLEVector`s can be created by specifying a vector to compress or the runvalues and run ends.

    x = RLEVector([1,1,2,2,3,3,4,4,4])
    x = RLEVector([4,5,6],[3,6,9])

## Describing `RLEVector` objects
`RLEVector`s implement the usual descriptive functions for an array as well as some that are
specific to the type.

* `length(x)` The full length of the vector, uncompressed
* `size(x)` Same as `length`, as for any other vector
* `size(x,dim)` Returns `(length(x),1) for dim == 1`
* `starts(x)` The index of the beginning of each run
* `widths(x)` The width of each run
* `ends(x)` The index of the end of each run
* `values(x)` The data value for each run
* `isempty(x)` Returns boolean, as for any other vector
* `nrun(x)` Returns the number of runs represented in the array
* `eltype(x)` Returns the element type of the runs
* `endtype(x)` Returns the element type of the run ends

"""
struct RLEVector{T1,T2 <: Integer} <: AbstractArray{T1,1}
  runvalues::Vector{T1}
  runends::Vector{T2}
  RLEVector{T1,T2}(runvalues,runends) where {T1,T2<:Integer} = new(runvalues,runends)
end

function RLEVector(runvalues::Vector{T1}, runends::Vector{T2}) where {T1, T2 <: Integer}
    runvalues, runends = ree!(runvalues,runends)
    RLEVector{T1,T2}(runvalues, runends)
end

function RLEVector(runvalues::BitVector, runends::Vector{T2}) where T2 <: Integer
    runvalues, runends = ree!(runvalues,runends)
    RLEVector{Bool,T2}(runvalues, runends)
end

function RLEVector(vec::Vector)
  runvalues, runends = ree(vec)
  RLEVector(runvalues, runends)
end

function RLEVector(x::UnitRange)
    RLEVector{eltype(x),Int64}(collect(x),collect(1:length(x)))
end

function RLEVector(x)
    RLEVector([x])
end

function RLEVector(runvalues::T1, runends::T2) where {T1,T2 <: Integer}
    RLEVector{T1,T2}([runvalues], [runends])
end

#  Having specific types of Rle would be useful for lists of the same type, but Julia does a good job noticing that
#  Could also be useful for method definitions
const FloatRle = RLEVector{Float64,UInt32}
const IntegerRle = RLEVector{Int64,UInt32}
const BoolRle = RLEVector{Bool,UInt32}
const StringRle = RLEVector{String,UInt32}
const RLEVectorList{T1,T2} = Vector{ RLEVector{T1,T2} }
@doc (@doc RLEVector) FloatRle,  IntegerRle, BoolRle, StringRle, RLEVectorList

# copy
Base.copy(x::RLEVector) = RLEVector(copy(x.runvalues), copy(x.runends))

# similar
function Base.similar(a::RLEVector{T1,T2}, ::Type{T}, dims::Tuple{Int}) where {T1,T2,T}
    RLEVector{T,T2}(Vector{T}(undef,1), [dims[1]])
end

# show
function Base.show(io::IO, ::MIME"text/plain", x::RLEVector)
    t = typeof(x)::DataType
    show(io, t)
    n = nrun(x)
    write(io,"\n Run values: ")
    Base.show_vector(io,x.runvalues,"[", "]")
    write(io,"\n Run ends: ")
    Base.show_vector(io,x.runends,"[", "]")
end

function Base.show(io::IO, x::RLEVector)
    write(io,"Values: ")
    Base.show_vector(io,values(x),"[", "]")
    write(io," Ends: ")
    Base.show_vector(io,ends(x),"[", "]")
end

function ree!(x::RLEVector)
    ree!(x.runvalues,x.runends)
end

# conversions
convert(::Type{Vector}, x::RLEVector) = collect(x)
convert(::Type{Set}, x::RLEVector) = Set(values(x))
convert(::Type{RLEVector}, x::Vector) = RLEVector(x)
promote_rule(::Type{Set}, ::Type{RLEVector}) = Set
convert(::Type{RLEVector{T1,T2}}, x::RLEVector) where {T1,T2} = RLEVector(convert(Vector{T1},values(x)),convert(Vector{T2},ends(x)))

# the basics
function collect(x::RLEVector)
  inverse_ree(x.runvalues,x.runends)
end

function isequal(x::RLEVector, y::RLEVector)
    isequal(x.runends,y.runends) && isequal(x.runvalues, y.runvalues)
end

Base.hash(a::RLEVector) = hash(a.runvalues, hash(a.runlengths, hash(:RLEVector)))

"""
    growat!(x::AbstractVector, i, insert_length)
    growat!(x::RLEVector, i, insert_length)
Increases length of vector by `insert_length` at index `i` while
moving values `i:n` down to accomodate the new spaces. For the
`RLEVector` method `insert_length` uninitialized *runs* are added.
"""
function growat!(x::AbstractVector, i, insert_length)
    len = length(x)
    resize!(x, len + insert_length)
    ind = len
    @inbounds for ind in len:-1:i
        x[ind + insert_length] = x[ind]
    end
    x
end

function growat!(x::RLEVector, i, insert_length)
    growat!(x.runvalues, i, insert_length)
    growat!(x.runends, i, insert_length)
    x
end


text:
################################################################################
# Quadruped Trajectory Generator
################################################################################

# Trajectory template
function low_foot_trajectory(r::T, N::Int) where T
	# r = radius of the foot traj
	# low foot trajectory decomposed into N steps (N+1 waypoints)
	X = Vector(range(-r, r, length = N+1))
	low_traj = [[x, 0, 0.0] for x in X]
	return low_traj
end

function high_foot_trajectory(r::T, β::T, N::Int) where T
	# high foot trajectory decomposed into N steps (N+1 waypoints)
	α = range(0, π, length = N+1)
	θ = π * (1 .- cos.(α)) / 2
	X = [r*cos(θi) for θi in θ]
	Z = [β*r*sin(θi) for θi in θ] # the facto 2 flatten the foo trajectory
	high_traj = [[X[i], 0.0, Z[i]] for i = 1:N+1]
	return high_traj
end

function foot_trajectory(r::T, β::T, N::Int) where T
	low_traj = low_foot_trajectory(r, N)
	high_traj = high_foot_trajectory(r, β, N)
	traj = [low_traj[1:end-1]; high_traj[1:end-1]]
	return traj
end

# Inverse kinematics: given the position of the first link, finds the knee and hip angles that will put the foot
# at the p_foot location
function IKatlas(mechanism::Mechanism, p_base, p_foot; leg::Symbol = :r)
	# starting point of the local search
	θ = [0.5, 1.0] # θhip, θknee
	for k = 1:10
		err = AtlasIKerror(mechanism, p_base, p_foot, θ; leg = leg)
		norm(err, Inf) < 1e-10 && continue
		∇ = FiniteDiff.finite_difference_jacobian(θ -> AtlasIKerror(mechanism, p_base, p_foot, θ; leg = leg), θ)
		θ -= ∇ \ err
	end
	return θ
end

function AtlasIKerror(mechanism::Mechanism, p_base, p_foot, θ; leg::Symbol = :r)
	set_minimal_coordinates!(mechanism, get_joint(mechanism, :floating_base), [p_base; zeros(3)])
	set_minimal_coordinates!(mechanism, get_joint(mechanism, Symbol(leg, :_leg_hpxyz)), [0.0, -θ[1], 0.0])
	set_minimal_coordinates!(mechanism, get_joint(mechanism, Symbol(leg, :_leg_kny)), [θ[2]])
	set_minimal_coordinates!(mechanism, get_joint(mechanism, Symbol(leg, :_leg_akxy)), [θ[1]-θ[2], 0.0])

	foot = get_body(mechanism, Symbol(leg, :_foot))
	foot_contacts = mechanism.contacts[findall(x -> x.parent_id == foot.id, mechanism.contacts)]
	p = mean([contact_location(contact, foot) for contact in foot_contacts]) # average of all contact locations for one foot
	err = p - p_foot
	return err[[1,3]]
end

function atlas_trajectory(mechanism::Mechanism{T}; timestep=0.05, β=0.5,
		αtorso=0.15, Δx=0.0, r=0.10, x=0.00, z=0.85, N=12, Ncycles=1) where T
	pL = [Δx, + 0.1145, 0]
	pR = [Δx, - 0.1145, 0]

	t = reverse(foot_trajectory(r, β, N))
	t_dephased = [t[N+1:2N]; t[1:N]]
	# Foot positions
	tL = [pL + ti for ti in t_dephased]
	tR = [pR + ti for ti in t]

	# Leg angles
	p_base = [x,0,z]
	θL = [IKatlas(mechanism, p_base, tL[i], leg = :l) for i = 1:2N]
	θR = [IKatlas(mechanism, p_base, tR[i], leg = :r) for i = 1:2N]

	# adding angular velocities
	ωL = [(θL[i%(2N)+1] - θL[i]) / timestep for i = 1:2N]
	ωR = [(θR[i%(2N)+1] - θR[i]) / timestep for i = 1:2N]

	# Minimal Coordinates
	nx = minimal_dimension(mechanism)
	X = [zeros(nx) for i = 1:2N]
	for i = 1:2N

		set_minimal_coordinates!(mechanism, get_joint(mechanism, :floating_base), [p_base; 0.0; αtorso; 0.0])
		set_minimal_coordinates!(mechanism, get_joint(mechanism, :back_bkxyz), [0.0, -αtorso, 0.0])
		set_minimal_coordinates!(mechanism, get_joint(mechanism, :l_leg_hpxyz), [0.0, -θL[i][1], 0.0])
		set_minimal_coordinates!(mechanism, get_joint(mechanism, :l_leg_kny), [θL[i][2]])
		set_minimal_coordinates!(mechanism, get_joint(mechanism, :l_leg_akxy), [θL[i][1]-θL[i][2], 0.0])

		set_minimal_coordinates!(mechanism, get_joint(mechanism, :r_leg_hpxyz), [0.0, -θR[i][1], 0.0])
		set_minimal_coordinates!(mechanism, get_joint(mechanism, :r_leg_kny), [θR[i][2]])
		set_minimal_coordinates!(mechanism, get_joint(mechanism, :r_leg_akxy), [θR[i][1]-θR[i][2], 0.0])

		set_minimal_velocities!(mechanism, get_joint(mechanism, :floating_base), [zeros(3); zeros(3)])
		set_minimal_velocities!(mechanism, get_joint(mechanism, :l_leg_hpxyz), [0.0, -ωL[i][1], 0.0])
		set_minimal_velocities!(mechanism, get_joint(mechanism, :l_leg_kny), [ωL[i][2]])
		set_minimal_velocities!(mechanism, get_joint(mechanism, :l_leg_akxy), [ωL[i][1]-ωL[i][2], 0.0])

		set_minimal_velocities!(mechanism, get_joint(mechanism, :r_leg_hpxyz), [0.0, -ωR[i][1], 0.0])
		set_minimal_velocities!(mechanism, get_joint(mechanism, :r_leg_kny), [ωR[i][2]])
		set_minimal_velocities!(mechanism, get_joint(mechanism, :r_leg_akxy), [ωR[i][1]-ωR[i][2], 0.0])
		X[i] .= get_minimal_state(mechanism)
	end

	X = vcat([deepcopy(X) for i = 1:Ncycles]...)
	# adding forward displacement and velocity
	for i = 1:2N*Ncycles
		X[i][1] += (i-1) * 2r / N
		X[i][7] += 2r / N / timestep
	end

	return X
end

function interpolate(parent_vertex, p2, N::Int)
	return [parent_vertex .+ α*(p2 - parent_vertex) for α in range(0, stop=1.0, length=N)]
end

function atlas_jump(mechanism::Mechanism{T}; timestep=0.05, αtorso=0.0, z=0.95, Δx=0.0, x = 0.0,
		Ns=[5,5,5,5,5], foot_high=0.20, base_high=0.20, base_low=-0.20) where T
	N = sum(Ns)
	# Foot position
	fL = [Δx, + 0.1145, 0]
	fH = [Δx, + 0.1145, foot_high]
	# Base position
	bL = [x, 0, z+base_low]
	bM = [x, 0, z]
	bH = [x, 0, z+base_high]

	pfoot = [
		interpolate(fL, fL, Ns[1]);
		interpolate(fL, fL, Ns[2]);
		interpolate(fL, fH, Ns[3]);
		interpolate(fH, fL, Ns[4]);
		interpolate(fL, fL, Ns[5]);
		]
	pbase = [
		interpolate(bM, bL, Ns[1]);
		interpolate(bL, bM, Ns[2]);
		interpolate(bM, bH, Ns[3]);
		interpolate(bH, bL, Ns[4]);
		interpolate(bL, bM, Ns[5]);
		]

	# Leg angles
	θL = [IKatlas(mechanism, pbase[i], pfoot[i], leg=:l) for i=1:N]
	θR = [IKatlas(mechanism, pbase[i], pfoot[i], leg=:r) for i=1:N]

	# adding angular velocities
	ωL = [(θL[i%N+1] - θL[i]) / timestep for i = 1:N]
	ωR = [(θR[i%N+1] - θR[i]) / timestep for i = 1:N]

	# Minimal Coordinates
	nx = minimal_dimension(mechanism)
	X = [zeros(nx) for i = 1:N]
	for i = 1:N
		set_minimal_coordinates!(mechanism, get_joint(mechanism, :floating_base), [pbase[i]; 0.0; αtorso; 0.0])
		set_minimal_coordinates!(mechanism, get_joint(mechanism, :back_bkxyz), [0.0, -αtorso, 0.0])
		set_minimal_coordinates!(mechanism, get_joint(mechanism, :l_leg_hpxyz), [0.0, -θL[i][1], 0.0])
		set_minimal_coordinates!(mechanism, get_joint(mechanism, :l_leg_kny), [θL[i][2]])
		set_minimal_coordinates!(mechanism, get_joint(mechanism, :l_leg_akxy), [θL[i][1]-θL[i][2], 0.0])

		set_minimal_coordinates!(mechanism, get_joint(mechanism, :r_leg_hpxyz), [0.0, -θR[i][1], 0.0])
		set_minimal_coordinates!(mechanism, get_joint(mechanism, :r_leg_kny), [θR[i][2]])
		set_minimal_coordinates!(mechanism, get_joint(mechanism, :r_leg_akxy), [θR[i][1]-θR[i][2], 0.0])

		# set_minimal_velocities!(mechanism, get_joint(mechanism, :floating_base), [zeros(3); zeros(3)])
		# set_minimal_velocities!(mechanism, get_joint(mechanism, :l_leg_hpxyz), [0.0, -ωL[i][1], 0.0])
		# set_minimal_velocities!(mechanism, get_joint(mechanism, :l_leg_kny), [ωL[i][2]])
		# set_minimal_velocities!(mechanism, get_joint(mechanism, :l_leg_akxy), [ωL[i][1]-ωL[i][2], 0.0])
		#
		# set_minimal_velocities!(mechanism, get_joint(mechanism, :r_leg_hpxyz), [0.0, -ωR[i][1], 0.0])
		# set_minimal_velocities!(mechanism, get_joint(mechanism, :r_leg_kny), [ωR[i][2]])
		# set_minimal_velocities!(mechanism, get_joint(mechanism, :r_leg_akxy), [ωR[i][1]-ωR[i][2], 0.0])
		X[i] .= get_minimal_state(mechanism)
	end

	# X = vcat([deepcopy(X) for i = 1:Ncycles]...)
	# # adding forward displacement and velocity
	# for i = 1:2N*Ncycles
	# 	X[i][1] += (i-1) * 2r / N
	# 	X[i][7] += 2r / N / timestep
	# end
	return X
end

################################################################################
# Compute trajectory
################################################################################
# vis= Visualizer()
# open(vis)

# mech = get_mechanism(:atlas, timestep=0.05, model_type=:armless, damper=1000.0)
# initialize!(mech, :atlas, tran=[1,0,0.0], rot=[0,0,0.], αhip=0.5, αknee=1.0)
#
# @elapsed storage = simulate!(mech, 0.05, record=true, verbose=true)
# visualize(mech, storage, vis=vis)
#
# X = atlas_jump(mech; Ns=[7,4,4,8,8], timestep=0.025, Δx=0.0, x=0.0, z=1.10, foot_high=0.20, base_high=+0.20, base_low=-0.20)
# zref = [minimal_to_maximal(mech, x) for x in X]
# storage = generate_storage(mech, zref)
# visualize(mech, storage, vis=vis)

#
# bodies = mech.bodies
# joints = mech.joints
# contacts = collect(mech.contacts)
# getfield.(contacts, :parent_id)
# get_body(mech, 15)
# nx = minimal_dimension(mech)
#
#
# p_base = [0, 0, 0.9385]
# p_base = [0, 0, 0.88]
# p_foot = [0, 0, 0.00]
# θ = [0.4, 0.8]
# IKerror(mech, p_base, p_foot, θ; leg=:r)
# IKatlas(mech, p_base, p_foot; leg=:r)


text:
using Oceananigans

include("utils.jl")

Nx, Ny, Nz = 1, 1, 512
Lx, Ly, Lz = 1, 1, 5000
Nt, Δt = 10000, 60

model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=4e-2, κ=4e-2)

@. model.tracers.T.data[:, :, 1:Int(Nz/2)]   -= 0.01
@. model.tracers.T.data[:, :, Int(Nz/2):end] += 0.01

# Add a NetCDF output writer that saves NetCDF files to the current directory
# "." with a filename prefix of "vertical_profile_" every 200 iterations.
nc_writer = NetCDFOutputWriter(dir=".", prefix="column_model_", frequency=200)
push!(model.output_writers, nc_writer)

time_step!(model, Nt, Δt)

make_vertical_profile_movie(model, nc_writer, "T", Nt, Δt)


text:
export  transFluxPlot

function transFluxPlot( mpSim::MPsim, timeGrid::Vector{T},
    transitions::MP.TransitionType...; showPlot::Bool=true,
    savePlot::Bool=false, filename::AbstractString="",
    timeFactor::Real=12.0 ) where T <: Real

    savePlot |= filename != ""

    # Don't do anything if the plot mustn't be shown or saved.
    if !( showPlot || savePlot )
        return
    end  # if !( showPlot || savePlot )

    # Generate the report.
    fluxReport = transitionFluxReport( mpSim, timeGrid, transitions... )

    if isempty( fluxReport )
        return
    end  # if isempty( fluxReport )

    generateTransFluxPlot( fluxReport, timeFactor, showPlot, savePlot,
        filename )

end  # transFluxPlot( mpSim, timeGrid, transitions, showPlot, savePlot,
     #   filename, timeFactor )

transFluxPlot( mpSim::MPsim, timeRes::Real,
    transitions::MP.TransitionType...; showPlot::Bool=true,
    savePlot::Bool=false, filename::AbstractString="",
    timeFactor::Real=12.0 ) =
    transFluxPlot( mpSim, MP.generateTimeGrid( mpSim, timeRes ),
    transitions..., showPlot=showPlot, savePlot=savePlot,
    filename=filename, timeFactor=timeFactor )

include( "private/transitionflux.jl" )


text:
@testset "Gen&Opt" begin
    isdir("basic") && rm("basic", recursive=true, force=true)
    newsite("basic")

    serve(single=true)
    # ---------------
    @test all(isdir, ("assets", "css", "libs", "pub", "src"))
    @test all(isfile, ("index.html",
                 map(e->joinpath("pub", "menu$e.html"), 1:3)...,
                 map(e->joinpath("css", e), ("basic.css", "judoc.css"))...,
                 )
              )
    # ---------------
    if JuDoc.JD_CAN_MINIFY
        presize1 = stat(joinpath("css", "basic.css")).size
        presize2 = stat("index.html").size
        optimize(prerender=false)
        @test stat(joinpath("css", "basic.css")).size < presize1
        @test stat("index.html").size < presize2
    end
    # ---------------
    # change the prepath
    index = read("index.html", String)
    @test occursin("=\"/css/basic.css", index)
    @test occursin("=\"/css/judoc.css", index)
    @test occursin("=\"/libs/highlight/github.min.css", index)
    @test occursin("=\"/libs/katex/katex.min.css", index)

    optimize(minify=false, prerender=false, prepath="prependme")
    index = read("index.html", String)
    @test occursin("=\"/prependme/css/basic.css", index)
    @test occursin("=\"/prependme/css/judoc.css", index)
    @test occursin("=\"/prependme/libs/highlight/github.min.css", index)
    @test occursin("=\"/prependme/libs/katex/katex.min.css", index)
end

if J.JD_CAN_PRERENDER && J.JD_CAN_HIGHLIGHT
@testset "Prerender" begin
  hs = raw"""
    <!doctype html>

    <meta charset=UTF-8>

Title

Blah

Consider an invertible matrix \(M\) made of blocks \(A\), \(B\), \(C\) and \(D\) with

      \[ M \quad\!\! =\quad\!\! \begin{pmatrix} A & B \\ C & D \end{pmatrix} \]

using Test
      # Woodbury formula
      b = 2
      println("hello $b")

    """
  jskx = J.js_prerender_katex(hs)
  # conversion of `\(M\)` (inline)
  @test occursin("""<math xmlns=\"http://www.w3.org/1998/Math/MathML\"><semantics><mrow><mi>M</mi></mrow>""", jskx)
  # conversion of the equation (display)
  @test occursin("""<math xmlns=\"http://www.w3.org/1998/Math/MathML\"><semantics><mrow><mi>M</mi>""", jskx)
  jshl = J.js_prerender_highlight(hs)
  # conversion of the code
  @test occursin("""
using""", jshl)
  @test occursin(raw""""hello $b"""", jshl)
end
end # if can prerender


text:
using Pkg; Pkg.activate(pwd())
using Revise
using Test, LinearAlgebra, Random, Statistics
using FileIO, JLD2
dir = pwd()
include(joinpath(dir, "src\\ops_utility.jl"))
include(joinpath(dir, "src\\ops_build.jl"))
include(joinpath(dir, "src\\ops_methods.jl"))
include(joinpath(dir, "src\\cfrops_solve.jl"))
include(joinpath(dir, "src\\dcfrops_solve.jl"))
include(joinpath(dir, "src\\optops_solve.jl"))
include(joinpath(dir, "src\\ops_results.jl"))

##########################################################
# build sequence-form reward
# CAUTION: takes 1 hr+ to build sequence-form reward for 10 city default
ncity = 10
g = AirDefenseGame(ncity, 6, 3, 4, 4, 2, 2)
A, An, na_stage = build_action(g)
ni, ni1, ni2, ni_stage = build_info(g)
ns1, ns2, ns1_stage, ns2_stage = build_nseq(g, na_stage, ni_stage)
Σ, seqn, seqactions = build_Σset(g, A, ns1, ns2)
(U, z), uh_time = @timed build_utility_hist(g, A, An)
reward_exp2, runtime2, nalloc, gc_time, misc = @timed build_utility_seq(g, gs, seqn, seqactions, expected = true)

dir = pwd()
fn = joinpath(dir, "data\\vars_opt_rewardfunction_expected_15city_temp.jld2")
# @save fn g U z reward_exp Σ seqn seqactions runtime nalloc
@show runtime
# # load(fn)

gs = GameSet(g)
gos = GameOptSet(gs, reward_exp)
(status, u_ne, r_ne, t_solve), t_total = @timed lp_nash(gos)
lp_best_response(gos, r_ne, fixedplayer = 2)

##########################################################
# Data-biased response and robust best response
g = AirDefenseGame(10,
        reshape([0.03674972063074833,0.7471235124528937,0.8006182686839263,0.1134119509104532,0.3720810993450727,0.6452465195265953,0.9927983579441313,0.8545575129265128,0.29473362663716407,0.9808676508834953,0.6242283949611096,0.6173536676330247,0.4739652533798455,0.6677615345281047,0.6338233451643072,0.8231268712815414,0.011816562466650193,0.048514560425303443,0.7974455234366102,0.22065505430418475], 10, 2),
        [0.47973468500756566,0.1596483573101528,0.7116043145363229,0.969005005082471,0.7139100481289631,0.4747849042596852,0.5066187533113924,0.7074385618615113,0.21009348198668265,0.29520436695319363],
        0.3,
        6,
        3,
        4,
        4,
        2,
        2,
        0.9,
        0.8,
        0.7,
        0.2,
        [1, 3, 4, 5, 8, 9])
gs = GameSet(g)
A, An, na_stage = build_action(g)
ni, ni1, ni2, ni_stage = build_info(g)
ns1, ns2, ns1_stage, ns2_stage = build_nseq(g, na_stage, ni_stage)
Σ, seqn, seqactions = build_Σset(g, A, ns1, ns2)
(U, z), uh_time = @timed build_utility_hist(g, A, An)
(reward_exp, reward_complete), runtime, _, _, _ = @timed build_utility_seq(g, gs, (ns1, ns2), seqn, seqactions, expected = true)

dir = pwd()
fn = joinpath(dir, "data\\vars_opt_rewardfunction_expected_15city_flipped_temp.jld2")
# @save fn g U z reward_exp Σ seqn seqactions runtime
@show runtime
# load(fn)

dir = pwd()
include(joinpath(dir, "src\\dbr_cfrops_solve.jl"))
using Plots; gr()
utemp, rtemp, stemp, σtemp = cfr(1_000, g, gs)
σfix = copy(σtemp)
pdbr = vcat(fill(0.95, gs.ni_stage[3]), fill(0.0, gs.ni_stage[5]), fill(0.85, gs.ni_stage[6]))
fn = "data\\vars_cfr_dbr_attackerfixed_temp.jld2"
# @save joinpath(dir, fn) σfix pdbr

T = 3_000

# NE for g
ucfr, _, _, σcfr = cfr(T, g, gs)
mean(ucfr)
plot(cumsum(ucfr) ./ collect(1:T))

# NE for g using cfr_dbr (PConf = 0 and σfix = arbitrary)
u2, _, _, σ2 = cfr_dbr(T, g, gs, zeros(length(pdbr)), σcfr .* 100)
mean(u2)
all(sum(i) ≈ 1 for i in σ2)

# NE for g using cfr_dbr (PConf = 1 and σfix = σNE)
u3, _, _, _ = cfr_dbr(T, g, gs, ones(length(pdbr)), σcfr)
mean(u3)

# calculate defender's robust CFR BR to σfix at 100% confidence
# this matches lp_best_response (at least sometimes...)
T = 3000
u1, _, _, _ = cfr_dbr(T, g, gs, ones(length(pdbr)), σfix)
mean(u1)
plot(cumsum(u1) ./ collect(1:T))

# calculate defender's BR to σfix
gos = GameOptSet(g, gs, reward_exp)
status_br, u_br, r_br = lp_best_response(gos, real_strat(σfix, gs, gos)[2], fixedplayer = 2)
u_br

status_brne, u_brne, r_brne = lp_best_response(gos, real_strat(σcfr, gs, gos)[2], fixedplayer = 2)
u_brne

st_ne, u_ne, r_ne, t_ne = lp_nash(gos, 5 * 60)
u_ne

status_brne, u_brne, r_brne = lp_best_response(gos, r_ne, fixedplayer = 2)
u_brne

T = 20_000
ucfr, _, _, σcfr = cfr(T, g, gs)
mean(ucfr)
status_brne, u_brne_1, r_brne = lp_best_response(gos, real_strat(σcfr, gs, gos)[2], fixedplayer = 2)
u_brne_1
status_brne, u_brne_2, r_brne = lp_best_response(gos, real_strat(σcfr, gs, gos)[1], fixedplayer = 1)
u_brne_2

# calculate defender's robust CFR BR to σfix at pdbr% confidence
u, _, _, _ = cfr_dbr(T, g, gs, pdbr, σfix)
mean(u)
plot(cumsum(u) ./ collect(1:T))

# calculate defender's exact robust BR to σfix
# TBD - need to code up robust BR...from overleaf

##########################################################
# deprecate once GameOptSet constructor works
# build AirDefenseGame, GameSet, and GameOptSet for both optimization and CFR
g = AirDefenseGame()
gs = GameSet(g)
A, An, na_stage = build_action(g)
ni, ni1, ni2, ni_stage = build_info(g)
ns1, ns2, ns1_stage, ns2_stage = build_nseq(g, na_stage, ni_stage)
Σ, seqn, seqactions = build_Σset(g, A, ns1, ns2)
Iset, (I1set, I2set), (seqI1, seqI2), (nextseqI1, nextseqI2) = build_Iset(g, A, na_stage, ns1_stage, ns2_stage, ni1, ni2, ns1, ns2)
# load expensive reward function
fn = joinpath(dir, "data\\vars_opt_rewardfunction.jld2")
@load fn reward
fn2 = joinpath(dir, "data\\vars_opt_rewardfunction_expected.jld2")
@load fn2 reward_exp
# Next line defaults to expected reward
gos = GameOptSet(ns1, ns2, ni1, ni2, ns1_stage, ns2_stage, reward_exp,
    gs.na_stage, gs.sensors, seqI1, seqI2, nextseqI1, nextseqI2)

##########################################################
# CFR
Random.seed!(579843)
# g = AirDefenseGame(0.95)
# gamesize(AirDefenseGame(15,8,5,8,8,3,3))
# g = AirDefenseGame(13,6,4,6,6,3,3)
# gamesize(g)
# gamesize(AirDefenseGame())
# g = AirDefenseGame(10,8,4,8,8,3,3)
g = AirDefenseGame()

g = AirDefenseGame(10,
        reshape([0.03674972063074833,0.7471235124528937,0.8006182686839263,0.1134119509104532,0.3720810993450727,0.6452465195265953,0.9927983579441313,0.8545575129265128,0.29473362663716407,0.9808676508834953,0.6242283949611096,0.6173536676330247,0.4739652533798455,0.6677615345281047,0.6338233451643072,0.8231268712815414,0.011816562466650193,0.048514560425303443,0.7974455234366102,0.22065505430418475], 10, 2),
        [0.47973468500756566,0.1596483573101528,0.7116043145363229,0.969005005082471,0.7139100481289631,0.4747849042596852,0.5066187533113924,0.7074385618615113,0.21009348198668265,0.29520436695319363],
        0.3,
        6,
        3,
        4,
        4,
        2,
        2,
        0.9,
        0.8,
        0.7,
        0.2,
        [1, 3, 4, 5, 8, 9])
gs = GameSet(g)
# gos = GameOptSet(gs, reward_exp)
T = 5_000
@time u, r, s, σ = cfr(T, g, gs) # 45 ms => 45 sec per 1k iterations
mean(u)

@time u1, r, s, σ, σ1, converged, iter = dcfr_full(g, gs, iterlimit = T)
mean(u1)

# psc = 0.55, pac = 0.8, T = 5_000: utility = -1.2374797
# psc = 0.95, pac = 0.8, T = 5_000: utility = -1.2249

##########################################################
# CFR Stopping Rule and Step size
using Plots; gr()
using Distributed
nprocs() == 1 && addprocs()
@everywhere using Pkg
@everywhere Pkg.activate(pwd())
@everywhere(using Revise, BenchmarkTools, Printf, ProgressMeter, FileIO, JLD2)
@everywhere(using Random, LinearAlgebra, Statistics, Distances, IterTools, DataFrames, DataFramesMeta)
@everywhere(using StaticArrays, SparseArrays)
@everywhere(using JuMP, Clp)
@everywhere(using Clp: ClpCInterface)
@everywhere dir = pwd()
@everywhere include(joinpath(dir, "src\\ops_utility.jl"))
@everywhere include(joinpath(dir, "src\\ops_build.jl"))
@everywhere include(joinpath(dir, "src\\ops_methods.jl"))
@everywhere include(joinpath(dir, "src\\cfrops_solve.jl"))
@everywhere include(joinpath(dir, "src\\dcfrops_solve.jl"))
@everywhere include(joinpath(dir, "src\\optops_solve.jl"))
@everywhere include(joinpath(dir, "src\\ops_results.jl"))

Random.seed!(579843)
g = AirDefenseGame(10,
        reshape([0.03674972063074833,0.7471235124528937,0.8006182686839263,0.1134119509104532,0.3720810993450727,0.6452465195265953,0.9927983579441313,0.8545575129265128,0.29473362663716407,0.9808676508834953,0.6242283949611096,0.6173536676330247,0.4739652533798455,0.6677615345281047,0.6338233451643072,0.8231268712815414,0.011816562466650193,0.048514560425303443,0.7974455234366102,0.22065505430418475], 10, 2),
        [0.47973468500756566,0.1596483573101528,0.7116043145363229,0.969005005082471,0.7139100481289631,0.4747849042596852,0.5066187533113924,0.7074385618615113,0.21009348198668265,0.29520436695319363],
        0.3,
        6,
        3,
        4,
        4,
        2,
        2,
        0.9,
        0.8,
        0.7,
        0.2,
        [1, 3, 4, 5, 8, 9])
gs = GameSet(g)
# @load joinpath(dir, "data\\vars_opt_rewardfunction_expected.jld2") reward_exp
# gos = GameOptSet(gs, reward_exp)

# @time u, r, s, σ, σ1, converged = cfr_stoprule(5, g, gs, tolerance = 5e-5)
# mean(u)
# um = cumsum(u) ./ collect(1:length(u))
# status, u_br, r_br = lp_best_response(gos, real_strat(σ, gs, gos)[1], fixedplayer = 1)
# plot(1:length(um), um, xlim = (0,1000))
# ud, rd, sd, σd, σ1d, converged_d = dcfr_full(g, gs, timelimit = 2, tol = 5e-5,
#                                                 α = 1.5, β = 0.5, γ = 2.0, discounted = true)
# non-discounted: mean(ud) = -1.245077
# @show mean(ud)
# umd = cumsum(ud) ./ collect(1:length(ud))
# plot!(1:length(umd), umd)
# plot(1:length(ud), ud)
# statusd, u_brd, r_brd = lp_best_response(gos, real_strat(σd, gs, gos)[1], fixedplayer = 1)
# rerr(-u_brd, -1.2595387377)

#
# alphas = [0.0, 0.0, 0.0, 0.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 1.5]
# betas = [-8.0, -8.0, 8.0, 8.0, -8.0, -8.0, 8.0, 8.0, 8.0, -8.0, 0.5]
# gammas = [0.0, 8.0, 0.0, 8.0, 0.0, 8.0, 0.0, 8.0, 0.0, 2.0, 2.0]
reps = 50
alphas = repeat([8.0, 8.0, 1.5], inner = reps)
betas = repeat([8.0, -8.0, 0.5], inner = reps)
gammas = repeat([0.0, 2.0, 2.0], inner = reps)
doe_discount = DataFrame(alpha = alphas, beta = betas, gamma = gammas)
results = pmap((a,b,c) -> dcfr_full(g, gs, timelimit = 10, tol = 5e-5,
                            α = a, β = b, γ = c, discounted = true), alphas, betas, gammas)
nsteps = 12_000
ums = [cumsum(r[1]) ./ collect(1:length(r[1])) for r in results]
ums_means = [[mean(ums[i][j] for i = (k * reps + 1):((k + 1) * reps)) for j in 1:nsteps] for k in 0:2]
ums_stds = [[std(ums[i][j] for i = (k * reps + 1):((k + 1) * reps)) for j in 1:nsteps] for k in 0:2]
# us = [r[1] for r in results]
# us_means = [[mean(us[i][j] for i = (k * reps + 1):((k + 1) * reps)) for j in 1:nsteps] for k in 0:2]
# us_stds = [[std(us[i][j] for i = (k * reps + 1):((k + 1) * reps)) for j in 1:nsteps] for k in 0:2]
pcolors = [:steelblue, :tomato, :mediumseagreen]
pstyles = [:solid, :dash, :dashdot]
fig = plot(xlabel = "Iterations", ylabel = "Defender Utility",
    xlims = (0, nsteps), ylims = (-1.27, -1.2))
# fig = plot(xlabel = "Iterations", ylabel = "Defender Utility",
#     xlims = (0, 2_000), ylims = (-1.5, -0.5))
for i in 1:3
    plot!(fig, 1:nsteps, ums_means[i], ribbon = 2.009 * ums_stds[i] ./ sqrt(reps), fillalpha = 0.2, linecolor = pcolors[i],
        label = string(alphas[i*reps], ", ", betas[i*reps], ", ", gammas[i*reps]),
        linestyle = pstyles[i])
    # plot!(fig, 1:nsteps, us_means[i], alpha = 0.5, seriestype = :line,
    #     color = pcolors[i])
    # plot!(fig, 1:nsteps, ums_means[i] .- ums_stds[i], linecolor = pcolors[i], linestyle = :dot)
    # plot!(fig, 1:nsteps, ums_means[i] .+ ums_stds[i], linecolor = pcolors[i], linestyle = :dot)
end
plot!(fig, 1:nsteps, fill(-1.2595387377, nsteps), label = "Nash Equilibrium",
    linestyle = :dot, color = :black)

dir = pwd()
fn = joinpath(dir, "data\\plot_stepsize_v4_temp.pdf")
# savefig(fn)
fn = joinpath(dir, "data\\vars_stepsize_fig_temp.jld2")
# @save fn results reps g doe_discount

##########################################################
# CFR DOE
using Pkg; Pkg.activate(pwd())
using Distributed, SharedArrays, Statistics, Dates
using CSV
# using Plots; gr()
nprocs() == 1 && addprocs()
@everywhere using Pkg
@everywhere Pkg.activate(pwd())
@everywhere(using Revise, BenchmarkTools, Printf, ProgressMeter, FileIO, JLD2)
@everywhere(using Random, LinearAlgebra, Statistics, Distances, IterTools, DataFrames, DataFramesMeta)
@everywhere(using StaticArrays, SparseArrays)
@everywhere(using JuMP, Clp)
@everywhere(using Clp: ClpCInterface)
@everywhere dir = pwd()
@everywhere include(joinpath(dir, "src\\ops_utility.jl"))
@everywhere include(joinpath(dir, "src\\ops_build.jl"))
@everywhere include(joinpath(dir, "src\\ops_methods.jl"))
@everywhere include(joinpath(dir, "src\\cfrops_solve.jl"))
@everywhere include(joinpath(dir, "src\\dcfrops_solve.jl"))
@everywhere include(joinpath(dir, "src\\optops_solve.jl"))
@everywhere include(joinpath(dir, "src\\ops_results.jl"))

# see scratch for screening design and RSM design
@everywhere function build_design(fn::String)
    factor_names = ["pp", "psc", "pdc",
                    "alpha", "beta", "gamma"]
    col_types = [Float64, Float64, Float64,
                    Float64, Float64, Float64]
    const_names = ["radius", "pac",
                    "ncity", "ndp", "nap",
                    "ndpdc", "ndpac", "ndc", "nac"]
    const_types = [Float64, Float64,
                    Int64, Int64, Int64,
                    Int64, Int64, Int64, Int64]
    const_values = [0.3, 0.2,
                    13, 6, 4,
                    6, 6, 3, 3]
    design = CSV.read(joinpath(dir, fn), types = col_types)
    rename!(design, old => new for (old, new) = zip(names(design), Symbol.(factor_names)))
    l, w = size(design)
    for i in 1:2
        design[Symbol(const_names[i])] = fill(const_values[i], l)
    end
    for i in 3:length(const_names)
        design[Symbol(const_names[i])] = fill(Int64(const_values[i]), l)
    end
    return design
end

@everywhere function build_constants(design)
    nruns = size(design, 1)
    Xmax = rand(maximum(design.ncity), 2)
    vmax = rand(maximum(design.ncity))
    Xs = [Xmax[1:nc, :] for nc in design.ncity]
    vs = [vmax[1:nc] for nc in design.ncity]
    iadstable = Vector{Any}[]
    for nc in unique(design.ncity), nndp in unique(design.ndp)
        push!(iadstable, [nc, nndp, sort!(sample(1:nc, nndp, replace = false))])
    end
    iadss = [iadstable[findfirst(x -> x[1] == design.ncity[i] && x[2] == design.ndp[i],
                                    iadstable)][3] for i in 1:nruns]
    return Xs, vs, iadss
end

@everywhere function build_doegames(design, Xs, vs, iadss)
    nruns = size(design, 1)
    gvec = Array{AirDefenseGame}(undef, nruns)
    gsvec = Array{GameSet}(undef, nruns)
    nnodevec = Array{String}(undef, nruns)
    for i in 1:nruns
        ncity = design.ncity[i][1]
        radius = design.radius[i][1]
        ndp = design.ndp[i][1]
        nap = design.nap[i][1]
        ndpdc = design.ndpdc[i][1]
        ndpac = design.ndpac[i][1]
        ndc = design.ndc[i][1]
        nac = design.nac[i][1]
        pp = design.pp[i][1]
        psc = design.psc[i][1]
        pdc = design.pdc[i][1]
        pac = design.pac[i][1]
        a = design.alpha[i][1]
        b = design.beta[i][1]
        c = design.gamma[i][1]
        gvec[i] = AirDefenseGame(ncity, Xs[i], vs[i], radius, ndp, nap, ndpdc, ndpac, ndc, nac,
                            pp, psc, pdc, pac, iadss[i])
        gsvec[i] = GameSet(gvec[i])
        nnodevec[i] = gamesize(gvec[i])[4]
    end
    return gvec, gsvec, nnodevec
end

@everywhere function calc_run(df_row, g, gs, runnum)
    a = df_row.alpha[1]
    b = df_row.beta[1]
    c = df_row.gamma[1]
    println("----------- Start Run $runnum -----------")
    (u, _, _, _, _, converged, iter), runtime = @timed dcfr_full(g, gs,
                iterlimit = 100_000, timelimit = 240, tol = 5e-5,
                α = a, β = b, γ = c, discounted = true)
    println("------##### End Run $runnum #####------")
    return mean(u), converged, iter, runtime
end

# design = build_design("data\\design_frac_207.csv")
design = build_design("data\\design_frac_207_augment.csv")
nruns = size(design, 1)
Xs, vs, iadss = build_constants(design)
gvec, gsvec, nnodevec = build_doegames(design, Xs, vs, iadss)
parse.(Float64, nnodevec) |> findmax

results_doe, t_doe = @timed pmap(x -> calc_run(design[x, :], gvec[x], gsvec[x], x), 1:nruns)
fn = joinpath(dir, "data\\results_doe_frac207augment_temp.jld2")
# @save fn results_doe, t_doe, gvec

dfr = hcat(design, DataFrame(utility = [ri[1] for ri in results_doe],
                converged = [ri[2] for ri in results_doe],
                iterations = [ri[3] for ri in results_doe],
                time = [ri[4] for ri in results_doe]))

# CSV.write(joinpath(dir, "data\\results_doe_frac207augment_temp.csv"), dfr)

##########################################################
# CFR Large

using Pkg; Pkg.activate(pwd())
using Revise
using Distributed, SharedArrays, Statistics, Dates
using CSV
# using Plots; gr()
nprocs() == 1 && addprocs()
@everywhere using Pkg
@everywhere Pkg.activate(pwd())
@everywhere(using Revise, BenchmarkTools, Printf, ProgressMeter, FileIO, JLD2)
@everywhere(using Random, LinearAlgebra, StatsBase, Distances, IterTools, DataFrames, DataFramesMeta)
@everywhere(using StaticArrays, SparseArrays)
@everywhere(using JuMP, Clp)
@everywhere(using Clp: ClpCInterface)
@everywhere dir = pwd()
@everywhere include(joinpath(dir, "src\\ops_utility.jl"))
@everywhere include(joinpath(dir, "src\\ops_build.jl"))
@everywhere include(joinpath(dir, "src\\ops_methods.jl"))
@everywhere include(joinpath(dir, "src\\cfrops_solve.jl"))
@everywhere include(joinpath(dir, "src\\dcfrops_solve.jl"))
@everywhere include(joinpath(dir, "src\\optops_solve.jl"))
@everywhere include(joinpath(dir, "src\\ops_results.jl"))

# collect and load reward function data
@everywhere function collect_sizedata_seqreward(ncity, tlim)
    g = AirDefenseGame(ncity, 6, 3, 4, 4, 2, 2)
    gs = GameSet(g)
    A, An, na_stage = build_action(g)
    ni, ni1, ni2, ni_stage = build_info(g)
    ns1, ns2, ns1_stage, ns2_stage = build_nseq(g, na_stage, ni_stage)
    Σ, seqn, seqactions = build_Σset(g, A, ns1, ns2)
    U, z = build_utility_hist(g, A, An)
    (reward_exp, completed), runtime = @timed build_utility_seq(g, gs, (ns1, ns2), seqn, seqactions,
                                                    expected = true, timelimit = tlim)
    return g, U, z, reward_exp, completed, Σ, seqn, seqactions, runtime
end
Random.seed!(579843)
ncities = [7, 8, 9, 10, 11, 12, 13, 14, 15] # number of cities (all other game params fixed)
n = length(ncities)
println("Start Time: $(now())")
results = pmap((a, b) -> collect_sizedata_seqreward(a, b), ncities, fill(12 * 60, length(ncities)))
dirlocal = "C:\\Users\\AKeith\\JuliaProjectsLocal\\StrategyGamesLocal"
fnlocal = joinpath(dirlocal, "data\\vars_opt_rewardfunction_expected_7to15city_temp.jld2")
# # @save fnlocal results
# fn = joinpath(dir, "data\\vars_opt_rewardfunction_expected_7to15city_temp.jld2") # 12 hr time limit
# # @save fn results
# @load joinpath(dirlocal, "data\\vars_opt_rewardfunction_expected_7to15city.jld2") results # approx 5 min
gvec = [ri[1] for ri in results]
gsvec = GameSet.(gvec)
gosvec = [GameOptSet(gvec[i], gsvec[i], rexps[i]) for i in 1:n]
rexps = [ri[4] for ri in results]
gsizes = [gamesize(gi)[4] for gi in gvec]
seqsizes = [size(ri[7][1])[1] * size(ri[7][2])[1] for ri in results]
statuses = [ri[5] for ri in results]
runtimes = [ri[end] for ri in results] ./ (60 * 60) # in hours

# collect and load LP
unes = SharedArray{Float64}(n)
tnes = SharedArray{Float64}(n)
statusnes = SharedArray{Bool}(n)
println("LP NE start time for $n city sizes: $(now())")
@distributed for i in 1:n
    g = gvec[i]
    gs = gsvec[i]
    gos = gosvec[i]
    maxseconds = 10 * 60
    (status_ne, u_ne, r_ne, t_solve), t_ne = @timed lp_nash(gos, maxseconds)
    unes[i] = u_ne
    tnes[i] = t_ne
    (status_ne == :Optimal) && (statusnes[i] = true)
end
fn = joinpath(dir, "results_size_lpne_7to15city_10min_temp.jld2")
# @save fn unes tnes statusnes
# @load joinpath(dir, "results_size_lpne_15city_10min.jld2")  unes tnes statusnes

# collect and load CFR
@everywhere function collect_sizedata_cfr(g::AirDefenseGame, gs::GameSet)
    (u, r, s, σ, σs, converged), runtime = @timed dcfr_full(g, gs, timelimit = 10, tol = 5e-5, α = 1.5, β = 0.0, γ = 2.0, discounted = false)
    return u, r, s, σ, σs, converged, runtime
end
Random.seed!(748685)
cfr_params = Dict("timelimit" => 10, "tol" => 5e-5, "α" => 1.5, "β" => 0.0, "γ" => 2.0, "discounted" => false)
println("Start Time: $(now())")
results_cfr = pmap((a, b) -> collect_sizedata_cfr(a, b), gvec, gsvec)
fn = joinpath(dir, "data\\results_size_cfr_7to15city_10min_temp.jld2")
# @save fn results_cfr gvec cfr_params
# @load joinpath(dir, "data\\results_size_cfr_15city_5min.jld2") results ncities T
ucfrs = [ri[1] for ri in results_cfr]
umcfrs = [cumsum(ri[1]) ./ collect(1:length(ri[1])) for ri in results_cfr]
umeancfrs = [mean(ri[1]) for ri in results_cfr]
tcfrs = [ri[7] for ri in results_cfr]
σcfrs = [ri[4] for ri in results_cfr]
convergecfrs = [ri[6] for ri in results_cfr]

# collect and load best responses
@everywhere fbr(g, gs, gos, σ) = lp_best_response(gos, real_strat(σ, gs, gos)[1], fixedplayer = 1)
results_br = pmap(fbr, gvec, gsvec, gosvec, σcfrs) # status_br, u_br, r_br
fn = joinpath(dir, "data\\results_size_br_7to15city_10min_temp.jld2")
# @save fn results_br gvec
ubrs = [-1 * ri[2] for ri in results_br]
uexploits = rerr.(ubrs, unes)

# make size summary data frame
df_size = DataFrame(ncity = ncities, nnodes = gsizes, nseqcombos = seqsizes,
                        rewardtime = runtimes, rewardcomplete = statuses,
                        lptime = tnes, lputility = unes, lpstatus = statusnes,
                        cfrtime = tcfrs, cfrutility = umeancfrs, cfrconverged = convergecfrs,
                        brutility = ubrs, exploitability = uexploits)
# CSV.write(joinpath(dir, "data\\results_size_7to15_temp.csv"), df_size)

# plots
using DataFrames, DataFramesMeta

dfr = DataFrame(ncity = Int64[], iter = Int64[], utility = Float64[], um = Float64[], converged = Bool[], runtime = Float64[])
for i in eachindex(us)
    for j in eachindex(us[i])
        push!(dfr, [ncities[i], j, us[i][j], ums[i][j]])
    end
end

dfx = @linq dfr |>
    where(:ncity .== 10) |>
    select(:iter, :um)

fig = plot(dfx[:iter], dfx[:um], xlabel = "Iterations", ylabel = "Defender Utility")
for city in ncities
    dfi = @linq dfr |>
        where(:ncity .== city) |>
        select(:iter, :um)
    plot!(fig, dfi[:iter], dfi[:um], label = "|M|=$city")
end
fig

######################################################################################
# Analyze default scenario
# Output: Defender utility by iteration plot, relative exploitability by iteration plot
include(joinpath(dir, "src\\optops_solve.jl"))
using Plots; gr()
Random.seed!(579843)
g = AirDefenseGame()
gs = GameSet(g)
@load joinpath(dir, "data\\vars_opt_rewardfunction_expected.jld2") reward_exp
gos = GameOptSet(g, gs, reward_exp)
T = 1_000
@time u, r, s, σ, u_br = cfr_exploit(T, g, gs, gos)
(status, u_ne, r1_ne, r2_ne, t_solve), t_total = @timed lp_nash(gos)
println("Status: $status, Obj Value: $u_ne, Solve Time: $t_solve, Build Time: $(t_total - t_solve)")

@show u_ne # -1.2595387377
um = cumsum(u) ./ collect(1:T) # incremental utility
u_rerr = rerr.(u_br, u_ne) # relative error calc of br to ne

# title = "Defensive utility raw"
plot(1:T, um, legend = false, ylims = (-1.26, -1.2),
    ylabel = "Defender Utility", xlabel = "Iterations")

# title = "CFR approximate relative error"
plot(1:T, [rerr.(um, u_true)] * 100, legend = false, ylims = (0.0,15),
    ylabel = "Relative Error (%)", xlabel = "Iterations")

# title = "CFR exploitability"
plot(T ÷ 10:T ÷ 10:T, u_rerr * 100, legend = false, ylims = (0.0,2),
    seriestype = :line, markerstrokealpha = 0,
    ylabel = "Exploitability (% Nash Equilibrium Utility)", xlabel = "Iterations")
plot!(T ÷ 10:T ÷ 10:T, u_rerr * 100, seriestype = :scatter, markerstrokealpha = 0,
    markercolor = :steelblue)
plot!(1:T, fill(0.5, T), linecolor = :red)
annotate!(T, 0.6, text("0.5% Exploitability", 10, :red, :right))
dir = pwd()
fn = joinpath(dir, "data\\plot_default_temp.pdf")
# savefig(fn)

function h_mle(σ, g, gs; chance_ind = [1,1,0,1,0,0])
    h = SVector(0, 0, 0, 0, 0, 0)
    depth = 1
    for i = 1:6
        if getplayer(depth) == 3
            h = setindex(h, chance_ind[depth], depth)
            depth += 1
        else
            maxval, maxind = findmax(σ[infoset(h, depth, g, gs)])
            h = setindex(h, maxind, depth)
            depth += 1
        end
    end
    h
end

function getutility(h, U, z)
    U[findfirst(x -> x == h, z)]
end

g = AirDefenseGame()
gs = GameSet(g)
U, z = build_utility_hist(g, gs.A, gs.An)
h1 = h_mle(σ, g, gs, chance_ind = [1,1,0,4,0,0])
h3 = SVector(1, 1, 1, 4, 1, 1)
dp3 = getlaydown(h3, g, coverage(g), gs.A)
plot_laydown(dp3)
dp1 = getlaydown(h1, g, coverage(g), gs.A)
plot_laydown(dp1)
dir = pwd()
fn = joinpath(dir, "data\\plot_default_[1,1,1,1,6,39].pdf")
savefig(fn)
h2 = h_mle(σ, g, gs, chance_ind = [1,1,0,4,0,0])
dp2 = getlaydown(h2, g, coverage(g), gs.A)
plot_laydown(dp2)
getutility(h1, U, z)
getutility(h2, U, z) # this should be higher since we have cyber sa

######################################################
# Distributed Pdetect investigation

using Distributed, Statistics, Plots
gr()
addprocs(10); nprocs()
@everywhere using Pkg
@everywhere Pkg.activate(pwd())
@everywhere(using Revise, BenchmarkTools, Printf, ProgressMeter, FileIO, JLD2)
@everywhere(using Random, LinearAlgebra, StatsBase, Distances, IterTools, DataFrames, DataFramesMeta)
@everywhere(using StaticArrays, SparseArrays)
@everywhere(using JuMP, Clp)
@everywhere(using Clp: ClpCInterface)
@everywhere dir = pwd()
@everywhere include(joinpath(dir, "src\\ops_utility.jl"))
@everywhere include(joinpath(dir, "src\\ops_build.jl"))
@everywhere include(joinpath(dir, "src\\ops_methods.jl"))
@everywhere include(joinpath(dir, "src\\cfrops_solve.jl"))
@everywhere include(joinpath(dir, "src\\dcfrops_solve.jl"))
@everywhere include(joinpath(dir, "src\\optops_solve.jl"))
@everywhere include(joinpath(dir, "src\\ops_results.jl"))

@everywhere function sensor_eval(p::Float64)
    Random.seed!(579843)
    g = AirDefenseGame(p)
    gs = GameSet(g)
    (u, r, s, σ, σs, converged), runtime = @timed dcfr_full(g, gs, timelimit = 10,
        tol = 5e-5, discounted = false)
    return u, σ, converged, runtime
end

ps = 0.5:0.05:1.0
results = pmap(sensor_eval, ps)
us = [mean(r[1]) for r in results]
plot(ps, us, xlabel = "Probability of Detection", ylabel = "Defender Utility",
    legend = false, seriestype = :line)
plot!(ps, us, seriestype = :scatter, markerstrokealpha = 0, markercolor = :steelblue)
dir = pwd()
fn = joinpath(dir, "data\\plot_pdetect_stoprule_temp.pdf")
# savefig(fn)

rmprocs(getindex(procs(), procs() .> 1))

gtoy = AirDefenseGame(2, [1.0 1;2 1], [1.0, 1], 0.2, 2, 1, 2, 2, 1, 1, 0.9, 0.99, 0.7, 0.8, [1,2])
gstoy = GameSet(gtoy)
T = 15_000
(utoy, r, s, σtoy), t_cfr = @timed cfr(T, gtoy, gstoy)
mean(utoy)
c1 = [chance(4, i, gtoy, gstoy) for i in 1:2]
s1 = σtoy[infoset(SVector(1,1,1,1,0,0), 5, gtoy, gstoy)]
s2 = σtoy[infoset(SVector(1,1,1,2,0,0), 5, gtoy, gstoy)]
plot_laydown(getlaydown(SVector(1,1,1,1,2,1), gtoy, coverage(gtoy), gstoy.A))
leafutility(SVector(1,1,1,1,1,1), 1, 1.0, 1.0, gtoy, gstoy)
leafutility(SVector(1,1,1,1,2,1), 1, 1.0, 1.0, gtoy, gstoy)
sh1 = h_mle(σtoy, gtoy, gstoy, chance_ind = [1,1,0,1,0,0])
sh2 = h_mle(σtoy, gtoy, gstoy, chance_ind = [1,1,0,2,0,0])

gtoy = AirDefenseGame(2, [1.0 1;2 1], [1.0, 1], 0.2, 2, 1, 2, 2, 1, 1, 0.9, 0.5, 0.7, 0.8, [1,2])
gstoy = GameSet(gtoy)
T = 15_000
(utoy, r, s, σtoy), t_cfr = @timed cfr(T, gtoy, gstoy)
mean(utoy)
c1 = [chance(4, i, gtoy, gstoy) for i in 1:2]
s1 = σtoy[infoset(SVector(1,1,1,1,0,0), 5, gtoy, gstoy)]
s2 = σtoy[infoset(SVector(1,1,1,2,0,0), 5, gtoy, gstoy)]
plot_laydown(getlaydown(SVector(1,1,1,1,2,1), gtoy, coverage(gtoy), gstoy.A))
leafutility(SVector(1,1,1,1,1,1), 1, 1.0, 1.0, gtoy, gstoy)
leafutility(SVector(1,1,1,1,2,1), 1, 1.0, 1.0, gtoy, gstoy)
sh1 = h_mle(σtoy, gtoy, gstoy, chance_ind = [1,1,0,1,0,0])
sh2 = h_mle(σtoy, gtoy, gstoy, chance_ind = [1,1,0,2,0,0])

g1 = AirDefenseGame(0.5)
gs1 = GameSet(g1)
@show c1 = [chance(4, i, g1, gs1) for i in 1:4]
leafutility(SVector(1,1,1,4,1,1), 1, 1.0, 1.0, g1, gs1)
leafutility(SVector(1,1,1,4,6,1), 1, 1.0, 1.0, g1, gs1)
T = 15_000
(u1, r, s, σ1), t_cfr = @timed cfr(T, g1, gs1)
mean(u1)
s11 = σ1[infoset(SVector(1,1,1,1,0,0), 5, g1, gs9)]
s12 = σ1[infoset(SVector(1,1,1,4,0,0), 5, g1, gs9)]
s1h1 = h_mle(σ1, g1, gs1, chance_ind = [1,1,0,1,0,0])
s1h2 = h_mle(σ1, g1, gs1, chance_ind = [1,1,0,4,0,0])

g9 = AirDefenseGame(0.99)
gs9 = GameSet(g9)
@show c9 = [chance(4, i, g9, gs9) for i in 1:4]
leafutility(SVector(1,1,1,4,1,1), 1, 1.0, 1.0, g9, gs9)
leafutility(SVector(1,1,1,4,6,1), 1, 1.0, 1.0, g9, gs9)
T = 15_000
(u9, r, s, σ9), t_cfr = @timed cfr(T, g9, gs9)
mean(u9)
s9h1 = h_mle(σ9, g9, gs9, chance_ind = [1,1,0,1,0,0])
s9h2 = h_mle(σ9, g9, gs9, chance_ind = [1,1,0,4,0,0])

sensor_eval(0.1, 10)
sensor_eval(0.9, 10)

function sensor_eval_toy(p::Float64, T::Int64)
    Random.seed!(579843)
    g = AirDefenseGame(2, [1.0 1;2 1], [1.0, 1], 0.2, 2, 1, 2, 2, 1, 1, 0.9, p, 0.7, 0.8, [1,2])
    @show g.v
    gs = GameSet(g)
    (u, r, s, σ), t_cfr = @timed cfr(T, g, gs)
    return u, r, s, σ, t_cfr
end

ps = 0.55:0.05:0.95
evals = [mean(sensor_eval_toy(p, 15000)[1]) for p in ps]

T = 15_000
results = pmap((a, b) -> sensor_eval(a, b), ps, fill(T, length(ps)))
us = [mean(r[1]) for r in results]
plot(ps, us, xlabel = "Probability of Detection", ylabel = "Defender Utility")

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense