text:
# This file was generated, do not modify it.

using MLJ, StatsBase, Random, PyPlot, CategoricalArrays, PrettyPrinting, DataFrames
X, y = @load_crabs
X = DataFrame(X)
@show size(X)
@show y[1:3]
first(X, 3) |> pretty

levels(y) |> pprint

train, test = partition(eachindex(y), 0.70, shuffle=true, rng=52)
@load XGBoostClassifier
xgb_model = XGBoostClassifier()

countmap(y[train]) |> pprint

xgb  = XGBoostClassifier()
xgbm = machine(xgb, X, y)

r = range(xgb, :num_round, lower=10, upper=500)
curve = learning_curve!(xgbm, resampling=CV(),
                        range=r, resolution=25,
                        measure=cross_entropy)

figure(figsize=(8,6))
plot(curve.parameter_values, curve.measurements)
xlabel("Number of rounds", fontsize=14)
ylabel("Cross entropy", fontsize=14)
xticks([10, 100, 250, 500], fontsize=12)
yticks(0.8:0.05:1, fontsize=12)

savefig("assets/literate/EX-crabs-xgb-curve1.svg") # hide

xgb.num_round = 100;

r1 = range(xgb, :max_depth, lower=3, upper=10)
r2 = range(xgb, :min_child_weight, lower=0, upper=5)

tm = TunedModel(model=xgb, tuning=Grid(resolution=8),
                resampling=CV(rng=11), ranges=[r1,r2],
                measure=cross_entropy)
mtm = machine(tm, X, y)
fit!(mtm, rows=train)

r = report(mtm)

md = r.parameter_values[:,1]
mcw = r.parameter_values[:,2]

figure(figsize=(8,6))
tricontourf(md, mcw, r.measurements)

xlabel("Maximum tree depth", fontsize=14)
ylabel("Minimum child weight", fontsize=14)
xticks(3:2:10, fontsize=12)
yticks(fontsize=12)

savefig("assets/literate/EX-crabs-xgb-heatmap.svg") # hide

xgb = fitted_params(mtm).best_model
@show xgb.max_depth
@show xgb.min_child_weight

xgbm = machine(xgb, X, y)
r = range(xgb, :gamma, lower=0, upper=10)
curve = learning_curve!(xgbm, resampling=CV(),
                        range=r, resolution=30,
                        measure=cross_entropy);

@show round(minimum(curve.measurements), sigdigits=3)
@show round(maximum(curve.measurements), sigdigits=3)

r1 = range(xgb, :subsample, lower=0.6, upper=1.0)
r2 = range(xgb, :colsample_bytree, lower=0.6, upper=1.0)
tm = TunedModel(model=xgb, tuning=Grid(resolution=8),
                resampling=CV(rng=234), ranges=[r1,r2],
                measure=cross_entropy)
mtm = machine(tm, X, y)
fit!(mtm, rows=train)

r = report(mtm)
ss = r.parameter_values[:,1]
cbt = r.parameter_values[:,2]

figure(figsize=(8,6))
tricontourf(ss, cbt, r.measurements)

xlabel("Sub sample", fontsize=14)
ylabel("Col sample by tree", fontsize=14)
xticks(fontsize=12)
yticks(fontsize=12)

savefig("assets/literate/EX-crabs-xgb-heatmap2.svg") # hide

xgb = fitted_params(mtm).best_model
@show xgb.subsample
@show xgb.colsample_bytree

ŷ = predict_mode(mtm, rows=test)
round(misclassification_rate(ŷ, y[test]), sigdigits=3)

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
int_rules_1_2_1_6 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.1.6*(g+h*x)^m*(a+b*x+c*x^2)^p*(d+e*x+f*x^2)^q =#
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * (~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q, ~x) => (~c / ~f) ^ ~p * Antiderivative((~g + ~h * ~x) ^ ~m * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~p + ~q), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~p, ~q], ~x) && (EqQ(~c * ~d - ~a * ~f, 0) && (EqQ(~b * ~d - ~a * ~e, 0) && ((IntegerQ(~p) || GtQ(~c / ~f, 0)) && (Not(IntegerQ(~q)) || LeafCount(~d + ~e * ~x + ~f * (~x) ^ 2) <= LeafCount(~a + ~b * ~x + ~c * (~x) ^ 2)))))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * (~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q, ~x) => (((~a) ^ IntPart(~p) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ FracPart(~p)) / ((~d) ^ IntPart(~p) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ FracPart(~p))) * Antiderivative((~g + ~h * ~x) ^ ~m * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~p + ~q), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~p, ~q], ~x) && (EqQ(~c * ~d - ~a * ~f, 0) && (EqQ(~b * ~d - ~a * ~e, 0) && (Not(IntegerQ(~p)) && (Not(IntegerQ(~q)) && Not(GtQ(~c / ~f, 0))))))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * (~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q, ~x) => ((~a + ~b * ~x + ~c * (~x) ^ 2) ^ FracPart(~p) / ((4 * ~c) ^ IntPart(~p) * (~b + 2 * ~c * ~x) ^ (2 * FracPart(~p)))) * Antiderivative((~g + ~h * ~x) ^ ~m * (~b + 2 * ~c * ~x) ^ (2 * ~p) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~p, ~q], ~x) && EqQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * (~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(f') * (~x) ^ 2) ^ ~q, ~x) => ((~a + ~b * ~x + ~c * (~x) ^ 2) ^ FracPart(~p) / ((4 * ~c) ^ IntPart(~p) * (~b + 2 * ~c * ~x) ^ (2 * FracPart(~p)))) * Antiderivative((~g + ~h * ~x) ^ ~m * (~b + 2 * ~c * ~x) ^ (2 * ~p) * (~d + ~f * (~x) ^ 2) ^ ~q, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~g, ~h, ~m, ~p, ~q], ~x) && EqQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~g + ~(h') * ~x) ^ ~(m') * (~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~(m'), ~x) => Antiderivative(((~d * ~g) / ~a + (~f * ~h * ~x) / ~c) ^ ~m * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~m + ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~p], ~x) && (EqQ((~c * (~g) ^ 2 - ~b * ~g * ~h) + ~a * (~h) ^ 2, 0) && (EqQ(((~c) ^ 2 * ~d * (~g) ^ 2 - ~a * ~c * ~e * ~g * ~h) + (~a) ^ 2 * ~f * (~h) ^ 2, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~g + ~(h') * ~x) ^ ~(m') * (~a + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~(m'), ~x) => Antiderivative(((~d * ~g) / ~a + (~f * ~h * ~x) / ~c) ^ ~m * (~a + ~c * (~x) ^ 2) ^ (~m + ~p), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~h, ~p], ~x) && (EqQ(~c * (~g) ^ 2 + ~a * (~h) ^ 2, 0) && (EqQ(((~c) ^ 2 * ~d * (~g) ^ 2 - ~a * ~c * ~e * ~g * ~h) + (~a) ^ 2 * ~f * (~h) ^ 2, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~g + ~(h') * ~x) ^ ~(m') * (~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(f') * (~x) ^ 2) ^ ~(m'), ~x) => Antiderivative(((~d * ~g) / ~a + (~f * ~h * ~x) / ~c) ^ ~m * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~m + ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~g, ~h, ~p], ~x) && (EqQ((~c * (~g) ^ 2 - ~b * ~g * ~h) + ~a * (~h) ^ 2, 0) && (EqQ((~c) ^ 2 * ~d * (~g) ^ 2 + (~a) ^ 2 * ~f * (~h) ^ 2, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~g + ~(h') * ~x) ^ ~(m') * (~a + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(f') * (~x) ^ 2) ^ ~(m'), ~x) => Antiderivative(((~d * ~g) / ~a + (~f * ~h * ~x) / ~c) ^ ~m * (~a + ~c * (~x) ^ 2) ^ (~m + ~p), ~x) <-- FreeQ([~a, ~c, ~d, ~f, ~g, ~h, ~p], ~x) && (EqQ(~c * (~g) ^ 2 + ~a * (~h) ^ 2, 0) && (EqQ((~c) ^ 2 * ~d * (~g) ^ 2 + (~a) ^ 2 * ~f * (~h) ^ 2, 0) && IntegerQ(~m)))
    #= Antiderivative(((~g)+(~h')*(~x))^(~m')*((~a')+(~b')*(~x)+(~c')*(~x)^2)^(~p)*((~d')+(~e')*(~x)+(~f')*(~x)^2)^ (~q),~x) := Antiderivative((g+h*x)^(m+p)*(a/g+c/h*x)^p*(d+e*x+f*x^2)^q,x) <--  FreeQ([a,b,c,d,e,f,g,h,m,q],x) && EqQ(c*g^2-b*g*h+a*h^2,0) &&  IntegerQ(p) =#
    #= Antiderivative(((~g)+(~h')*(~x))^(~m')*((~a)+(~c')*(~x)^2)^(~p)*((~d')+(~e')*(~x)+(~f')*(~x)^2)^(~q),(~x)*Symbol) := Antiderivative((g+h*x)^(m+p)*(a/g+c/h*x)^p*(d+e*x+f*x^2)^q,x) <--  FreeQ([a,c,d,e,f,g,h,m,q],x) && NeQ(e^2-4*d*f,0) &&  EqQ(c*g^2+a*h^2,0) && IntegerQ(p) =#
    #= Antiderivative(((~g)+(~h')*(~x))^(~m')*((~a')+(~b')*(~x)+(~c')*(~x)^2)^(~p)*((~d')+(~f')*(~x)^2)^(~q),(~x)*Symbol) := Antiderivative((g+h*x)^(m+p)*(a/g+c/h*x)^p*(d+f*x^2)^q,x) <--  FreeQ([a,b,c,d,f,g,h,m,q],x) && NeQ(b^2-4*a*c,0) &&  EqQ(c*g^2-b*g*h+a*h^2,0) && IntegerQ(p) =#
    #= Antiderivative(((~g)+(~h')*(~x))^(~m')*((~a)+(~c')*(~x)^2)^(~p)*((~d')+(~f')*(~x)^2)^(~q),~x) :=   Antiderivative((g+h*x)^(m+p)*(a/g+c/h*x)^p*(d+f*x^2)^q,x) <--  FreeQ([a,c,d,f,g,h,m,q],x) && EqQ(c*g^2+a*h^2,0) && IntegerQ(p) =#
    #= Antiderivative(((~g)+(~h')*(~x))^(~m')*((~a')+(~b')*(~x)+(~c')*(~x)^2)^(~p)*((~d')+(~e')*(~x)+(~f')*(~x)^2)^ (~q),~x) :=  (a+b*x+c*x^2)^FracPart(p)/((g+h*x)^FracPart(p)*(a/g+(c*x)/h)^FracPart( p))*Antiderivative((g+h*x)^(m+p)*(a/g+c/h*x)^p*(d+e*x+f*x^2)^q,x) <--  FreeQ([a,b,c,d,e,f,g,h,m,q],x) && EqQ(c*g^2-b*g*h+a*h^2,0) &&  Not(IntegerQ(p)) =#
    #= Antiderivative(((~g)+(~h')*(~x))^(~m')*((~a)+(~c')*(~x)^2)^(~p)*((~d')+(~e')*(~x)+(~f')*(~x)^2)^(~q),(~x)*Symbol) :=  (a+c*x^2)^FracPart(p)/((g+h*x)^FracPart(p)*(a/g+(c*x)/h)^FracPart(p))* Antiderivative((g+h*x)^(m+p)*(a/g+c/h*x)^p*(d+e*x+f*x^2)^q,x) <--  FreeQ([a,c,d,e,f,g,h,m,q],x) && NeQ(e^2-4*d*f,0) &&  EqQ(c*g^2+a*h^2,0) && Not(IntegerQ(p)) =#
    #= Antiderivative(((~g)+(~h')*(~x))^(~m')*((~a')+(~b')*(~x)+(~c')*(~x)^2)^(~p)*((~d')+(~f')*(~x)^2)^(~q),(~x)*Symbol) :=  (a+b*x+c*x^2)^FracPart(p)/((g+h*x)^FracPart(p)*(a/g+(c*x)/h)^FracPart( p))*Antiderivative((g+h*x)^(m+p)*(a/g+c/h*x)^p*(d+f*x^2)^q,x) <--  FreeQ([a,b,c,d,f,g,h,m,q],x) && NeQ(b^2-4*a*c,0) &&  EqQ(c*g^2-b*g*h+a*h^2,0) && Not(IntegerQ(p)) =#
    #= Antiderivative(((~g)+(~h')*(~x))^(~m')*((~a)+(~c')*(~x)^2)^(~p)*((~d')+(~f')*(~x)^2)^(~q),~x) :=   (a+c*x^2)^FracPart(p)/((g+h*x)^FracPart(p)*(a/g+(c*x)/h)^FracPart(p) )*Antiderivative((g+h*x)^(m+p)*(a/g+c/h*x)^p*(d+f*x^2)^q,x) <--  FreeQ([a,c,d,f,g,h,m,q],x) && EqQ(c*g^2+a*h^2,0) && Not(IntegerQ(p)) =#
    @apply_utils Antiderivative((~x) ^ ~p * (~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q, ~x) => Antiderivative((~a / ~e + (~c / ~f) * ~x) ^ ~p * (~e * ~x + ~f * (~x) ^ 2) ^ (~p + ~q), ~x) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~q], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ((~c * (~e) ^ 2 - ~b * ~e * ~f) + ~a * (~f) ^ 2, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative((~x) ^ ~p * (~a + ~(c') * (~x) ^ 2) ^ ~p * (~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q, ~x) => Antiderivative((~a / ~e + (~c / ~f) * ~x) ^ ~p * (~e * ~x + ~f * (~x) ^ 2) ^ (~p + ~q), ~x) <-- FreeQ([~a, ~c, ~e, ~f, ~q], ~x) && (EqQ(~c * (~e) ^ 2 + ~a * (~f) ^ 2, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~g + ~(h') * ~x) / ((~a + ~(c') * (~x) ^ 2) ^ (1 / 3) * (~d + ~(f') * (~x) ^ 2)), ~x) => ((sqrt(3) * ~h * atan(1 / sqrt(3) - (2 ^ (2 / 3) * (1 - (3 * ~h * ~x) / ~g) ^ (2 / 3)) / (sqrt(3) * (1 + (3 * ~h * ~x) / ~g) ^ (1 / 3)))) / (2 ^ (2 / 3) * (~a) ^ (1 / 3) * ~f) + (~h * log(~d + ~f * (~x) ^ 2)) / (2 ^ (5 / 3) * (~a) ^ (1 / 3) * ~f)) - (3 * ~h * log((1 - (3 * ~h * ~x) / ~g) ^ (2 / 3) + 2 ^ (1 / 3) * (1 + (3 * ~h * ~x) / ~g) ^ (1 / 3))) / (2 ^ (5 / 3) * (~a) ^ (1 / 3) * ~f) <-- FreeQ([~a, ~c, ~d, ~f, ~g, ~h], ~x) && (EqQ(~c * ~d + 3 * ~a * ~f, 0) && (EqQ(~c * (~g) ^ 2 + 9 * ~a * (~h) ^ 2, 0) && GtQ(~a, 0)))
    @apply_utils Antiderivative((~g + ~(h') * ~x) / ((~a + ~(c') * (~x) ^ 2) ^ (1 / 3) * (~d + ~(f') * (~x) ^ 2)), ~x) => ((1 + (~c * (~x) ^ 2) / ~a) ^ (1 / 3) / (~a + ~c * (~x) ^ 2) ^ (1 / 3)) * Antiderivative((~g + ~h * ~x) / ((1 + (~c * (~x) ^ 2) / ~a) ^ (1 / 3) * (~d + ~f * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~c, ~d, ~f, ~g, ~h], ~x) && (EqQ(~c * ~d + 3 * ~a * ~f, 0) && (EqQ(~c * (~g) ^ 2 + 9 * ~a * (~h) ^ 2, 0) && Not(GtQ(~a, 0))))
    @apply_utils Antiderivative((~g + ~(h') * ~x) * (~a + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(f') * (~x) ^ 2) ^ ~q, ~x) => ~g * Antiderivative((~a + ~c * (~x) ^ 2) ^ ~p * (~d + ~f * (~x) ^ 2) ^ ~q, ~x) + ~h * Antiderivative(~x * (~a + ~c * (~x) ^ 2) ^ ~p * (~d + ~f * (~x) ^ 2) ^ ~q, ~x) <-- FreeQ([~a, ~c, ~d, ~f, ~g, ~h, ~p, ~q], ~x)
    @apply_utils Antiderivative((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q * (~(g') + ~(h') * ~x), ~x) => Antiderivative(ExpandIntegrand((~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q * (~g + ~h * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (IGtQ(~p, 0) && IntegerQ(~q))))
    @apply_utils Antiderivative((~a + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q * (~(g') + ~(h') * ~x), ~x) => Antiderivative(ExpandIntegrand((~a + ~c * (~x) ^ 2) ^ ~p * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q * (~g + ~h * ~x), ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (IntegersQ(~p, ~q) && (GtQ(~p, 0) || GtQ(~q, 0))))
    @apply_utils Antiderivative((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q * (~(g') + ~(h') * ~x), ~x) => (((~g * ~b - 2 * ~a * ~h) - (~b * ~h - 2 * ~g * ~c) * ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q) / (((~b) ^ 2 - 4 * ~a * ~c) * (~p + 1)) - (1 / (((~b) ^ 2 - 4 * ~a * ~c) * (~p + 1))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~q - 1) * Simp(((~e * ~q * (~g * ~b - 2 * ~a * ~h) - ~d * (~b * ~h - 2 * ~g * ~c) * (2 * ~p + 3)) + (2 * ~f * ~q * (~g * ~b - 2 * ~a * ~h) - ~e * (~b * ~h - 2 * ~g * ~c) * (2 * ~p + ~q + 3)) * ~x) - ~f * (~b * ~h - 2 * ~g * ~c) * (2 * ~p + 2 * ~q + 3) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (LtQ(~p, -1) && GtQ(~q, 0))))
    @apply_utils Antiderivative((~a + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q * (~(g') + ~(h') * ~x), ~x) => ((~a * ~h - ~g * ~c * ~x) * (~a + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q) / (2 * ~a * ~c * (~p + 1)) + (2 / (4 * ~a * ~c * (~p + 1))) * Antiderivative((~a + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~q - 1) * Simp((~g * ~c * ~d * (2 * ~p + 3) - ~a * (~h * ~e * ~q)) + (~g * ~c * ~e * (2 * ~p + ~q + 3) - ~a * (2 * ~h * ~f * ~q)) * ~x + ~g * ~c * ~f * (2 * ~p + 2 * ~q + 3) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (LtQ(~p, -1) && GtQ(~q, 0)))
    @apply_utils Antiderivative((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(f') * (~x) ^ 2) ^ ~q * (~(g') + ~(h') * ~x), ~x) => (((~g * ~b - 2 * ~a * ~h) - (~b * ~h - 2 * ~g * ~c) * ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~f * (~x) ^ 2) ^ ~q) / (((~b) ^ 2 - 4 * ~a * ~c) * (~p + 1)) - (1 / (((~b) ^ 2 - 4 * ~a * ~c) * (~p + 1))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~f * (~x) ^ 2) ^ (~q - 1) * Simp((-(~d) * (~b * ~h - 2 * ~g * ~c) * (2 * ~p + 3) + (2 * ~f * ~q * (~g * ~b - 2 * ~a * ~h)) * ~x) - ~f * (~b * ~h - 2 * ~g * ~c) * (2 * ~p + 2 * ~q + 3) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (LtQ(~p, -1) && GtQ(~q, 0)))
    @apply_utils Antiderivative((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q * (~(g') + ~(h') * ~x), ~x) => (((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~q + 1)) / (((~b) ^ 2 - 4 * ~a * ~c) * ((~c * ~d - ~a * ~f) ^ 2 - (~b * ~d - ~a * ~e) * (~c * ~e - ~b * ~f)) * (~p + 1))) * (~g * ~c * (2 * ~a * ~c * ~e - ~b * (~c * ~d + ~a * ~f)) + (~g * ~b - ~a * ~h) * ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (~b * ~e + 2 * ~a * ~f)) + ~c * (~g * ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (~b * ~e + 2 * ~a * ~f)) - ~h * ((~b * ~c * ~d - 2 * ~a * ~c * ~e) + ~a * ~b * ~f)) * ~x) + (1 / (((~b) ^ 2 - 4 * ~a * ~c) * ((~c * ~d - ~a * ~f) ^ 2 - (~b * ~d - ~a * ~e) * (~c * ~e - ~b * ~f)) * (~p + 1))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q * Simp(((((~b * ~h - 2 * ~g * ~c) * ((~c * ~d - ~a * ~f) ^ 2 - (~b * ~d - ~a * ~e) * (~c * ~e - ~b * ~f)) * (~p + 1) + (((~b) ^ 2 * (~g * ~f) - ~b * (~h * ~c * ~d + ~g * ~c * ~e + ~a * ~h * ~f)) + 2 * (~g * ~c * (~c * ~d - ~a * ~f) - ~a * (-(~h) * ~c * ~e))) * (~a * ~f * (~p + 1) - ~c * ~d * (~p + 2))) - ~e * ((~g * ~c) * (2 * ~a * ~c * ~e - ~b * (~c * ~d + ~a * ~f)) + (~g * ~b - ~a * ~h) * ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (~b * ~e + 2 * ~a * ~f))) * (~p + ~q + 2)) - (2 * ~f * ((~g * ~c) * (2 * ~a * ~c * ~e - ~b * (~c * ~d + ~a * ~f)) + (~g * ~b - ~a * ~h) * ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (~b * ~e + 2 * ~a * ~f))) * (~p + ~q + 2) - (((~b) ^ 2 * ~g * ~f - ~b * (~h * ~c * ~d + ~g * ~c * ~e + ~a * ~h * ~f)) + 2 * (~g * ~c * (~c * ~d - ~a * ~f) - ~a * (-(~h) * ~c * ~e))) * (~b * ~f * (~p + 1) - ~c * ~e * (2 * ~p + ~q + 4))) * ~x) - ~c * ~f * (((~b) ^ 2 * (~g * ~f) - ~b * (~h * ~c * ~d + ~g * ~c * ~e + ~a * ~h * ~f)) + 2 * (~g * ~c * (~c * ~d - ~a * ~f) + ~a * ~h * ~c * ~e)) * (2 * ~p + 2 * ~q + 5) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~q], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (LtQ(~p, -1) && (NeQ((~c * ~d - ~a * ~f) ^ 2 - (~b * ~d - ~a * ~e) * (~c * ~e - ~b * ~f), 0) && Not(Not(IntegerQ(~p)) && ILtQ(~q, -1))))))
    @apply_utils Antiderivative((~a + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q * (~(g') + ~(h') * ~x), ~x) => (((~a + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~q + 1)) / ((-4 * ~a * ~c) * (~a * ~c * (~e) ^ 2 + (~c * ~d - ~a * ~f) ^ 2) * (~p + 1))) * (~g * ~c * (2 * ~a * ~c * ~e) + (-(~a) * ~h) * (2 * (~c) ^ 2 * ~d - ~c * (2 * ~a * ~f)) + ~c * (~g * (2 * (~c) ^ 2 * ~d - ~c * (2 * ~a * ~f)) - ~h * (-2 * ~a * ~c * ~e)) * ~x) + (1 / ((-4 * ~a * ~c) * (~a * ~c * (~e) ^ 2 + (~c * ~d - ~a * ~f) ^ 2) * (~p + 1))) * Antiderivative((~a + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q * Simp(((((-2 * ~g * ~c) * ((~c * ~d - ~a * ~f) ^ 2 - (-(~a) * ~e) * (~c * ~e)) * (~p + 1) + (2 * (~g * ~c * (~c * ~d - ~a * ~f) - ~a * (-(~h) * ~c * ~e))) * (~a * ~f * (~p + 1) - ~c * ~d * (~p + 2))) - ~e * ((~g * ~c) * (2 * ~a * ~c * ~e) + (-(~a) * ~h) * (2 * (~c) ^ 2 * ~d - ~c * (2 * ~a * ~f))) * (~p + ~q + 2)) - (2 * ~f * ((~g * ~c) * (2 * ~a * ~c * ~e) + (-(~a) * ~h) * (2 * (~c) ^ 2 * ~d + -(~c) * (2 * ~a * ~f))) * (~p + ~q + 2) - (2 * (~g * ~c * (~c * ~d - ~a * ~f) - ~a * (-(~h) * ~c * ~e))) * (-(~c) * ~e * (2 * ~p + ~q + 4))) * ~x) - ~c * ~f * (2 * (~g * ~c * (~c * ~d - ~a * ~f) - ~a * (-(~h) * ~c * ~e))) * (2 * ~p + 2 * ~q + 5) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~h, ~q], ~x) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (LtQ(~p, -1) && (NeQ(~a * ~c * (~e) ^ 2 + (~c * ~d - ~a * ~f) ^ 2, 0) && Not(Not(IntegerQ(~p)) && ILtQ(~q, -1)))))
    @apply_utils Antiderivative((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(f') * (~x) ^ 2) ^ ~q * (~(g') + ~(h') * ~x), ~x) => (((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~f * (~x) ^ 2) ^ (~q + 1)) / (((~b) ^ 2 - 4 * ~a * ~c) * ((~b) ^ 2 * ~d * ~f + (~c * ~d - ~a * ~f) ^ 2) * (~p + 1))) * ((~g * ~c) * (-(~b) * (~c * ~d + ~a * ~f)) + (~g * ~b - ~a * ~h) * ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (2 * ~a * ~f)) + ~c * (~g * ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (2 * ~a * ~f)) - ~h * (~b * ~c * ~d + ~a * ~b * ~f)) * ~x) + (1 / (((~b) ^ 2 - 4 * ~a * ~c) * ((~b) ^ 2 * ~d * ~f + (~c * ~d - ~a * ~f) ^ 2) * (~p + 1))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p + 1) * (~d + ~f * (~x) ^ 2) ^ ~q * Simp((((~b * ~h - 2 * ~g * ~c) * ((~c * ~d - ~a * ~f) ^ 2 - (~b * ~d) * (-(~b) * ~f)) * (~p + 1) + (((~b) ^ 2 * (~g * ~f) - ~b * (~h * ~c * ~d + ~a * ~h * ~f)) + 2 * (~g * ~c * (~c * ~d - ~a * ~f))) * (~a * ~f * (~p + 1) - ~c * ~d * (~p + 2))) - (2 * ~f * ((~g * ~c) * (-(~b) * (~c * ~d + ~a * ~f)) + (~g * ~b - ~a * ~h) * ((2 * (~c) ^ 2 * ~d + (~b) ^ 2 * ~f) - ~c * (2 * ~a * ~f))) * (~p + ~q + 2) - (((~b) ^ 2 * (~g * ~f) - ~b * (~h * ~c * ~d + ~a * ~h * ~f)) + 2 * (~g * ~c * (~c * ~d - ~a * ~f))) * (~b * ~f * (~p + 1))) * ~x) - ~c * ~f * (((~b) ^ 2 * (~g * ~f) - ~b * (~h * ~c * ~d + ~a * ~h * ~f)) + 2 * (~g * ~c * (~c * ~d - ~a * ~f))) * (2 * ~p + 2 * ~q + 5) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~g, ~h, ~q], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (LtQ(~p, -1) && (NeQ((~b) ^ 2 * ~d * ~f + (~c * ~d - ~a * ~f) ^ 2, 0) && Not(Not(IntegerQ(~p)) && ILtQ(~q, -1)))))
    @apply_utils Antiderivative((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q * (~(g') + ~(h') * ~x), ~x) => (~h * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~q + 1)) / (2 * ~f * (~p + ~q + 1)) - (1 / (2 * ~f * (~p + ~q + 1))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p - 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q * Simp(~h * ~p * (~b * ~d - ~a * ~e) + ~a * (~h * ~e - 2 * ~g * ~f) * (~p + ~q + 1) + (2 * ~h * ~p * (~c * ~d - ~a * ~f) + ~b * (~h * ~e - 2 * ~g * ~f) * (~p + ~q + 1)) * ~x + (~h * ~p * (~c * ~e - ~b * ~f) + ~c * (~h * ~e - 2 * ~g * ~f) * (~p + ~q + 1)) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~q], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (GtQ(~p, 0) && NeQ(~p + ~q + 1, 0))))
    @apply_utils Antiderivative((~a + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q * (~(g') + ~(h') * ~x), ~x) => (~h * (~a + ~c * (~x) ^ 2) ^ ~p * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ (~q + 1)) / (2 * ~f * (~p + ~q + 1)) + (1 / (2 * ~f * (~p + ~q + 1))) * Antiderivative((~a + ~c * (~x) ^ 2) ^ (~p - 1) * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q * Simp(((~a * ~h * ~e * ~p - ~a * (~h * ~e - 2 * ~g * ~f) * (~p + ~q + 1)) - 2 * ~h * ~p * (~c * ~d - ~a * ~f) * ~x) - (~h * ~c * ~e * ~p + ~c * (~h * ~e - 2 * ~g * ~f) * (~p + ~q + 1)) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~h, ~q], ~x) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (GtQ(~p, 0) && NeQ(~p + ~q + 1, 0)))
    @apply_utils Antiderivative((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~d + ~(f') * (~x) ^ 2) ^ ~q * (~(g') + ~(h') * ~x), ~x) => (~h * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p * (~d + ~f * (~x) ^ 2) ^ (~q + 1)) / (2 * ~f * (~p + ~q + 1)) - (1 / (2 * ~f * (~p + ~q + 1))) * Antiderivative((~a + ~b * ~x + ~c * (~x) ^ 2) ^ (~p - 1) * (~d + ~f * (~x) ^ 2) ^ ~q * Simp(~h * ~p * (~b * ~d) + ~a * (-2 * ~g * ~f) * (~p + ~q + 1) + (2 * ~h * ~p * (~c * ~d - ~a * ~f) + ~b * (-2 * ~g * ~f) * (~p + ~q + 1)) * ~x + (~h * ~p * (-(~b) * ~f) + ~c * (-2 * ~g * ~f) * (~p + ~q + 1)) * (~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~g, ~h, ~q], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (GtQ(~p, 0) && NeQ(~p + ~q + 1, 0)))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / ((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * (~d + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([q = Simplify((((((~c) ^ 2 * (~d) ^ 2 - ~b * ~c * ~d * ~e) + ~a * ~c * (~e) ^ 2 + (~b) ^ 2 * ~d * ~f) - 2 * ~a * ~c * ~d * ~f) - ~a * ~b * ~e * ~f) + (~a) ^ 2 * (~f) ^ 2)], (1 / q) * Antiderivative(Simp(((((~g * (~c) ^ 2 * ~d - ~g * ~b * ~c * ~e) + ~a * ~h * ~c * ~e + ~g * (~b) ^ 2 * ~f) - ~a * ~b * ~h * ~f) - ~a * ~g * ~c * ~f) + ~c * (((~h * ~c * ~d - ~g * ~c * ~e) + ~g * ~b * ~f) - ~a * ~h * ~f) * ~x, ~x) / (~a + ~b * ~x + ~c * (~x) ^ 2), ~x) + (1 / q) * Antiderivative(Simp(((((-(~h) * ~c * ~d * ~e + ~g * ~c * (~e) ^ 2 + ~b * ~h * ~d * ~f) - ~g * ~c * ~d * ~f) - ~g * ~b * ~e * ~f) + ~a * ~g * (~f) ^ 2) - ~f * (((~h * ~c * ~d - ~g * ~c * ~e) + ~g * ~b * ~f) - ~a * ~h * ~f) * ~x, ~x) / (~d + ~e * ~x + ~f * (~x) ^ 2), ~x) <-- NeQ(q, 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NeQ((~e) ^ 2 - 4 * ~d * ~f, 0))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / ((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * (~d + ~(f') * (~x) ^ 2)), ~x) => With([q = Simplify((((~c) ^ 2 * (~d) ^ 2 + (~b) ^ 2 * ~d * ~f) - 2 * ~a * ~c * ~d * ~f) + (~a) ^ 2 * (~f) ^ 2)], (1 / q) * Antiderivative(Simp((((~g * (~c) ^ 2 * ~d + ~g * (~b) ^ 2 * ~f) - ~a * ~b * ~h * ~f) - ~a * ~g * ~c * ~f) + ~c * ((~h * ~c * ~d + ~g * ~b * ~f) - ~a * ~h * ~f) * ~x, ~x) / (~a + ~b * ~x + ~c * (~x) ^ 2), ~x) + (1 / q) * Antiderivative(Simp(((~b * ~h * ~d * ~f - ~g * ~c * ~d * ~f) + ~a * ~g * (~f) ^ 2) - ~f * ((~h * ~c * ~d + ~g * ~b * ~f) - ~a * ~h * ~f) * ~x, ~x) / (~d + ~f * (~x) ^ 2), ~x) <-- NeQ(q, 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~g, ~h], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~g + ~(h') * ~x) / ((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => -2 * ~g * Subst(Antiderivative(1 / ((~b * ~d - ~a * ~e) - ~b * (~x) ^ 2), ~x), ~x, sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (EqQ(~c * ~e - ~b * ~f, 0) && EqQ(~h * ~e - 2 * ~g * ~f, 0))))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / ((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => (-((~h * ~e - 2 * ~g * ~f)) / (2 * ~f)) * Antiderivative(1 / ((~a + ~b * ~x + ~c * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) + (~h / (2 * ~f)) * Antiderivative((~e + 2 * ~f * ~x) / ((~a + ~b * ~x + ~c * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (EqQ(~c * ~e - ~b * ~f, 0) && NeQ(~h * ~e - 2 * ~g * ~f, 0))))
    @apply_utils Antiderivative(~x / ((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~d + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => -2 * ~e * Subst(Antiderivative((1 - ~d * (~x) ^ 2) / (((~c * ~e - ~b * ~f) - ~e * ((2 * ~c * ~d - ~b * ~e) + 2 * ~a * ~f) * (~x) ^ 2) + (~d) ^ 2 * (~c * ~e - ~b * ~f) * (~x) ^ 4), ~x), ~x, (1 + ((~e + sqrt((~e) ^ 2 - 4 * ~d * ~f)) * ~x) / (2 * ~d)) / sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && EqQ(~b * ~d - ~a * ~e, 0)))
    @apply_utils Antiderivative((~g + ~(h') * ~x) / ((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~d + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => ~g * Subst(Antiderivative(1 / (~a + (~c * ~d - ~a * ~f) * (~x) ^ 2), ~x), ~x, ~x / sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (EqQ(~b * ~d - ~a * ~e, 0) && EqQ(2 * ~h * ~d - ~g * ~e, 0))))
    @apply_utils Antiderivative((~g + ~(h') * ~x) / ((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~d + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => (-((2 * ~h * ~d - ~g * ~e)) / ~e) * Antiderivative(1 / ((~a + ~b * ~x + ~c * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) + (~h / ~e) * Antiderivative((2 * ~d + ~e * ~x) / ((~a + ~b * ~x + ~c * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (EqQ(~b * ~d - ~a * ~e, 0) && NeQ(2 * ~h * ~d - ~g * ~e, 0))))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / ((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => -2 * ~g * (~g * ~b - 2 * ~a * ~h) * Subst(Antiderivative(1 / Simp(~g * (~g * ~b - 2 * ~a * ~h) * ((~b) ^ 2 - 4 * ~a * ~c) - (~b * ~d - ~a * ~e) * (~x) ^ 2, ~x), ~x), ~x, Simp((~g * ~b - 2 * ~a * ~h) - (~b * ~h - 2 * ~g * ~c) * ~x, ~x) / sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (NeQ(~b * ~d - ~a * ~e, 0) && EqQ(((~h) ^ 2 * (~b * ~d - ~a * ~e) - 2 * ~g * ~h * (~c * ~d - ~a * ~f)) + (~g) ^ 2 * (~c * ~e - ~b * ~f), 0))))
    @apply_utils Antiderivative((~g + ~(h') * ~x) / ((~a + ~(c') * (~x) ^ 2) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => -2 * ~a * ~g * ~h * Subst(Antiderivative(1 / Simp(2 * (~a) ^ 2 * ~g * ~h * ~c + ~a * ~e * (~x) ^ 2, ~x), ~x), ~x, Simp(~a * ~h - ~g * ~c * ~x, ~x) / sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && EqQ((~a * (~h) ^ 2 * ~e + 2 * ~g * ~h * (~c * ~d - ~a * ~f)) - (~g) ^ 2 * ~c * ~e, 0)
    @apply_utils Antiderivative((~g + ~(h') * ~x) / ((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~d + ~(f') * (~x) ^ 2)), ~x) => -2 * ~g * (~g * ~b - 2 * ~a * ~h) * Subst(Antiderivative(1 / Simp(~g * (~g * ~b - 2 * ~a * ~h) * ((~b) ^ 2 - 4 * ~a * ~c) - ~b * ~d * (~x) ^ 2, ~x), ~x), ~x, Simp((~g * ~b - 2 * ~a * ~h) - (~b * ~h - 2 * ~g * ~c) * ~x, ~x) / sqrt(~d + ~f * (~x) ^ 2)) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && EqQ((~b * (~h) ^ 2 * ~d - 2 * ~g * ~h * (~c * ~d - ~a * ~f)) - (~g) ^ 2 * ~b * ~f, 0))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / ((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c * ~g - ~h * (~b - q)) / q) * Antiderivative(1 / (((~b - q) + 2 * ~c * ~x) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) - ((2 * ~c * ~g - ~h * (~b + q)) / q) * Antiderivative(1 / ((~b + q + 2 * ~c * ~x) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && PosQ((~b) ^ 2 - 4 * ~a * ~c)))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / ((~a + ~(c') * (~x) ^ 2) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([q = Rt(-(~a) * ~c, 2)], (~h / 2 + (~c * ~g) / (2q)) * Antiderivative(1 / ((-q + ~c * ~x) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) + (~h / 2 - (~c * ~g) / (2q)) * Antiderivative(1 / ((q + ~c * ~x) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && PosQ(-(~a) * ~c))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / ((~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~d + ~(f') * (~x) ^ 2)), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c * ~g - ~h * (~b - q)) / q) * Antiderivative(1 / (((~b - q) + 2 * ~c * ~x) * sqrt(~d + ~f * (~x) ^ 2)), ~x) - ((2 * ~c * ~g - ~h * (~b + q)) / q) * Antiderivative(1 / ((~b + q + 2 * ~c * ~x) * sqrt(~d + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && PosQ((~b) ^ 2 - 4 * ~a * ~c))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / ((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([q = Rt((~c * ~d - ~a * ~f) ^ 2 - (~b * ~d - ~a * ~e) * (~c * ~e - ~b * ~f), 2)], (1 / (2q)) * Antiderivative(Simp((~h * (~b * ~d - ~a * ~e) - ~g * ((~c * ~d - ~a * ~f) - q)) - (~g * (~c * ~e - ~b * ~f) - ~h * ((~c * ~d - ~a * ~f) + q)) * ~x, ~x) / ((~a + ~b * ~x + ~c * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) - (1 / (2q)) * Antiderivative(Simp((~h * (~b * ~d - ~a * ~e) - ~g * ((~c * ~d - ~a * ~f) + q)) - (~g * (~c * ~e - ~b * ~f) - ~h * ((~c * ~d - ~a * ~f) - q)) * ~x, ~x) / ((~a + ~b * ~x + ~c * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && (NeQ(~b * ~d - ~a * ~e, 0) && NegQ((~b) ^ 2 - 4 * ~a * ~c))))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / ((~a + ~(c') * (~x) ^ 2) * sqrt(~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([q = Rt((~c * ~d - ~a * ~f) ^ 2 + ~a * ~c * (~e) ^ 2, 2)], (1 / (2q)) * Antiderivative(Simp((-(~a) * ~h * ~e - ~g * ((~c * ~d - ~a * ~f) - q)) + (~h * ((~c * ~d - ~a * ~f) + q) - ~g * ~c * ~e) * ~x, ~x) / ((~a + ~c * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x) - (1 / (2q)) * Antiderivative(Simp((-(~a) * ~h * ~e - ~g * ((~c * ~d - ~a * ~f) + q)) + (~h * ((~c * ~d - ~a * ~f) - q) - ~g * ~c * ~e) * ~x, ~x) / ((~a + ~c * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~e) ^ 2 - 4 * ~d * ~f, 0) && NegQ(-(~a) * ~c))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / ((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~d + ~(f') * (~x) ^ 2)), ~x) => With([q = Rt((~c * ~d - ~a * ~f) ^ 2 + (~b) ^ 2 * ~d * ~f, 2)], (1 / (2q)) * Antiderivative(Simp((~h * ~b * ~d - ~g * ((~c * ~d - ~a * ~f) - q)) + (~h * ((~c * ~d - ~a * ~f) + q) + ~g * ~b * ~f) * ~x, ~x) / ((~a + ~b * ~x + ~c * (~x) ^ 2) * sqrt(~d + ~f * (~x) ^ 2)), ~x) - (1 / (2q)) * Antiderivative(Simp((~h * ~b * ~d - ~g * ((~c * ~d - ~a * ~f) + q)) + (~h * ((~c * ~d - ~a * ~f) - q) + ~g * ~b * ~f) * ~x, ~x) / ((~a + ~b * ~x + ~c * (~x) ^ 2) * sqrt(~d + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NegQ((~b) ^ 2 - 4 * ~a * ~c))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / (sqrt(~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~d + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([s = Rt((~b) ^ 2 - 4 * ~a * ~c, 2), t = Rt((~e) ^ 2 - 4 * ~d * ~f, 2)], ((sqrt(~b + s + 2 * ~c * ~x) * sqrt(2 * ~a + (~b + s) * ~x) * sqrt(~e + t + 2 * ~f * ~x) * sqrt(2 * ~d + (~e + t) * ~x)) / (sqrt(~a + ~b * ~x + ~c * (~x) ^ 2) * sqrt(~d + ~e * ~x + ~f * (~x) ^ 2))) * Antiderivative((~g + ~h * ~x) / (sqrt(~b + s + 2 * ~c * ~x) * sqrt(2 * ~a + (~b + s) * ~x) * sqrt(~e + t + 2 * ~f * ~x) * sqrt(2 * ~d + (~e + t) * ~x)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NeQ((~e) ^ 2 - 4 * ~d * ~f, 0))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / (sqrt(~a + ~(b') * ~x + ~(c') * (~x) ^ 2) * sqrt(~d + ~(f') * (~x) ^ 2)), ~x) => With([s = Rt((~b) ^ 2 - 4 * ~a * ~c, 2), t = Rt(-4 * ~d * ~f, 2)], ((sqrt(~b + s + 2 * ~c * ~x) * sqrt(2 * ~a + (~b + s) * ~x) * sqrt(t + 2 * ~f * ~x) * sqrt(2 * ~d + t * ~x)) / (sqrt(~a + ~b * ~x + ~c * (~x) ^ 2) * sqrt(~d + ~f * (~x) ^ 2))) * Antiderivative((~g + ~h * ~x) / (sqrt(~b + s + 2 * ~c * ~x) * sqrt(2 * ~a + (~b + s) * ~x) * sqrt(t + 2 * ~f * ~x) * sqrt(2 * ~d + t * ~x)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~g, ~h], ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / ((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ (1 / 3) * (~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([q = ((-9 * ~c * (~h) ^ 2) / (2 * ~c * ~g - ~b * ~h) ^ 2) ^ (1 / 3)], ((sqrt(3) * ~h * q * atan(1 / sqrt(3) - (2 ^ (2 / 3) * (1 - (3 * ~h * (~b + 2 * ~c * ~x)) / (2 * ~c * ~g - ~b * ~h)) ^ (2 / 3)) / (sqrt(3) * (1 + (3 * ~h * (~b + 2 * ~c * ~x)) / (2 * ~c * ~g - ~b * ~h)) ^ (1 / 3)))) / ~f + (~h * q * log(~d + ~e * ~x + ~f * (~x) ^ 2)) / (2 * ~f)) - (3 * ~h * q * log((1 - (3 * ~h * (~b + 2 * ~c * ~x)) / (2 * ~c * ~g - ~b * ~h)) ^ (2 / 3) + 2 ^ (1 / 3) * (1 + (3 * ~h * (~b + 2 * ~c * ~x)) / (2 * ~c * ~g - ~b * ~h)) ^ (1 / 3))) / (2 * ~f)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (EqQ(~c * ~e - ~b * ~f, 0) && (EqQ((~c) ^ 2 * ~d - ~f * ((~b) ^ 2 - 3 * ~a * ~c), 0) && (EqQ((((~c) ^ 2 * (~g) ^ 2 - ~b * ~c * ~g * ~h) - 2 * (~b) ^ 2 * (~h) ^ 2) + 9 * ~a * ~c * (~h) ^ 2, 0) && GtQ((-9 * ~c * (~h) ^ 2) / (2 * ~c * ~g - ~b * ~h) ^ 2, 0))))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) / ((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ (1 / 3) * (~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2)), ~x) => With([q = -(~c) / ((~b) ^ 2 - 4 * ~a * ~c)], ((q * (~a + ~b * ~x + ~c * (~x) ^ 2)) ^ (1 / 3) / (~a + ~b * ~x + ~c * (~x) ^ 2) ^ (1 / 3)) * Antiderivative((~g + ~h * ~x) / ((q * ~a + ~b * q * ~x + ~c * q * (~x) ^ 2) ^ (1 / 3) * (~d + ~e * ~x + ~f * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (EqQ(~c * ~e - ~b * ~f, 0) && (EqQ((~c) ^ 2 * ~d - ~f * ((~b) ^ 2 - 3 * ~a * ~c), 0) && (EqQ((((~c) ^ 2 * (~g) ^ 2 - ~b * ~c * ~g * ~h) - 2 * (~b) ^ 2 * (~h) ^ 2) + 9 * ~a * ~c * (~h) ^ 2, 0) && Not(GtQ(4 * ~a - (~b) ^ 2 / ~c, 0)))))
    @apply_utils Antiderivative((~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q * (~(g') + ~(h') * ~x), ~x) => Unintegrable((~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q * (~g + ~h * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~p, ~q], ~x)
    @apply_utils Antiderivative((~(a') + ~(c') * (~x) ^ 2) ^ ~p * (~(d') + ~(e') * ~x + ~(f') * (~x) ^ 2) ^ ~q * (~(g') + ~(h') * ~x), ~x) => Unintegrable((~a + ~c * (~x) ^ 2) ^ ~p * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q * (~g + ~h * ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~h, ~p, ~q], ~x)
    @apply_utils Antiderivative((~(g') + ~(h') * ~u) ^ ~(m') * (~(a') + ~(b') * ~u + ~(c') * (~u) ^ 2) ^ ~(p') * (~(d') + ~(e') * ~u + ~(f') * (~u) ^ 2) ^ ~(q'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~g + ~h * ~x) ^ ~m * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~p, ~q], ~x) && (LinearQ(~u, ~x) && NeQ(~u, ~x))
    @apply_utils Antiderivative((~(g') + ~(h') * ~u) ^ ~(m') * (~(a') + ~(c') * (~u) ^ 2) ^ ~(p') * (~(d') + ~(e') * ~u + ~(f') * (~u) ^ 2) ^ ~(q'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~g + ~h * ~x) ^ ~m * (~a + ~c * (~x) ^ 2) ^ ~p * (~d + ~e * ~x + ~f * (~x) ^ 2) ^ ~q, ~x), ~x, ~u) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~p, ~q], ~x) && (LinearQ(~u, ~x) && NeQ(~u, ~x))
end

ext: jl

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


text:
export Term,
    Formula,
    PredicateSymbol,
    FunctionSymbol,
    Variable,
    Function,
    AtomicFormula,
    Negation,
    EFormula,
    AFormula,
    Conjunction,
    Disjunction,
    Literal,
    Clause,
    CNF,
    Substitution,
    Unifiable,
    OverMaxSearchDepthError

abstract type Term end
abstract type Formula end

@auto_hash_equals struct PredicateSymbol
    symbol::String
end

@auto_hash_equals struct FunctionSymbol
    symbol::String
end

@auto_hash_equals struct Variable <: Term
    variable::String
end

@auto_hash_equals struct Function <: Term
    name::FunctionSymbol
    arguments::Vector{Term}
end

@auto_hash_equals struct AtomicFormula <: Formula
    name::PredicateSymbol
    arguments::Vector{Term}
end

@auto_hash_equals struct Negation <: Formula
    formula::Formula
end

@auto_hash_equals struct EFormula <: Formula
    variable::Variable
    formula::Formula
end

@auto_hash_equals struct AFormula <: Formula
    variable::Variable
    formula::Formula
end

@auto_hash_equals struct Conjunction <: Formula
    formula1::Formula
    formula2::Formula
end

@auto_hash_equals struct Disjunction <: Formula
    formula1::Formula
    formula2::Formula
end

@auto_hash_equals struct Literal
    negative::Bool
    formula::AtomicFormula
end

const Clause = Set{Literal}
const CNF = Set{Clause}  # Conjunctive normal form
const Substitution = Vector{Pair{Variable, Term}}
const Unifiable = Union{Literal, Term}

String(x::PredicateSymbol) = x.symbol
String(x::FunctionSymbol) = x.symbol
String(x::Variable) = x.variable
String(x::Function) = "$(String(x.name))($(join(String.(x.arguments), ", ")))"
String(x::AtomicFormula) = "$(String(x.name))($(join(String.(x.arguments), ", ")))"
String(x::Negation) = "~$(String(x.formula))"
String(x::EFormula) = "?{$(String(x.variable))}$(String(x.formula))"
String(x::AFormula) = "*{$(String(x.variable))}$(String(x.formula))"
String(x::Conjunction) = "($(String(x.formula1)) & $(String(x.formula2)))"
String(x::Disjunction) = "($(String(x.formula1)) | $(String(x.formula2)))"
String(x::Literal) = "$(x.negative ? "~" : "")$(String(x.formula))"
String(x::Clause) = "{$(join(String.(collect(x)), ", "))}"
String(x::CNF) = "{$(join(String.(collect(x)), ", "))}"

Formula(literal::Literal) = begin
    literal.negative ? Negative(literal.formula) : literal.formula
end

Base.convert(::Type{Formula}, cnf::Set{Set{Literal}}) = begin
    if length(cnf) == 1
        Formula(first(cnf))
    else
        head, tail = headtail(cnf)
        @match cnf begin
        x::CNF => Conjunction(Formula(head), Formula(tail))
        x::Clause => Disjunction(Formula(head), Formula(tail))
        end
    end
end

type OverMaxSearchDepthError <: Exception
end

type NotUnifiableError <: Exception
end


text:
using Compat
try
    using Test
catch
    using Base.Test
end

using Signals
using Signals: state

include("push_pull.jl")
include("benchmark.jl")
include("operators.jl")
include("binding.jl")
include("time.jl")


text:
@testset "Test Dof handlers" begin

using jFEMTools
using SparseArrays

mesh = unitSquareMesh(TriangleCell, (2,2));
dim = 2
element = PoissonVirtualElement(dim,1);
u = TrialFunction(VEMFunctionSpace(mesh,element))
dh = DofHandler(mesh, [u]);

@test dh.cell_dofs == [1, 2, 3, 2, 4, 3, 2, 5, 4, 5, 6, 4, 3, 4, 7, 4, 8, 7, 4, 6, 8, 6, 9, 8]
@test dh.cell_dofs_offset == [1, 4, 7, 10, 13, 16, 19, 22, 25]
@test jFEMTools.dof_range(dh, u, 1) == 1:3
@test jFEMTools.vertexdofs(dh, u) == [1,2,5,3,4,6,7,8,9]
K = jFEMTools.create_sparsity_pattern(dh)
@test size(K) == (9,9)
@test nnz(K) == 41
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
# Pawel Zielinski
module matrixgen

using LinearAlgebra

export blockmat

function matcond(n::Int, c::Float64)
	# Function generates a random square matrix A of size n with
	# a given condition number c.
	# Inputs:
	#	n: size of matrix A, n>1
	#	c: condition of matrix A, c>= 1.0
	#
	# Usage: matcond(10, 100.0)
	#
	# Pawel Zielinski
    if n < 2
        error("size n should be > 1")
    end
    if c < 1.0
        error("condition number  c of a matrix  should be >= 1.0")
    end
    (U, S, V) = svd(rand(n, n))
    return U * diagm(0 => [LinRange(1.0, c, n);]) * V'
end

  function blockmat(n::Int, l::Int, ck::Float64, outputfile::String)
		# Function generates a random block sparse matrix A of size n with
		# a given condition number ck of inner block Ak and it save the output
		# matrix in a text file.
		# Inputs:
		#	n: size of block matrix A, n>1
		# l: size of inner matrices Ak, n mod l =0 (n is  divisible by l)
		#	ck: condition of inner matrix Ak, ck>= 1.0
		# outputfile: name of the output text file
		#
		# Usage: blockmat(100, 4 ,10.0, "A.txt")
		#
		#
		#  the output file format
	  #  n  l              <--- the size of block matrix A, the size of inner matrices Ak
		#  i1  j1   A[i1,j1] <--- a non-zero element of block matrix A
		#  i2  j2   A[i2,j2] <--- a non-zero element of block matrix A
		#  i3  j3   A[i3,j3] <--- a non-zero element of block matrix A
		#  ...
		#  ...
		#  EOF
		#

    if n < 2
        error("size n should be > 1")
    end
		if n % l != 0
			error("n is not divisible by l")
		end

		nb = div(n, l)
		Ak = Matrix{Float64}(undef, l, l)
		open(outputfile, "w") do f
			println(f, n, " ", l)
			for k in 1:nb
				Ak = matcond(l, ck)
				for i in 1:l, j in 1:l
					println(f, (k - 1) * l + i, " ", (k - 1) * l + j, " ", Ak[i,j])
				end
				if k < nb
			  	 for i in 1:l
						println(f, (k - 1) * l + i, " ", k * l + i, " ", 0.3 * rand())
				 	 end
				end
				if k > 1
			   for i in 1:l
					 println(f, (k - 1) * l + i, " ", (k - 1) * l, " ", 0.3 * rand())
				 end
				end
			end
		end	 # do
	end # blockmat

end # matrixgen

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
module URITemplate

	function is_valid_literal(ch)
		if ch > 0x80
			return true
		end
		return ch != '\'' && ch != '\<' && ch != '>' && ch != '<' && ch != '"' &&
			   ch != '^' && ch != '`' && ch != '{' && ch != '|' && ch != '}' && ch != '%'
	end

	ops = ['+','#','.','/',';','?','&','=',',','!','@','|']

	isunreserved(ch) = isalnum(ch) || ch == '-' || ch == '.' || ch == '_' || ch == '~'

    isreserved(ch) = ch == ':' || ch == '/' || ch == '?' || ch == '#' || ch == '[' ||
    				 ch == ']' || ch == '@' || ch == '!' || ch == '$' || ch == '&' ||
    				 ch == '\''|| ch == '(' || ch == ')' || ch == '*' || ch == '+' ||
    				 ch == ',' || ch == ';' || ch == '='

	pctencode(s::IO,c::UInt8) = write(s,'%',base(16,c,2))

	function pctencode(s::IO, c::Char, allowR = false)
	    if c < 0x80
	    	if isunreserved(c)
	    		write(s,c)
	    	elseif allowR && isreserved(c)
	    		write(s,c)
	    	else
	        	pctencode(s, uint8(c))
	        end
	        return 1
	    elseif c < 0x800
	        pctencode(s, uint8(( c >> 6          ) | 0xC0))
	        pctencode(s, uint8(( c        & 0x3F ) | 0x80))
	        return 2
	    elseif c < 0x10000
	        pctencode(s, uint8(( c >> 12         ) | 0xE0))
	        pctencode(s, uint8(((c >> 6)  & 0x3F ) | 0x80))
	        pctencode(s, uint8(( c        & 0x3F ) | 0x80))
	        return 3
	    elseif c < 0x110000
	        pctencode(s, uint8(( c >> 18         ) | 0xF0))
	        pctencode(s, uint8(((c >> 12) & 0x3F ) | 0x80))
	        pctencode(s, uint8(((c >> 6)  & 0x3F ) | 0x80))
	        pctencode(s, uint8(( c        & 0x3F ) | 0x80))
	        return 4
	    else
	        return pctencode(s, '\ufffd')
	    end
	end

	function pctencode(s::IO, string::AbstractString, allowR = false)
		for c in string
			pctencode(s,c,allowR)
		end
	end

	function is_valid_pctencoding(string,i)
		done(string,i) && return false
		(ch,i) = next(string,i)
		(!ishex(ch) || done(string,i)) && return false
		(ch,i) = next(string,i)
		!ishex(ch) && return false
		return true
	end

	if VERSION < v"0.4.0"
		# eltype of Dict is now pair and not a Tuple
		function keytype( dict )
			return eltype(dict)[1]
		end
		# Size hint is replaced by sizehint!
		sizehint! = sizehint
	end

	function expand(template::AbstractString,variables)
		if !( keytype( variables ) <: AbstractString)
			variables = [string(k) => v for (k,v) in variables]
		end

		out = IOBuffer()

		# As a heuristic the result will probably be about as long as the template
		# in either case it's probably not much shorter, so we can avoid spurious
		# allocation in the early phases, without too much overhead.
		sizehint!(out.data,sizeof(template))

		i = start(template)
		while !done(template,i)
			(ch,i) = next(template,i)
			if ch == '%' && !is_valid_pctencoding(template,i)
				error("'%' encountered but percent encoding sequence invalid")
			end
			if ch != '{'
				if !is_valid_literal(ch)
					error("Non-literal character '$ch' encountered!")
				end
				write(out,ch)
				continue
			end
			# Expression
			j = i
			while !done(template,i)
				(ch,i) = next(template,i)
				if ch == '}'
					break
				end
			end
			done(template,prevind(template,i)) && error("Template ended while scanning expression")
			#The expression excluding '{' and '}'
			ex = SubString(template,j,prevind(template,prevind(template,i)))
			isempty(ex) && error("Expression may not be empty!")
			(ch,j) = next(ex,start(ex))
			op = :NUL
			#   .------------------------------------------------------------------.
   			#	|          NUL     +      .       /       ;      ?      &      #   |
   			#	|------------------------------------------------------------------|
   			#	| first |  ""     ""     "."     "/"     ";"    "?"    "&"    "#"  |
   			#   | sep   |  ","    ","    "."     "/"     ";"    "&"    "&"    ","  |
   			#   | named | false  false  false   false   true   true   true   false |
   			#   | ifemp |  ""     ""     ""      ""      ""     "="    "="    ""   |
   			#   | allow |   U     U+R     U       U       U      U      U     U+R  |
   			#   `------------------------------------------------------------------'
			first = ""
			sep = ","
			named = false
			allowR = false
			ifemp = ""
			if ch == '+'
				op = :+
				allowR = true
			elseif ch == '.'
				first = "."
				sep = "."
				op = :.
			elseif ch == '/'
				first = "/"
				sep = "/"
				op = :/
			elseif ch == ';'
				first = ";"
				sep = ";"
				named = true
				op = :semicolon
			elseif ch == '?'
				first = "?"
				sep = "&"
				named = true
				ifemp = "="
				op = :?
			elseif ch == '&'
				first = "&"
				sep = "&"
				named = true
				ifemp = "="
				op = :&
			elseif ch == '#'
				first = "#"
				sep = ","
				op = :hash
				allowR = true
			elseif ch in ops
				error("Unimplemented template operator")
			end
			k = j
			if op != :NUL
				(ch,k) = next(ex,j)
				#k = nextind(ex,k)
			else
				j = prevind(ex,k)
			end

			first_defined = true
			#Process the list of variable names
			while true
				if done(ex,k) || (!isalnum(ch) && ch != '_' && ch != '.' && !(ch == '%' && is_valid_pctencoding(ex,j))) #End of variable name
					explode = false
					limitlength = false
					limit = 0
					if k > sizeof(ex)
						varend = sizeof(ex)
					else
						varend = prevind(ex,nextind(ex,k))
					end
					if ch == '*'
						explode = true
						if !done(ex,k)
							(ch,k) = next(ex,k)
							varend = prevind(ex,varend)
						end
						varend = prevind(ex,varend)
					elseif ch == ','
						varend = prevind(ex,prevind(ex,varend))
					elseif ch == ':'
						limitlength = true
						kl,l = k,0
						varend = prevind(ex,prevind(ex,varend))
						while !done(ex,k)
							l>3 && break
							(ch,k) = next(ex,k)
							if !('0' <= ch <= '9')
								break
							end
							l+=1
						end
						l == 0 && error("Zero-length : postfix not allowed")
						if done(ex,k)
							limit = parseint(ex[kl:sizeof(ex)],10)
						else
							limit = parseint(ex[kl:prevind(ex,prevind(ex,k))],10)
						end
					elseif !done(ex,k)
						error("Spurious characters past the end of expression.")
					end
					varname = ex[prevind(ex,nextind(ex,j)):varend]
					if !haskey(variables,varname)
						if done(ex,k)
							break
						end
						(ch,k) = next(ex,k)
					end
					if first_defined
						write(out,first)
						first_defined = false
					else
						write(out,sep)
					end
					val = variables[varname]
					if isa(val,AbstractString)
						if named
							write(out,varname)
							if isempty(val)
								write(out,ifemp)
							else
								write(out,'=')
							end
						end
						if limitlength
							l,m = 0,start(val)
							while l < limit && !done(val,m)
								(ch,m) = next(val,m)
								pctencode(out,ch,allowR)
								l+=1
							end
						else
							pctencode(out,val,allowR)
						end
					elseif !explode
						if named
							write(out,varname)
							if isempty(val)
								write(out,ifemp)
							else
								write(out,'=')
							end
						end
						f = true
						if isa(val,Associative)
							for (key,v) in val
								if !f
									print(out,',')
								else
									f = false
								end
								pctencode(out,key,allowR)
								write(out,',')
								pctencode(out,v,allowR)
							end
						else
							for v in val
								if !f
									print(out,',')
								else
									f = false
								end
								pctencode(out,v,allowR)
							end
						end
					else
						if named
							firstentry = true
							if isa(val,Associative)
								for (key,v) in val
									if !firstentry
										write(out,sep)
									else
										firstentry = false
									end
									write(out,key)
									if isempty(v)
										write(out,ifemp)
									else
										write(out,'=')
										pctencode(out,v)
									end
								end
							else
								for v in val
									if !firstentry
										write(out,sep)
									else
										firstentry = false
									end
									write(out,varname)
									if isempty(v)
										write(out,ifemp)
									else
										write(out,'=')
										pctencode(out,v)
									end
								end
							end
						elseif isa(val,Associative)
							f = true
							for (key,v) in val
								if !f
									print(out,sep)
								else
									f = false
								end
								pctencode(out,key,allowR)
								write(out,'=')
								pctencode(out,v,allowR)
							end
						else
							f = true
							for v in val
								if !f
									print(out,sep)
								else
									f = false
								end
								pctencode(out,v,allowR)
							end
						end
					end
					j=k
				end
				if done(ex,k)
					break
				end
				(ch,k) = next(ex,k)
			end
		end
		takebuf_string(out)
	end

end


text:
module How2Mix

# Unassuming code in an unassuming module...
module SubFoo

function f()
    x = rand()
    y = rand()
    return x + y
end

end

using Mixtape
import Mixtape: CompilationContext, transform, allow
using CodeInfoTools
using MacroTools

# 101: How2Mix
struct MyMix <: CompilationContext end

# A few little utility functions for working with Expr instances.
swap(e) = e
function swap(e::Expr)
    new = MacroTools.postwalk(e) do s
        isexpr(s, :call) || return s
        s.args[1] == Base.rand || return s
        return 4
    end
    return new
end

# This is pre-inference - you get to see a CodeInfoTools.Builder instance.
function transform(::MyMix, src)
    b = CodeInfoTools.Builder(src)
    for (v, st) in b
        b[v] = swap(st)
    end
    return CodeInfoTools.finish(b)
end

# MyMix will only transform functions which you explicitly allow.
# You can also greenlight modules.
allow(ctx::MyMix, m::Module) = m == SubFoo

# This loads up a call interface which will cache the result of the pipeline.
Mixtape.@load_abi()
@assert(call(SubFoo.f; ctx = MyMix()) == 8)
@assert(SubFoo.f() != 8)

end # module


text:
# copyright Lénaïc Chizat
# the code shows the training dynamics in parameter and predictor space for a 2-layer relu NN on a two class classification task
# we optimize the exponential loss (the behavior is sensibly the same with the logistic loss)
# we use a specific step-size schedule motivated by the theory in this ref https://arxiv.org/abs/2002.04486 (sections 4 & 5)

# HOW TO RUN THE CODE
# in a prompt first run:
# include("exp_bothlayers.jl")
# Then for fast prototyping, run for instance the following code (takes 1 min)
# illustration(4, 60, 200, 0.4, 100_000, 10, 0.05);
# These parameters are explained in the header of the function "illustration"
# A good illustration is obtained as follows (takes 20 min):
# Random.seed!(9); # 9 generates a nice random problem
# illustration(4, 60, 400, 0.4, 1000000, 400, 0.005);

##
using LinearAlgebra, Random
using PyPlot, ProgressMeter # these packages need to be installed (via "] add NamePackage" )
using Colors
using3D()

##

Base.@kwdef mutable struct TwoLayer
    d::Int  # input dim
    m::Int  # output dim
    W::Array{Float64, 2}  # the input and output weighst
    TwoLayer(d::Int, m::Int) = new(d, m, randn(m, d+1))
end

"""
Initialise a model
"""
function init!(l::TwoLayer, alpha1::Real=1.0, alpha2::Real=1.0, sym::Bool=false)
    d, m = l.d, l.m
    # input weights initialized uniformly on the unit sphere
    l.W[:,1:d] ./= sqrt.(sum(l.W[:,1:d].^2, dims=2))
    l.W[:, 1:d] .*= alpha1
    # output weights initialized by -1 or 1 with equal proba
    l.W[1:div(m,2),end] .= alpha2
    l.W[div(m,2)+1:end,end] .= -alpha2
    if sym
        @assert iseven(m)
        W[div(m,2)+1:end,1:d]  .=  W[1:div(m,2),1:d]
    end
end

mutable struct Result
    niter::Int
    m::Int
    d::Int
    Ws::Array{Float64, 3}     # store optimization path
    loss::Vector{Float64}  # loss is -log of the empirical risk
    margins::Vector{Float64}
    betas::Vector{Float64}

    Result(niter::Int, model::TwoLayer) = new(niter, model.m, model.d,
                                              zeros(model.m, model.d+1, niter),
                                              zeros(niter),
                                              zeros(niter),
                                              zeros(niter)
                                             )

end
##
"""
Gradient descent to train a 2-layers ReLU neural network for the exponential loss
We use the step-size schedule from in https://arxiv.org/abs/2002.04486 (sections 4 & 5)
INPUT: X (training input), Y (training output), m (nb neurons), both: training both layers or just the output layer
OUTPUT: Ws (training trajectory)
"""
function twonet(model::TwoLayer, X, Y, m, stepsize, niter)::Result

    (n,d) = size(X) # n samples in R^d
    # initialize
    result = Result(niter, model)

    @showprogress 1 "Training neural network..." for iter = 1:niter
        result.Ws[:,:,iter] = model.W
        act  =  max.( model.W[:,1:end-1] * X', 0.0) # activations
        out  =  (1/m) * sum( model.W[:,end] .* act , dims=1) # predictions of the network
        perf = Y .* out[:]
        margin = minimum(perf)
        temp = exp.(margin .- perf) # numerical stabilization of exp
        gradR = temp .* Y ./ sum(temp)' # gradient of the loss
        grad_w1 = (model.W[:,end] .* float.(act .> 0) * ( X .* gradR  ))  # gradient for input weights
        grad_w2 = act * gradR  # gradient for output weights

        grad = cat(grad_w1, grad_w2, dims=2) # size (m × d+1)
        result.betas[iter] = sum(model.W.^2)/m
        result.loss[iter] = margin - log(sum(exp.(margin .- perf))/n)
        result.margins[iter] = margin/result.betas[iter] # margin in F_1 norm
        model.W = model.W + stepsize * grad/(sqrt(iter+1))
    end

    return result
end

##

"Coordinates of the 2d cluster centers, p is k^2 the number of clusters"
function cluster_center(p,k, Delta)
    p1 = mod1.(p,k)
    p2 = fld1.(p, k)
    x1 =  Delta*3*(p1 .- 1) .- 1/2
    x2 =  Delta*3*(p2 .- 1) .- 1/2
    return x1, x2
end

mutable struct Data
    d::Int  # input dimension
    p::Int  # output dimension
    n::Int  # number of data points
    X::Array{Float64}
    Y::Array{Float64}
end

function Data(k::Int, n::Int)
    sd = 0 # number of spurious dimensions with random noise
    Delta = 1/(3k-1) # interclass distance
    A = ones(k^2) # cluster affectation
    A[randperm(k^2)[1:div(k^2,2)]] .= -1  # random -1 label

    # sample from this data distribution
    P = rand(1:k^2,n) # cluster label
    T = 2π*rand(n)  # shift angle
    R = Delta*rand(n) # shift magnitude
    X = #[ 2(rand(n, 2).-.5) ones(n) ]
    [ cluster_center(P,k, Delta)[1] + R .* cos.(T) cluster_center(P,k, Delta)[2] + R .* sin.(T) ones(n) ]
        # dims=2)  #
    # X = cat(
        # ones(n),
        # cluster_center(P,k)[1] + R .* cos.(T),
        # cluster_center(P,k)[2] + R .* sin.(T),
        # rand(n,sd) .- 1/2,
        # dims=2)  #
    Y = A[P]
    Data(2, 1, n, X, Y)
end

##
"""
Plot the classifier for a test case, comparing training both  layer
k: number of clusters per dimension in the data set (choose 3, 4 or 5)
n: number of training samples
m: number of neurons
stepsize, niter: parameters for training
"""
function create_test_case(k=4, n=60)
    # data distribution (itself random)
    # plot training set
    data = Data(k, n)
    X, Y = data.X, data.Y
    X1 = X[(Y .== 1),:]
    X2 = X[(Y .== -1),:]
    fig = figure(figsize=[3,3])
    plot(X1[:,1],X1[:,2],"+k")
    plot(X2[:,1],X2[:,2],"_k")
    axis("equal");axis("off");
    display(fig)
    data
end

function train_test_case(data::Data; m=200, stepsize=0.05, niter=100_000)

    X, Y = data.X, data.Y
    (n, d) = size(X)
    model = TwoLayer(d, m)
    init!(model)
    # train the neural network
    twonet(model, X, Y, m, stepsize, niter)
end

"""
nframes, resolution: parameters for the plots
"""
function plot_results(result::Result, nframes::Int=20, resolution::Real=0.05)

    # define the sequence of time steps ts to be plotted, with a power law
    niter = result.niter
    a = (niter-1)/(nframes-1)^4
    ts = setdiff(Int.(floor.(a*(0:nframes-1).^4)) .+ 1)  # in order to get unique frames
    Ws, margins, betas, loss = result.Ws, result.margins, result.betas, result.loss
    Ws = Ws[:,:,ts]

    m = size(Ws, 1)
    X, Y = data.X, data.Y
    # these are the positions we plot for the particles
    # Wproj of size (n,d,nframes)

    # Wproj = result.Ws[:,1:end-1,:] .* abs.(result.Ws[:,end:end,:])
    Wcut = Ws[:, 1:end-1, :]

    # the norm, size n*1*nframes
    Wout = result.Ws[:, end, :]  # will color the points
    vmin, vmax = minimum(Wout), maximum(Wout)
    Ncolors = 100
    cmapidx = floor.(Int, (Ncolors-1) .* (Wout .- vmin) ./ (vmax-vmin) .+ 1) # between 1 and N
    to_rgb(c::RGB) = [red(c), green(c), blue(c)]
    cmap = colormap("RdBu", Ncolors)
    cmap = to_rgb.(cmap)
    # WN   = sqrt.(sum(Wproj.^2, dims = 2))
    # Wdir = Wproj ./ WN
    # Wlog = tanh.(0.5*WN) .* Wdir

    @showprogress 1 "Plotting images..." for k = 1:length(ts)
        ioff() # turns off interactive plotting
        fig = figure(figsize=[7,4])
        ax1 = subplot(121, projection="3d")
        ax1.set_position([0,0.1,0.5,0.8])

        indt = k<11 ? (1:k) : (k-10):k  # take the 10 most recent k for the tail

        for i = 1:size(Wcut,1)  #for all the  units
            # x, y, z is 2, 3, 1 ?
            plot3D(Wcut[i,1,indt],Wcut[i,2,indt],Wcut[i,3,indt], color="k", linewidth=0.2) # tail
            # plot3D(Wlog[i,2,indt],Wlog[i,3,indt],Wlog[i,1,indt], color="k", linewidth=0.2) # tail
        end
        # plot the positive in blue and the negative in red
        scatter3D(Wcut[:,1,k],Wcut[:,2,k],Wcut[:,3,k],"o",color=cmap[cmapidx[:, k]])

        # ax1.set_xlim3d(-1, 1)
        # ax1.set_ylim3d(-1, 1)
        # ax1.set_zlim3d(-1, 1)
        # ax1.set_xticks([-1/2, 0, 1/2])
        # ax1.set_yticks([-1/2, 0, 1/2])
        # ax1.set_zticks([-1/2, 0, 1/2])
        # ax1.view_init(25-20*sinpi(k/length(ts)),45-45*cospi(k/length(ts)))

        ax2 = subplot(122)
        # ax2.set_position([0.45,0.25,0.5,0.5])

        f(x1,x2,k) = (1/m) * sum( Ws[:,end,k] .* max.( Ws[:,1:3,k] * [x1;x2;1], 0.0)) # prediction function
        # with relu

        xs = -1:resolution:1
        tab = [f(xs[i],xs[j],k) for i=1:length(xs), j=1:length(xs)]
        pcolormesh(xs', xs, tanh.(tab'), cmap="coolwarm", shading="gouraud", vmin=-1.0, vmax=1.0, edgecolor="face")

        contour(xs', xs, tanh.(tab'), levels =0, colors="k", antialiased = true, linewidths=2)

        # plot training set
        X1 = X[(Y .== 1),:]
        X2 = X[(Y .== -1),:]
        plot(X1[:,1],X1[:,2],"+k")
        plot(X2[:,1],X2[:,2],"_k")
        axis("equal");axis("off");
        ax2.set_xticks([-1/2, 0, 1/2])
        ax2.set_yticks([-1/2, 0, 1/2])

        PyPlot.savefig("dynamics_$(k).png",bbox_inches="tight", dpi=300)
        close(fig)
    end

end
##

max_stars_repo_name: brechetp/ImplicitBias

max_issues_repo_name: brechetp/ImplicitBias

max_forks_repo_name: brechetp/ImplicitBias


text:
s = BenchmarkGroup()
suite["dismantling"] = s

Random.seed!(1)
for G in GLIST
    for n in [100]
        l = 2
        g = G(n , 4n, seed=1)
        s["ci l=$l","$g"] = @benchmarkable dismantle_ci($g, $l, $(n÷4))
    end
end


text:
module SDatasets

export data_ball

function say(sth)
  println(sth)
end

function data_ball(N=1, center=[0 0], radius=[1];
          class=NaN, distribution="uniform")
  dim = size(center, 2)
  n_center = size(center, 1)
  n_radius = size(radius, 1)

  width = isnan(class) ? dim : dim + 1
  data = Array(Float64, N, width)
  for i = 1:N
    pt = Array(Float64, 1, dim)
    r = radius[(i - 1) % n_radius + 1]
    while true
      pt[:] = (rand(dim) - .5) * (r * 2)
      if (pt*pt')[1] <= r * r
        break
      end
    end

    data[i,1:dim] = pt + center[(i - 1) % n_radius + 1, :]
    if !isnan(class)
      data[i, width] = class
    end
  end

  data
end

end # module


text:
abstract type MarkableFloat <: AbstractFloat end

struct MarkableFloat64{M} <: MarkableFloat
   value::Float64
end

struct MarkableFloat32{M} <: MarkableFloat
   value::Float32
end

value(x::MarkableFloat64{M}) where {M} = x.value
value(x::MarkableFloat32{M}) where {M} = x.value

MarkableFloat64(x::Float64) = MarkableFloat64{0}(x)
MarkableFloat32(x::Float32) = MarkableFloat32{0}(x)

ismarked(x::MarkableFloat64{0}) = false
ismarked(x::MarkableFloat32{0}) = false
ismarked(x::MarkableFloat64{M}) where {M} = true
ismarked(x::MarkableFloat32{M}) where {M} = true

isunmarked(x::MarkableFloat64{0}) = true
isunmarked(x::MarkableFloat32{0}) = true
isunmarked(x::MarkableFloat64{M}) where {M} = false
isunmarked(x::MarkableFloat32{M}) where {M} = false

marking(x::MarkableFloat64{0}) = nothing
marking(x::MarkableFloat32{0}) = nothing
marking(x::MarkableFloat64{M}) = M
marking(x::MarkableFloat32{M}) = M

Markable64s = Union{MarkableFloat64{0}, MarkableFloat64{1}, MarkableFloat64{2}}
Markable32s = Union{MarkableFloat32{0}, MarkableFloat32{1}, MarkableFloat32{2}}


text:
## analytic.jl :special cases where eigenfunctions are known analytically

# (separable) exponential covariance function with p = 0.5
function compute_analytic(cov::CovarianceFunction{1,Exponential{T}} where {T},n::N where {N<:Integer},pts::V where {V<:AbstractVector})
	λ = cov.cov.λ
    ω = findroots(λ, n)
    ev = @. 2*λ/(λ^2*ω^2+1)
    n = @. sqrt(2)/2*sqrt(1/ω*(λ^2*ω^2*cos(ω)*sin(ω)+λ^2*ω^3-2*λ*ω*cos.(ω)^2-cos(ω)*sin(ω)+ω)+2*λ)
	ef = diagm(0 => 1.0/n)*( sin.(ω*pts') + λ*diagm(0 => ω)*cos.(ω*pts') )

	SpectralData(sqrt.(ev),ef')
end

# find all positive (>0) zeros of the transcendental function tan(ω) = 2*λ*ω/(λ^2*ω^2-1)
function findroots(λ::T, n::N) where {T <: AbstractFloat, N <: Integer}

	# define the transcendental function
	f(ω) = (λ^2*ω^2-1)*sin(ω)-2*λ*ω*cos(ω)

	# find range around singularity 1/λ
	left_point_of_range = (2*floor(1/(π*λ)-1/2))*π/2 # left odd multiple of π/2
	right_point_of_range = (2*ceil(1/(π*λ)-1/2)+1)*π/2 # right odd multiple of π/2

	# find roots before 1/λ, if any
	if left_point_of_range ≠ π/2
		roots = zeros(min(n,round(UInt64,floor(abs(1/λ/π-1/2)))))
		left_point = π/2
		right_point = 1*π
		for i = 1:length(roots)
			@inbounds roots[i] = bisect_root(f,left_point,right_point)[1]
			right_point = roots[i] + π
			left_point = left_point + π
		end
	else
		roots = zeros(0)
	end

	# find roots inside range around 1/λ
	( length(roots) ≥ n || floor(1/(π*λ)-1/2) < 0 ) ||
	push!(roots,bisect_root(f,left_point_of_range+eps(T),1/λ)[1]) # first intersection point
	( length(roots) ≥ n || ceil(1/(π*λ)-1/2) < 0 ) ||
	push!(roots,bisect_root(f,1/λ,right_point_of_range)[1]) # second intersection point

	# if the first root is zero, cut it off
	roots[1] == 0 ? shift!(roots) : [] # empty expression

	# find roots after 1/λ
	startindex = 1 + length(roots)
	if n-length(roots) > 0
		roots = [roots; zeros(n-length(roots))]
		left_point = (2*ceil(1/(π*λ)-1/2)+2)*π/2
		right_point = (2*ceil(1/(π*λ)-1/2)+3)*π/2
		for i = startindex:length(roots)
			@inbounds roots[i] = bisect_root(f,left_point,right_point)[1]
			right_point = roots[i] + π
			left_point = left_point + π
		end
	end
	return roots
end

# bissection method to find the zeros of a function in a particular interval [x1,x2]
function bisect_root(fn::Function, x1::T, x2::T) where {T<:Real}
	xm = middle(x1, x2)
	s1 = sign(fn(x1))
	s2 = sign(fn(x2))
	while x1 < xm < x2
		sm =  sign(fn(xm))

		if s1 != sm
			x2 = xm
			s2 = sm
		else
			x1 = xm
			s1 = sm
		end

		xm = middle(x1, x2)
	end

	return x1, x2
end

# helper function to find the real "mid point" of two given floating point numbers
function middle(x1::Float64, x2::Float64)
	# use the usual float rules for combining non-finite numbers
	if !isfinite(x1) || !isfinite(x2)
		return x1 + x2
	end

	# always return 0.0 when inputs have opposite sign
	if sign(x1) != sign(x2) && x1 != 0.0 && x2 != 0.0
		return 0.0
	end

	negate = x1 < 0.0 || x2 < 0.0

	x1_int = reinterpret(UInt64, abs(x1))
	x2_int = reinterpret(UInt64, abs(x2))
	unsigned = reinterpret(Float64, (x1_int + x2_int) >> 1)

	negate ? -unsigned : unsigned
end


text:
# physical constants
μ₀ = 4π * 1e-7
ħ = 1.055e-34
γH = 2.675e8
γC = 6.726e7
rCH = 1.117e-10  # Tugarinov 2004
rHH = sqrt(3) * rCH * sin(110.4*π/180)  # Tugarinov 2004
P2cosβ = -1/3.  # 109 degree tetrahedral geometry

max_stars_repo_name: chriswaudby/qq-hahnecho

max_issues_repo_name: chriswaudby/qq-hahnecho

max_forks_repo_name: chriswaudby/qq-hahnecho


text:
using Plots
using ConstraintSolver
dir = pkgdir(ConstraintSolver)

x = [
    "niallsudoku_5500",
    "niallsudoku_5501",
    "niallsudoku_5502",
    "niallsudoku_5503",
    "niallsudoku_6417",
    "niallsudoku_6249",
];

plot(; xaxis = ("Problem"), yaxis = ("Time in s"), title = "Killer sudoku special")
cs_010 = [0.73, 3.69, 0.16, 0.066, 0.13, 0.06]
plot!(x, cs_010, label = "CS v0.1.0", color = :red, seriestype = :scatter)
cs_017 = [0.607, 1.602, 0.173, 0.08, 0.27, 0.09]
plot!(x, cs_017, label = "CS v0.1.7", color = :orange, seriestype = :scatter)
or = [2.82, 1.67, 0.46, 0.74, 0.17, 0.17];
plot!(x, or, label = "OR-Tools", color = :black, seriestype = :scatter)

savefig(joinpath(dir, "benchmark/results/killer_sudoku_special/plots/current.png"))

plot(; xaxis = ("Problem"), yaxis = ("Time in s"), title = "Killer sudoku normal rules")
cs_010 = [0.16, 0.39, 0.09, 0.39, 0.30, 0.05]
plot!(x, cs_010, label = "CS v0.1.0", color = :red, seriestype = :scatter)
cs_017 = [0.096, 0.50, 0.095, 0.50, 0.48, 0.06]
plot!(x, cs_017, label = "CS v0.1.7", color = :orange, seriestype = :scatter)
or = [1.18, 1.64, 0.37, 1.07, 0.04, 0.04];
plot!(x, or, label = "OR-Tools", color = :black, seriestype = :scatter)

savefig(joinpath(dir, "benchmark/results/killer_sudoku/plots/current.png"))


text:
abstract type Extension end

struct ExtensionNone <: Extension end
struct ExtensionNothing <: Extension end
struct ExtensionReplicate <: Extension end
struct ExtensionCircular <: Extension end
struct ExtensionSymmetric <: Extension end
struct ExtensionConstant{T} <: Extension
    value::T
end

"""
    Kernel{X,F}

A kernel object with axes `X` that is wrapping a kernel function of type `F`.
"""
struct Kernel{X,F}
    f::F

    """
        Kernel{X}(f)

    Create a kernel with axes `X` wrapping a kernel function `f`.

    The kernel function `f` defines a reduction of values within an `X`-sized
    view. When the kernel is applied to data of an array `a` the kernel function
    gets called with one argument `w::Window` that provides a local view on the
    data.

    ```@example
    # Laplacian 3x3 Kernel (i.e. axes (-1:1, -1:1))
    function f(w)
        return w[0,-1] + w[-1,0] - 4*w[0,0] + w[1,0] + w[0,1]
    end
    Kernel{(-1:1,-1:1)}(f)
    ```

    For best performance you should annotate the kernel function `f` with `@inline`
    and index accesses within using `@inbounds`.
    """
    function Kernel{X}(f::F) where {X,F<:Function}
        X isa NTuple{<:Any,UnitRange{Int}} ||
            throw(ArgumentError("invalid axes: $X"))
        return new{X,F}(f)
    end
end

"""
    Window{T,N,X,K,A}

A stack-allocated array view fit for a kernel type `K` with indexing on axes
`X` relative to some position in the parent array.
"""
# FIXME: inheriting from AbstractArray prevents us from constant-propagating
# getindex(w, i). See also https://github.com/JuliaLang/julia/issues/35531
#
#struct Window{T,N,X,K,A} <: AbstractArray{T,N}
struct Window{T,N,X,K,A}
    position::CartesianIndex{N}
    parent::A
    kernel::K

    function Window{X}(k::Kernel, a::AbstractArray{T,N}, pos::CartesianIndex{N}) where {T,N,X}
        return new{T,N,X,typeof(k),typeof(a)}(pos, a, k)
    end
end

"""
    ExtensionArray{T,N,A,Ext} <: AbstractArray{T,N}

An array-wrapper mimicking its parent but additonally specifying an extension
for determining how to index out-of-bounds.
"""
struct ExtensionArray{T,N,A,Ext} <: AbstractArray{T,N}
    parent::A
    extension::Ext

    function ExtensionArray(parent::AbstractArray, extension::Extension)
        return new{eltype(parent),ndims(parent),typeof(parent),typeof(extension)}(parent, extension)
    end
end


text:
##x import Base.sort

function printList(levels)
  for l in levels    println(l.energy)   end
end

"""
`Basics.sort()`  ... sorts various object due to different criteria

  + `("multiplet: by energy", multiplet::Multiplet)`
    ... to sort all levels in the multiplet into a sequence of increasing energy; a (new) multiplet::Multiplet is returned.
"""
function Basics.sort(sa::String, multiplet::Multiplet)
    !(sa == "multiplet: by energy")   &&   error("Unsupported keystring = $sa")

    sortedLevels = Base.sort( multiplet.levels , lt=isless)
    newLevels = Level[];   index = 0
    for lev in sortedLevels
        index = index + 1
        push!(newLevels, Level(lev.J, lev.M, lev.parity, index, lev.energy, lev.relativeOcc, lev.hasStateRep, lev.basis, lev.mc) )
    end

    newMultiplet = JAC.Multiplet(multiplet.name, newLevels)

    return( newMultiplet )
end

"""
  + `("multiplet: by energy", multiplet::JAC.Hfs.IJF_Multiplet)`
    ... to sort all hyperfine levels in the multiplet into a sequence of increasing energy; a (new) multiplet::JAC.Hfs.IJF_Multiplet
        is returned.
"""
function Basics.sort(sa::String, multiplet::JAC.Hfs.IJF_Multiplet)
    !(sa == "multiplet: by energy")   &&   error("Unsupported keystring = $sa")

    sortedLevels = Base.sort( multiplet.levelFs , lt=isless)
    newLevels = JAC.Hfs.IJF_Level[];   index = 0
    for lev in sortedLevels
        index = index + 1
        push!(newLevels, JAC.Hfs.IJF_Level(lev.I, lev.F, lev.M, lev.parity, lev.energy, lev.hasStateRep, lev.basis, lev.mc) )
    end

    newMultiplet = JAC.Hfs.IJF_Multiplet(multiplet.name, newLevels)

    return( newMultiplet )
end


text:
import Base: iterate

closedform_systems(rt::RecurrenceTemplate, iter) =
    (ClosedFormTemplate(rt, part) for part in (isnothing(iter) ? partitions(size(rt)) : iter))

# TODO: How to deal with templates where a permutation is useless?
recurrence_systems(vars::Vector{Symbol}, shape::MatrixShape, iter=permutations(vars); kwargs...) =
    (RecurrenceTemplate(copy(perm), shape; kwargs...) for perm in (isnothing(iter) ? permutations(vars) : iter))

synthesis_problems(inv::Invariant, recsys_iter, part_iter; kwargs...) =
    (SynthesisProblem(inv, r, c) for r in recsys_iter for c in closedform_systems(r, part_iter))

synthesis_problems(inv::Invariant, vars::Vector{Symbol}, shape::MatrixShape, perm_iter, part_iter; kwargs...) =
    synthesis_problems(inv, recurrence_systems(vars, shape, perm_iter; kwargs...), part_iter)

# ------------------------------------------------------------------------------

function strategy_all(inv::Invariant, vars::Vector{Symbol}, shape::MatrixShape; kwargs...)
    @assert shape == UpperTriangular || shape == UnitUpperTriangular || shape == UserSpecific
    synthesis_problems(inv, vars, shape, nothing, nothing; kwargs...)
end

function strategy_permutations(inv::Invariant, vars::Vector{Symbol}, shape::MatrixShape, roots::Vector{Int}; kwargs...)
    @assert shape == UpperTriangular || shape == UnitUpperTriangular || shape == UserSpecific
    synthesis_problems(inv, vars, shape, nothing, [roots]; kwargs...)
end

function strategy_partitions(inv::Invariant, vars::Vector{Symbol}, shape::MatrixShape; kwargs...)
    rec = RecurrenceTemplate(copy(vars), shape; kwargs...)
    synthesis_problems(inv, [rec], nothing; kwargs...)
end

function strategy_fixed(inv::Invariant, vars::Vector{Symbol}, shape::MatrixShape, roots::Vector{Int}; kwargs...)
    rec = RecurrenceTemplate(vars, shape; kwargs...)
    synthesis_problems(inv, [rec], [roots]; kwargs...)
end

strategy_mixed(inv::Invariant, vars::Vector{Symbol}; kwargs...) =
    Iterators.flatten([
        strategy_all(inv, vars, UnitUpperTriangular; kwargs...),
        strategy_all(inv, vars, UpperTriangular; kwargs...),
        strategy_partitions(inv, vars, FullSymbolic; kwargs...)
    ])

# ------------------------------------------------------------------------------

abstract type AbstractPermutations end
struct AllPermutations <: AbstractPermutations
    vars::Vector{Symbol}
end
struct FixedPermutations <: AbstractPermutations
    perms::Vector{Vector{Symbol}}
end

FixedPermutations(ps::Vector{Symbol}...) = FixedPermutations(collect(ps))

perms(p::AllPermutations) = Combinatorics.permutations(p.vars)
perms(p::FixedPermutations) = p.perms

function permutations(r::RecTemplate, p::AbstractPermutations)
    if has_constant_one(r)
        return ([vs; CONST_ONE_SYM] for vs in perms(p))
    end
    return perms(p)
end

Base.summary(io::IO, p::AllPermutations; uppercase=true) = print(io, uppercase ? "A" : "a", "ll permutations of $(p.vars)")
Base.summary(io::IO, p::FixedPermutations; uppercase=true) = print(io, uppercase ? "F" : "f", "ixed permutations $(p.perms)")
Base.show(io::IO, p::AbstractPermutations) = summary(io, p)

abstract type AbstractPartitions end
struct AllPartitions <: AbstractPartitions
    n::Int
end
struct FixedPartitions <: AbstractPartitions
    parts::Vector{Vector{Int}}
end

FixedPartitions(ps::Vector{Int}...) = FixedPartitions(collect(ps))

partitions(p::AllPartitions) = Combinatorics.partitions(p.n)
partitions(p::FixedPartitions) = p.parts

Base.summary(io::IO, p::AllPartitions; uppercase=true) = print(io, uppercase ? "A" : "a", "ll integer partitions of $(p.n)")
Base.summary(io::IO, p::FixedPartitions; uppercase=true) = print(io, uppercase ? "F" : "f", "ixed partitions $(p.parts)")
Base.show(io::IO, p::AbstractPartitions) = summary(io, p)

# ------------------------------------------------------------------------------

abstract type AbstractStrategy end

solutions(s::AbstractStrategy; kwargs...) = solutions(synthesis_problems(s); kwargs...)
models(s::AbstractStrategy; kwargs...) = models(synthesis_problems(s); kwargs...)

struct Strategy <: AbstractStrategy
    inv::Invariant
    rec::RecTemplate
    perms::AbstractPermutations
    parts::AbstractPartitions

    function Strategy(inv::Invariant, rec::RecTemplate; perms::AbstractPermutations, parts::AbstractPartitions)
        inv_vars = program_variables(inv)
        rec_vars = variables(rec)
        @assert issubset(inv_vars, rec_vars)
        new(inv, rec, perms, parts)
    end
end

synthesis_problems(s::Strategy) =
    synthesis_problems(s.inv, recurrence_systems(s.rec, permutations(s.rec, s.perms)), partitions(s.parts))

recurrence_systems(r::RecTemplate, perms) =
    (RecurrenceTemplate(r, copy(p)) for p in perms)

function Base.show(io::IO, s::Strategy)
    compact = get(io, :compact, false)
    # if compact
        print(io, "$(typeof(s)) with ")
        summary(io, s.perms, uppercase=false)
        print(io, " and ")
        summary(io, s.parts, uppercase=false)
    # else
        # print(io, "$(typeof(r)) for some permutation σ of $(vars)")
    # end
end

struct CombinedStrategy <: AbstractStrategy
    ls::Vector{Strategy}

    CombinedStrategy(ls::Strategy...) = new(collect(ls))
end

synthesis_problems(s::CombinedStrategy) =
    Iterators.flatten(synthesis_problems(t) for t in s.ls)

function IterativeStrategy(inv::Invariant; constant_one::Bool=true)
    vars = program_variables(inv)
    size = constant_one ? length(vars)+1 : length(vars)
    rec1 = RecTemplate(copy(vars), constant_one=constant_one, shape=UnitUpperTriangular)
    str1 = Strategy(inv, rec1, perms=AllPermutations(vars), parts=FixedPartitions([size]))
    rec2 = RecTemplate(copy(vars), constant_one=constant_one, shape=UpperTriangular)
    str2 = Strategy(inv, rec2, perms=AllPermutations(vars), parts=AllPartitions(size))
    rec3 = RecTemplate(copy(vars), constant_one=constant_one, shape=FullSymbolic)
    str3 = Strategy(inv, rec3, perms=FixedPermutations(vars), parts=AllPartitions(size))
    return CombinedStrategy(str1, str2, str3)
end

# ------------------------------------------------------------------------------

struct Solutions
    solver::NLSolver
    problem::SynthesisProblem
    maxsol::Number
    timeout::Int

    function Solutions(p::SynthesisProblem; maxsol::Number=1, timeout::Int=2, solver::Type{<:NLSolver}=Z3Solver, constraints::ClauseSet=ClauseSet())
        @assert maxsol >= 1
        s = create_solver(p, solver)
        add!(s, constraints)
        new(s, p, maxsol, timeout)
    end
end

iterate(ms::Solutions) = iterate(ms, 0)

function iterate(ms::Solutions, state)
    (state === nothing || state >= ms.maxsol) && return
    model, result = solve(ms.problem, ms.solver, timeout=ms.timeout)
    model === nothing && return result, nothing
    next_constraints!(ms.solver, ms.problem, model)
    result, state+1
end

function next_constraints!(s::NLSolver, sp::SynthesisProblem, m::NLModel)
    vs = [Symbol(string(v)) for v in body(sp) if !isconstant(v)]
    cs = Clause([Constraint{NEQ}(:($v - ($(numerator(m[v])) / $(denominator(m[v]))))) for v in vs])
    add!(s, ClauseSet(cs))
end

# ------------------------------------------------------------------------------

solutions(strategy; kwargs...) =
    Iterators.flatten(Solutions(problem; kwargs...) for problem in strategy)

models(strategy; kwargs...) =
    Iterators.filter(Absynth.issat, solutions(strategy; kwargs...))

function synth(inv::Invariant; timeout=2, solver=Z3Solver, constant_one=true, kwargs...)
    progress = ProgressUnknown("I'm trying hard:")
    strategy = IterativeStrategy(inv; constant_one=constant_one)
    for s in solutions(strategy; maxsol=1, timeout=timeout, solver=solver)
        if issat(s)
            finish!(progress)
            return s.recsystem
        end
        next!(progress)
    end
    finish!(progress)
end


text:
#
# Set-up
#

# Check availability NVCC
nvcc = get(ENV, "NVCC", "nvcc")
nvcc_path = try
    chomp(readall(`which $nvcc` .> DevNull))
catch
    error("Could not find NVCC -- consider specifying with NVCC environment variable")
end
nvcc_ver = Nullable{VersionNumber}()
for line in readlines(`$(nvcc_path) --version`)
    m = match(r"release ([0-9.]+)", line)
    if m != nothing
        nvcc_ver = Nullable(VersionNumber(m.captures[1]))
    end
end
if isnull(nvcc_ver)
    error("Could not parse NVCC version info")
end

# Determine host compiler version requirements
const hostcc_support = [
    (v"5.0", v"4.6.4"),
    (v"5.5", v"4.7.2"),
    (v"6.0", v"4.8.1") ]
if get(nvcc_ver) < hostcc_support[1][1]
    error("No support for CUDA < $(hostcc_req[1][1])")
end
hostcc_maxver = Nullable{VersionNumber}()
for i = 1:length(hostcc_support)
    if get(nvcc_ver) == hostcc_support[i][1]
        hostcc_maxver = Nullable(hostcc_support[i][2])
        break
    end
end
if isnull(hostcc_maxver)
    error("Unknown NVCC version $(get(nvcc_ver))")
end

# Collect possible hostcc executable names
hostcc_names = [ "gcc" ]
for ver in [v"4.4" v"4.5" v"4.6" v"4.7" v"4.8" v"4.9"]
    push!(hostcc_names, "gcc-$(ver.major).$(ver.minor)")
    push!(hostcc_names, "gcc$(ver.major)$(ver.minor)")
end

# Check availability host compiler
hostcc_possibilities = []
for hostcc in hostcc_names
    hostcc_path = try
        chomp(readall(`which $hostcc` .> DevNull))
    catch
        continue
    end

    verstring = chomp(readlines(`$hostcc_path --version`)[1])
    m = match(Regex("^$hostcc \\(.*\\) ([0-9.]+)"), verstring)
    if m == nothing
        warn("Could not parse GCC version info (\"$verstring\")")
        continue
    end
    hostcc_ver = VersionNumber(m.captures[1])

    if hostcc_ver <= get(hostcc_maxver)
        push!(hostcc_possibilities, (hostcc_path, hostcc_ver))
    end
end
if length(hostcc_possibilities) == 0
    error("Could not find a suitable host compiler")
end
sort!(hostcc_possibilities; rev=true, lt=(a, b) -> a[2]<b[2])
hostcc = hostcc_possibilities[1]

# Determine compilation options
flags = [ "--compiler-bindir", hostcc[1] ]
if haskey(ENV, "ARCH")
    append!(flags, [ "--gpu-architecture", ENV["ARCH"] ])
end

#
# Compilation
#

using CUDA

scriptdir = dirname(Base.source_path())

# Generate a temporary file with specific suffix
function mkstemps(suffix::String)
    b = joinpath(tempdir(), "tmpXXXXXX$suffix")
    p = ccall(:mkstemps, Int32, (Ptr{UInt8}, Cint), b, length(suffix)) # modifies b
    return (b, fdio(p, true))
end

# Macro for C++ compilation
macro compile(kernel, code)
    kernel_name = string(kernel)
    containing_file = Base.source_path()
    :( _compile($kernel_name, $code, $containing_file) )
end
function _compile(kernel, code, containing_file)
    # Write the source into a compilable file
    (source, io) = mkstemps(".cu")
    write(io, """
extern "C"
{
$code
}
""")
    close(io)

    # Manage the build directory
    builddir = "$scriptdir/.build"
    if !isdir(builddir)
        mkdir(builddir)
    end

    # Check if we need to compile
    output = "$builddir/$kernel.ptx"
    if isfile(output)
        need_compile = (stat(containing_file).mtime > stat(output).mtime)
    else
        need_compile = true
    end

    # Compile the source, if necessary
    if need_compile
        try
            run(`$(nvcc_path) $flags -ptx -o $output $source`)
        catch
            error("compilation of kernel $kernel failed (typo in C++ source?)")
        end

        if !isfile(output)
            error("compilation of kernel $kernel failed (no output generated)")
        end
    end

    # Pass the module to the CUDA driver
    mod = try
        CuModule(output)
    catch
        error("loading of kernel $kernel failed (invalid CUDA code?)")
    end

    # Load the function pointer
    func = try
        CuFunction(mod, kernel)
    catch
        error("could not find kernel $kernel in the compiled binary (wrong function name?)")
    end

    # Generate a function for accessing the function pointer
    fname = symbol("$kernel")
    @eval begin
        function $fname()
            $func
        end
    end
end

# Include all files
scriptdir = dirname(Base.source_path())
for file in readdir(scriptdir)
    m = match(r"^(.*)\.jl$", file)
    if m != nothing && m.captures[1] != "load"
        include(file)
    end
end


text:
function edge_hab(dom, p = 0.5)
    lower = ceil(Int, 100 * (1 - p))
    upper = floor(Int, 100 * p)
    hab = zeros(Int, size(dom))
    hab[lower:upper, lower:upper] .= 1
    src = LinearIndices(size(dom))[findall(==(1), hab)]
    dgr = DistanceGradient(src)
    -rand(dgr, 100, 100)
end

function general_hab_distr(dom, mean = zeros(size(dom)))
    cov_kernel = Matérn32Cov(1.0, 20.0)
    cov_mat = cov(cov_kernel, dom)
    DomainDistribution(MvNormal(vec(mean), cov_mat), dom)
end

function general_hab_pref(gh)
    cdf(Normal(), gh)
end

function rocky_hab_distr(dom; wts = [0.75, 0.25])
    cov_kernel = Matérn32Cov(1.0, 5.0)
    cov_mat = cov(cov_kernel, dom)
    distr = MvNormal(cov_mat)
    ClassifiedDomainDistribution(distr, dom, wts)
end

function rocky_hab_pref_gen(pref = 1.0)
    function pref_fn(rh)
        rh ? pref : one(pref)
    end
    pref_fn
end

function generate_habitat(dom)
    gen_hab = rand(general_hab_distr(dom))
    rocky_hab = rand(rocky_hab_distr(dom))
    Habitat(gen_hab, rocky_hab)
end

#-------------------------------------------------------------------------------

function move_rate()
    mov_fn(d) = Matérn32Cov(1.0, 2.5)(d)
    MovementRate(mov_fn)
end

function hab_pref(rocky_pref)
    # HabitatPreference(general_hab_pref, rocky_hab_pref_gen(rocky_pref))
    HabitatPreference(gh -> 1, rocky_hab_pref_gen(rocky_pref))
end

function make_moveop(hab::Habitat, rocky_pref, dom::AbstractFisheryDomain)
    MovementModel(hab, hab_pref(rocky_pref), move_rate(), dom)
end

#-------------------------------------------------------------------------------

struct HabQSpec{T} <: SpatQSimSpec
    realization::Int
    rocky_pref::T
    prep_file::String

    function HabQSpec(realization::Int,
                      rocky_pref::T,
                      prep_file::String) where T<:Real
        new{T}(realization, rocky_pref, prep_file)
    end
end

function HabQSpec(realization::Integer,
                  rocky_pref::T;
                  base_dir = ".") where T<:Real
    prep_fn = prep_path(HabQSpec, realization; base_dir = base_dir)
    HabQSpec(realization, rocky_pref, prep_fn)
end

sim_value(spec::HabQSpec) = spec.rocky_pref

#-------------------------------------------------------------------------------

struct HabQPrep{D, H, M, P}
    domain::D
    habitat::H
    movement::M
    init_pop::P

    function HabQPrep(domain::D, habitat::H, movement::M, init_pop::P) where {D<:AbstractFisheryDomain,
                                                                              H<:Habitat,
                                                                              M<:MovementModel,
                                                                              P<:PopState}
        new{D, H, M, P}(domain, habitat, movement, init_pop)
    end
end

domain(prep::HabQPrep) = prep.domain
habitat(prep::HabQPrep) = prep.habitat
movement(prep::HabQPrep) = prep.movement
init_pop(prep::HabQPrep) = prep.init_pop

function HabQPrep(spec::HabQSpec; K = 100, save = false, base_dir = ".")
    dom = GriddedFisheryDomain()
    hab_fn = prep_file(spec)
    if !isfile(hab_fn)
        hab = prepare_habitat(spec, save = save, base_dir = base_dir)
    else
        hab = get_habitat(spec)
    end
    move = make_moveop(hab, sim_value(spec), dom)
    init_pop = eqdist(move, K)

    HabQPrep(dom, hab, move, init_pop)
end

function prepare_habitat(simtype::Type{<:HabQSpec}, realization::Integer; save = true, base_dir = ".")
    dom = GriddedFisheryDomain()
    hab = generate_habitat(dom)
    if save
        fn = prep_path(simtype, realization)
        h5open(fn, "w") do fid
            dom_size = size(hab)
            ghab = create_dataset(fid, "general_hab", datatype(Float64), dataspace(dom_size...))
            ghab[:, :] = hab[1]
            rhab = create_dataset(fid, "rocky_hab", datatype(Bool), dataspace(dom_size...))
            rhab[:, :] = hab[2]
        end
    end
    hab
end
function prepare_habitat(spec::HabQSpec; save = true, base_dir = ".")
    prepare_habitat(typeof(spec), realization(spec); save = save, base_dir = base_dir)
end

function get_habitat(spec::HabQSpec)
    prep_fn = prep_file(spec)
    fid = h5open(prep_fn, "r")
    ghab = read_dataset(fid, "general_hab")
    rhab = read_dataset(fid, "rocky_hab")
    close(fid)
    Habitat(ghab, rhab)
end

#-------------------------------------------------------------------------------

function survey_vessel(spec::HabQSpec, prep::HabQPrep)
    dom = domain(prep)
    hab = habitat(prep)
    q = HabitatCatchability(hab, base_catchability().catchability, gh -> 1.0, rh -> rh ? 0.1 : 1.0)
    q_real = getfield.(getindex.(Ref(q), eachindex(dom)), :catchability)
    survey_q = Catchability(reshape(q_real, size(dom)...))
    Vessel(survey_targeting(spec), survey_q, tweedie_shape(spec), tweedie_dispersion(spec))
end

function comm_vessel(spec::HabQSpec, prep::HabQPrep)
    dom = domain(prep)
    hab = habitat(prep)
    q = HabitatCatchability(hab, base_catchability().catchability, gh -> 1.0, rh -> rh ? 0.9 : 1.0)
    q_real = getfield.(getindex.(Ref(q), eachindex(dom)), :catchability)
    comm_q = Catchability(reshape(q_real, size(dom)...))
    comm_target = DynamicPreferentialTargeting(init_pop(prep).P,
                                               p -> comm_q.catchability .* p)
    Vessel(comm_target, comm_q, tweedie_shape(spec), tweedie_dispersion(spec))
end

#-------------------------------------------------------------------------------

function run_sims(simtype::Type{HabQSpec},
                  repl_range::AbstractRange{<:Integer};
                  base_dir = ".",
                  checkpoint = true,
                  csv = false,
                  feather = true)
    for rlz in repl_range
        run_sims(simtype,
                 rlz;
                 base_dir = base_dir,
                 checkpoint = checkpoint,
                 csv = csv,
                 feather = feather)
    end
    nothing
end

function complete_results(spec; csv = false, feather = true)
    files = file_paths(spec)
    complete = isfile(files[1])
    if csv
        complete &= all(isfile.(files[2:3]))
    end
    if feather
        complete &= all(isfile.(files[4:5]))
    end
    complete
end

function run_sims(simtype::Type{HabQSpec},
                  rlz::Integer;
                  base_dir = ".",
                  checkpoint = true,
                  csv = false,
                  feather = true)
    study_range = sim_values(simtype)
    for val in study_range
        spec = simtype(rlz, val, prep_path(simtype, rlz))
        files = file_paths(spec)

        # If checkpointing, skip the rest of the inner for loop body if the
        # results files are present.
        if (checkpoint && complete_results(spec; csv = csv, feather = feather))
            continue
        end

        prep = HabQPrep(spec; K = 100, save = true, base_dir = base_dir)
        setup = SpatQSimSetup(spec, prep)
        result = simulate(setup)
        save(result; csv = csv, feather = feather)
    end
    nothing
end


text:
@testset "periodic hard-spheres" begin

ϵ   = 0.1

q = [[0.5,  0.8], [0.8,  0.5]]
v = [[0.0, -1.0], [-1.0, 0.0]]

p  = HardSpheres( q, v, ϵ )

pc = PeriodicCollisions( p )

@show step!( p, pc )
@show p.v
@show step!( p, pc )
@show p.q
@show p.v

end