text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_sastre_m6_opt_rho2_22(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_sastre_m6_opt_rho2_22!(A_copy)
end

@inline function exp_sastre_m6_opt_rho2_22!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=9
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 T2k9
    # Computing Bb2 = x*I+x*A
    coeff1=0.002597101166550988
    coeff2=0.9907226316794635
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.002597101166550988
    coeff2=0.9907226316794635
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-0.005681896024400657
    coeff2=0.023301056977275802
    coeff3=0.999729943278336
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=-0.011734927258548565
    coeff2=0.000831101913142186
    coeff3=9.825887307216484e-5
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.0371131949784736
    coeff2=0.4122623848215059
    coeff3=0.00982117020658891
    coeff4=0.9691387395826522
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.0993170088559773
    coeff2=0.1075843791249758
    coeff3=0.03444785392998929
    coeff4=1.0058330089499217
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.000852548961418871
    coeff2=0.012977435119527533
    coeff3=0.9913750471409082
    coeff4=0.007013193540886666
    coeff5=-0.21625812151699
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.008832765141592277
    coeff2=0.0022901357546005413
    coeff3=0.0003677844696610227
    coeff4=0.008680829603225283
    coeff5=0.0006838601561276326
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.07225676735259205
    coeff2=0.40637944910599555
    coeff3=0.044422614636220305
    coeff4=0.09025415645822
    coeff5=0.1540082683372147
    coeff6=1.0657896264887965
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.06348900975427792
    coeff2=0.10321951688742986
    coeff3=0.04045126277415917
    coeff4=0.028287777007324065
    coeff5=0.13721732464684955
    coeff6=1.0702928736641044
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=3.144403235434427
    coeff2=0.5115991466872888
    coeff3=0.1652954236470456
    coeff4=7.924897961202367
    coeff5=0.9723882198958138
    coeff6=0.0692074219098335
    coeff7=0.009097568257407438
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.0055528225167482785
    coeff2=0.3634011920588302
    coeff3=0.09318615672033989
    coeff4=-0.06811866746595016
    coeff5=0.16962731509782766
    coeff6=1.85596241477658
    coeff7=1.0800725044535462
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing T2k9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.9658865558277218
    coeff2=0.021151469113934267
    coeff3=-0.017762277056814803
    coeff4=-0.0005931703373916548
    coeff5=0.09923354639251344
    coeff6=0.005656023271494485
    coeff7=0.0161502730596768
    coeff8=0.9681060677762616
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    return memslots1 # Returning T2k9
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
### A Pluto.jl notebook ###
# v0.12.17

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 85d5a830-4116-11eb-1479-630b53b50e94
begin
	using Pkg
    Pkg.activate(pwd())
end

# ╔═╡ 4862c2d0-4116-11eb-3f26-1bdefcb9ce3d
begin
	using UCIData
	using PlutoUI
	using DataFrames
	using Statistics
	using LinearAlgebra
	using Plots
	using Clustering
	using Distances
end

# ╔═╡ c73d5420-4126-11eb-19b3-bf6c04aef086
md"""
# Data Mining
## Assignment 2 - K-means Clustering
### Aadam & Obaidullah (CS1945 & CS1947)
"""

# ╔═╡ 360ab320-411d-11eb-2d11-c9fe7b3137d1
plotly()

# ╔═╡ 73936e9e-4116-11eb-0594-1d202bd073cd
function Kmeans(X, k; max_iters = 300, tol = 1e-5)
    # Reshape 2D array to a 1D array with length of all training examples
    X_array_list = collect(eachrow(X))

    # Save some info on the incoming data
    N = length(X_array_list)  # Length of all training examples
    n = length(X_array_list[1])  # Length of a single training example
    distances = zeros(N)  # Empty vector for all training examples.

    # Step 1: Random initialization
    reps_centroids = [zeros(n) for grp = 1:k]  # Initiate centroids for each
    labels = rand(1:k, N)  # Randomly assign labels to all training examples

    J_previous = Inf

    for iter = 1:max_iters

        # Step 2: Update the representative centroids for each group
        for j = 1:k
            # get group indices for each group
            group_idx = [i for i = 1:N if labels[i] == j]

            # use group indices to locate each group
            reps_centroids[j] = mean(X_array_list[group_idx]);
        end;

        # Step 3: Update the group labels
        for i = 1:N
            # compute the distance between each example and the updated centroid
            nearest_rep_distance =
			[norm(X_array_list[i] - reps_centroids[x]) for x = 1:k]

            # update distances and label arrays
            # findmin returns the min value & index location
            distances[i], labels[i] = findmin(nearest_rep_distance)
        end;

        # Step 4: Compute the clustering cost
        J = (norm(distances)^ 2) / N

        # Show progress and terminate if J stopped decreasing.
        println("Iteration ", iter, ": Jclust = ", J, ".")

        # Final Step 5: Check for convergence
        if iter > 1 && abs(J - J_previous) < (tol * J)
            # Terminate algorithm with the assumption that K-means has converged
            return labels, reps_centroids
        end

        J_previous = J
    end
	return labels, reps_centroids
end

# ╔═╡ aac40900-41cc-11eb-1103-97d9373bdd65
function calculate_average_distance(i, X, result)
	dist = 0
	center = result.centers[i]
	for instance in X[:, result.assignments .== i]
		dist += Euclidean()(instance, center)
	end
	return dist/result.counts[i]
end

# ╔═╡ 12615c80-41c7-11eb-2aed-c1eb96003ff1
function dbi(X, result)
	num_of_clusters = length(result.counts)

	#average distance of all points in a cluster to its center, for each cluster
    average_centroid_distances = [
		calculate_average_distance(i, X, result) for i in 1:num_of_clusters]

	sum_Dij = 0

	#calculate the sum Dij
    for i in 1:num_of_clusters
        max_Dij = 0

        #calculate the max Dij
        for j in 1:num_of_clusters
            if i != j

				#calculate di + dj
				numerator =
				average_centroid_distances[i] + average_centroid_distances[j]

				#calculate dij
				denominator =
				Euclidean()( result.centers[i], result.centers[j] )

				current_Dij = numerator/denominator

				if max_Dij < current_Dij
					max_Dij = current_Dij
				end
			end

		end
		sum_Dij += max_Dij
	end

    return sum_Dij / num_of_clusters
end

# ╔═╡ c7001000-42c6-11eb-23d1-314932816f0a
function big_delta(maskci, distances)
	maximum(distances[maskci, :][:, maskci])
end

# ╔═╡ c6ca5a00-42c6-11eb-386e-2f13a99d75b6
function delta(maskci, maskcj, distances)
	minimum(distances[maskci, :][:, maskcj])
end

# ╔═╡ c69d2f80-42c6-11eb-10cf-e70c07e40e97
function dunn(X, result)
	n = length(result.counts) # number of clusters

	distances = pairwise(Euclidean(), X, dims=2)

	deltas = ones(n, n) * Inf
	big_deltas = zeros(n)

	for k in 1:n
		for l in 1:n
			if k == l
				continue
			end
            deltas[k, l] = delta(
				result.assignments .== k,
				result.assignments .== l,
				distances)
		end

        big_deltas[k] = big_delta(result.assignments .== k, distances)

	end
	minimum(deltas)/maximum(big_deltas)
end

# ╔═╡ 9accd320-4117-11eb-0fc4-15e8e85040a8
datasets = ["iris" => "iris", "glass-identification" => "glass"];

# ╔═╡ 8b43a9a0-4118-11eb-127d-d57dcea189e1
md"""
Dataset: $(@bind selected_dataset Select(datasets))

"""

# ╔═╡ 73c57b20-4116-11eb-3dba-01c9f7b91fa7
data = UCIData.dataset(selected_dataset);

# ╔═╡ 5cddc6a0-411c-11eb-2380-1b9442938a9b
features = collect(Matrix(data[:, 2:end-1])');

# ╔═╡ bec93fa0-4119-11eb-026e-0be14d6a5779
md"""
Vlaue of **k**: $(@bind k Slider(2:15; show_value=true))
"""
# size(UCIData.dataset(selected_dataset))[1]

# ╔═╡ 737b2bb0-4116-11eb-124b-d981032a9638
result = kmeans(features, k);

# ╔═╡ 725d2e10-42c2-11eb-3ee7-4ba089ebb440
md"""
### Davies Bouldin Index: $(round(dbi(features, result), digits=5))
### Dunn Index : $(round(dunn(features, result), digits=5))
"""

# ╔═╡ ce4134a0-4124-11eb-2a7c-9761e81b2064
md"""
## Plot Clusters
Feature $(@bind feature1 Select(string.(collect(1:size(features)[1])))) vs.
Feature $(@bind feature2 Select(string.(collect(1:size(features)[1])));)
"""

# ╔═╡ 73454ea0-4116-11eb-0938-f9a06f3dbaf8
begin
	theme(:dark)
	scatter(
		features[parse(Int,feature1), :],
		features[parse(Int, feature2), :],
		marker_z = result.assignments,
		markershape = :diamond,
		markersize = 3,
		xlabel="Feature $(feature1)",
		ylabel="Feature $(feature2)",
		color = :seaborn_bright,
		legend = false)
	scatter!(
		result.centers[parse(Int,feature1), :],
		result.centers[parse(Int, feature2), :],
		markersize = 6,
		marker_z = result.centers',
		color = :temperaturemap
	)
end

# ╔═╡ 6cd3fe00-41c4-11eb-2e84-376a14af7398
md"""
### Cluster Centers
$(
df = DataFrame(hcat(result.centers', result.counts));
last_col_name = names(df)[end];
df[!,last_col_name] = convert.(Int,df[:,last_col_name]);
rename!(df, last_col_name => :NodeCount)
)
"""

# ╔═╡ Cell order:
# ╟─c73d5420-4126-11eb-19b3-bf6c04aef086
# ╟─85d5a830-4116-11eb-1479-630b53b50e94
# ╠═4862c2d0-4116-11eb-3f26-1bdefcb9ce3d
# ╟─360ab320-411d-11eb-2d11-c9fe7b3137d1
# ╠═73936e9e-4116-11eb-0594-1d202bd073cd
# ╠═aac40900-41cc-11eb-1103-97d9373bdd65
# ╠═12615c80-41c7-11eb-2aed-c1eb96003ff1
# ╠═c7001000-42c6-11eb-23d1-314932816f0a
# ╠═c6ca5a00-42c6-11eb-386e-2f13a99d75b6
# ╠═c69d2f80-42c6-11eb-10cf-e70c07e40e97
# ╟─9accd320-4117-11eb-0fc4-15e8e85040a8
# ╟─8b43a9a0-4118-11eb-127d-d57dcea189e1
# ╠═73c57b20-4116-11eb-3dba-01c9f7b91fa7
# ╠═5cddc6a0-411c-11eb-2380-1b9442938a9b
# ╟─bec93fa0-4119-11eb-026e-0be14d6a5779
# ╠═737b2bb0-4116-11eb-124b-d981032a9638
# ╠═725d2e10-42c2-11eb-3ee7-4ba089ebb440
# ╟─ce4134a0-4124-11eb-2a7c-9761e81b2064
# ╟─73454ea0-4116-11eb-0938-f9a06f3dbaf8
# ╟─6cd3fe00-41c4-11eb-2e84-376a14af7398

ext: jl


text:
using ReachabilityAnalysis

const μ = 3.986e14 * 60^2
const r = 42164.0e3
const r² = r^2
const mc = 500.0
const n² = μ / r^3
const n = sqrt(n²)

const two_n = 2*n
const μ_r² = μ/r²;

const K₁ = [-28.8287 0.1005 -1449.9754 0.0046 0.0;
            -0.087 -33.2562 0.00462 -1451.5013 0.0]
const K₂ = [-288.0288 0.1312 -9614.9898 0.0 0.0;
            -0.1312 -288.0 0.0 -9614.9883   0.0]

const K₁mc = K₁/mc
const K₂mc = K₂/mc

function mymul!(v, A, x)
    @inbounds for ind = 1:length(v)
        v[ind] = zero(x[1])
        for jind = 1:length(x)
            v[ind] += A[ind, jind] * x[jind]
        end
    end
    return nothing
end

@taylorize function spacecraft_approaching!(du, u, p, t)
    x, y, vx, vy, t = u

    rx = r + x
    rx² = rx^2
    y² = y^2
    rc = sqrt(rx² + y²)
    rc³ = rc^3
    μ_rc³ = μ / rc³

    uxy = Vector{typeof(x)}(undef, 2)
    mymul!(uxy, K₁mc, u)

    # x' = vx
    du[1] = vx

    # y' = vy
    du[2] = vy

    # vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x) + ux/mc
    du[3] = (n²*x + two_n*vy) + ((μ_r² - μ_rc³*rx) + uxy[1])

    # vy' = n²y - 2n*vx - μ/(rc^3)y + uy/mc
    du[4] = (n²*y - two_n*vx) - (μ_rc³*y - uxy[2])

    # t' = 1
    du[5] = one(x)

    return du
end

@taylorize function spacecraft_attempt!(du, u, p, t)
    x, y, vx, vy, t = u

    rx = r + x
    rx² = rx^2
    y² = y^2
    rc = sqrt(rx² + y²)
    rc³ = rc^3
    μ_rc³ = μ / rc³

    uxy = Vector{typeof(x)}(undef, 2)
    mymul!(uxy, K₂mc, u)

    # x' = vx
    du[1] = vx

    # y' = vy
    du[2] = vy

    # vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x) + ux/mc
    du[3] = (n²*x + two_n*vy) + ((μ_r² - μ_rc³*rx) + uxy[1])

    # vy' = n²y - 2n*vx - μ/(rc^3)y + uy/mc
    du[4] = (n²*y - two_n*vx) - (μ_rc³*y - uxy[2])

    # t' = 1
    du[5] = one(x)

    return du
end

@taylorize function spacecraft_aborting!(du, u, p, t)
    x, y, vx, vy, t = u

    rx = r + x
    rx² = rx^2
    y² = y^2
    rc = sqrt(rx² + y²)
    rc³ = rc^3
    μ_rc³ = μ / rc³

    # x' = vx
    du[1] = vx

    # y' = vy
    du[2] = vy

    # vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x)
    du[3] = (n²*x + two_n*vy) + (μ_r² - μ_rc³*rx)

    # vy' = n²y - 2n*vx - μ/(rc^3)y
    du[4] = (n²*y - two_n*vx) - μ_rc³*y

    # t' = 1
    du[5] = one(x)

    return du
end

using Symbolics

const var = @variables x y vx vy t

function spacecraft(; abort_time=(120.0, 150.0))

    n = 4 + 1  # number of variables
    t_abort_lower, t_abort_upper = abort_time

    automaton = LightAutomaton(3)

    # mode 1 "approaching"
    invariant = HalfSpace(x <= -100, var)
    approaching = @system(x' = spacecraft_approaching!(x), dim:5, x ∈ invariant)

    # mode 2 "rendezvous attempt"
    invariant = HalfSpace(x >= -100, var)
    attempt = @system(x' = spacecraft_attempt!(x), dim:5, x ∈ invariant)

    # mode 3 "aborting"
    invariant = Universe(n)
    aborting = @system(x' = spacecraft_aborting!(x), dim: 5, x ∈ invariant)

    # transition "approach" -> "attempt"
    add_transition!(automaton, 1, 2, 1)
    guard = HalfSpace(x >= -100, var)
    t1 = @map(x -> x, dim: n, x ∈ guard)

    # transition "approach" -> "abort"
    add_transition!(automaton, 1, 3, 2)
    guard_time = HPolyhedron([t >= t_abort_lower, t <= t_abort_upper], var)
    t2 = @map(x -> x, dim: n, x ∈ guard_time)

    # transition "attempt" -> "abort"
    add_transition!(automaton, 2, 3, 3)
    t3 = @map(x -> x, dim: n, x ∈ guard_time)

    H = HybridSystem(automaton=automaton,
                     modes=[approaching, attempt, aborting],
                     resetmaps=[t1, t2, t3])

     # initial condition in mode 1
     X0 = Hyperrectangle([-900., -400., 0., 0., 0.], [25., 25., 0., 0., 0.])
     init = [(1, X0)]

     return InitialValueProblem(H, init)
end

const tan30 = tand(30)

LineOfSightCone = HPolyhedron([x >= -100, y >= x*tan30, -y >= x*tan30], var)

Target = Hyperrectangle(zeros(2), [0.2, 0.2]);

function line_of_sight(sol)
    all_idx = findall(x -> x == 2, location.(sol))  # attempt
    for idx in all_idx
        verif = all(set(R) ⊆ LineOfSightCone for R in sol[idx][2:end])
        !verif && return false
    end
    return true
end

function absolute_velocity(R)
    vx, vy = 3, 4
    vx2 = set(overapproximate(project(R, vars=(vx,)), Interval))
    vy2 = set(overapproximate(project(R, vars=(vy,)), Interval))
    return sqrt(vx2.dat + vy2.dat)
end

function velocity_constraint(sol)
    all_idx = findall(x -> x == 2, location.(sol)) # attempt mode
    for idx in all_idx
        # maximum velocity measured in m/min
        verif = all(absolute_velocity(R) < 0.055 * 60. for R in sol[idx])
        !verif && return false
    end
    return true
end

function target_avoidance(sol)
    all_idx = findall(x -> x == 3, location.(sol)) # aborting mode
    for idx in all_idx
        verif = all(is_intersection_empty(set(Projection(R, [x, y])), Target) for R in sol[idx])
        !verif && return false
    end
    return true
end

using Plots

function solve_spacecraft(prob; k=25, s=missing)

    # template directions
    boxdirs = BoxDirections(5)

    # transition from mode 1 to mode 2
    sol12 = solve(prob,
                tspan=(0.0, 200.0),
                alg=TMJets20(abstol=1e-5, maxsteps=10_000, orderT=5, orderQ=1, disjointness=BoxEnclosure()),
                max_jumps=1,
                intersect_source_invariant=false,
                intersection_method=TemplateHullIntersection(boxdirs),
                clustering_method=LazyClustering(1),
                disjointness_method=BoxEnclosure())

    sol12jump = overapproximate(sol12[2](120 .. 150), Zonotope)
    t0 = tstart(sol12jump[1])
    sol12jump_c = cluster(sol12jump, 1:length(sol12jump), BoxClustering(k, s))

    # transition from mode 2 to mode 3
    H = system(prob)
    sol3 = solve(IVP(mode(H, 3), [set(X) for X in sol12jump_c]),
                 tspan=(t0, 200.0),
                 alg=TMJets20(abstol=1e-10, orderT=7, orderQ=1, disjointness=BoxEnclosure()))
    d = Dict{Symbol, Any}(:loc_id => 3)

    F12 = [fp for fp in sol12.F]
    F23 = [Flowpipe(fp.Xk, d) for fp in sol3.F]
    return HybridFlowpipe(vcat(F12, F23))
end

prob = spacecraft()
sol = solve_spacecraft(prob; k=25, s=missing)
solz = overapproximate(sol, Zonotope);

plot(solz, vars=(1, 2), xlab="x", ylab="y")

idx_approaching = findall(x -> x == 1, location.(solz))
idx_attempt = findall(x -> x == 2, location.(solz))
idx_aborting = findall(x -> x == 3, location.(solz))

using Plots, Plots.PlotMeasures, LaTeXStrings

fig = plot(legend=:bottomright, tickfont=font(10, "Times"), guidefontsize=15,
           xlab=L"x", ylab=L"y", lw=0.0,
           xtick=[-750, -500, -250, 0, 250.], ytick=[-400, -300, -200, -100, 0.],
           xlims=(-1000.0, 400.0), ylims=(-450.0, 0.0), size=(600, 600),
           bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm)

plot!(fig, solz[idx_approaching[1]], vars=(1, 2), lw=0.0, color=:lightgreen, alpha=1., lab="Approaching")
for k in 2:length(idx_approaching)
    plot!(fig, solz[idx_approaching[k]], vars=(1, 2), lw=0.0, color=:lightgreen, alpha=1.)
end

plot!(fig, solz[idx_attempt[1]], vars=(1, 2), lw=0.0, color=:red, alpha=1., lab="Attempt")
for k in 2:length(idx_attempt)
    plot!(fig, solz[idx_attempt[k]], vars=(1, 2), lw=0.0, color=:red, alpha=1.)
end

plot!(fig, solz[idx_aborting[1]], vars=(1, 2), lw=0.0, color=:cyan, alpha=1., lab="Aborting")
for k in 2:length(idx_aborting)
    plot!(fig, solz[idx_aborting[k]], vars=(1, 2), lw=0.0, color=:cyan, alpha=1.)
end

fig


text:
using .Curve: BN
using .Config: ORDER, @ID

struct LagrangeCoeffGenerator
    coeffs::Dict{@ID,BN}
    order::BN
    LagrangeCoeffGenerator() = new(Dict{@ID,BN}(), ORDER)
    # calculate the weights ahead of time
    function LagrangeCoeffGenerator(keys)
        lc = LagrangeCoeffGenerator()
        for x in keys
            push!(lc, x);
        end
        return lc
    end
end

Base.length(lc::LagrangeCoeffGenerator) = length(lc.coeffs)
Base.getindex(lc::LagrangeCoeffGenerator, i::@ID) = getindex(lc.coeffs, i)

# NOTE: not thread safe
function Base.push!(lc::LagrangeCoeffGenerator, i::@ID)
    coeff = get(lc.coeffs, i, nothing)
    if coeff === nothing
        # update existing coeffs
        bi = BN(i) # move to the field OUTSIDE of the loop
        for (k, c) in lc.coeffs
            lc.coeffs[k] = mod(c * invmod(i - k, lc.order) * bi, lc.order)
        end
        # update new coeff
        coeff = one(BN)
        for (k, _) in lc.coeffs
            coeff = mod(coeff * invmod(k - i, lc.order) * k, lc.order)
        end
        lc.coeffs[i] = coeff
    end
    return coeff
end

function Base.delete!(lc::LagrangeCoeffGenerator, i::@ID)
    if length(lc.coeffs) != length(delete!(lc.coeffs, i))
        # update existing coeffs
        invi = invmod(i, lc.order)
        for (k, c) in lc.coeffs
            lc.coeffs[k] = mod(c * invi * (i - k), lc.order)
        end
    end
end


text:
"""
Along same lines as  `@test_throws` but to test if a macro throw an exception when it is
expanded.
"""
macro test_macro_throws(err_expr, expr)
    quote
        err = nothing
        try
            @macroexpand($(esc(expr)))
        catch load_err
            # all errors thrown at macro expansion time are LoadErrors, we need to unwrap
            @assert load_err isa LoadError
            err = load_err.error
        end
        # Reuse `@test_throws` logic
        if err!==nothing
            @test_throws $(esc(err_expr)) ($(Meta.quot(expr)); throw(err))
        else
            @test_throws $(esc(err_expr)) $(Meta.quot(expr))
        end
    end
end

# struct need to be defined outside of tests for julia 1.0 compat
struct NonDiffExample
    x
end

struct NonDiffCounterExample
    x
end

@testset "rule_definition_tools.jl" begin
    @testset "@non_differentiable" begin
        @testset "two input one output function" begin
            nondiff_2_1(x, y) = fill(7.5, 100)[x + y]
            @non_differentiable nondiff_2_1(::Any, ::Any)
            @test frule((Zero(), 1.2, 2.3), nondiff_2_1, 3, 2) == (7.5, DoesNotExist())
            res, pullback = rrule(nondiff_2_1, 3, 2)
            @test res == 7.5
            @test pullback(4.5) == (NO_FIELDS, DoesNotExist(), DoesNotExist())
        end

        @testset "one input, 2-tuple output function" begin
            nondiff_1_2(x) = (5.0, 3.0)
            @non_differentiable nondiff_1_2(::Any)
            @test frule((Zero(), 1.2), nondiff_1_2, 3.1) == ((5.0, 3.0), DoesNotExist())
            res, pullback = rrule(nondiff_1_2, 3.1)
            @test res == (5.0, 3.0)
            @test isequal(
                pullback(Composite{Tuple{Float64, Float64}}(1.2, 3.2)),
                (NO_FIELDS, DoesNotExist()),
            )
        end

        @testset "constrained signature" begin
            nonembed_identity(x) = x
            @non_differentiable nonembed_identity(::Integer)

            @test frule((Zero(), 1.2), nonembed_identity, 2) == (2, DoesNotExist())
            @test frule((Zero(), 1.2), nonembed_identity, 2.0) == nothing

            res, pullback = rrule(nonembed_identity, 2)
            @test res == 2
            @test pullback(1.2) == (NO_FIELDS, DoesNotExist())

            @test rrule(nonembed_identity, 2.0) == nothing
        end

        @testset "Pointy UnionAll constraints" begin
            pointy_identity(x) = x
            @non_differentiable pointy_identity(::Vector{<:AbstractString})

            @test frule((Zero(), 1.2), pointy_identity, ["2"]) == (["2"], DoesNotExist())
            @test frule((Zero(), 1.2), pointy_identity, 2.0) == nothing

            res, pullback = rrule(pointy_identity, ["2"])
            @test res == ["2"]
            @test pullback(1.2) == (NO_FIELDS, DoesNotExist())

            @test rrule(pointy_identity, 2.0) == nothing
        end

        @testset "kwargs" begin
            kw_demo(x; kw=2.0) = x + kw
            @non_differentiable kw_demo(::Any)

            @testset "not setting kw" begin
                @assert kw_demo(1.5) == 3.5

                res, pullback = rrule(kw_demo, 1.5)
                @test res == 3.5
                @test pullback(4.1) == (NO_FIELDS, DoesNotExist())

                @test frule((Zero(), 11.1), kw_demo, 1.5) == (3.5, DoesNotExist())
            end

            @testset "setting kw" begin
                @assert kw_demo(1.5; kw=3.0) == 4.5

                res, pullback = rrule(kw_demo, 1.5; kw=3.0)
                @test res == 4.5
                @test pullback(1.1) == (NO_FIELDS, DoesNotExist())

                @test frule((Zero(), 11.1), kw_demo, 1.5; kw=3.0) == (4.5, DoesNotExist())
            end
        end

        @testset "Constructors" begin
            @non_differentiable NonDiffExample(::Any)

            @test isequal(
                frule((Zero(), 1.2), NonDiffExample, 2.0),
                (NonDiffExample(2.0), DoesNotExist())
            )

            res, pullback = rrule(NonDiffExample, 2.0)
            @test res == NonDiffExample(2.0)
            @test pullback(1.2) == (NO_FIELDS, DoesNotExist())

            # https://github.com/JuliaDiff/ChainRulesCore.jl/issues/213
            # problem was that `@nondiff Foo(x)` was also defining rules for other types.
            # make sure that isn't happenning
            @test frule((Zero(), 1.2), NonDiffCounterExample, 2.0) === nothing
            @test rrule(NonDiffCounterExample, 2.0) === nothing
        end

        @testset "Varargs" begin
            fvarargs(a, xs...) = sum((a, xs...))
            @testset "xs::Float64..." begin
                @non_differentiable fvarargs(a, xs::Float64...)

                y, pb = rrule(fvarargs, 1)
                @test y == fvarargs(1)
                @test pb(1) == (Zero(), DoesNotExist())

                y, pb = rrule(fvarargs, 1, 2.0, 3.0)
                @test y == fvarargs(1, 2.0, 3.0)
                @test pb(1) == (Zero(), DoesNotExist(), DoesNotExist(), DoesNotExist())

                @test frule((1, 1), fvarargs, 1, 2.0) == (fvarargs(1, 2.0), DoesNotExist())

                @test frule((1, 1), fvarargs, 1, 2) == nothing
                    @test rrule(fvarargs, 1, 2) == nothing
            end

            @testset "::Float64..." begin
                @non_differentiable fvarargs(a, ::Float64...)

                y, pb = rrule(fvarargs, 1, 2.0, 3.0)
                @test y == fvarargs(1, 2.0, 3.0)
                @test pb(1) == (Zero(), DoesNotExist(), DoesNotExist(), DoesNotExist())

                @test frule((1, 1), fvarargs, 1, 2.0) == (fvarargs(1, 2.0), DoesNotExist())
            end

            @testset "::Vararg{Float64}" begin
                @non_differentiable fvarargs(a, ::Vararg{Float64})

                y, pb = rrule(fvarargs, 1, 2.0, 3.0)
                @test y == fvarargs(1, 2.0, 3.0)
                @test pb(1) == (Zero(), DoesNotExist(), DoesNotExist(), DoesNotExist())

                @test frule((1, 1), fvarargs, 1, 2.0) == (fvarargs(1, 2.0), DoesNotExist())
            end

            @testset "::Vararg" begin
                @non_differentiable fvarargs(a, ::Vararg)
                @test frule((1, 1), fvarargs, 1, 2) == (fvarargs(1, 2), DoesNotExist())

                y, pb = rrule(fvarargs, 1, 1)
                @test y == fvarargs(1, 1)
                @test pb(1) == (Zero(), DoesNotExist(), DoesNotExist())
            end

            @testset "xs..." begin
                @non_differentiable fvarargs(a, xs...)
                @test frule((1, 1), fvarargs, 1, 2) == (fvarargs(1, 2), DoesNotExist())

                y, pb = rrule(fvarargs, 1, 1)
                @test y == fvarargs(1, 1)
                @test pb(1) == (Zero(), DoesNotExist(), DoesNotExist())
            end
        end

        @testset "Not supported (Yet)" begin
            # Where clauses are not supported.
            @test_macro_throws(
                ErrorException,
                (@non_differentiable where_identity(::Vector{T}) where T<:AbstractString)
            )
        end
    end

    @testset "@scalar_rule" begin
        @testset "@scalar_rule with multiple output" begin
            simo(x) = (x, 2x)
            @scalar_rule(simo(x), 1f0, 2f0)

            y, simo_pb = rrule(simo, π)

            @test simo_pb((10f0, 20f0)) == (NO_FIELDS, 50f0)

            y, ẏ = frule((NO_FIELDS, 50f0), simo, π)
            @test y == (π, 2π)
            @test ẏ == Composite{typeof(y)}(50f0, 100f0)
            # make sure type is exactly as expected:
            @test ẏ isa Composite{Tuple{Irrational{:π}, Float64}, Tuple{Float32, Float32}}
        end
    end
end


text:
@testset "InputShapeVertex" begin
    import NaiveNASflux: FluxDense, FluxConv, FluxRnn
    import NaiveNASflux: actrank, actdim
    v = inputvertex("in", 3, FluxDense())
    @test layertype(v) == FluxDense()
    @test name(v) == "in"
    @test nout(v) == 3

    @testset "Copy with $label" for (label, cfun) in (
        (deepcopy, deepcopy),
        ("fmap", g -> fmap(identity, g))
    )
        c = cfun(v)
        @test layertype(c) == FluxDense()
        @test name(c) == "in"
        @test nout(c) == 3
    end

    @test nout(inputvertex("in", 4, FluxDense()))   == nout(denseinputvertex("in", 4))
    @test nout(inputvertex("in", 3, FluxConv{1}())) == nout(conv1dinputvertex("in", 3))
    @test nout(inputvertex("in", 3, FluxConv{2}())) == nout(conv2dinputvertex("in", 3))
    @test nout(inputvertex("in", 3, FluxConv{3}())) == nout(conv3dinputvertex("in", 3))
    @test nout(inputvertex("in", 2, FluxRnn()))     == nout(rnninputvertex("in", 2))

    @test actdim(inputvertex("in", 4, FluxDense()))   == actdim(denseinputvertex("in", 4))
    @test actdim(inputvertex("in", 3, FluxConv{2}())) == actdim(convinputvertex("in", 3, 2))
    @test actdim(inputvertex("in", 3, FluxConv{1}())) == actdim(conv1dinputvertex("in", 3))
    @test actdim(inputvertex("in", 3, FluxConv{2}())) == actdim(conv2dinputvertex("in", 3))
    @test actdim(inputvertex("in", 3, FluxConv{3}())) == actdim(conv3dinputvertex("in", 3))
    @test actdim(inputvertex("in", 2, FluxRnn()))     == actdim(rnninputvertex("in", 2))

    @test actrank(inputvertex("in", 4, FluxDense()))   == actrank(denseinputvertex("in", 4))
    @test actrank(inputvertex("in", 3, FluxConv{2}())) == actrank(convinputvertex("in", 3, 2))
    @test actrank(inputvertex("in", 3, FluxConv{1}())) == actrank(conv1dinputvertex("in", 3))
    @test actrank(inputvertex("in", 3, FluxConv{2}())) == actrank(conv2dinputvertex("in", 3))
    @test actrank(inputvertex("in", 3, FluxConv{3}())) == actrank(conv3dinputvertex("in", 3))
    @test actrank(inputvertex("in", 2, FluxRnn()))     == actrank(rnninputvertex("in", 2))
end

@testset "$pt Vertex" for (pt, arg) in (
    (MaxPool, ((4,4),)),
    (MeanPool, ((4,4),)),
    (GlobalMaxPool, ()),
    (GlobalMeanPool, ()))
    using Random
    l = pt(arg...)
    v = fluxvertex(l, conv2dinputvertex("in", 3))

    @test NaiveNASflux.layertype(v) isa NaiveNASflux.FluxPoolLayer
    indata = randn(MersenneTwister(1), 6,6,3,2)

    @test l(indata) == v(indata)
end

@testset "$dt Vertex" for dt in (Dropout, AlphaDropout)
    using Random
    l = dt(0.4)
    v = fluxvertex(l, conv2dinputvertex("in", 3))

    @test NaiveNASflux.layertype(v) isa NaiveNASflux.FluxDropOut
    indata = randn(MersenneTwister(1), 3,2)

    @test l(indata) == v(indata)
end

@testset "Size mutations" begin
    import NaiveNASflux: weights, bias, layertype

    @testset "Dense to Dense" begin
        inpt = inputvertex("in", 4)
        dl1 = Dense(4, 5)
        dl2 = Dense(5, 3)
        bias(dl1)[1:end] = 1:5
        bias(dl2)[1:end] = 1:3
        dense1 = fluxvertex(dl1, inpt)
        dense2 = fluxvertex(dl2, dense1)

        @test inputs(dense2) == [dense1]
        @test outputs(dense1) == [dense2]
        @test inputs(dense1) == [inpt]
        @test outputs(inpt) == [dense1]

        @test layer(dense1) == dl1
        @test layer(dense2) == dl2
        @test layertype(dense1) == FluxDense()
        @test layertype(dense2) == FluxDense()

        @test [nout(inpt)] == nin(dense1) == [4]
        @test [nout(dense1)] == nin(dense2) == [5]
        @test nout(dense2) == 3

        inds = Bool[1, 1, 0, 1, 1]
        NaiveNASlib.Δsize!(dense1) do v
            v == dense1 || return 1
            return inds .- 0.5
        end

        @test [nout(inpt)] == nin(dense1) == [4]
        @test [nout(dense1)] == nin(dense2) == [4]
        @test nout(dense2) == 3

        W2exp, b2exp = weights(dl2)[:, inds], bias(dl2)
        W1exp, b1exp = weights(dl1)[inds, :], bias(dl1)[inds]

        @test size(CompGraph([inpt], [dense2])(collect(Float32, 1:nout(inpt)))) == (3,)

        assertlayer(layer(dense2), W2exp, b2exp)
        assertlayer(layer(dense1), W1exp, b1exp)
    end

    @testset "Invariant parametric layer" begin
        inpt = inputvertex("in", 3)
        cv = fluxvertex(Conv((1,1), nout(inpt) => 4), inpt)
        # layerfun = identity or else LazyMutable will take case or the sizes for us
        bv = fluxvertex(BatchNorm(nout(cv)), cv; layerfun=identity)

        @test nin(bv) == [nout(cv)] == [4]

        Δnin!(v -> 1, bv, -1)

        @test nin(bv) == [nout(cv)] == [3]

        # Lets go out of our way to make it size inconsistent to verify that nin/nout is the size of the BatchNorm parameters.
        NaiveNASlib.applyΔsize!(NaiveNASlib.NeuronIndices(), bv, [1:2], 1:2)

        @test [nout(bv)] == nin(bv) == [2]
        @test nout(cv) == 3
    end

    @testset "Invariant non-parametric layer" begin
        inpt = inputvertex("in", 3)
        cv = fluxvertex(Conv((1,1), nout(inpt) => 4), inpt)
        bv = fluxvertex(MeanPool((2,2)), cv)

        @test nin(bv) == [nout(cv)] == [4]

        Δnin!(v -> 1, bv, -1)

        @test nin(bv) == [nout(cv)] == [3]
    end

    @testset "DepthwiseConv" begin
        import NaiveNASflux: outdim, wrapped, FluxConv
        lazymutable(v::AbstractVertex) = lazymutable(base(v))
        lazymutable(v::CompVertex) = lazymutable(v.computation)
        lazymutable(m::NaiveNASflux.AbstractMutableComp) = lazymutable(wrapped(m))
        lazymutable(m::LazyMutable) = m
        lazyouts(v) = lazymutable(v).outputs
        lazyins(v) = lazymutable(v).inputs

        @testset "Depthwise single layer" begin
            import NaiveNASflux: neuronutility
            # just to check that I have understood the wiring of the weight
            @testset "4 inputs times 2" begin
                inpt = inputvertex("in", 4, FluxConv{2}())
                dc = fluxvertex("dc", DepthwiseConv(reshape(Float32[10 10 10 10;20 20 20 20], 1, 1, 2, 4), Float32[0,0,0,0,1,1,1,1]), inpt)
                @test neuronutility(dc) == [20, 40, 20, 40, 21, 41, 21, 41]
                @test reshape(dc(fill(1f0, (1,1,4,1))), :) == [10, 20, 10, 20, 11, 21, 11, 21]
                @test Δnout!( dc => -4)
                @test lazyouts(dc) == [2, 4, 6, 8]
                @test reshape(dc(fill(1f0, (1,1,4,1))), :) == [20, 20, 21, 21]
                # pass insert function so we get zeros instead of randn
                @test Δsize!(WithKwargs(ΔNout(dc, 4); insert =(args...) -> (args...) -> 0))
                @test lazyouts(dc) ==                         [ 1,-1,  2,-1,  3,-1,  4,-1]
                @test reshape(dc(fill(1f0, (1,1,4,1))), :) == [20, 0, 20, 0, 21, 0, 21, 0]
            end

            @testset "2 inputs times 3" begin
                inpt = inputvertex("in", 2, FluxConv{2}())
                dc = fluxvertex("dc", DepthwiseConv(reshape(Float32[10 10;20 20;30 30], 1, 1, 3, 2), Float32[0,0,1,1,2,2]), inpt)
                @test reshape(dc(fill(1f0, (1,1,2,1))), :) == [10, 20, 31, 11, 22, 32]
                @test Δnout!(dc => -2)
                @test lazyouts(dc) == [2,3,5,6]
                @test reshape(dc(fill(1f0, (1,1,2,1))), :) == [20, 31, 22, 32]
                # pass insert function so we get zeros instead of randn
                @test Δsize!(WithKwargs(ΔNout(dc, 4); insert =(args...) -> (args...) -> 0))
                @test lazyouts(dc) ==                         [ 1,  2, -1,-1,  3,  4, -1,-1]
                @test reshape(dc(fill(1f0, (1,1,2,1))), :) == [20, 31,  0, 0, 22, 32,  0, 0]
            end

            @testset "1 input times 5" begin
                inpt = inputvertex("in", 1, FluxConv{2}())
                dc = fluxvertex("dc", DepthwiseConv(reshape(Float32.(10:10:50), 1, 1, 5, 1), Float32.(1:5)), inpt)
                @test reshape(dc(fill(1f0, (1,1,1,1))), :) == [11, 22, 33, 44, 55]
                @test Δnout!(dc=>-2)
                @test lazyouts(dc) == 3:5
                @test reshape(dc(fill(1f0, (1,1,1,1))), :) == [33, 44, 55]
                # pass insert function so we get zeros instead of randn
                @test Δsize!(WithKwargs(ΔNout(dc, 3); insert =(args...) -> (args...) -> 0))
                @test lazyouts(dc) ==                     vcat(1:3,        fill(-1, 3))
                @test reshape(dc(fill(1f0, (1,1,1,1))), :) == [33, 44, 55, 0, 0, 0]
            end

            @testset "3 inputs times 7" begin
                inpt = inputvertex("in", 3, FluxConv{2}())
                dc = fluxvertex("dc", DepthwiseConv(reshape(repeat(Float32.(10:10:70), 3), 1,1,7,3), Float32.(1:21)), inpt)
                @test reshape(dc(fill(10f0, (1,1,3,1))), :) == repeat(100:100:700, 3) .+ (1:21)
                @test Δnout!(dc => -9) do v
                    v == dc || return 1
                    val = ones(nout(v))
                    val[[2,13,14]] .= -10
                    val[1] = 10
                    return val
                end
                # TODO: Examine failure and expected result more carefully when not tired af
                @test lazyouts(dc) == [1,3,4,5,8,10,11,12,15,17,18,19]
                @test reshape(dc(fill(10f0, (1,1,3,1))), :) == [101,303,404,505,108,310,411,512,115,317,418,519]
                # pass insert function so we get zeros instead of randn
                @test Δsize!(WithKwargs(ΔNout(dc, 6); insert =(args...) -> (args...) -> 0))
                @test lazyouts(dc) ==                      vcat(1:4,           -1,-1, 5:8,           -1,-1, 9:12,          -1,-1)
                @test reshape(dc(fill(10f0, (1,1,3,1))), :) == [101,303,404,505,0, 0, 108,310,411,512,0 ,0, 115,317,418,519,0, 0]
            end

            @testset "DepthwiseConvAllowNinChangeStrategy" begin
                import NaiveNASflux: DepthwiseConvAllowNinChangeStrategy
                import NaiveNASlib: ΔNout
                inpt = inputvertex("in", 2, FluxConv{2}())
                dc = fluxvertex("dc", DepthwiseConv((1,1), nout(inpt) => 3*nout(inpt)), inpt)

                # Get output multiplier == 4 (nout = 4 * nin) by adding one more outgroup (4 = 3 + 1)
                okstrat = DepthwiseConvAllowNinChangeStrategy([1], [4], ΔNout(dc => 2))
                @test Δsize!(okstrat, dc)
                @test nout(dc) == 8
                @test nin(dc) == [2]

                failstrat = DepthwiseConvAllowNinChangeStrategy([10], [0], ΔNout(dc => 2))
                @test @test_logs (:warn, r"Could not change nout of dc") match_mode=:any Δsize!(failstrat, dc) == false
            end

            @testset "DepthwiseConvSimpleΔSizeStrategy" begin
                using NaiveNASflux: DepthwiseConvSimpleΔSizeStrategy
                using NaiveNASlib: ΔNout
                inpt = inputvertex("in", 2, FluxConv{2}())
                dc = fluxvertex("dc", DepthwiseConv((1,1), nout(inpt) => 3*nout(inpt)), inpt)

                okstrat = DepthwiseConvSimpleΔSizeStrategy(4, ΔNout(dc => 2))
                @test Δsize!(okstrat, dc)
                @test nout(dc) == 8
                @test nin(dc) == [2]

                # We tested complete failure above, so lets make the relaxation work here
                failstrat = DepthwiseConvSimpleΔSizeStrategy(5, ΔNout(dc => 3))
                @test_logs (:warn, r"Could not change nout of dc") Δsize!(failstrat, dc)
                @test nout(dc) == 10
                @test nin(dc) == [2]
            end
        end

        @testset "DepthwiseConv groupsize 2 into groupsize 1" begin
            inpt = inputvertex("in", 4, FluxConv{2}())
            dc1 = fluxvertex("dc1", DepthwiseConv((2,2), nout(inpt) => 2 * nout(inpt)), inpt)
            dc2 = fluxvertex("dc2", DepthwiseConv((2,2), nout(dc1) => nout(dc1)), dc1)

            @test @test_logs (:warn, r"Could not change nout of") Δnout!(v -> 1, dc1 => 2)
            @test [nout(dc1)] == nin(dc2) == [nout(dc2) ]== [12]
            @test lazyins(dc1) == [1:nout(inpt)]
            @test [lazyouts(dc1)] == lazyins(dc2) == [[1, 2, -1, 3, 4, -1, 5, 6, -1, 7, 8, -1]]

            # Add deterministic utilityfunction which wants to do non-contiguous selection across groups
            @test  @test_logs (:warn, r"Could not change nout of") Δnout!(dc1 => -2) do v
                util = repeat([1, 2], nout(v) ÷ 2)
                util[3] = -10 # don't want group 3 to make testcase deterministic
                return util
            end
            @test [nout(dc1)] == nin(dc2) == [nout(dc2)] == [8]

            @test lazyins(dc1) == [1:nout(inpt)]
            @test [lazyouts(dc1)] == lazyins(dc2) == [1:8]

            # Test that we actually succeeded in making a valid model
            y1 = dc1(ones(Float32, 3,3, nout(inpt), 2))
            @test size(y1, outdim(dc1)) == nout(dc1)
            y2 = dc2(y1)
            @test size(y2, outdim(dc2)) == nout(dc2)
        end

        @testset "DepthwiseConv groupsize 3 into groupsize 5" begin
            inpt = inputvertex("in", 4, FluxConv{2}())
            dc1 = fluxvertex("dc1", DepthwiseConv((2,2), nout(inpt) => 3 * nout(inpt)), inpt)
            dc2 = fluxvertex("dc2", DepthwiseConv((2,2), nout(dc1) => 5 * nout(dc1)), dc1)
            dc3 = fluxvertex("dc3", DepthwiseConv((2,2), nout(dc2) => nout(dc2)), dc2)

            # TODO: Check compgraph output pre and post?
            # Need to insert zeros then?

            @test @test_logs (:warn, r"Could not change nout of") Δnout!(v -> 1, dc1 => 2)
            @test [nout(dc1)] == nin(dc2) == [16]
            @test [nout(dc2)] == nin(dc3) == [96] # TODO: Why so big??

           # Add deterministic utilityfunction which wants to do non-contiguous selection across groups
            @test @test_logs (:warn, r"Could not change nout of") Δnout!(dc1 => -2) do v
                util = repeat([1, 2], nout(v) ÷ 2)
                util[1] = -10 # don't want group 1 to make testcase deterministic
                return util
            end
            @test [nout(dc1)] == nin(dc2) == [12]
            @test [nout(dc2)] == nin(dc3) == [96]

            @test lazyins(dc1) == [1:nout(inpt)]
            # NaiveNASlib might not pick non-new indices (i.e not -1) due to our artificial weight function above
            # default neuron utility function would give zero utility to new neurons
            @test [lazyouts(dc1)] == lazyins(dc2) == [[2, 3, -1, 5, 6, -1, 8, 9, -1, 11, 12, -1]]

            # All neurons had a positive utility, so NaiveNASlib should inrease to next valid size
            @test [lazyouts(dc2)] == lazyins(dc3)

            # Test that we actually succeeded in making a valid model
            y1 = dc1(ones(Float32,5,5, nout(inpt), 2))
            @test size(y1, outdim(dc1)) == nout(dc1)
            y2 = dc2(y1)
            @test size(y2, outdim(dc2)) == nout(dc2)
            y3 = dc3(y2)
            @test size(y3, outdim(dc3)) == nout(dc3)
        end

        @testset "Depthwise conv change input size from Conv" begin
            import NaiveNASflux: weights
            inpt = inputvertex("in", 4, FluxConv{2}())
            v1 = fluxvertex("v1", Conv((1,1), nout(inpt) => 3), inpt)
            v2 = fluxvertex("v2", DepthwiseConv((1,1), nout(v1) => 2 * nout(v1)), v1)
            v3 = fluxvertex("v3", Conv((1,1), nout(v2) => 3), v2)

            graph = CompGraph(inpt, v3)
            indata = randn(Float32, 1,1,4,1)
            expout = graph(indata)

            @test Δnout!(v -> 1, v1 => 1, v2 => 2)
            @test [nout(v1)] == nin(v2) == [4]
            @test [nout(v2)] == nin(v3) == [8]
            @test lazyouts(v2) == [1, 2, 3, 4, 5, 6, -1, -1]

            @test graph(indata) == expout

            @test Δnout!(v2 => -2) do v
                v == v2 || return 1
                return vcat(1:6, -10, -10) # Prefer to drop the last two inds
            end
            @test [nout(v1)] == nin(v2) == [3]
            @test [nout(v2)] == nin(v3) == [6]
            @test lazyouts(v2) == 1:6

            @test graph(indata) == expout

            @test size(v2(ones(Float32, 1,1,nout(v1),1))) == (1,1,nin(v3)[],1)
        end

        @testset "Depthwise conv change output size" begin
            import NaiveNASflux: weights
            inpt = inputvertex("in", 4, FluxConv{2}())
            v1 = fluxvertex("v1", Conv((1,1), nout(inpt) => 3), inpt)
            v2 = fluxvertex("v2", DepthwiseConv((1,1), nout(v1) => 2 * nout(v1)), v1)
            v3 = fluxvertex("v3", Conv((1,1), nout(v2) => 3), v2)

            graph = CompGraph(inpt, v3)
            indata = randn(Float32, 1,1,4,1)
            expout = graph(indata)

            @test Δnout!(v -> 1, v2 => 3 * nin(v1)[])
            @test [nout(v1)] == nin(v2) == [3]
            @test [nout(v2)] == nin(v3) == [18]
            @test lazyouts(v2) == [1, 2, -1, -1, -1, -1, 3, 4, -1, -1, -1, -1, 5, 6, -1, -1, -1, -1]

            @test graph(indata) == expout

            @test Δnout!(v2 => -6) do v
                v == v2 || return 1
                repeat(vcat(ones(4), -1, -1), 3)
            end
            @test [nout(v1)] == nin(v2) == [3]
            @test [nout(v2)] == nin(v3) == [12]
            @test lazyouts(v2) ==  [1, 2, 3, 4, 7, 8, 9, 10, 13, 14, 15, 16]

            @test graph(indata) == expout
        end
    end

    @testset "Concatenate activations" begin

        function testgraph(layerfun, nin1, nin2)
            vfun = (v, s) -> fluxvertex(layerfun(nout(v), s), v)
            return testgraph_vfun(vfun, nin1, nin2)
        end

        function testgraph_vfun(vertexfun, nin1::Integer, nin2::Integer)
            in1 = inputvertex("in1", nin1)
            in2 = inputvertex("in2", nin2)
            return testgraph_vfun(vertexfun, in1, in2)
        end

        function testgraph_vfun(vertexfun, in1, in2)
            l1 = vertexfun(in1, 3)
            l2 = vertexfun(in2, 7)

            joined = concat(l1, l2)
            l3 = vertexfun(joined, 9)

            @test nout(joined) == nout(l1) + nout(l2) == 10
            return CompGraph([in1, in2], l3)
        end

        @testset "Concatenate Dense" begin
            nin1 = 2
            nin2 = 5
            @test size(testgraph(Dense, nin1, nin2)(ones(nin1), ones(nin2))) == (9,)
        end

        @testset "Concatenate $rnntype" for rnntype in (RNN, GRU, LSTM)
            nin1 = 2
            nin2 = 5
            indata1 = reshape(collect(Float32, 1:nin1*4), nin1, 4)
            indata2 = reshape(collect(Float32, 1:nin2*4), nin2, 4)

            @test size(testgraph(rnntype, nin1, nin2)(indata1, indata2)) == (9,4)
        end

        @testset "Concatenate $convtype" for convtype in (Conv, ConvTranspose, CrossCor)
            nin1 = 2
            nin2 = 5
            indata1 = reshape(collect(Float32, 1:nin1*4*4), 4, 4, nin1, 1)
            indata2 = reshape(collect(Float32, 1:nin2*4*4), 4, 4, nin2, 1)

            convfun = (nin,nout) -> convtype((3,3), nin=>nout, pad = (1,1))
            @test size(testgraph(convfun, nin1, nin2)(indata1, indata2)) == (4,4,9,1)
        end

        @testset "Concatenate Pooled Conv" begin
            nin1 = 2
            nin2 = 5
            indata1 = reshape(collect(Float32, 1:nin1*4*4), 4, 4, nin1, 1)
            indata2 = reshape(collect(Float32, 1:nin2*4*4), 4, 4, nin2, 1)
            function vfun(v, s)
                cv = fluxvertex(Conv((3,3), nout(v)=>s, pad = (1,1)), v)
                return fluxvertex(MaxPool((3,3), pad=(1,1), stride=(1,1)), cv)
            end
            @test size(testgraph_vfun(vfun, nin1, nin2)(indata1, indata2)) == (4,4,9,1)
        end

        @testset "Concatenate InputShapeVertex" begin
            nin1 = 6
            nin2 = 4
            in1 = inputvertex("in1", nin1, FluxDense())
            in2 = inputvertex("in2", nin2, FluxDense())
            @test size(testgraph_vfun((v,s) -> v, in1, in2)(ones(nin1), ones(nin2))) == (10,)
        end

        @testset "Concatenate BatchNorm only" begin
            nin1 = 3
            nin2 = 7
            indata1 = reshape(collect(Float32, 1:nin1*4*4), 4, 4, nin1, 1)
            indata2 = reshape(collect(Float32, 1:nin2*4*4), 4, 4, nin2, 1)
            in1 = conv2dinputvertex("in1", nin1)
            in2 = conv2dinputvertex("in2", nin2)
            vfun(v,s) = fluxvertex(BatchNorm(nout(v)), v)

            @test size(testgraph_vfun(vfun, in1, in2)(indata1, indata2)) == (4,4,10,1)
        end

        @testset "Concatenate elementwise and concatenated" begin
            import NaiveNASflux: FluxDense
            in1 = denseinputvertex("in1", 3)
            in2 = denseinputvertex("in2", 3)

            v1 = "v1" >> in1 + in2
            v2 = concat("v2", in1, in2)
            v3 = concat("v3", v1, v1)
            v4 = "v4" >> v2 + v3
            v5 = "v5" >> v2 + v3
            v6 = concat("v6", v4, v5, in1)

            @test nin(v6) == [nout(v4), nout(v5), nout(in1)] == [6, 6, 3]
            @test nout(v6) == 15
        end

        @testset "Concatentate dimension mismatch fail" begin
            d1 = fluxvertex(Dense(2,3), inputvertex("in1", 2))
            c1 = fluxvertex(Conv((3,3), 4=>5), inputvertex("in2", 4))
            r1 = fluxvertex(RNN(6,7), inputvertex("in3", 6))
            @test_throws DimensionMismatch concat(d1, c1)
            @test_throws DimensionMismatch concat(r1, c1)
            @test_throws DimensionMismatch concat(d1, r1)
        end

        @testset "Concat with name" begin
            d1 = fluxvertex(Dense(2,3), inputvertex("in1", 2))
            d2 = fluxvertex(Dense(2,5), inputvertex("in1", 2))
            c = concat("c", d1, d2)

            @test name(c) == "c"
            @test nout(c) == 8
            @test nin(c) == [3, 5]
        end
    end

    @testset "Tricky structures" begin
        import NaiveNASflux: NoParams
        mutable struct Probe
            activation
        end
        function (p::Probe)(x)
            p.activation = x
            return x
        end

        function probe(in)
            p = Probe(nothing)
            return invariantvertex(NoParams(p), in), p
        end

        conv3x3(inpt::AbstractVertex, nch::Integer) = fluxvertex(Conv((3,3), nout(inpt)=>nch, pad=(1,1)), inpt)
        batchnorm(inpt) = fluxvertex(BatchNorm(nout(inpt)), inpt)
        mmaxpool(inpt) = fluxvertex(MaxPool((2,2)), inpt)

        @testset "Residual Conv block" begin
            inpt = inputvertex("in", 3)
            conv1 = conv3x3(inpt, 5)
            pv1, p1 = probe(conv1)
            bn1 = batchnorm(pv1)
            conv2 = conv3x3(bn1, 5)
            pv2, p2 = probe(conv2)
            bn2 = batchnorm(pv2)
            add = bn2 + bn1
            mp = mmaxpool(add)
            out = conv3x3(mp, 3)

            graph = CompGraph([inpt], [out])

            # Test that the graph works to begin with
            indata = reshape(collect(Float32, 1:3*4*4), 4, 4, 3, 1)
            @test size(graph(indata)) == (2, 2, 3, 1)
            @test size(p1.activation) == (4, 4, 5, 1)
            @test size(p2.activation) == (4, 4, 5, 1)

            Δnin!(v -> 1:nout(v), out => -1)

            @test size(graph(indata)) == (2, 2, 3, 1)
            @test size(p1.activation) == (4, 4, 4, 1)
            @test size(p2.activation) == (4, 4, 4, 1)
        end

        @testset "Residual fork Conv block" begin
            inpt = inputvertex("in", 3)
            conv1 = conv3x3(inpt, 5)
            pv1, p1 = probe(conv1)
            bn1 = batchnorm(pv1)
            conv1a = conv3x3(bn1, 3)
            pv1a, p1a = probe(conv1a)
            bn1a = batchnorm(pv1a)
            conv1b = conv3x3(bn1, 2)
            pv1b, p1b = probe(conv1b)
            bn1b = batchnorm(pv1b)
            mv = concat(bn1a, bn1b)
            add = mv + bn1
            mp = mmaxpool(add)
            out = conv3x3(mp, 3)

            graph = CompGraph([inpt], [out])

            # Test that the graph works to begin with
            indata = reshape(collect(Float32, 1:3*4*4), 4, 4, 3, 1)
            @test size(graph(indata)) == (2, 2, 3, 1)
            @test size(p1.activation) == (4, 4, 5, 1)
            @test size(p1a.activation) == (4, 4, 3, 1)
            @test size(p1b.activation) == (4, 4, 2, 1)

            Δnin!(v -> 1:nout(v), out => -1)

            @test size(graph(indata)) == (2, 2, 3, 1)
            @test size(p1.activation) == (4, 4, 4, 1)
            @test size(p1a.activation) == (4, 4, 2, 1)
            @test size(p1b.activation) == (4, 4, 2, 1)
        end

        rnnvertex(inpt, outsize) = fluxvertex("rnn", RNN(nout(inpt), outsize), inpt)
        densevertex(inpt, outsize) = fluxvertex("dense", Dense(nout(inpt), outsize), inpt)

        @testset "RNN to Dense" begin
            inpt = inputvertex("in", 4)
            rnn = rnnvertex(inpt, 5)
            pv, p = probe(rnn)
            dnn = densevertex(pv, 3)

            graph = CompGraph([inpt], [dnn])

            indata = [collect(Float32, 1:nin(rnn)[]) for i =1:10]
            @test size(hcat(graph.(indata)...)) == (3,10)
            @test size(p.activation) == (5,)

            Δnin!(v -> 1:nout(v), dnn, 1)

            @test size(hcat(graph.(indata)...)) == (3,10)
            @test size(p.activation) == (6,)

            Δnout!(v -> 1:nout(v), rnn, -2)

            @test size(hcat(graph.(indata)...)) == (3,10)
            @test size(p.activation) == (4,)
        end
    end
end

@testset "functor" begin
    import Functors: functor, fmap
    import Flux: params
    import NaiveNASflux: weights, bias, FluxDense
    inpt = inputvertex("in", 2, FluxDense())
    v1 = fluxvertex(Dense(2, 3), inpt)
    v2 = fluxvertex(Dense(3, 4), v1)
    v3 = concat(v2, v1)
    v4 = fluxvertex(Dense(nout(v3), 2), v3)
    g1 = CompGraph(inpt, v4)

    pars1 = params(g1).order
    @test pars1[1] == weights(layer(v1))
    @test pars1[2] == bias(layer(v1))
    @test pars1[3] == weights(layer(v2))
    @test pars1[4] == bias(layer(v2))
    @test pars1[5] == weights(layer(v4))
    @test pars1[6] == bias(layer(v4))

    # Basically what Flux.gpu does except function is CuArrays.cu(x) instead of 2 .* x
    testfun(x) = x
    testfun(x::AbstractArray) = 2 .* x
    g2 = fmap(testfun, g1)

    pars2 = params(g2).order.data
    @test pars2 == 2 .* pars1

    indata = randn(nout(inpt), 1)

    g3 = Flux.f64(g1)
    @test g3(indata) ≈ g1(indata) rtol=1e-7
end

@testset "Trainable insert values" begin
    using Random
    import NaiveNASflux: weights, bias
    import Flux: params, Descent, train!, mse

    @testset "Dense-Dense-Dense" begin
        Random.seed!(0)
        iv = denseinputvertex("in", 3)
        v1 = fluxvertex("v1", Dense(3,3), iv)
        v2 = fluxvertex("v2", Dense(3,4), v1)
        v3 = fluxvertex("v3", Dense(4,2), v2)

        g = CompGraph(iv, v3)

        indata = randn(3,4)
        expectedout = g(indata)

        Δnout!(v -> 1, v1 => 2)
        NaiveNASflux.forcemutation(g)

        @test g(indata) ≈ expectedout

        train!((x,y) -> mse(g(x), y), params(g), [(randn(nin(v1)[],8), randn(nout(v3) ,8))], Descent(0.5))

        @test minimum(abs.(weights(layer(v1)))) > 0
        @test minimum(abs.(weights(layer(v2)))) > 0
        @test minimum(abs.(weights(layer(v3)))) > 0
    end

    @testset "Conv-Bn-Conv" begin
        import NaiveNASflux: FluxConv
        Random.seed!(0)
        iv = inputvertex("in", 2, FluxConv{2}())
        v1 = fluxvertex("v1", Conv((1,1), 2 => 2), iv)
        v2 = fluxvertex("v2", BatchNorm(2), v1)
        v3 = fluxvertex("v3", Conv((1,1), 2 => 2), v2)

        g = CompGraph(iv, v3)

        indata = randn(Float32, 2,2,2,8)
        expectedout = g(indata)

        Δnout!(v -> 1, v1 => 2)
        NaiveNASflux.forcemutation(g)

        @test g(indata) == expectedout

        train!((x,y) -> mse(g(x), y), params(g), [(randn(Float32,2,2,2,8), randn(Float32,2,2,2,8))], Descent(0.5))

        @test minimum(abs.(weights(layer(v1)))) > 0
        @test minimum(abs.(weights(layer(v3)))) > 0
    end

    @testset "Conv-Conv-Conv" begin
        import NaiveNASflux: FluxConv
        Random.seed!(0)
        iv = inputvertex("in", 2, FluxConv{2}())
        v1 = fluxvertex("v1", Conv((1,1), 2 => 2), iv; layerfun=ActivationContribution ∘ LazyMutable)
        v2 = fluxvertex("v2", Conv((1,1), 2 => 2), v1; layerfun=ActivationContribution ∘ LazyMutable)
        v3 = fluxvertex("v3", Conv((1,1), 2 => 2), v2; layerfun=ActivationContribution ∘ LazyMutable)

        g = CompGraph(iv, v3)

        indata = randn(Float32, 2,2,2,8)
        expectedout = g(indata)

        Δnout!(v->1, v1 => 2, v2 => 1)
        NaiveNASflux.forcemutation(g)

        @test g(indata) == expectedout

        train!((x,y) -> mse(g(x), y), params(g), [(randn(Float32,2,2,2,8), randn(Float32,2,2,2,8))], Descent(0.5))

        @test minimum(abs.(weights(layer(v1)))) > 0
        @test minimum(abs.(weights(layer(v2)))) > 0
        @test minimum(abs.(weights(layer(v3)))) > 0
    end
end

@testset "setlayer" begin
    v = fluxvertex(Dense(3, 4, Flux.relu), inputvertex("in", 3))
    NaiveNASflux.setlayer!(v, (;σ=tanh))
    @test layer(v).σ == tanh
end


text:
using MPI
using p4est_wrapper
using Test

# Initialize MPI if not initialized yet
if !MPI.Initialized()
    MPI.Init()
end

#############################################################################
# Definition of constants, data structures, functions and function callbacks
#############################################################################

## Global variable which is updated across calls to init_fn_callback_2d
current_quadrant_index = Cint(0)
## Global variable which is updated across calls to refine_replace_callback_2d
num_calls = Cint(0)

## Better to use a C-enum. But I did not use it in order to keep the Julia
## version of this C example as simple as possible
const nothing_flag = Cint(0)
const refine_flag  = Cint(1)
const coarsen_flag = Cint(2)

# This C callback function is called once per quadtree quadrant. Here we are assuming
# that p4est->user_pointer has been set prior to the first call to this call
# back function to an array of ints with as many entries as forest quadrants. This call back function
# initializes the quadrant->p.user_data void * pointer of all quadrants such that it
# points to the corresponding entry in the global array mentioned in the previous sentence.
function init_fn_callback_2d(forest_ptr::Ptr{p4est_t},
                            which_tree::p4est_topidx_t,
                            quadrant_ptr::Ptr{p4est_quadrant_t})
    @assert which_tree == 0
    # Extract a reference to the first (and uniquely allowed) tree
    forest = forest_ptr[]
    tree = p4est_tree_array_index(forest.trees, 0)[]
    quadrant = quadrant_ptr[]
    q = p4est_quadrant_array_index(tree.quadrants, current_quadrant_index)
    @assert p4est_quadrant_compare(q,quadrant_ptr) == 0
    user_data  = unsafe_wrap(Array, Ptr{Cint}(forest.user_pointer), current_quadrant_index+1)[current_quadrant_index+1]
    unsafe_store!(Ptr{Cint}(quadrant.p.user_data), user_data, 1)
    global current_quadrant_index = (current_quadrant_index+1) % (tree.quadrants.elem_count)
    return nothing
end

const init_fn_callback_2d_c = @cfunction(init_fn_callback_2d, Cvoid, (Ptr{p4est_t}, p4est_topidx_t, Ptr{p4est_quadrant_t}))

function refine_callback_2d(::Ptr{p4est_t},
                            which_tree::p4est_topidx_t,
                            quadrant_ptr::Ptr{p4est_quadrant_t})
    @assert which_tree == 0
    quadrant = quadrant_ptr[]
    return Cint(unsafe_wrap(Array, Ptr{Cint}(quadrant.p.user_data),1)[] == refine_flag)
end

const refine_callback_2d_c = @cfunction(refine_callback_2d, Cint, (Ptr{p4est_t}, p4est_topidx_t, Ptr{p4est_quadrant_t}))

# For each refined quadrant:
#   num_calls is even -> leave untouched all children
#   num_calls is odd  -> mark all children for coarsening
function refine_replace_callback_2d(forest_ptr::Ptr{p4est_t},
                                    which_tree::p4est_topidx_t,
                                    num_outgoing::Cint,
                                    outgoing_ptr::Ptr{Ptr{p4est_quadrant_t}},
                                    num_incoming::Cint,
                                    incoming_ptr::Ptr{Ptr{p4est_quadrant_t}})
    @assert which_tree   == 0
    @assert num_outgoing == 1
    @assert num_incoming == p4est_wrapper.P4EST_CHILDREN
    incoming = unsafe_wrap(Array, incoming_ptr, p4est_wrapper.P4EST_CHILDREN)
    for quadrant_index = 1:p4est_wrapper.P4EST_CHILDREN
        quadrant = incoming[quadrant_index][]
        if (num_calls % 2 == 0)
            unsafe_store!(Ptr{Cint}(quadrant.p.user_data), nothing_flag, 1)
        else
            unsafe_store!(Ptr{Cint}(quadrant.p.user_data), coarsen_flag, 1)
        end
    end
    global num_calls += 1
    return nothing
end

const refine_replace_callback_2d_c = @cfunction(refine_replace_callback_2d, Cvoid, (Ptr{p4est_t}, p4est_topidx_t, Cint, Ptr{Ptr{p4est_quadrant_t}}, Cint, Ptr{Ptr{p4est_quadrant_t}}))

function coarsen_callback_2d(forest_ptr::Ptr{p4est_t},
                            which_tree::p4est_topidx_t,
                            quadrants_ptr::Ptr{Ptr{p4est_quadrant_t}})
    @assert which_tree == 0
    coarsen = Cint(1)
    quadrants = unsafe_wrap(Array, quadrants_ptr, p4est_wrapper.P4EST_CHILDREN)
    for quadrant_index = 1:p4est_wrapper.P4EST_CHILDREN
        quadrant = quadrants[quadrant_index][]
        coarsen = Cint(unsafe_wrap(Array, Ptr{Cint}(quadrant.p.user_data), 1)[] ==  coarsen_flag)
        if (!Bool(coarsen)) return coarsen end
    end
    return coarsen
end

const coarsen_callback_2d_c = @cfunction(coarsen_callback_2d, Cint, (Ptr{p4est_t}, p4est_topidx_t, Ptr{Ptr{p4est_quadrant_t}}))

## Refine those cells with even identifier    (0,2,4,6,8,...)
## Leave untouched cells with odd identifier  (1,3,5,7,9,...)
function allocate_and_set_refinement_and_coarsening_flags(forest_ptr::Ptr{p4est_t})
    forest = forest_ptr[]
    tree = p4est_tree_array_index(forest.trees, 0)[]
    return [i%2 == 0 ? nothing_flag : refine_flag for i = 1:tree.quadrants.elem_count]
end

function perform_single_mesh_adaptation_step(forest_ptr::Ptr{p4est_t}, step::Int)
    user_data = allocate_and_set_refinement_and_coarsening_flags(forest_ptr)
    p4est_reset_data(forest_ptr, Cint(sizeof(Cint)), init_fn_callback_2d_c, pointer(user_data))
    p4est_refine_ext(forest_ptr, 0, -1, refine_callback_2d_c, C_NULL, refine_replace_callback_2d_c)
    p4est_coarsen(forest_ptr, 0, coarsen_callback_2d_c, C_NULL)
    p4est_partition(forest_ptr, 1, C_NULL)
    p4est_vtk_write_file(forest_ptr, C_NULL, string("step_", step))
end

#############################################################################
# Main program
#############################################################################

mpicomm = p4est_wrapper.P4EST_ENABLE_MPI ? MPI.COMM_WORLD : Cint(0)
sc_init(mpicomm, Cint(true), Cint(true), C_NULL, p4est_wrapper.SC_LP_DEFAULT)
p4est_init(C_NULL, p4est_wrapper.SC_LP_DEFAULT)

unitsquare_connectivity = p4est_connectivity_new_unitsquare()
unitsquare_forest = p4est_new(mpicomm, unitsquare_connectivity, 0, C_NULL, C_NULL)

for i=1:7
    perform_single_mesh_adaptation_step(unitsquare_forest,i)
end

p4est_destroy(unitsquare_forest);
p4est_connectivity_destroy(unitsquare_connectivity)
sc_finalize()

# Finalize MPI if initialized and session is not interactive
if (MPI.Initialized() && !isinteractive())
    MPI.Finalize()
end


text:
function empty_regression_deviation(F1 = Float64)
    return RegressionDeviation{F1}(name = :empty)
end

"""
	$(SIGNATURES)

Set regressors to exclude from scalar deviation.
"""
exclude_regressors!(d :: RegressionDeviation{F1}, exclV) where F1 =
    d.excludeV = exclV;

is_excluded(d :: RegressionDeviation{F1}, name :: Symbol) where F1 =
    name ∈ d.excludeV;

function drop_excluded(d :: RegressionDeviation{F1}, nameV) where F1
    if isempty(d.excludeV)
        outV = copy(nameV);
    else
        outV = filter(n -> !is_excluded(d, n), nameV);
    end
    return outV
end

"""
	$(SIGNATURES)

Returns vector of names, coefficients and std errors. Always sorted by name. Otherwise, the order is unpredictable.
"""
function get_unpacked_data_values(d :: RegressionDeviation{F1};
    dropExcluded :: Bool = false) where F1
    dataRt = get_data_values(d);
    nameV = sort(EconometricsLH.get_names(d.dataV));
    if dropExcluded
        nameV = drop_excluded(d, nameV);
    end
    coeffV, seV = get_coeff_se_multiple(d.dataV, nameV);
    return nameV, coeffV, seV
end

# Return coefficients and std errors in same order as for data
function get_unpacked_model_values(d :: RegressionDeviation{F1};
    dropExcluded :: Bool = false) where F1

    nameV = sort(EconometricsLH.get_names(d.modelV));
    if dropExcluded
        nameV = drop_excluded(d, nameV);
    end
    coeffV, seV = get_coeff_se_multiple(d.modelV, nameV);
    return nameV, coeffV, seV
end

"""
	$(SIGNATURES)

Set model values for a RegressionDeviation.

If model regression is missing regressors, the option of setting these regressors to 0 is provided (e.g. for dummies without occurrences).
"""
function set_model_values(d :: RegressionDeviation, modelV :: RegressionTable);
    # setMissingRegressors :: Bool = false)

    # if setMissingRegressors
    #     set_missing_regressors!(modelV, get_names(d.dataV));
    # end

    if !have_same_regressors([d.dataV, modelV])
        rModelV = EconometricsLH.get_names(modelV);
        rDataV = EconometricsLH.get_names(d.dataV);
        error("""Regressors do not match: $d
            Model: $rModelV
            Data:  $rDataV""")
    end
    d.modelV = deepcopy(modelV);
    return nothing
end

function get_weights(d :: RegressionDeviation)
    return 1.0
end

# For RegressionDeviation: sum of (model - data) / se
# se2coeffLb governs the scaling of the std errors. s.e./beta >= se2coeffLb
# Excludes regressors if this is set in the `RegressionDeviation`.
function scalar_dev(d :: RegressionDeviation; se2coeffLb :: Float64 = 0.1,
    inclScalarWt :: Bool = true)
    nameV, coeffV, seV = get_unpacked_data_values(d; dropExcluded = true);
    mNameV, mCoeffV, _ = get_unpacked_model_values(d; dropExcluded = true);
    # Names and sort order must match
    @assert isequal(nameV, mNameV)  "Names do not match: \n $nameV \n $mNameV"

    seV = max.(seV, se2coeffLb .* abs.(coeffV));
    scalarDev = scalar_deviation(mCoeffV, coeffV, 1.0 ./ seV; p = norm_p(d));
    if inclScalarWt
        scalarDev *= d.scalarWt;
    end
    scalarStr = sprintf1(d.fmtStr, scalarDev);
    return scalarDev, scalarStr
end

"""
	$(SIGNATURES)

Show a RegressionDeviation
"""
function regression_show_fct(d :: RegressionDeviation;
    showModel :: Bool = true, fPath :: String = "")

    nameV, coeffV, seV = get_unpacked_data_values(d);
    dataM = hcat(nameV, round.(coeffV, digits = 3), round.(seV, digits = 3));
    headerV = ["Regressor", "Data", "s.e."];

    if showModel
        _, mCoeffV, mSeV = get_unpacked_model_values(d);
        dataM = hcat(dataM,  round.(mCoeffV, digits = 3), round.(mSeV, digits = 3));
        headerV = vcat(headerV, ["Model", "s.e."])
    end

    io = open_show_path(d, fPath = fPath);
    pretty_table(io, dataM,  headerV);
    close_show_path(d, io);
end

Base.show(io :: IO, d :: RegressionDeviation{F1}) where F1 =
    Base.print(io, "$(name(d)):  ", short_description(d));

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


text:
"""
    clone(package_or_url::String, to::String=pwd(); force::Bool=false)

Clone a package or remote URL to a local directory.

# Intro

This command helps you clone a remote package repo from github.
The corresponding package or repo name will automatically remove
the `.jl` extension if it has one, e.g `https://github.com/JuliaLang/Example.jl`
will be downloaded as `Example`. And when cloning a repo that has
no write access, it will prompt to ask for a fork.

# Arguments

- `package_or_url`: package name or URL.
- `to`: local directory to download the git repo to.

# Flags

- `-f,--force`: force overwrite existing path.
- `-y,--yes`: skip prompt with answer yes.
"""
@cast function clone(package_or_url::String, to::String=pwd(); force::Bool=false, yes::Bool=false)
    ispath(to) || mkpath(to)
    try
        if isurl(package_or_url)
            clone_url(package_or_url, to, force, yes)
        else
            clone_package(package_or_url, to, force, yes)
        end
    catch e
        cmd_error("fail to clone $package_or_url")
    end
    return
end

function clone_url(url::String, to::String, force::Bool, yes::Bool)
    if endswith(url, "jl.git")
        _clone(url, joinpath(to, basename(url)[1:end-7]), force, yes)
    else
        _clone(url, joinpath(to, basename(url)), force, yes)
    end
    return
end

function clone_package(package::String, to::String, force::Bool, yes::Bool)
    info = find_package(package)
    isnothing(info) && cmd_error("cannot find $package in local registries")
    pkg_toml = get_registry_file(info.reg, joinpath(info.path, "Package.toml"))
    _clone(pkg_toml["repo"], joinpath(to, pkg_toml["name"]), force, yes)
    return
end

function _clone(url::String, to::String, force::Bool, yes::Bool)
    force && ispath(to) && rm(to; force=true, recursive=true)
    username = readchomp(`git config user.name`)
    auth = GitHub.authenticate(read_github_auth())
    rp = fetch_repo_from_url(url; auth=auth)

    local has_access
    try
        has_access = iscollaborator(rp, username; auth=auth)
    catch e
        has_access = false
    end

    if has_access
        git_clone(url, to)
    elseif prompt("do not have access to $url, fork?"; yes)
        @info "fork upstream repo: $(rp.full_name)"
        owned_repo = create_fork(rp; auth)
        git_clone(owned_repo.clone_url.uri, to)
        @info "setting upstream to $url"
        redirect_stdio(;stdout=devnull, stderr=devnull) do
            git_set_upstream(to; url)
        end
    else # no access only clone
        git_clone(url, to)
    end
    return
end


text:
#==============================================================================#
# EC2/test/runtests.jl
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

using AWSEC2
using AWSCore
using Test

AWSCore.set_debug_level(1)

#-------------------------------------------------------------------------------
# Load credentials...
#-------------------------------------------------------------------------------

aws = AWSCore.aws_config()

#-------------------------------------------------------------------------------
# EC2 tests
#-------------------------------------------------------------------------------

@test ec2_id(aws, "Not a real server name!!") == nothing

r = ec2(aws, Dict("Action" => "DescribeImages",
                  "Filter.1.Name" => "owner-alias",
                  "Filter.1.Value" => "amazon",
                  "Filter.2.Name" => "name",
                  "Filter.2.Value" => "amzn-ami-hvm-2015.09.1.x86_64-gp2"))

@test r["imagesSet"]["item"]["description"] ==
      "Amazon Linux AMI 2015.09.1 x86_64 HVM GP2"

#==============================================================================#
# End of file.
#==============================================================================#


text:
__precompile__(true)

module Grabbit

using Compat
using Compat: @debug
using Compat.Printf
using JSON
using DataStructures
using ArgCheck
using Missings
using EnglishText

# grab files with structure directory/filenames
# two steps:
# 1. parse config into entities, and then read in layout from filesystem based
#    on config
# 2. search layout given query ("get(; kw...)")

# grabbit manages this with the Layouts, Domains, Entities, and Files.  Files
# are tagged with Entities + Domains, Domains have a root path and list of
# entities (and files?).  domains apply in order of proximity (so when entities
# conflict, closer domains can "mask" top-level ones).
#
# domains are puzzling me...there's somethign weird about how the root is
# computed: https://github.com/grabbles/grabbit/issues/57.  Okay that's a bug.
# They should apply to the subdirectory their config is found in.

# The logic is:
#
# walk directory tree, from top down.
#
# if you see a config file, create a new domain and compose it with the current one.
#
# for each file: extract matching entities, generate tags, and add tags to file
#
#
# Querying:
#
# get files matching tags.
# get files matching tags, restricted by domain
# get files matching tags, restricted by extension
# get entity ids/dirs matching tags.
#
# design considerations:
#
# what's the most natural data structure for representing all this information?
# one-direction nesting (Domain [-> Domain] -> Entity -> File), or a
# many-to-many mapping (Domain -> { Entities, Files }, Entity -> { Domain,
# Files}, File -> { Domains, Entities }).  The first feels more...elegant
# somehow but the second is how the python implenentation works and is a little
# easier to reason about for querying: you just search the vector of things you
# have to find the matching ones.  But it feels like you're missing out on some
# of the benefits of having the _structure_: you don't have to search below
# something where you know there's a mismatch...but then again, the way the
# matching happens in the python one there could be a mismatch (if there's two
# `sub-NN` where there are different `NN`s...)

# python implementation does it this way:
#
# keep a list of domains. for every file, check whether domain applies by
# matching root of domain against file path.  extract entities from applicable
# domains.  match entitites against file path.  construct File object with
# applicable domains and matching entities, and add the file to the relevant
# domains/entitites
#
# layout: entities are stored with _ID_ key ("domain.name").  when getting
# unique values for entity, or count of files, searches by name and by id.
#
# file: entities are stored with name key
#
# query by entity name, optionally with domain qualification (as domains kw)

mutable struct Entity
    name::String
    pattern::Regex
    mandatory::Bool
    values::OrderedSet
end

mutable struct Domain
    name::String
    config::Dict
    root::String
    parent::Union{Domain,Compat.Nothing}
    entities::Dict{String,Entity}
    include                     # predicte function for whether to include
end

Base.show(io::IO, d::Domain) = print(io, "Domain $(d.name) ($(d.root))")

function Domain(config::Dict, parent=nothing)
    @argcheck(!all(haskey.(Ref(config), ["include", "exclude"])),
              "Cannot specify both include and exclude regex")

    d = Domain(config["name"],
               config,
               config["root"],
               parent,
               Dict{String,Entity}(),
               make_include_predicate(config))

    for e in config["entities"]
        e = Entity(e)
        d.entities[e.name] = e
    end

    return d
end

make_include_predicate(config) =
    if haskey(config, "include")
        f -> Compat.occursin(Regex(make_pattern(config["include"])), f)
    elseif haskey(config, "exclude")
        f -> Compat.occursin(Regex(make_pattern(config["exclude"])), f) ? false : missing
    else
        f -> missing
    end

_include(d::Domain, f::AbstractString) = d.include(f) & _include(d.parent, f)
_include(::Compat.Nothing, f::AbstractString) = missing
include(d::Domain, f::AbstractString) = (incl = _include(d, f); incl === missing ? true : incl)
# curry
include(d::Domain) = f -> include(d, f)

function Entity(config::Dict)
    name = config["name"]
    pattern = Regex(config["pattern"])
    mandatory = convert(Bool, get(config, "mandatory", false))
    e = Entity(name, pattern, mandatory, OrderedSet())
    return e
end

Base.show(io::IO, e::Entity) = println(io, "Entity $(e.name)")

function match!(entity::Entity, fn::AbstractString)
    m = match(entity.pattern, fn)
    if m !== nothing
        value = m.captures[1]
        push!(entity.values, value)
        return value
    elseif entity.mandatory
        error("Mandatory entity $(entity.name) failed to match file $(f.path).")
    else
        return nothing
    end
end

mutable struct File
    path::String
    filename::String
    dirname::String
    tags::Dict
end

Base.show(io::IO, f::File) = show(io, f.path)

function File(fn::AbstractString, domain::Domain, entities::Dict{String,Entity})
    f = File(fn, basename(fn), dirname(fn), Dict())
    @debug "  Parsing file $fn"

    for entity in values(entities)
        m = match!(entity, fn)
        if m !== nothing
            f.tags[entity.name] = m
            @debug "    ✔ $(entity.name): $m"
        else
            @debug "    ✘ $(entity.name)"
        end
    end
    return f
end

tags(f::File) = f.tags
path(f::File) = f.path
Base.basename(f::File) = f.filename
Base.dirname(f::File) = f.dirname

mutable struct Layout
    root::String
    entities::OrderedDict{String,Entity}
    mandatory::Set{Entity}
    domains::OrderedDict{String,Domain}
    files::Vector{File}
    config_filenames::Set{String}
end

function Base.show(io::IO, layout::Layout)
    files = ItemQuantity( length(layout.files), "file")
    print(io, "Layout of $files in $(layout.root)")
    for (domain_name, domain) in layout.domains
        print(io, "\n  Domain \"$domain_name\":")
        maxlen = maximum(length(k) for k in keys(domain.entities))
        for (k, e) in domain.entities
            print(io, "\n    ")
            print(io, rpad("\"$k\"", maxlen+3))
            chars_left = displaysize(io)[2] - (maxlen+3) - 4
            reprs = []
            for v in e.values
                r = repr(v)
                if chars_left < length(r)+5
                    push!(reprs, "…")
                    break
                else
                    push!(reprs, repr(v))
                    chars_left -= length(r)+2
                end
            end
            print(io, "[")
            join(io, reprs, ", ")
            print(io, "]")
        end
    end
end

function Base.showcompact(io::IO, layout::Layout)
    files = ItemQuantity( length(layout.files), "file")
    domains = ItemQuantity(length(layout.domains), "domain")
    entities = ItemQuantity(length(layout.entities), "entity")
    print(io, "Layout of $files in $(layout.root) ($entities in $domains)")
end

# placeholder
is_config(l::Layout, filename::AbstractString) = false
is_config(l::Layout) = f -> is_config(l, f)

"""
    parse_config(root, config)

Read in config JSON file (relative to root), adding a "root" key if it's missing
"""
parse_config(root::AbstractString, config::AbstractString) =
    merge(Dict("root"=>root), JSON.parsefile(joinpath(root, config)))

Layout(root::AbstractString, config::AbstractString) = Layout(parse_config(root, config))

const DEFAULT_ENTITIES =
    OrderedDict("extension" => Entity("extension",
                               r"\.([^/]*?)$",
                               false,
                               OrderedSet()))

function Layout(config::Dict)
    root = config["root"]
    domain = Domain(config)
    entities = merge(Dict(DEFAULT_ENTITIES), domain.entities)

    # filenames to look for inside tree to create new domains
    config_files = get(config, "config_filename", [])

    l = Layout(root,
               merge(DEFAULT_ENTITIES,
                     OrderedDict{String,Entity}("$(domain.name).$k"=>v for (k,v) in entities)),
               Set{Entity}(e for (k, e) in entities if e.mandatory),
               OrderedDict(domain.name=>domain),
               File[],
               config_files)

    # walk dir top down
    parsedir!(l, root, domain, entities)

    return l
end

Layout(r, e, m, d, f, config_files::AbstractString) =
    Layout(r, e, m, d, f, Set{String}((config_files, )))
Layout(r, e, m, d, f, config_files::AbstractArray) =
    Layout(r, e, m, d, f, Set{String}(string(f) for f in config_files))

function parsedir!(layout::Layout, current, domain::Domain, entities::Dict{String,Entity})
    @debug "Parsing directory $current"
    contents = [joinpath(current, x) for x in readdir(current) if include(domain, x)]
    dirs = filter(isdir, contents)
    files = filter(!isdir, contents)

    configs = filter(is_config(layout), files)
    filter!(!is_config(layout), files)

    for config in filter(is_config(layout), files)
        error("Additional config files not supported yet.")
    end

    for file in files
        push!(layout.files, File(file, domain, entities))
    end

    for dir in dirs
        parsedir!(layout, dir, domain, entities)
    end

    return layout
end

################################################################################
# Querying a layout
################################################################################

make_pattern(val::Number) = "0*$val"
make_pattern(val::AbstractArray) = "(" * join(make_pattern.(val), "|") * ")"
make_pattern(val) = "$val"

function make_regex(val, regex_search=false)
    pattern = make_pattern(val)
    regex_search ? Regex(pattern) : Regex("^$pattern\$")
end

function make_query(layout, filters)
    filters = Dict(string(k)=>make_regex(v) for (k,v) in filters)
    function(f::File)
        @debug "Querying file $(f.filename):"
        for (name, filter) in filters
            if !haskey(f.tags, name)
                @debug "  ✘ no $name"
                return false
            elseif !Compat.occursin(filter, f.tags[name])
                @debug "  ✘ $name ($(f.tags[name])) ≠ $(filter.pattern)"
                return false
            else
                @debug "  ✔ $name = $(filter.pattern)"
            end
        end
        return true
    end
end

function Base.get(layout::Layout, queries::Dict=Dict(); kw...)
    queries = merge!(Dict{Any,Any}(kw), queries)
    query = make_query(layout, queries)
    filter(query, layout.files)
end

Base.get(layout::Layout, r::Regex) = [f for f in layout.files if ismatch(r, path(f))]

function tags(layout)
    # get all the tags, in teh form of a Dict(entity=>values)
    ts = Dict()
    for (k, e) in layout.entities
        vals = get!(ts, e.name, OrderedSet())
        union!(vals, e.values)
    end
    return ts
end

function Base.get(layout::Layout, target::AbstractString)
    reduce(union!, OrderedSet(), e.values for (k,e) in layout.entities if e.name==target)
end

Base.get(file::File, target::AbstractString) = file.tags[target]

export
    Domain,
    Entity,
    Layout,
    File,
    basename,
    dirname,
    path,
    get,
    tags

end # module


text:
using PyPlot, PyCall
@pyimport mpl_toolkits.basemap as basemap

include("general_config.jl")
include("regions.jl")

lat = collect(range(-90, stop=90, length=50))
lon = collect(range(0, stop=360, length=100))

plt[:figure](figsize=(12,10))
m = basemap.Basemap(projection="merc",llcrnrlat=-80,urcrnrlat=80,
            llcrnrlon=0,urcrnrlon=360,lat_ts=20,resolution="c",lon_0=180.0)

lab_lons = collect(0:30:360)
lab_lats = collect(-90:30:90)

m[:drawcoastlines](linewidth=1.5)
m[:fillcontinents](color="#888888", lake_color="aqua")
m[:drawmeridians](lab_lons, labels=repeat([true, ], outer=(length(lab_lons),)))
m[:drawparallels](lab_lats, labels=repeat([true, ], outer=(length(lab_lats),)))

#m[:drawmapboundary](fill_color="aqua")
#m[:scatter](190, 0, color="r", s=200, marker="o", latlon=true)

for k in keys(regions)
    println("Region: ", k)
    v = regions[k]
    r = []
    push!(r, m(v[1][1], v[2][1]))
    push!(r, m(v[1][1], v[2][2]))
    push!(r, m(v[1][2], v[2][2]))
    push!(r, m(v[1][2], v[2][1]))

    m[:plot]([r[1][1], r[2][1], r[3][1], r[4][1], r[1][1]], [r[1][2], r[2][2], r[3][2], r[4][2], r[1][2]], linewidth=1.0)
    plt[:text]((r[1][1] + r[3][1])/2.0, (r[1][2] + r[2][2])/2.0, k, va="center", ha="center")

end

plt[:savefig](joinpath(img_path, "regions.png"), dpi=101)
plt[:show]()

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
using JuMP, EAGO, Gurobi

m = Model(with_optimizer(EAGO.Optimizer, relaxed_optimizer = Gurobi.Optimizer(OutputFlag=0)))

# ----- Variables ----- #
x_Idx = Any[2, 3, 4, 5, 6, 7, 8, 9]
@variable(m, x[x_Idx])
JuMP.set_lower_bound(x[5], 0.0)
JuMP.set_lower_bound(x[9], 0.0)
JuMP.set_lower_bound(x[8], 0.0)
JuMP.set_lower_bound(x[7], 0.0)
JuMP.set_lower_bound(x[4], 0.0)
JuMP.set_lower_bound(x[6], 0.0)
JuMP.set_lower_bound(x[3], 0.0)
JuMP.set_lower_bound(x[2], 0.115)
JuMP.set_upper_bound(x[2], 0.115)
JuMP.set_upper_bound(x[3], 1.0)
JuMP.set_upper_bound(x[4], 1.0)
JuMP.set_upper_bound(x[5], 1.0)
JuMP.set_upper_bound(x[6], 1.0)
JuMP.set_upper_bound(x[7], 1.0)
JuMP.set_upper_bound(x[8], 1.0)
JuMP.set_upper_bound(x[9], 1.0)

# ----- Constraints ----- #
@constraint(m, e2, x[2]-0.1287*x[3]-0.1096*x[4]-0.0501*x[5]-0.1524*x[6]-0.0763*x[7]-0.1854*x[8]-0.062*x[9] == 0.0)
@constraint(m, e3, x[3]+x[4]+x[5]+x[6]+x[7]+x[8]+x[9] == 1.0)

# ----- Objective ----- #
@NLobjective(m, Min, 0.5*(42.18*x[3]*x[3]+20.18*x[3]*x[4]+10.88*x[3]*x[5]+5.3*x[3]*x[6]+12.32*x[3]*x[7]+23.84*x[3]*x[8]+17.41*x[3]*x[9]+20.18*x[4]*x[3]+70.89*x[4]*x[4]+21.58*x[4]*x[5]+15.41*x[4]*x[6]+23.24*x[4]*x[7]+23.8*x[4]*x[8]+12.62*x[4]*x[9]+10.88*x[5]*x[3]+21.58*x[5]*x[4]+25.51*x[5]*x[5]+9.6*x[5]*x[6]+22.63*x[5]*x[7]+13.22*x[5]*x[8]+4.7*x[5]*x[9]+5.3*x[6]*x[3]+15.41*x[6]*x[4]+9.6*x[6]*x[5]+22.33*x[6]*x[6]+10.32*x[6]*x[7]+10.46*x[6]*x[8]+x[6]*x[9]+12.32*x[7]*x[3]+23.24*x[7]*x[4]+22.63*x[7]*x[5]+10.32*x[7]*x[6]+30.01*x[7]*x[7]+16.36*x[7]*x[8]+7.2*x[7]*x[9]+23.84*x[8]*x[3]+23.8*x[8]*x[4]+13.22*x[8]*x[5]+10.46*x[8]*x[6]+16.36*x[8]*x[7]+42.23*x[8]*x[8]+9.9*x[8]*x[9]+17.41*x[9]*x[3]+12.62*x[9]*x[4]+4.7*x[9]*x[5]+x[9]*x[6]+7.2*x[9]*x[7]+9.9*x[9]*x[8]+16.42*x[9]*x[9]))

JuMP.optimize!(m)
run_time = backend(m).optimizer.model.optimizer._run_time
objval = objective_value(m)
println("run time: $run_time")


text:
module TakingImmutabilitySeriously

include("ImmutableArrays.jl")
using .ImmutableArrays

export freeze, thaw

include("cassette.jl")

end # module


text:
################################################################################
#
#          EllCrv/Models.jl : Different models of elliptic curves
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# (C) 2016 Tommy Hofmann
# (C) 2016 Robin Ammon
# (C) 2016 Sofia Brenner
# (C) 2022 Jeroen Hanselman
#
################################################################################

export short_weierstrass_model, is_short_weierstrass_model, simplified_model,
is_simplified_model, integral_model, is_integral_model

################################################################################
#
# Short Weierstrass model
#
################################################################################

@doc Markdown.doc"""
    short_weierstrass_model(E::EllCrv{fmpq}) ->
      (EE::EllCrv, EllCrvIso, EllCrvIso)

Transform a curve given in long Weierstrass form over QQ to short Weierstrass
form. Return short form and both transformations for points on the curve;
first transformation from E (long form) to EE (short form),
second transformation is the inverse of this map.
"""
function short_weierstrass_model(E::EllCrv)

  R = base_field(E)
  p = characteristic(R)

  if p == 2 || (p == 3 && j_invariant(E) == 0)
      error("Conversion to short weierstrass form is not possible.")
  end
  return simplified_model(E)
  #return _short_weierstrass_model(E)
end

function _short_weierstrass_model(E::EllCrv{T}) where T
  R = base_field(E)
  p = characteristic(R)

  if (p == 2) || (p == 3)
      error("Converting to short form not possible in characteristic 2 and 3")
  end

  a1, _, a3= a_invars(E)

  b2, b4, b6, b8 = b_invars(E)

  c4, c6 = c_invars(E)

  Anew = -divexact(c4, 48)
  Bnew = -divexact(c6, 864)

  EE = EllipticCurve([Anew, Bnew])::EllCrv{T}

  # we are hitting https://github.com/JuliaLang/julia/issues/15276

  _b2 = deepcopy(b2)
  _a1 = deepcopy(a1)
  _a3 = deepcopy(a3)

  # transforms a point on E (long form) to a point on EE (short form)
  trafo = function(P::EllCrvPt)

    if P.is_infinite
      return infinity(EE)
    end

    xnew = P.coordx + divexact(_b2, 12)
    ynew = P.coordy + divexact(_a1*P.coordx + _a3, 2)
    Q = EE([xnew, ynew])::EllCrvPt{T}
    return Q
  end

  # transforms a point on EE (short form) back to a point on E (long form)
  ruecktrafo = function(R::EllCrvPt)
    if R.is_infinite
        return infinity(E)
    end

    xnew = R.coordx - divexact(_b2, 12)
    ynew = R.coordy - divexact(_a1*xnew + _a3, 2)
    S = E([xnew, ynew])
    return S::EllCrvPt{T}
  end

  # type annotation necessary due to #15276
  return EE::EllCrv{T}, trafo, ruecktrafo
end

@doc Markdown.doc"""
    is_short_weierstrass_model(E::EllCrv) -> Bool

Return true if E is in short Weierstrass form.
"""
function is_short_weierstrass_model(E::EllCrv)
  return E.short
end

################################################################################
#
#  Simplified models
#
################################################################################

@doc Markdown.doc"""
    simplified_model(E::EllCrv) ->
      (EE::EllCrv, function(EllCrvPt), function(EllCrvPt))

Transform an elliptic curve to simplified Weierstrass form as defined in Connell.
Return simplified form and both transformations for points on the curve;
first transformation from E (original) to EE (simplified),
second transformation is the inverse of this map.

Returns short Weierstrass form if $\char K \neq 2, 3$,
Returns equation of the form $y^2 + xy = x^3 + a2x^2 + a6$
if $\char K = 2$  and $j(E) \neq 0$,
Returns equation of the form $y^2 + a3y = x^3 + a4x + a6$
if $\char K = 2$ and $j(E) = 0$.
Returns equation of the form $y^2 = x^3 + a2x^2 + a6$
if $\char K = 3$ and $j(E) \neq 0$,
Returns equation of the form $y^2 = x^3 + a4x + a6$
if $\char K = 3$ and $j(E) = 0$.
"""
#Magma returns minimal model if base field is QQ. Not sure if we want the same.
function simplified_model(E::EllCrv)
  K = base_field(E)
  a1, a2, a3, a4, a6 = a_invars(E)
  if characteristic(K) == 2
    if j_invariant(E) == 0
      return transform_rstu(E, [a2, 0, 0, 1])
    else
      return transform_rstu(E, [a3//a1, 0, a4//a1 + a3^2//a1^3, a1])
    end
  end

  if characteristic(K)==3
    if j_invariant(E) == 0
      return transform_rstu(E, [0, a1, a3, 1])
    else
      b2, b4 = b_invars(E)
      return transform_rstu(E, [-b4//b2, a1, a3 - a1*b4//b2, 1])
    end
  end

  b2, b4 = b_invars(E)

  return transform_rstu(E, [-b2//12, -a1//2, -a3//2 + a1*b2//24, 1])
end

@doc Markdown.doc"""
    is_simplified_model(E::EllCrv) -> Bool

Return true if E is a simplified model.
"""
function is_simplified_model(E::EllCrv)
  K = base_field(E)
  a1, a2, a3, a4, a6 = a_invars(E)
  if characteristic(K) == 2
    if j_invariant(E) == 0
      return (a1, a2) == (0, 0)
    else
      return (a1, a3, a4) == (1, 0, 0)
    end
  end

  if characteristic(K) == 3
    if j_invariant(E) == 0
      return (a1, a2, a3) == (0, 0, 0)
    else
      return (a1, a3, a4) == (0, 0, 0)
    end
  end

  return is_short_weierstrass_model(E)
end

################################################################################
#
#  Integral models
#
################################################################################

@doc Markdown.doc"""
    integral_model(E::EllCrv{fmpq}) -> (F::EllCrv{fmpz}, function, function)

Given an elliptic curve $E$ over $\mathbf Q$ in short form, returns an
isomorphic curve $F$ with model over $\mathbf Z$. The second and third
return values are the isomorpisms $E \to F$ and $F \to E$.
"""
function integral_model_old(E::EllCrv{fmpq})
  _, _, _, A, B = a_invars(E)

  mue = lcm(denominator(A), denominator(B))
  Anew = mue^4 * A
  Bnew = mue^6 * B
  E_int = EllipticCurve([Anew, Bnew])

  trafo_int = function(P) # transformes a point on E into a point on E_int
    if P.is_infinite
      return infinity(E_int)
    end

    xnew = mue^2 * P.coordx
    ynew = mue^3 * P.coordy
    Q = E_int([xnew, ynew])
    return Q
  end

  trafo_rat = function(R) # transformes a point on E_int back into a point on E
    if R.is_infinite
      return infinity(E)
    end

    xnew = divexact(R.coordx, mue^2)
    ynew = divexact(R.coordy, mue^3)
    S = E([xnew, ynew])
    return S
  end

  return E_int::EllCrv{fmpq}, trafo_int, trafo_rat
end

@doc Markdown.doc"""
    integral_model(E::EllCrv{T}) -> (F::EllCrv{T}, EllCrvIso, EllCrvIso)
      where T<:Union{fmpq, nf_elem}

Given an elliptic curve $E$ over QQ or a number field $K$, returns an
isomorphic curve $F$ with model over $\mathcal{O}_K$. The second and third
return values are the isomorpisms $E \to F$ and $F \to E$.
"""
function integral_model(E::EllCrv{T}) where T<:Union{fmpq, nf_elem}

  a1, a2, a3, a4, a6 = map(denominator, a_invars(E))
  mu = lcm(a1, a2, a3, a4, a6)
  return transform_rstu(E, [0, 0, 0, 1//mu])
end

@doc Markdown.doc"""
    is_integral_model(E::EllCrv{T}) -> Bool where T<:Union{fmpq, nf_elem}

Given an elliptic curve $E$ over QQ or a number field $K$, return
true if $E$ is an integral model of $E$.
"""
function is_integral_model(E::EllCrv{T}) where T<:Union{fmpq, nf_elem}

  a1, a2, a3, a4, a6 = map(denominator, a_invars(E))
  mu = lcm(a1, a2, a3, a4, a6)
  if mu == 1
    return true
  end

  return false
end