text:
Lar = LinearAlgebraicRepresentation
#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - PART 1 METHODS
#    frag_edge
#    intersect_edges
#    merge_vertices!
#-------------------------------------------------------------------------------

"""
	frag_edge(
		model::CAGD.Model,
		edge_idx::Int,
		bigPI::Array{Array{Int64,1},1}
	)::Tuple{Lar.Points, Lar.ChainOp}

Splits the `edge_idx`-th edge of the Edge Topology `model.T[1]`.

This method splits the `edge_idx`-th edge between the Edges in the Topology
`EV = model.T[1]` into several parts by confronting it with the others that
intersect its bounding box (see also [`CAGD.intersect_edges`](@ref)).
Returns a set made of the detected vertices over the edge (with redundancies)
and the associated cochain (without redundancies).

See also:
 - [`CAGD.planar_arrangement_1`](@ref)
 - [`CAGD.frag_edge_channel`](@ref)
It uses: [`CAGD.intersect_edges`](@ref)
---
# Examples
```jldoctest
julia> model = CAGD.Model(hcat([
	[1.0, 0.0], [0.0, 1.0], [0.0, 0.5], [0.5, 1.0], [1.0, 1.0]
]...));
julia> CAGD.addModelCells!(model, 1, [[1, 2], [2, 5], [3, 4], [4, 5]])
julia> bigPI = CAGD.spaceIndex(model);
julia> CAGD.frag_edge(model, 1, bigPI)[1]
4×2 Array{Float64,2}:
 1.0   0.0
 0.0   1.0
 0.25  0.75
 0.0   1.0
julia> Matrix(CAGD.frag_edge(model, 1, bigPI)[2])
2×4 Array{Int8,2}:
 1  0  1  0
 0  1  1  0
```
"""
function frag_edge(
		model::CAGD.Model,
		edge_idx::Int,
		bigPI::Array{Array{Int64,1},1}
	)::Tuple{Lar.Points, Lar.ChainOp}
	V = convert(Lar.Points, model.G')
	EV = model.T[1]
	alphas = Dict{Float64, Int}()
	edge = EV[edge_idx, :]
	verts = V[edge.nzind, :]
	for i in bigPI[edge_idx]
		if i != edge_idx # && edge_idx in bigPI[i]
			intersection = CAGD.intersect_edges(V, edge, EV[i, :])
			for (point, alpha) in intersection
				verts = [verts; point]
				alphas[alpha] = size(verts, 1)
			end
		end
	end
	alphas[0.0], alphas[1.0] = [1, 2]
	alphas_keys = sort(collect(keys(alphas)))
	edge_num = length(alphas_keys)-1
	verts_num = size(verts, 1)
	ev = SparseArrays.spzeros(Int8, edge_num, verts_num)
	for i in 1:edge_num
		ev[i, alphas[alphas_keys[i]]] = 1
		ev[i, alphas[alphas_keys[i+1]]] = 1
	end
	return verts, ev
end

"""
	intersect_edges(
		V::Lar.Points,
		edge1::Lar.Cell,
		edge2::Lar.Cell
	)::Array{Tuple{Lar.Points, Float64}, 1}

Finds the intersection points (if there exist) between the two given edges.

Compute the new points over `edge1` given by the intersection with `edge2`.
For each intersection point it evaluates the scalar corresponding to the
relative distance from the first vertex of the edge.
Note that vertices of `edge1` are therefore never detected (see example 2).

See also: [`CAGD.frag_edge`](@ref)
---
# Examples
```jldoctest
# Cross
julia> V = [1.0 0.0; 0.0 1.0; 0.0 0.5; 0.5 1.0];
julia> copEV = Lar.coboundary_0([[1, 2], [3, 4]]);
julia> CAGD.intersect_edges(V, copEV[1, :], copEV[2, :])
1-element Array{Tuple{Array{T,2} where T,Float64},1}:
 ([0.25 0.75], 0.75)
```

```jldoctest
# Collinear
julia> V = [1.0 0.0; 0.0 1.0; 0.75 0.25; 0.5 0.5];
julia> copEV = Lar.coboundary_0([[1, 2], [3, 4]]);
julia> CAGD.intersect_edges(V, copEV[1, :], copEV[2, :])
2-element Array{Tuple{Array{T,2} where T,Float64},1}:
 ([0.75 0.25], 0.25)
 ([0.5 0.5], 0.5)
julia> CAGD.intersect_edges(V, copEV[2, :], copEV[1, :])
0-element Array{Tuple{Array{T,2} where T,Float64},1}
```
"""
function intersect_edges(
		V::Lar.Points,
		edge1::Lar.Cell,
		edge2::Lar.Cell
	)::Array{Tuple{Lar.Points, Float64}, 1}

	err = 10e-8

	x1, y1, x2, y2 = vcat(map(c->V[c, :], edge1.nzind)...)
	x3, y3, x4, y4 = vcat(map(c->V[c, :], edge2.nzind)...)
	ret = Array{Tuple{Lar.Points, Float64}, 1}()

	v1 = [x2-x1, y2-y1]
	v2 = [x4-x3, y4-y3]
	v3 = [x3-x1, y3-y1]
	ang1 = dot(LinearAlgebra.normalize(v1), LinearAlgebra.normalize(v2))
	ang2 = dot(LinearAlgebra.normalize(v1), LinearAlgebra.normalize(v3))
	parallel = 1-err < abs(ang1) < 1+err
	colinear = parallel && (1-err < abs(ang2) < 1+err || -err < norm(v3) < err)
	if colinear
		o = [x1 y1]
		v = [x2 y2] - o
		alpha = 1/dot(v,v')
		ps = [x3 y3; x4 y4]
		for i in 1:2
			a = alpha*dot(v',(reshape(ps[i, :], 1, 2)-o))
			if 0 < a < 1
				push!(ret, (ps[i:i, :], a))
			end
		end
	elseif !parallel
		denom = (v2[2])*(v1[1]) - (v2[1])*(v1[2])
		a = ((v2[1])*(-v3[2]) - (v2[2])*(-v3[1])) / denom
		b = ((v1[1])*(-v3[2]) - (v1[2])*(-v3[1])) / denom

		if -err < a < 1+err && -err <= b <= 1+err
			p = [(x1 + a*(x2-x1))  (y1 + a*(y2-y1))]
			push!(ret, (p, a))
		end
	end
	return ret
end

"""
merge_vertices!(
	model::CAGD.Model,
	[ edge_map::Array{Array{Int64,1},1} = Array{Array{Int64,1},1}() ],
	[ err::Float64 = 1e-4 ]
)::Nothing

Modify `model` by compacting its points closer than `err` in a single one.

This method check one at time each vertex ``v`` in `model.G` and identifies
each other vertex within `err` with ``v`` itself.
The edges cochain `model.T[1]` is coherently modified (multiple edges between
two vertices are not allowed).
If an `edge_map` is given in input (this could be usefull during the planar
arrangement), then also the map is coherently modified and given back in output.

See also: [`CAGD.planar_arrangement_1`](@ref)
---
# Examples
```jldoctest
julia> model = CAGD.Model([
	0.5 0.0 0.5 1.0 0.5 1.0
	0.5 0.0 0.5 1.0 0.5 1.0
]);
julia> CAGD.addModelCells!(model, 1, [[1, 4], [3, 2], [5, 6], [1, 6], [5, 3]]);
julia> CAGD.merge_vertices!(model)

julia> model.G
2×3 Array{Float64,2}:
 0.5  0.0  1.0
 0.5  0.0  1.0

julia> Matrix(model.T[1])
2×3 Array{Int8,2}:
 1  0  1
 1  1  0
```
"""
function merge_vertices!(
		model::CAGD.Model,
		edge_map::Array{Array{Int64,1},1} = Array{Array{Int64,1},1}(),
		err::Float64 = 1e-4
	)::Nothing
	V = convert(Lar.Points, model.G')
	EV = model.T[1]
	vertsnum = size(V, 1)
	edgenum = size(EV, 1)
	newverts = zeros(Int, vertsnum)
	# KDTree constructor needs an explicit array of Float64
	V = Array{Float64,2}(V)
	kdtree = KDTree(permutedims(V))

	# merge congruent vertices
	todelete = []
	i = 1
	for vi in 1:vertsnum
		if !(vi in todelete)
			nearvs = Lar.inrange(kdtree, V[vi, :], err)
			newverts[nearvs] .= i
			nearvs = setdiff(nearvs, vi)
			todelete = union(todelete, nearvs)
			i = i + 1
		end
	end
	nV = V[setdiff(collect(1:vertsnum), todelete), :]

	# merge congruent edges
	edges = Array{Tuple{Int, Int}, 1}(undef, edgenum)
	oedges = Array{Tuple{Int, Int}, 1}(undef, edgenum)
	for ei in 1:edgenum
		v1, v2 = EV[ei, :].nzind
		edges[ei] = Tuple{Int, Int}(sort([newverts[v1], newverts[v2]]))
		oedges[ei] = Tuple{Int, Int}(sort([v1, v2]))
	end
	nedges = union(edges)
	nedges = filter(t->t[1]!=t[2], nedges)
	nedgenum = length(nedges)
	nEV = spzeros(Int8, nedgenum, size(nV, 1))
	# maps pairs of vertex indices to edge index
	etuple2idx = Dict{Tuple{Int, Int}, Int}()
	# builds `edge_map`
	for ei in 1:nedgenum
		nEV[ei, collect(nedges[ei])] .= 1
		etuple2idx[nedges[ei]] = ei
	end
	if !isempty(edge_map)
		for i in 1:length(edge_map)
			row = edge_map[i]
			row = map(x->edges[x], row)
			row = filter(t->t[1]!=t[2], row)
			row = map(x->etuple2idx[x], row)
			edge_map[i] = row
		end
	end
	# return new vertices and new edges
	model.G = convert(Lar.Points, nV')
	model.T[1] = nEV
	model.T[2] = SparseArrays.spzeros(Int8, 0, size(model, 1, 1))
	return
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - CLEAN DECOMPOSITION
#-------------------------------------------------------------------------------

"""
	cleandecomposition(
		model::CAGD.Model,
		sigma::Union{Lar.Chain, Array{Int,1}},
		[ edge_map::Array{Array{Int64,1},1} = Array{Array{Int64,1},1}() ]
	)

This function clears the `model` from all edges outside ``σ``.

This function take an arranged-edge `model` and a `Lar.Chain` ``σ`` made of
edges indices of the 1-chains `model.T[1]`.
The method drops from the Topology all the edges that are outside ``σ`` and
gives back the dropped model.
If an `edge_map` is given as an input it is also purged from the edges outside
w.r.t. ``σ`` and it is given as an output as well.

See also: [`CAGD.planar_arrangement`](@ref)
---
# Examples
```jldoctest
# Nested Triangles
julia> model = CAGD.Model([
	0.0  2.0  4.0  1.0  3.0  2.0  2.0 -2.0  6.0
	0.0  0.0  0.0  1.5  1.5  3.0 -3.0  3.0  3.0
]);
julia> CAGD.addModelCells!(model, 1, [
	[1, 3], [1, 6], [3, 6], [2, 4], [2, 5], [4, 5], [7, 8], [7, 9], [8, 9]
]);

julia> σ = SparseArrays.sparse([1; 1; 1; 0; 0; 0; 0; 0; 0]);

julia> model = CAGD.cleandecomposition(model, convert(Lar.Chain, σ))

julia> model.G
2×6 Array{Float64,2}:
 0.0  2.0  4.0  1.0  3.0  2.0
 0.0  0.0  0.0  1.5  1.5  3.0

julia> Matrix(model.T[1])
6×6 Array{Int8,2}:
 1  0  1  0  0  0
 1  0  0  0  0  1
 0  0  1  0  0  1
 0  1  0  1  0  0
 0  1  0  0  1  0
 0  0  0  1  1  0
```
"""
function cleandecomposition(
		model::CAGD.Model,
		sigma::Union{Lar.Chain, Array{Int,1}},
		edge_map::Array{Array{Int64,1},1} = Array{Array{Int64,1},1}()
	)::Union{CAGD.Model, Tuple{CAGD.Model, Array{Array{Int64,1},1}}}

	model = deepcopy(model)
	edge_map = deepcopy(edge_map)
	# Model point extraction to use Lar.point_in_face                           ##
	V = convert(Lar.Points, model.G')

	if issparse(sigma)
		sigma = sigma.nzind
	end
	todel = Array{Int,1}()
	sigma_edges = model.T[1][sigma, :]
	for e in 1 : model.T[1].m
		# Since the model is already arranged an edge can only be a sigma-edge,
		#  an intersection-free innner edge w.r.t. sigma face (not to purge),
		#  or an intersection-free outer edge (has to be deleted).
		if !(e in sigma)
			v1, v2 = CAGD.getModelCellVertices(model, 1, e)
			centroid = .5*(v1 + v2)

			if ! CAGD.point_in_face(centroid, V, sigma_edges)
				push!(todel, e)
			end
		end
	end

	CAGD.deleteModelCells!(model, 1, todel)
	CAGD.modelPurge!(model)

	if isempty(edge_map)
		return model
	end

	# Edges in edge_map must be updated too according to the same approach
	#  edges are deleted from models.
	for i in reverse(todel)
		for row in edge_map

			filter!(x->x!=i, row)

			for j in 1:length(row)
				if row[j] > i
					row[j] -= 1
				end
			end
		end
	end
	return model, edge_map
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - BICONNECTED COMPONENTS
#-------------------------------------------------------------------------------

"""
	biconnected_components(EV::Lar.ChainOp)

Compute sets of 2-cells on the same 3-cells biconnected components.

The method evaluates the 2-cells ``σ_i`` wich lies on the same 3-cells
biconnected component ``χ_j`` and gives back an array that contains an array
of ``σ_i`` for each ``j``.
Do note that if the cochain `EV` contains more copies of the same 2-cell then
it will be considered like a 3-cell.

Intermediate part of Planar Arrangement's algorithmic pipeline.

See also:
 - [`CAGD.planar_arrangement`](@ref) for the complete pipeline.
 - [`CAGD.planar_arrangement_2`](@ref).
---

# Examples
```jldoctest
julia> copEV = SparseArrays.sparse(Array{Int8, 2}([
		[1 1 0 0 0 0] #1 -> 1,2  |
		[1 0 1 0 0 0] #2 -> 1,3  |
		[1 0 0 1 0 0] #3 -> 1,4   |
		[1 0 0 0 1 0] #4 -> 1,5   |
		[1 0 0 0 0 1] #5 -> 1,6
		[0 1 1 0 0 0] #6 -> 2,3  |
		[0 0 0 1 1 0] #7 -> 4,5   |
	]));
julia> CAGD.biconnected_components(copEV)
2-element Array{Array{Int64,1},1}:
 [2, 6, 1]
 [4, 7, 3]
```

```jldoctest
julia> copEV = SparseArrays.sparse(Array{Int8, 2}([
		[1 1 0] #1 -> 1,2  |
		[1 1 0] #2 -> 1,2  |
		[1 0 1] #3 -> 1,2
	]));
julia> CAGD.biconnected_components(copEV)
1-element Array{Array{Int64,1},1}:
 [2, 1]
```
"""
function biconnected_components(EV::Lar.ChainOp)

	ps = Array{Tuple{Int, Int, Int}, 1}()
	es = Array{Tuple{Int, Int}, 1}()
	todel = Array{Int, 1}()
	visited = Array{Int, 1}()
	bicon_comps = Array{Array{Int, 1}, 1}()
	hivtx = 1

	function an_edge(point) # TODO: fix bug
		# error? : BoundsError: attempt to access 0×0 SparseMatrix ...
		edges = setdiff(EV[:, point].nzind, todel)
		if length(edges) == 0
			edges = [false]
		end
		edges[1]
	end

	function get_head(edge, tail)
		setdiff(EV[edge, :].nzind, [tail])[1]
	end

	function v_to_vi(v)
		i = findfirst(t->t[1]==v, ps)
		# seems findfirst changed from 0 to Nothing
		if typeof(i) == Nothing
			return false
		elseif i == 0
			return false
		else
			return ps[i][2]
		end
	end

	push!(ps, (1,1,1))
	push!(visited, 1)
	exit = false
	while !exit
		edge = an_edge(ps[end][1])
		if edge != false
			tail = ps[end][2]
			head = get_head(edge, ps[end][1])
			hi = v_to_vi(head)
			if hi == false
				hivtx += 1
				push!(ps, (head, hivtx, ps[end][2]))
				push!(visited, head)
			else
				if hi < ps[end][3]
					ps[end] = (ps[end][1], ps[end][2], hi)
				end
			end
			push!(es, (edge, tail))
			push!(todel, edge)
		else
			if length(ps) == 1
				found = false
				pop!(ps)
				for i in 1:size(EV,2)
					if !(i in visited)
						hivtx = 1
						push!(ps, (i, hivtx, 1))
						push!(visited, i)
						found = true
						break
					end
				end
				if !found
					exit = true
				end

			else
				if ps[end][3] == ps[end-1][2]
					edges = Array{Int, 1}()
					while true
						edge, tail = pop!(es)
						push!(edges, edge)
						if tail == ps[end][3]
							if length(edges) > 1
								push!(bicon_comps, edges)
							end
							break
						end
					end

				else
					if ps[end-1][3] > ps[end][3]
						ps[end-1] = (ps[end-1][1], ps[end-1][2], ps[end][3])
					end
				end
				pop!(ps)
			end
		end
	end
	bicon_comps = sort(bicon_comps, lt=(x,y)->length(x)>length(y))
	return bicon_comps
end

function monoconnected_components(EV::Lar.ChainOp)::Array{Array{Int,1},1}
	mono = Array{Array{Int,1},1}()
	tocheck = collect(1 : EV.m)
	while(!isempty(tocheck))
		vs = Array{Int,1}()
		vs_i = Array{Int,1}()
		vs_o = Array{Int,1}()
		evs = [tocheck[1]]
		evs_i = [tocheck[1]]
		evs_o = Array{Int,1}()

		while(!isempty(vs_i) || !isempty(evs_i))
			if !isempty(evs_i)
				e = evs_i[1]
				push!(evs_o, e);
				push!(vs, EV[e, :].nzind...)
			end
			vs_i = setdiff(vs, vs_o)
			if !isempty(vs_i)
				v = vs_i[1]
				push!(vs_o, v);
				push!(evs, EV[:, v].nzind...)
			end
			evs_i = setdiff(evs, evs_o)
		end

		push!(mono, unique(evs))
	end
	return mono
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - REMOVE MAPPING DANGLING EDGES
#-------------------------------------------------------------------------------

"""
	remove_mapping_dangling_edges(
		edge_map::Array{Array{Int64,1},1},
		bicon_comps::Array{Array{Int64,1},1}
	)::Array{Array{Int64,1},1}

Remove from the edge_map the edges not belonging to a biconnected component.

This utility function deletes from the `edge_map` all the edges that do not
belongs to a biconnected component by ordinatelly recompatting the indices of
the remaning edges.

See also: [`CAGD.planar_arrangement`](@ref).
"""
function remove_mapping_dangling_edges(
		edge_map::Array{Array{Int64,1},1},
		bicon_comps::Array{Array{Int64,1},1}
	)::Array{Array{Int64,1},1}

	edge_map = deepcopy(edge_map)

	edges = sort(union(bicon_comps...))
								  #size(model.T[1, 1]) == max(edge_map...)
	todel = sort(setdiff(collect(1:max(max(edge_map...)...)), edges))

	for i in reverse(todel)
		for row in edge_map

			filter!(x->x!=i, row)

			for j in 1:length(row)
				if row[j] > i
					row[j] -= 1
				end
			end
		end
	end

	return edge_map
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - PART 2 METHODS
#    componentgraph (aka TGW) from tgw2.jl
#    cell_merging
#-------------------------------------------------------------------------------

"""
	cell_merging(n, containment_graph, V, EVs, boundaries, shells, shell_bboxes)

Cells composing for the Topological Gift Wrapping algorithm.

This is the online part of the TGW algorithm.

See also: [`CAGD.planar_arrangement_2`](@ref).
"""
function cell_merging(
		n, containment_graph, V, EVs, boundaries, shells, shell_bboxes
	)
	function bboxes(V::Lar.Points, indexes::Lar.ChainOp)
		boxes = Array{Tuple{Any, Any}}(undef, indexes.n)
		for i in 1:indexes.n
			v_inds = indexes[:, i].nzind
			boxes[i] = Lar.bbox(V[v_inds, :])
		end
		boxes
	end
	# initialization
	sums = Array{Tuple{Int, Int, Int}}(undef, 0)
	# assembling child components with father components
	for father in 1:n
		if sum(containment_graph[:, father]) > 0
			father_bboxes = bboxes(
				V, abs.(EVs[father]')*abs.(boundaries[father]')
			)
			for child in 1:n
				if containment_graph[child, father] > 0
					child_bbox = shell_bboxes[child]
					for b in 1:length(father_bboxes)
						if Lar.bbox_contains(father_bboxes[b], child_bbox)
							push!(sums, (father, b, child))
							break
						end
					end
				end
			end
		end
	end
	# offset assembly initialization
	EV = vcat(EVs...)
	edgenum = size(EV, 1)
	facenum = sum(map(x->size(x,1), boundaries))
	FE = spzeros(Int8, facenum, edgenum)
	shells2 = spzeros(Int8, length(shells), edgenum)
	r_offsets = [1]
	c_offset = 1
	# submatrices construction
	for i in 1:n
		min_row = r_offsets[end]
		max_row = r_offsets[end] + size(boundaries[i], 1) - 1
		min_col = c_offset
		max_col = c_offset + size(boundaries[i], 2) - 1
		FE[min_row:max_row, min_col:max_col] = boundaries[i]
		shells2[i, min_col:max_col] = shells[i]
		push!(r_offsets, max_row + 1)
		c_offset = max_col + 1
	end
	# offsetting assembly of component submatrices
	for (f, r, c) in sums
		FE[r_offsets[f]+r-1, :] += shells2[c, :]
	end

	return EV, FE
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - PART 1
#-------------------------------------------------------------------------------

"""
	planar_arrangement_1(
		model::CAGD.Model,
		[sigma::Lar.Chain = spzeros(Int8, 0)],
		[return_edge_map::Bool = false]
	)

First part of arrangement's algorithmic pipeline.

This function computes the pairwise intersection between each edge of a given
2D cellular complex 1-skeleton. The computation is speeded up via the evaluation
of the Spatial Index. See [`CAGD.spaceindex`](@ref).

See also: [`CAGD.planar_arrangement`](@ref) for the complete pipeline.
It uses:
 - [`CAGD.frag_edges_channel`](@ref)
 - [`CAGD.frag_edges`](@ref)
 - [`CAGD.merge_vertices!`](@ref)
---
# Examples
```jldoctest
julia> model = CAGD.Model([
	0.0 0.5 0.0 0.5 0.3 1.0 0.3 1.0;
	0.0 0.0 1.0 1.0 0.5 0.5 1.0 1.0
]);
julia> CAGD.addModelCells!(model, 1, [
	[1, 2], [3, 4], [1, 3], [2, 4], [5, 6], [7, 8], [5, 7], [6, 8]
]);
julia> (model, edge_map) =
	CAGD.planar_arrangement_1(model, spzeros(Int8, 0), true);
julia> model.G
2×9 Array{Float64,2}:
 0.0  0.5  0.0  0.5  0.3  0.5  0.3  1.0  1.0
 0.0  0.0  1.0  1.0  1.0  0.5  0.5  0.5  1.0

julia> Matrix(model.T[1])
11×9 Array{Int8,2}:
 1  1  0  0  0  0  0  0  0
 0  0  1  1  0  0  0  0  0
 0  0  0  1  1  0  0  0  0
 1  0  1  0  0  0  0  0  0
 0  1  0  0  0  1  0  0  0
 0  0  0  0  1  1  0  0  0
 0  0  0  0  0  1  1  0  0
 0  0  0  0  0  1  0  1  0
 0  0  0  0  1  0  0  0  1
 0  0  0  1  0  0  1  0  0
 0  0  0  0  0  0  0  1  1

julia> edge_map
8-element Array{Array{Int64,1},1}:
 [1]
 [2, 3]
 [4]
 [5, 6]
 [7, 8]
 [3, 9]
 [10]
 [11]

```
"""
function planar_arrangement_1(
		model::CAGD.Model,
		sigma::Lar.Chain = spzeros(Int8, 0),
		return_edge_map::Bool = false
	)::Union{
		CAGD.Model,
		Tuple{CAGD.Model, Array{Array{Int64,1},1}},
		Tuple{CAGD.Model, Array{Int,1}},
		Tuple{CAGD.Model, Array{Int,1}, Array{Array{Int64,1},1}}
	}
	# data structures initialization
	# V = convert(Lar.Points, model.G')
	# copEV = model.T[1]
	edgenum = size(model, 1, 1)
	edge_map = Array{Array{Int, 1}, 1}(undef,edgenum)
	rV = Lar.Points(zeros(0, 2))
	rEV = SparseArrays.spzeros(Int8, 0, 0)
	finalcells_num = 0

	# spaceindex computation
	bigPI = CAGD.spaceIndex(model, 1)

	# sequential (iterative) processing of edge fragmentation
	for i in 1:edgenum
		v, ev = CAGD.frag_edge(model, i, bigPI)
		newedges_nums = map(x->x+finalcells_num, collect(1:size(ev, 1)))
		edge_map[i] = newedges_nums
		finalcells_num += size(ev, 1)
		rV = convert(Lar.Points, rV)
		rV, rEV = Lar.skel_merge(rV, rEV, v, ev)
	end

	# merging of close vertices and edges (2D congruence)
	model = CAGD.Model(convert(Lar.Points, rV'))
	CAGD.addModelCells!(model, 1, rEV)
	CAGD.merge_vertices!(model, edge_map)

	if isempty(sigma)
		if return_edge_map
			return model, edge_map
		else
			return model
		end
	end

	# Sigma update by mapping each of its edges into newly born edges
	new_sigma = Array{Int,1}()
	map(i->new_sigma=union(new_sigma, edge_map[i]), sigma.nzind)
	if return_edge_map
		return model, new_sigma, edge_map
	else
		return model, new_sigma
	end
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - PART 2
#-------------------------------------------------------------------------------

"""
	planar_arrangement_2(
		model::CAGD.Model,
		bicon_comps::{Array{Array,1},1}
	)::CAGD.Model

Second part of arrangement's algorithmic pipeline.

This function is the complete Topological Gift Wrapping (TGW) algorithm that
is firstly locally used in order to decompose the 2-cells and then globally
to generate the 3-cells of the arrangement of the ambient space ``E^3``.

During this process each dangling 1-cell is removed.
Do note that the isolated 2-cells are not removed by this procedure.

See also: [`CAGD.planar_arrangement`](@ref) for the complete pipeline.
It uses:
 - [`CAGD.biconnected_components`](@ref)
 - [`CAGD.componentgraph`](@ref)
 - [`CAGD.cell_merging`](@ref)
---

# Examples
```jldoctest
# Triforce
julia> model = CAGD.Model([
	0.0 0 2 0 0 1 4 4 3 2 4 4
	0.0 2 0 3 4 4 2 0 0 4 4 2
]);

julia> CAGD.addModelCells!(model, 1, [
	[ 1,  2], [ 2,  3], [ 3,  1],
	[ 4,  5], [ 5,  6], [ 6,  7], [ 7,  8], [ 8,  9], [ 9, 4],
	[10, 11], [11, 12], [12, 10]
])

julia> bicon_comps = CAGD.biconnected_components(model.T[1])
3-element Array{Array{Int64,1},1}:
 [9, 8, 7, 6, 5, 4]
 [3, 2, 1]
 [12, 11, 10]

julia> model = CAGD.planar_arrangement_2(model, bicon_comps);

julia> Matrix(model.T[2])
3×12 Array{Int8,2}:
 -1  -1  -1  -1  -1  1   0   0  0   0   0  0
  0   0   0   0   0  0  -1  -1  1   0   0  0
  0   0   0   0   0  0   0   0  0  -1  -1  1
```
"""
function planar_arrangement_2(
		model::CAGD.Model,
		bicon_comps::Array{Array{Int64,1},1}
	)::CAGD.Model

	#==
	test to remove
	@assert isequal(
		bicon_comps,
		CAGD.biconnected_components(model.T[1])
	)
	==#

	# component graph
	n, containment_graph, V, EVs, boundaries, shells, shell_bboxes =
		CAGD.componentgraph(
			convert(Lar.Points, model.G'), model.T[1], bicon_comps
		)

	copEV, FE = CAGD.cell_merging(
	   	n, containment_graph, V, EVs, boundaries, shells, shell_bboxes
	)

	model = CAGD.Model(convert(Lar.Points, V'))
	CAGD.addModelCells!(model, 1, copEV)
	CAGD.addModelCells!(model, 2, FE)
	CAGD.modelPurge!(model)
	return model
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE
#-------------------------------------------------------------------------------

"""
	planar_arrangement(
		model::CAGD.Model,
		[sigma::Lar.Chain],
		[return_edge_map::Bool]
	)::Union{CAGD.Model, Tuple{CAGD.Model, Array{Array{Int, 1}, 1}}}

Compute the arrangement on the given cellular complex 1-skeleton in 2D.

A cellular complex is arranged when the intersection of every possible pair
of cell of the complex is empty and the union of all the cells is the
whole Euclidean space.
The methods default input/output is a `CAGD.Model` with at least Geometry and
1-cell Topoly.
Mutiprocessing is by default disabled.

It uses:
 - [`CAGD.planar_arrangement_1`](@ref)
 - [`CAGD.cleandecomposition`](@ref)
 - [`CAGD.biconnected_components`](@ref)
 - [`CAGD.remove_mapping_dangling_edges`](@ref)
 - [`CAGD.planar_arrangement_2`](@ref)

## Additional arguments:
 - `sigma::Chain`: if specified, the method return only cells contained inside
		this cell. _Default_: empty cell.
 - `return_edge_map::Bool`: builds an `edge_map` which maps the input edges
		to the output's. _Defaults_: `false`.
"""
function planar_arrangement(
		model::CAGD.Model,
		sigma::Lar.Chain = spzeros(Int8, 0),                                    ## TEST
		return_edge_map::Bool = false
	)::Union{CAGD.Model, Tuple{CAGD.Model, Array{Array{Int, 1}, 1}}}

	# Check dimensional condition
	@assert length(model) == 2

	#planar_arrangement_1
	if isempty(sigma)
		model, edge_map =
			CAGD.planar_arrangement_1(model, sigma, true)
	else
		model, sigma, edge_map =
			CAGD.planar_arrangement_1(model, sigma, true)
	end

#   print("PLANAR_ARRANGEMENT_1 COMPLETE\n")

	# cleandecomposition
	if !isempty(sigma)
		model, edge_map = CAGD.cleandecomposition(model, sigma, edge_map)
#       print("DECOMPOSITION CLEANED\n")
	end

	# generation of biconnected components
	bicon_comps = CAGD.biconnected_components(model.T[1])
	#mocon_comps = CAGD.monoconnected_components(model.T[1])

#	print("BICONNECTED COMPONENTS FOUND\n")

	if isempty(bicon_comps)
		println("No biconnected components found.\n")
		if (return_edge_map)
			return (model(), nothing)
		else
			return model()
		end
	end

	# Removing Dangling edges from edge map
	if return_edge_map
		@assert size(model, 1, 1) == max((edge_map...)...)                           # CTR
		edge_map = CAGD.remove_mapping_dangling_edges(edge_map, bicon_comps)
	end

#	print("DANGLING EDGES REMOVED\n")

	# Planar_arrangement_2
	model = CAGD.planar_arrangement_2(model, bicon_comps)

#	print("PLANAR ARRANGEMENT 2 COMPLETE\n")

	if size(model, 0, 2) - size(model, 1, 1) + size(model, 2, 1) != 1
		@show size(model, 0, 2) - size(model, 1, 1) + size(model, 2, 1)
	end

	if return_edge_map
		 return model, edge_map
	else
		 return model
	end
end

"""
	point_in_face(point, V::Points, copEV::ChainOp)
Check if `point` is not outside the area of the face bounded by the edges in `copEV`
"""
function point_in_face(point, V::Lar.Points, copEV::Lar.ChainOp)

	function pointInPolygonClassification(V,EV)

		function crossingTest(new, old, status, count)
			if status == 0
				status = new
				return status, (count + 0.5)
			else
				if status == old
					return 0, (count + 0.5)
				else
					return 0, (count - 0.5)
				end
			end
		end

		# Outputs a tileCoder that takes as input a point and classify its
		#  octant position with respect to the supercrescent list N-S-E-W:
		#  1001 | 0001 | 0101
		#  -----+------+-----
		#  1000 | 0000 | 0100
		#  -----+------+-----
		#  0101 | 0010 | 0110
		function setTile(box)
			tiles = [[9,1,5],[8,0,4],[10,2,6]]
			b1,b2,b3,b4 = box
			function tileCode(point)
				x,y = point
				code = 0
				if y>b1 code=code|1 end
				if y<b2 code=code|2 end
				if x>b3 code=code|4 end
				if x<b4 code=code|8 end
				return code
			end
			return tileCode
		end

		function pointInPolygonClassification0(pnt)
			# Set the tile coder w.r.t. the point we are checking
			x,y = pnt
			xmin,xmax,ymin,ymax = x,x,y,y
			tilecode = setTile([ymax,ymin,xmax,xmin])
			count,status = 0,0

			for k in 1:EV.m
				edge = EV[k,:]
				p1, p2 = V[edge.nzind[1], :], V[edge.nzind[2], :]
				(x1,y1),(x2,y2) = p1,p2
				c1,c2 = tilecode(p1),tilecode(p2)
				# For each edge we identify via:
				#  - xor: the different sides
				#  -  or: all the direction involved
				#  - and: the cardinal direction (if unique)
				c_edge, c_un, c_int = xor(c1, c2), c1|c2, c1&c2

				if (c_edge == 0) & (c_un == 0) return "p_on"
				elseif (c_edge == 12) & (c_un == c_edge) return "p_on"
				elseif c_edge == 3
					if c_int == 0 return "p_on"
					elseif c_int == 4 count += 1 end
				elseif c_edge == 15
					x_int = ((y-y2)*(x1-x2)/(y1-y2))+x2
					if x_int > x count += 1
					elseif x_int == x return "p_on" end
				elseif (c_edge == 13) & ((c1==4) | (c2==4))
						status, count = crossingTest(1,2,status,count)
				elseif (c_edge == 14) & ((c1==4) | (c2==4))
						status, count = crossingTest(2,1,status,count)
				elseif c_edge == 7 count += 1
				elseif c_edge == 11 count = count
				elseif c_edge == 1
					if c_int == 0 return "p_on"
					elseif c_int == 4
						status, count = crossingTest(1,2,status,count)
					end
				elseif c_edge == 2
					if c_int == 0 return "p_on"
					elseif c_int == 4
						status, count = crossingTest(2,1,status,count)
					end
				elseif (c_edge == 4) & (c_un == c_edge) return "p_on"
				elseif (c_edge == 8) & (c_un == c_edge) return "p_on"
				elseif c_edge == 5
					if (c1==0) | (c2==0) return "p_on"
					else
						status, count = crossingTest(1,2,status,count)
					end
				elseif c_edge == 6
					if (c1==0) | (c2==0) return "p_on"
					else
						status, count = crossingTest(2,1,status,count)
					end
				elseif (c_edge == 9) & ((c1==0) | (c2==0)) return "p_on"
				elseif (c_edge == 10) & ((c1==0) | (c2==0)) return "p_on"
				end
			end

			if (round(count)%2)==1
				return "p_in"
			else
				return "p_out"
			end
		end
		return pointInPolygonClassification0
	end

	ret = pointInPolygonClassification(V, copEV)(point)
	return "p_in" == ret || ret == "p_on"
end


text:
rrule(::Type{<:Diagonal}, d::AbstractVector) = Diagonal(d), Rule(diag)
rrule(::typeof(diag), A::AbstractMatrix) = diag(A), Rule(Diagonal)

ext: jl


text:
using Embeddings
using Base.Test

"""
    tempdatadeps(fun)

Run the function and delete all created datadeps afterwards
"""
function tempdatadeps(fun)
    tempdir = mktempdir()
    try
        info("sending all datadeps to $tempdir")
        withenv("DATADEPS_LOAD_PATH"=>tempdir) do
            fun()
        end
    finally
        try
            info("removing $tempdir")
            rm(tempdir, recursive=true, force=true)
        catch err
            warn("Something went wrong with removing $tempdir")
            warn(err)
        end
    end
end

# uncomment the below to not use tempdatadeps (i.e. if not debugging)
# tempdatadeps(fun) = fun()

"""
@testset_nokeep_data

Use just like @testset,
but know that it deletes any downloaded data dependencies when it is done.
"""
macro testset_nokeep_data(name, expr)
    quote
        tempdatadeps() do
            @testset $name $expr
        end
    end |> esc
end

@testset_nokeep_data "Word2Vec" begin
    embs_full = load_embeddings(Word2Vec)

    @test size(embs_full.embeddings) == (300, length(embs_full.vocab))

    embs_mini = load_embeddings(Word2Vec; max_vocab_size=100)
    @test length(embs_mini.vocab)==100

    @test embs_mini.embeddings == embs_full.embeddings[:, 1:100]
    @test embs_mini.vocab == embs_full.vocab[1:100]

    @test "for" ∈ embs_mini.vocab

    embs_specific =  load_embeddings(Word2Vec; keep_words=Set(["red", "green", "blue"]))

    @test size(embs_specific.embeddings) == (300, 3)
    @test Set(embs_specific.vocab) == Set(["red", "green", "blue"])
end

@testset "FastText" begin
    @testset_nokeep_data "English 1" begin
        @testset "Basic" begin
            embs1 = load_embeddings(FastText_Text; max_vocab_size=100)
            @test length(embs1.vocab)==100
            @test size(embs1.embeddings) == (300, 100)
        end

        @testset "Specific" begin
            embs_specific =  load_embeddings(FastText_Text; keep_words=Set(["red", "green", "blue"]))
            @test size(embs_specific.embeddings) == (300, 3)
            @test Set(embs_specific.vocab) == Set(["red", "green", "blue"])
        end
    end

    @testset_nokeep_data "French" begin
        embs_fr = load_embeddings(FastText_Text{:fr}; max_vocab_size=100)
        @test length(embs_fr.vocab)==100
        @test size(embs_fr.embeddings) == (300, 100)
    end

    @testset_nokeep_data "English file number 2" begin
        embs_specific =  load_embeddings(FastText_Text, 2; keep_words=Set(["red", "green", "blue"]))
        @test size(embs_specific.embeddings) == (300, 3)
        @test Set(embs_specific.vocab) == Set(["red", "green", "blue"])
    end

end


text:
using Mustache
using Test

@testset " comments " begin

	## Comment blocks should be removed from the template.
tpl = """12345{{! Comment Block! }}67890"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """1234567890"""

	## Multiline comments should be permitted.
tpl = """12345{{!
  This is a
  multi-line comment...
}}67890
"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """1234567890
"""

	## All standalone comment lines should be removed.
tpl = """Begin.
{{! Comment Block! }}
End.
"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """Begin.
End.
"""

	## All standalone comment lines should be removed.
tpl = """Begin.
  {{! Indented Comment Block! }}
End.
"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """Begin.
End.
"""

	## "\r\n" should be considered a newline for standalone tags.
tpl = """|
{{! Standalone Comment }}
|"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """|
|"""

	## Standalone tags should not require a newline to precede them.
tpl = """  {{! I'm Still Standalone }}
!"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """!"""

	## Standalone tags should not require a newline to follow them.
tpl = """!
  {{! I'm Still Standalone }}"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """!
"""

	## All standalone comment lines should be removed.
tpl = """Begin.
{{!
Something's going on here...
}}
End.
"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """Begin.
End.
"""

	## All standalone comment lines should be removed.
tpl = """Begin.
  {{!s
    Something's going on here...
  }}
End.
"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """Begin.
End.
"""

	## Inline comments should not strip whitespace
tpl = """  12 {{! 34 }}
"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """  12
"""

	## Comment removal should preserve surrounding whitespace.
tpl = """12345 {{! Comment Block! }} 67890"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """12345  67890"""
end


text:
module MachO

using StructIO

# Bring in ObjectFile definitions
using ObjectFile
import ObjectFile: DynamicLink, DynamicLinks, RPath, ObjectHandle, Section, Sections, SectionRef,
                   Segment, Segments, SegmentRef, StrTab, Symbols, SymtabEntry, SymbolRef,
                   getindex, length, iterate, keys, lastindex, eltype, handle, header, path,
                   rpaths, canonical_rpaths, find_library, readmeta, seek, seekstart, skip,
                   iostream, position, read, readuntil, eof, endianness, is64bit, isrelocatable,
                   isexecutable, islibrary, isdynamic, mangle_section_name, mangle_symbol_name,
                   format_string, section_header_offset, section_header_size, section_header_type,
                   segment_header_offset, segment_header_size, segment_header_type, startaddr,
                   symtab_entry_offset, symtab_entry_size, symtab_entry_type, find_libraries,
                   findfirst, deref, section_name, section_size,
                   section_offset, section_address, section_number, segment_name, segment_offset,
                   segment_file_size, segment_memory_size, segment_address, strtab_lookup,
                   symbol_name, symbol_value, isundef, isglobal, islocal, isweak, symbol_number
# Load in imported C #define constants
include("constants.jl")

# Start to bring in concrete types, in the order they're needed
include("MachOHeader.jl")
include("MachOHandle.jl")
include("MachOLoadCmd.jl")
include("MachOSegment.jl")
include("MachOSection.jl")
include("MachODynamicLink.jl")
include("MachOStrTab.jl")
include("MachOSymbol.jl")

# We do not yet support Fat (Universal) MachO binaries, as I have yet to come
# up with a nice abstraction over them that fits in well with COFF/ELF.
include("MachOFat.jl")

end # module MachO


text:
@testset "Basic Monomial Functionalities" begin

    @testset "automatic sorting" begin
        x = semi_elementary_monomial((3,2,1),(1,0,0))
        @test x.sp_term == (0,1,2)
        x = semi_elementary_monomial((2,1,0),(1,0,0))
        @test x.sp_term == (0,1,2)
    end

    @testset "pull out common factors" begin
        x = semi_elementary_monomial((3,4,5),(0,0,0))
        @test x == semi_elementary_monomial((0,1,2),(0,0,3))
    end

    @testset "test <, 0 terms same" begin
        x = semi_elementary_monomial((3,5,5),(0,0,0))
        y = semi_elementary_monomial((3,4,6),(0,0,0))
        @test (x < y) == true
    end

    @testset "test <, 1 term same" begin
        x = semi_elementary_monomial((3,3,6),(0,0,0))
        y = semi_elementary_monomial((3,4,6),(0,0,0))
        @test (x < y) == true
    end

    @testset "test <, all sp terms same" begin
        x = semi_elementary_monomial((3,3,6),(0,0,0))
        y = semi_elementary_monomial((3,3,6),(1,0,0))
        @test (x < y) == true
    end

    @testset "test <, autosimplifying" begin
        x = semi_elementary_monomial((2,2,2),(1,0,0))
        y = semi_elementary_monomial((3,3,3),(0,0,0))
        @test (x < y) == false
    end

    @testset "test is_elementary" begin
        x = semi_elementary_monomial((3,3,3),(0,0,0))
        @test is_elementary(x) == true
        x = semi_elementary_monomial((1,0,0),(0,0,0))
        @test is_elementary(x) == true
        x = semi_elementary_monomial((0,0,1),(0,0,0))
        @test is_elementary(x) == true
    end

    @testset "test dim" begin
        x = semi_elementary_monomial((3,3,3),(0,0,0))
        @test dim(x) == 3
        x = semi_elementary_monomial((3,3,3,3),(0,0,0,0))
        @test dim(x) == 4
    end

    @testset "string conversion" begin
        x = semi_elementary_monomial((3,3,3),(0,0,0))
        @test to_string(x) == "e3^3 "
        x = semi_elementary_monomial((3,3,0),(0,0,0))
        @test to_string(x) == "S(0, 3, 3)"
        x = semi_elementary_monomial((3,3,1),(1,0,0))
        @test to_string(x) == "S(0, 2, 2)e1 e3 "
    end
end

max_stars_repo_name: yushoteke/SymmetricPolynomials

max_issues_repo_name: yushoteke/SymmetricPolynomials

max_forks_repo_name: yushoteke/SymmetricPolynomials


text:
    const MaybePointer = Union{Ptr, UndefInitializer}

    # Definition and constructors:
    """
    ```julia
    MallocArray{T,N} <: AbstractArray{T,N}
    ```
    `N`-dimensional dense heap-allocated array with elements of type `T`.

    Much like `Base.Array`, except (1) backed by memory that is not tracked by
    the Julia garbage collector (is directly allocated with `malloc`) so is
    StaticCompiler-safe, (2) should be `free`d when no longer in use, and
    (3) indexing returns views rather than copies.
    """
    struct MallocArray{T,N} <: DenseArray{T,N}
        pointer::Ptr{T}
        length::Int
        size::NTuple{N, Int}
    end
    """
    ```julia
    MallocMatrix{T} <: AbstractMatrix{T}
    ```
    Two-dimensional dense heap-allocated array with elements of type `T`.
    As `Base.Matrix` is to `Base.Array`, but with `MallocArray`.
    """
    const MallocMatrix{T} = MallocArray{T,2}
    """
    ```julia
    MallocVector{T} <: AbstractVector{T}
    ```
    Two-dimensional dense heap-allocated array with elements of type `T`.
    As `Base.Vector` is to `Base.Array`, but with `MallocArray`.
    """
    const MallocVector{T} = MallocArray{T,1}

    """
    ```julia
    MallocArray{T}(undef, dims)
    MallocArray{T,N}(undef, dims)
    ```
    Construct an uninitialized `N`-dimensional `MallocArray` containing elements
    of type `T`. `N` can either be supplied explicitly, as in `Array{T,N}(undef, dims)`,
    or be determined by the length or number of `dims`. `dims` may be a tuple or
    a series of integer arguments corresponding to the lengths in each dimension.
    If the rank `N` is supplied explicitly, then it must match the length or
    number of `dims`. Here `undef` is the `UndefInitializer`.

    ## Examples
    ```julia
    julia> A = MallocArray{Float64}(undef, 3,3) # implicit N
    3×3 MallocMatrix{Float64}:
     3.10504e231   6.95015e-310   2.12358e-314
     1.73061e-77   6.95015e-310   5.56271e-309
     6.95015e-310  0.0           -1.29074e-231

    julia> free(A)
    0

    julia> A = MallocArray{Float64, 3}(undef, 2,2,2) # explicit N
    2×2×2 MallocArray{Float64, 3}:
    [:, :, 1] =
     3.10504e231  2.0e-323
     2.32036e77   6.94996e-310

    [:, :, 2] =
     6.95322e-310  5.0e-324
     6.95322e-310  5.56271e-309

    julia> free(A)
    0
    ```
    """
    @inline function MallocArray{T,N}(::UndefInitializer, length::Int, dims::Dims{N}) where {T,N}
        @assert Base.allocatedinline(T)
        @assert length == prod(dims)
        p = Ptr{T}(malloc(length*sizeof(T)))
        MallocArray{T,N}(p, length, dims)
    end
    @inline MallocArray{T,N}(x::MaybePointer, dims::Dims{N}) where {T,N} = MallocArray{T,N}(x, prod(dims), dims)
    @inline MallocArray{T}(x::MaybePointer, dims::Dims{N}) where {T,N} = MallocArray{T,N}(x, prod(dims), dims)
    @inline MallocArray{T,N}(x::MaybePointer, dims::Vararg{Int}) where {T,N} = MallocArray{T,N}(x, prod(dims), dims)
    @inline MallocArray{T}(x::MaybePointer, dims::Vararg{Int}) where {T} = MallocArray{T}(x, dims)

    # Indirect constructors
    @inline Base.similar(a::MallocArray{T,N}) where {T,N} = MallocArray{T,N}(undef, size(a))
    @inline Base.similar(a::MallocArray{T}, dims::Dims{N}) where {T,N} = MallocArray{T,N}(undef, dims)
    @inline Base.similar(a::MallocArray{T}, dims::Vararg{Int}) where {T} = MallocArray{T}(undef, dims)
    @inline Base.similar(a::MallocArray, ::Type{T}, dims::Dims{N}) where {T,N} = MallocArray{T,N}(undef, dims)
    @inline Base.similar(a::MallocArray, ::Type{T}, dims::Vararg{Int}) where {T} = MallocArray{T}(undef, dims)
    @inline function Base.copy(a::MallocArray{T,N}) where {T,N}
        new_a = MallocArray{T,N}(undef, size(a))
        copyto!(new_a, a)
        return new_a
    end

    # Destructor:
    @inline free(a::MallocArray) = free(a.pointer)

    # Fundamentals
    @inline Base.unsafe_convert(::Type{Ptr{T}}, a::MallocArray) where {T} = Ptr{T}(a.pointer)
    @inline Base.pointer(a::MallocArray) = a.pointer
    @inline Base.length(a::MallocArray) = a.length
    @inline Base.sizeof(a::MallocArray{T}) where {T} = a.length * sizeof(T)
    @inline Base.size(a::MallocArray) = a.size

    # Some of the AbstractArray interface:
    @inline Base.IndexStyle(::MallocArray) = IndexLinear()
    @inline Base.stride(a::MallocArray, dim::Int) = (dim <= 1) ? 1 : stride(a, dim-1) * size(a, dim-1)
    @inline Base.firstindex(::MallocArray) = 1
    @inline Base.lastindex(a::MallocArray) = a.length
    @inline Base.getindex(a::MallocArray{T}, i::Int) where T = unsafe_load(pointer(a)+(i-1)*sizeof(T))
    @inline Base.getindex(a::MallocArray, ::Colon) = a
    @inline Base.getindex(a::MallocArray{T}, r::UnitRange{<:Integer}) where T = MallocArray(pointer(a)+(first(r)-1)*sizeof(T), length(r), (length(r),))

    @inline Base.getindex(a::MallocArray, r::UnitRange{<:Integer}, inds::Vararg{Int}) = getindex(a, r, inds)
    @inline function Base.getindex(a::MallocArray{T}, r::UnitRange{<:Integer}, inds::Dims{N}) where {T,N}
        i0 = 0
        for i=1:N
            i0 += (inds[i]-1) * stride(a, i+1)
        end
        return MallocArray{T,1}(pointer(a)+i0*sizeof(T), length(r), (length(r),))
    end

    @inline Base.getindex(a::MallocArray, ::Colon, inds::Vararg{Int}) = getindex(a, :, inds)
    @inline function Base.getindex(a::MallocArray{T}, ::Colon, inds::Dims{N}) where {T,N}
        i0 = 0
        for i=1:N
            i0 += (inds[i]-1) * stride(a, i+1)
        end
        return MallocArray{T,1}(pointer(a)+i0*sizeof(T), size(a,1), (size(a,1),))
    end
    @inline Base.getindex(a::MallocArray, ::Colon, ::Colon, inds::Vararg{Int}) = getindex(a, :, :, inds)
    @inline function Base.getindex(a::MallocArray{T}, ::Colon, ::Colon, inds::Dims{N}) where {T,N}
        i0 = 0
        for i=1:N
            i0 += (inds[i]-1) * stride(a, i+2)
        end
        return MallocArray{T,2}(pointer(a)+i0*sizeof(T), size(a,1)*size(a,2), (size(a,1), size(a,2)))
    end
    @inline Base.getindex(a::MallocArray, ::Colon, ::Colon, ::Colon, inds::Vararg{Int}) = getindex(a, :, :, :, inds)
    @inline function Base.getindex(a::MallocArray{T}, ::Colon, ::Colon, ::Colon, inds::Dims{N}) where {T,N}
        i0 = 0
        for i=1:N
            i0 += (inds[i]-1) * stride(a, i+3)
        end
        return MallocArray{T,3}(pointer(a)+i0*sizeof(T), size(a,1)*size(a,2)*size(a,3), (size(a,1), size(a,2), size(a,3)))
    end

    @inline Base.setindex!(a::MallocArray{T}, x::T, i::Int) where T = unsafe_store!(pointer(a)+(i-1)*sizeof(T), x)
    @inline Base.setindex!(a::MallocArray{T}, x, i::Int) where T = unsafe_store!(pointer(a)+(i-1)*sizeof(T), convert(T,x))
    @inline function Base.setindex!(a::MallocArray, x, r::UnitRange{Int})
        is₀ = first(r)-1
        ix₀ = firstindex(x)-1
        @inbounds for i = 1:length(r)
            setindex!(a, x[i+ix₀], i+is₀)
        end
    end
    @inline function Base.setindex!(a::MallocArray{T}, x::Union{T,<:Number}, r::UnitRange{Int}) where T
        xₜ = convert(T,x)
        is₀ = first(r)-1
        @inbounds for i = 1:length(r)
            setindex!(a, xₜ, i+is₀)
        end
    end
    @inline function Base.setindex!(a::MallocArray, x, ::Colon)
        ix₀ = firstindex(x)-1
        @inbounds for i = 1:length(a)
            setindex!(a, x[i+ix₀], i)
        end
    end
    @inline function Base.setindex!(a::MallocArray{T}, x::Union{T,<:Number}, ::Colon) where T
        xₜ = convert(T,x)
        @inbounds for i = 1:length(a)
            setindex!(a, xₜ, i)
        end
    end

    # Other nice functions
    @inline function Base.:(==)(a::MallocArray{A}, b::MallocArray{B}) where {A,B}
        (N = length(a)) == length(b) || return false
        pa, pb = pointer(a), pointer(b)
        for n in 0:N-1
            unsafe_load(pa + n*sizeof(A)) == unsafe_load(pb + n*sizeof(B)) || return false
        end
        return true
    end
    @inline function Base.:(==)(a::MallocArray, b::NTuple{N, <:Number}) where N
        N == length(a) || return false
        for n in 1:N
            a[n] == b[n] || return false
        end
        return true
    end
    @inline function Base.:(==)(a::NTuple{N, <:Number}, b::MallocArray) where N
        N == length(b) || return false
        for n in 1:N
            a[n] == b[n] || return false
        end
        return true
    end

    # Reshaping and Reinterpreting
    @inline function Base.reshape(a::MallocArray{T}, dims::Dims{N})  where {T,N}
        @assert prod(dims) == length(a)
        MallocArray{T,N}(pointer(a), dims)
    end
    @inline Base.reshape(a::MallocArray, dims::Vararg{Int}) = reshape(a, dims)

    @inline function Base.reinterpret(::Type{Tᵣ}, a::MallocArray{Tᵢ,N}) where {N,Tᵣ,Tᵢ}
        @assert Base.allocatedinline(Tᵣ)
        @assert length(a)*sizeof(Tᵢ) % sizeof(Tᵣ) == 0
        @assert size(a,1)*sizeof(Tᵢ) % sizeof(Tᵣ) == 0
        lengthᵣ = length(a)*sizeof(Tᵢ)÷sizeof(Tᵣ)
        sizeᵣ = ntuple(i -> i==1 ? size(a,i)*sizeof(Tᵢ)÷sizeof(Tᵣ) : size(a,i), Val(N))
        pointerᵣ = Ptr{Tᵣ}(pointer(a))
        MallocArray{Tᵣ,N}(pointerᵣ, lengthᵣ, sizeᵣ)
    end

    # Custom printing
    @inline Base.print(a::MallocArray) = printf(a)
    @inline Base.println(a::MallocArray) = (printf(a); newline())
    @inline Base.print(fp::Ptr{FILE}, a::MallocArray) = printf(fp, a)
    @inline Base.println(fp::Ptr{FILE}, a::MallocArray) = (printf(fp, a); newline(fp))


text:
import Base.ht_keyindex

# was thinking maybe if i dont try to hash everytime then it will be faster
# but it's now

function sumby_contig{S,T}(by::Vector{T}, val::Vector{S})
  l = length(by)
  l == length(val) || throw(DimensionMismatch())
  res = Dict{T, S}()

  @inbounds tmp_val = val[1]
  @inbounds last_id = by[1]

  for i in 2:l
    @inbounds id = by[i]
    if id == last_id # keep adding if the value is the same
      @inbounds tmp_val += val[i]
    else
      index = ht_keyindex(res, last_id)
      if index > 0
        @inbounds res.vals[index] += tmp_val
      else
        @inbounds res[last_id] = tmp_val
      end
      @inbounds tmp_val = val[i]
      last_id = id
    end
  end

  @inbounds id = by[l]
  index = ht_keyindex(res, id)
  if index > 0
    @inbounds res.vals[index] += tmp_val
  else
    @inbounds res[id] = val[l]
  end

  return Dict(k => res[k] for k in keys(res))
end

two_groups = rand(1:2,N)
@time sumby_contig(two_groups, v1); # 36 seconds so only marginally better
@time sumby_contig(two_groups, v1) #
@time sumby(two_groups,v1)
@time sumby_contig(id6, v1);


text:
new_var(Q, Nel) = zeros(Q,Q,Nel)

function deriv_x!(elems::Array{Element2d,1}, f::Array{Float64,3}, dx::Array{Float64,3})

    Q = size(f,1)
    ws = zeros(Q,Q)
    Nel = length(elems)
    for e = 1:Nel
        deriv_x!(elems[e], sub(f,:,:,e), sub(dx,:,:,e), ws)
    end
end

function deriv_y!(elems::Array{Element2d,1}, f::Array{Float64,3}, dy::Array{Float64,3})

    Q = size(f,1)
    ws = zeros(Q,Q)
    Nel = length(elems)
    for e = 1:Nel
        deriv_y!(elems[e], sub(f,:,:,e), sub(dy,:,:,e), ws)
    end
end


text:
#! Transformations of the meshes. Uses:
#! LinearAlgebra.UniformScaling for scaling
#! Rotation.jl for rotation
#! CoordinateTransformations.jl for translation
# Voir https://stackoverflow.com/questions/10546320/remove-rotation-from-a-4x4-homogeneous-transformation-matrix
# pour extraire la rotation et la translation depuis la matrice 4*4, puis transformer
# cette matrice 4*4 en CoordinateTransformations.jl transformations et après on ne gèrera
# que ces transformations là.

"""
    transform_mesh!(node::Node, transformation)

Add a new `CoordinateTransformations.jl` transformation to the node geometry
`transformation` field. The transformation is composed with the previous
transformation if any.

`transformation` must be a `CoordinateTransformations.jl` transformation.

It is also possible to invert a transformation using `inv` from
`CoordinateTransformations.jl`.

# Examples

```julia
using PlantGeom, MultiScaleTreeGraph, GLMakie, Rotations, CoordinateTransformations

file = joinpath(dirname(dirname(pathof(PlantGeom))), "test", "files", "simple_plant.opf")
opf = read_opf(file)

# Visualize the mesh as is:
viz(opf)

# Copy the OPF, and translate the whole plant by 15 in the y direction (this is in cm, the mesh comes from XPlo):
opf2 = deepcopy(opf)
transform!(opf2, x -> transform_mesh!(x, Translation(0, 15, 0)))
viz!(opf2) # Visualize it again in the same figure

# Same but rotate the whole plant around the X axis:
opf3 = deepcopy(opf)
transform!(opf3, x -> transform_mesh!(x, LinearMap(RotX(0.3))))
# NB: we use Rotations.jl's RotX here. Input in radian, use rad2deg and deg2rad if needed.
viz!(opf3)

# Same but rotate only the second leaf around the Z axis:
opf4 = deepcopy(opf)
# Build the meshes from the reference meshes (need it because we want the coordinates of the parent):
transform!(opf4, refmesh_to_mesh!)

# Get the second leaf in the OPF:
leaf_node = get_node(opf4, 8)

# Get the parent node (internode) Z coordinates:
parent_zmax = zmax(leaf_node.parent)

# Define a rotation of the mesh around the Z axis defined by the parent node max Z:
transformation = recenter(LinearMap(RotZ(1.0)), Point3(0.0, 0.0, parent_zmax))

# Update the transformation matrix of the leaf and its mesh:
transform_mesh!(leaf_node, transformation)

# Plot the result:
viz(opf)
viz!(opf4)
```
"""
function transform_mesh!(node::MultiScaleTreeGraph.Node, transformation)
    if node[:geometry] !== nothing
        node[:geometry].transformation = transformation ∘ node[:geometry].transformation
        # If the node has a computed mesh, update it:
        if node[:geometry].mesh !== nothing
            refmesh_to_mesh!(node)
        end
    end
end


text:
module Tests

using JuMP
import GLPK
const lp_solver = optimizer_with_attributes(GLPK.Optimizer, "presolve" => GLPK.GLP_ON)
import Polyhedra
const lib = Polyhedra.DefaultLibrary{Float64}(lp_solver)

include("utilities.jl")

include("square.jl")
include("invariant.jl")
include("controlled_invariant.jl")
include("switched_system.jl")

end


text:
#=
julia> ilog2(15),ilog2(-15),ilog2(16),ilog2(-16),ilog2(17),ilog2(-17)
(3,3,4,4,4,4)

julia> ilog2(15.),ilog2(-15.),ilog2(16.),ilog2(-16.),ilog2(17.),ilog2(-17.)
(3,3,4,4,4,4)

julia> ilog2(1/15.),ilog2(-1/15.),ilog2(1/16.),ilog2(-1/16.),ilog2(1/17.),ilog2(-1/17.)
(3,3,3,3,4,4)
=#

@inline ilog2_core{T<:Integer}(x::T) = trailing_zeros(prevpow2(abs(x)))
@inline ilog2_core{T<:AbstractFloat}(x::T) = abs(floor(Int, log2(abs(x))))

ilog2{T<:Integer}(x::T) =
    ifelse( x!=zero(T), ilog2_core(x), zero(T))

ilog2{T<:AbstractFloat}(x::T) =
    ifelse( x!=zero(T), ilog2_core(x)-(abs(x)<one(T)), zero(T))


text:
using Test
using MIPVerify
using MIPVerify: BatchRunParameters, UnrestrictedPerturbationFamily, mkpath_if_not_present, create_summary_file_if_not_present, verify_target_indices
using MIPVerify: run_on_sample_for_untargeted_attack, run_on_sample_for_targeted_attack
using DataFrames

@timed_testset "unit.jl" begin
    mnist = read_datasets("MNIST")

    @testset "BatchRunParameters" begin
        brp = BatchRunParameters(
            Sequential([], "name"),
            UnrestrictedPerturbationFamily(),
            1,
            0
        )
        @testset "Base.show" begin
            io = IOBuffer()
            Base.show(io, brp)
            @test String(take!(io)) == "name__unrestricted__1__0"
        end
    end

    @testset "mkpath_if_not_present" begin
        mktempdir() do dir
            path = joinpath(dir, "1")
            mkpath_if_not_present(path)
            @test ispath(path)
            mkpath_if_not_present(path)
            @test ispath(path)
        end
    end

    @testset "create_summary_file_if_not_present" begin
        mktempdir() do dir
            file_path = joinpath(dir, "summary.csv")
            create_summary_file_if_not_present(file_path)
            @test isfile(file_path)
        end
    end

    @testset "verify_target_indices" begin
        @test_throws AssertionError verify_target_indices([0], mnist.test)
        @test_throws AssertionError verify_target_indices([10001], mnist.test)
    end

    @testset "run_on_sample_for_untargeted_attack" begin
        sample_numbers = [101, 102, 103, 104]
        dt = DataFrame(
            SampleNumber = [101, 102, 102, 103, 103],
            SolveStatus = ["Infeasible", "Optimal", "UserObjLimit", "UserObjLimit", "UserLimit"],
            ObjectiveValue = [NaN, 0.8, 0.9, 0.99, NaN])
        expected_results = Dict(
            MIPVerify.never => [false, false, false, true],
            MIPVerify.always => [true, true, true, true],
            MIPVerify.resolve_ambiguous_cases => [false, false, true, true],
            MIPVerify.refine_insecure_cases => [false, true, false, true]
        )
        for solve_rerun_option in keys(expected_results)
            @testset "Solve rerun option: $solve_rerun_option" begin
                actual = map(x -> run_on_sample_for_untargeted_attack(x, dt, solve_rerun_option), sample_numbers)
                expected = expected_results[solve_rerun_option]
                @test actual==expected
            end
        end
    end

    @testset "run_on_sample_for_targeted_attack" begin
        sample_numbers = [(101, 1), (102, 1), (103, 1), (104, 1), (101, 2)]
        dt = DataFrame(
            SampleNumber = [101, 102, 102, 103, 103, 101],
            TargetIndexes = map(x -> "[$x]", [1, 1, 1, 1, 1, 2]),
            SolveStatus = ["Infeasible", "Optimal", "UserObjLimit", "UserObjLimit", "UserLimit", "Optimal"],
            ObjectiveValue = [NaN, 0.8, 0.9, 0.99, NaN, 0.1])
        expected_results = Dict(
            MIPVerify.never => [false, false, false, true, false],
            MIPVerify.always => [true, true, true, true, true],
            MIPVerify.resolve_ambiguous_cases => [false, false, true, true, false],
            MIPVerify.refine_insecure_cases => [false, true, false, true, false],
            MIPVerify.retarget_infeasible_cases => [true, false, false, true, false]
        )
        for solve_rerun_option in keys(expected_results)
            @testset "Solve rerun option: $solve_rerun_option" begin
                actual = map(x -> run_on_sample_for_targeted_attack(x..., dt, solve_rerun_option), sample_numbers)
                expected = expected_results[solve_rerun_option]
                @test actual==expected
            end
        end
    end
end


text:
#   Generating exports
export XML_SetElementDeclHandler
export XML_SetAttlistDeclHandler
export XML_SetXmlDeclHandler
export XML_ParserCreate
export XML_ParserCreateNS
export XML_ParserCreate_MM
export XML_ParserReset
export XML_SetEntityDeclHandler
export XML_SetElementHandler
export XML_SetStartElementHandler
export XML_SetEndElementHandler
export XML_SetCharacterDataHandler
export XML_SetProcessingInstructionHandler
export XML_SetCommentHandler
export XML_SetCdataSectionHandler
export XML_SetStartCdataSectionHandler
export XML_SetEndCdataSectionHandler
export XML_SetDefaultHandler
export XML_SetDefaultHandlerExpand
export XML_SetDoctypeDeclHandler
export XML_SetStartDoctypeDeclHandler
export XML_SetEndDoctypeDeclHandler
export XML_SetUnparsedEntityDeclHandler
export XML_SetNotationDeclHandler
export XML_SetNamespaceDeclHandler
export XML_SetStartNamespaceDeclHandler
export XML_SetEndNamespaceDeclHandler
export XML_SetNotStandaloneHandler
export XML_SetExternalEntityRefHandler
export XML_SetExternalEntityRefHandlerArg
export XML_SetSkippedEntityHandler
export XML_SetUnknownEncodingHandler
export XML_DefaultCurrent
export XML_SetReturnNSTriplet
export XML_SetUserData
export XML_SetEncoding
export XML_UseParserAsHandlerArg
export XML_UseForeignDTD
export XML_SetBase
export XML_GetBase
export XML_GetSpecifiedAttributeCount
export XML_GetIdAttributeIndex
export XML_Parse
export XML_GetBuffer
export XML_ParseBuffer
export XML_StopParser
export XML_ResumeParser
export XML_GetParsingStatus
export XML_ExternalEntityParserCreate
export XML_SetParamEntityParsing
export XML_SetHashSalt
export XML_GetErrorCode
export XML_GetCurrentLineNumber
export XML_GetCurrentColumnNumber
export XML_GetCurrentByteIndex
export XML_GetCurrentByteCount
export XML_GetInputContext
export XML_FreeContentModel
export XML_MemMalloc
export XML_MemRealloc
export XML_MemFree
export XML_ParserFree
export XML_ErrorString
export XML_ExpatVersion
export XML_ExpatVersionInfo
export XML_GetFeatureList
export XML_Char
export XML_LChar
export XML_Index
export XML_Size
export XML_Parser
export XML_Bool
export XML_STATUS_ERROR
export XML_STATUS_OK
export XML_STATUS_SUSPENDED
export XML_ERROR_NONE
export XML_ERROR_NO_MEMORY
export XML_ERROR_SYNTAX
export XML_ERROR_NO_ELEMENTS
export XML_ERROR_INVALID_TOKEN
export XML_ERROR_UNCLOSED_TOKEN
export XML_ERROR_PARTIAL_CHAR
export XML_ERROR_TAG_MISMATCH
export XML_ERROR_DUPLICATE_ATTRIBUTE
export XML_ERROR_JUNK_AFTER_DOC_ELEMENT
export XML_ERROR_PARAM_ENTITY_REF
export XML_ERROR_UNDEFINED_ENTITY
export XML_ERROR_RECURSIVE_ENTITY_REF
export XML_ERROR_ASYNC_ENTITY
export XML_ERROR_BAD_CHAR_REF
export XML_ERROR_BINARY_ENTITY_REF
export XML_ERROR_ATTRIBUTE_EXTERNAL_ENTITY_REF
export XML_ERROR_MISPLACED_XML_PI
export XML_ERROR_UNKNOWN_ENCODING
export XML_ERROR_INCORRECT_ENCODING
export XML_ERROR_UNCLOSED_CDATA_SECTION
export XML_ERROR_EXTERNAL_ENTITY_HANDLING
export XML_ERROR_NOT_STANDALONE
export XML_ERROR_UNEXPECTED_STATE
export XML_ERROR_ENTITY_DECLARED_IN_PE
export XML_ERROR_FEATURE_REQUIRES_XML_DTD
export XML_ERROR_CANT_CHANGE_FEATURE_ONCE_PARSING
export XML_ERROR_UNBOUND_PREFIX
export XML_ERROR_UNDECLARING_PREFIX
export XML_ERROR_INCOMPLETE_PE
export XML_ERROR_XML_DECL
export XML_ERROR_TEXT_DECL
export XML_ERROR_PUBLICID
export XML_ERROR_SUSPENDED
export XML_ERROR_NOT_SUSPENDED
export XML_ERROR_ABORTED
export XML_ERROR_FINISHED
export XML_ERROR_SUSPEND_PE
export XML_ERROR_RESERVED_PREFIX_XML
export XML_ERROR_RESERVED_PREFIX_XMLNS
export XML_ERROR_RESERVED_NAMESPACE_URI
export XML_CTYPE_EMPTY
export XML_CTYPE_ANY
export XML_CTYPE_MIXED
export XML_CTYPE_NAME
export XML_CTYPE_CHOICE
export XML_CTYPE_SEQ
export XML_CQUANT_NONE
export XML_CQUANT_OPT
export XML_CQUANT_REP
export XML_CQUANT_PLUS
export XML_Content
export XML_ElementDeclHandler
export XML_AttlistDeclHandler
export XML_XmlDeclHandler
export XML_StartElementHandler
export XML_EndElementHandler
export XML_CharacterDataHandler
export XML_ProcessingInstructionHandler
export XML_CommentHandler
export XML_StartCdataSectionHandler
export XML_EndCdataSectionHandler
export XML_DefaultHandler
export XML_StartDoctypeDeclHandler
export XML_EndDoctypeDeclHandler
export XML_EntityDeclHandler
export XML_UnparsedEntityDeclHandler
export XML_NotationDeclHandler
export XML_StartNamespaceDeclHandler
export XML_EndNamespaceDeclHandler
export XML_NotStandaloneHandler
export XML_ExternalEntityRefHandler
export XML_SkippedEntityHandler
export XML_UnknownEncodingHandler
export XML_INITIALIZED
export XML_PARSING
export XML_FINISHED
export XML_SUSPENDED
export XML_PARAM_ENTITY_PARSING_NEVER
export XML_PARAM_ENTITY_PARSING_UNLESS_STANDALONE
export XML_PARAM_ENTITY_PARSING_ALWAYS
export XML_FEATURE_END
export XML_FEATURE_UNICODE
export XML_FEATURE_UNICODE_WCHAR_T
export XML_FEATURE_DTD
export XML_FEATURE_CONTEXT_BYTES
export XML_FEATURE_MIN_SIZE
export XML_FEATURE_SIZEOF_XML_CHAR
export XML_FEATURE_SIZEOF_XML_LCHAR
export XML_FEATURE_NS
export XML_FEATURE_LARGE_SIZE
export XML_FEATURE_ATTR_INFO
export XML_TRUE
export XML_TRUE


text:
using FileIO
using Plots
include("ex3_base.jl")
import Example3
using LaTeXStrings
using CSV, DataFrames
pgfplots()

solutions = "solutions/"

function contour_plot(results; use_log=false)
    levels = 10
    if !use_log
        z = results["grid"]["cost"]
    else
        z = log10.(results["grid"]["cost"])
        if results["options"]["pem"][:sim_len] == 1024
            levels = 0:0.25:7
        elseif results["options"]["pem"][:sim_len] == 32
            levels = -4:0.25:2
        end
    end
    p_contour = contourf(results["grid"]["gl"], results["grid"]["ka"], z,
                        color=:viridis, levels=levels, xguide=L"g/l", yguide=L"k_a")
    scatter!(p_contour, [Example3.GL], [Example3.KA], color="red", legend=false, ms=8)
    if length(results["solutions"]) != 0
        g_l_list = [min(max(θ[1], results["options"]["grid_cost_funtion"][:gl][1]), results["options"]["grid_cost_funtion"][:gl][2]) for θ in results["solutions"]["theta"]]
        ka_list = [min(max(θ[2], results["options"]["grid_cost_funtion"][:ka][1]), results["options"]["grid_cost_funtion"][:ka][2]) for θ in results["solutions"]["theta"]]
        scatter!(p_contour, g_l_list', ka_list', color="blue", ms=8)
        plot!(p_contour, legend=false)
    end
    return p_contour
end

for root_dir in take!(walkdir(solutions))[2]
    # Plot contour
    first = true
    for dir_aux in take!(walkdir(joinpath(solutions, root_dir)))[2]
        dir = joinpath(solutions, root_dir, dir_aux)
        results = load(joinpath(dir, "results.jld2"))["results"]
        if first
            # Plot input
            p_input = plot(results["data"]["u"], legend=false, xguide=L"k", yguide=L"u", grid=false)
            savefig(p_input, joinpath(solutions, root_dir, "input_signal.pdf"))
            savefig(p_input, joinpath(solutions, root_dir, "input_signal.tex"))
            # Plot output
            p_output = plot(results["data"]["y"]/π, legend=false, xguide=L"k", yguide=L"$\frac{y}{\pi}$",
                           grid=false)
            savefig(p_output, joinpath(solutions, root_dir, "output_signal.pdf"))
            savefig(p_output, joinpath(solutions, root_dir, "output_signal.tex"))
        end
        first=false
        # Plot countour
        pyplot()
        p_contour = contour_plot(results, use_log=false)
        savefig(p_contour, joinpath(dir, "countour.pdf"))
        savefig(p_contour, joinpath(dir, "countour.png"))
        # Plot log contour
        p_contourl = contour_plot(results, use_log=true)
        savefig(p_contourl, joinpath(dir, "countour_log.pdf"))
        savefig(p_contourl, joinpath(dir, "countour_log.png"))
        pgfplots()
        p_contourl = contour_plot(results, use_log=true)
        savefig(p_contourl, joinpath(dir, "countour_log.tex"))
        # Save solutions
        CSV.write(joinpath(dir, "solutions.csv"),
                  DataFrame(hcat(results["solutions"]["theta"]...)'),
                  header=false)
       # Save grid
       ka = repeat(Array(results["grid"]["ka"])[1:2:end], outer=50)
       gl = repeat(Array(results["grid"]["gl"])[1:2:end], inner=50)
       cost = reshape(results["grid"]["cost"][1:2:end, 1:2:end]', (2500,))
       CSV.write(joinpath(dir, "grid.csv"),
                 DataFrame(hcat(gl, ka, cost)),
                 header=false)

    end
end


text:
module Kernels

export Standard, SE, Perioidc

using LinearAlgebra

abstract type Kernel end

nax = [CartesianIndex()]

struct Standard <: Kernel
    hyperparameters::Dict{String, Float64}
    data
    function Standard(hyperparameters)
        new(hyperparameters, nothing)
    end
end
"""
Xm : DxM
Xn : DxN
"""
function (ker::Standard)(Xm::Array{Float64, 2}, Xn::Array{Float64, 2};
                            orig=false,
                            noise=false,
                            hyperparameters=nothing,
                            gradient_only=false,
                            hessian_only=false,
                            potential_only=false)

    hyperparameters = hyperparameters == nothing ? ker.hyperparameters :
                                                    hyperparameters

    println("Size of Xm, Xn: ", size(Xm), size(Xn))
    ll = get(hyperparameters, "ll", 1.)
    σ_f = get(hyperparameters, "sigma_f", 1.)
    # Xn = Xn == nothing ? copy(Xm) : Xn
    D, M = size(Xm) # Xm : DxM
    N = size(Xn)[2] # Xn : DxN
    println(size(Xm), size(Xn))
    Xmn = Xm[:, :, nax] .- Xn[:, nax, :]               # DxMxN
    dists = dropdims(sum(Xmn.^2, dims=1), dims=1)      # MxN
    K = σ_f .* exp.(-dists / 2ll)      # M x N
    println(size(K))
    # for m=1:M
    #     for n=1:N
    #         ll = l[m] * l[n]
    #         K[m, n] = sigma_f * exp(-(Xm[m] - Xn[n])^2 / 2ll)
    #     end
    # end
    if orig
        noise_f = get(hyperparameters, "noise_f", 0)
        return @. K + noise_f * I(M)
    end

    # Derivative coefficient D x M x N
    dc_gd = -Xmn / ll
    # DxMxN x 1xMxN -> DxMxN -> MDxN
    Kgd = dc_gd .* K[nax, :, :]
    Kgd = vcat([Kgd[i, :, :] for i in 1:size(Kgd, 1)]...)
    println(size(Kgd))
    if potential_only
        return vcat(K, Kgd)                  # M(1+D) x N
    end
    # DxMxN * 1xMxN -> MxND
    Kdg = -dc_gd .* K[nax, :, :]
    Kdg = hcat([Kdg[i, :, :] for i in 1:size(Kdg, 1)]...)
    # DxMxN -> MxDxN
    Xnm = permutedims(Xmn, [2, 1, 3])
    # DxMx1xN  * 1xMxDxN  -> D x M x D x N
    dc_dd_glob = -Xmn[:, :, nax, :] .* Xnm[nax, :, :, :] / ll^2
    # dc_dd_glob = np.einsum('inm, jnm -> injm', dc_gd, -dc_gd)
    # ∂_mn exp(Xn - Xm)^2
    dc_dd_diag = I(D)[:, nax, :, nax] / ll
    # DxMxDxN - DxMxDxN
    Kdd = @. (dc_dd_glob + dc_dd_diag) * K[nax, :, nax, :]
    # DM x DN
    Kdd = reshape(Kdd, D * M, D * N)
    if gradient_only
        # (D+1)M x DN
        return vcat(Kdg, Kdd)
    end
    if hessian_only
        # Delta _ dd
        dnm = [1:D;]
        # MxDxN
        dc_dg = -Xnm / ll
        # DxMxN * MxDxN -> MxDxDxN
        dc_hg_glob = -Xmn[:, :, nax, :] .* Xnm[nax, :, :, :] / ll^2
        # dc_hg_glob = np.einsum('inm, jnm -> nijm', -dc_gd, -dc_gd)
        # DxNxM -> NxDxM -> NxDDxM
        dc_hg_diag = zeros(M, D, D, N)
        # ∂_mm K(Xm,Xn) MxDxDxN
        # dc_hg_diag[:, dnm, dnm, :] = -1 / ll
        for i in 1:D
            dc_hg_diag[:, i, i, :] .= -1 / ll
        end
        # MxDxDxN + MxDxDxN -> MxDxDxN
        Khg = (dc_hg_glob + dc_hg_diag)
        # Bacground term: ∂_mmn K(Xm,Xn) DxMx1x1xN * 1xMxDxDxN -> DxMxDxDxN
        dc_hd_back = dc_gd[:, :, nax, nax, :] * Khg[nax, :, :, :, :]
        # Diagonal term : DxMxDxDxN * 1xMxDxDxN -> DxMxDxDxN
        dc_hd_diag = zeros(D, M, D, D, N)
        # Global term :
        dc_hd_glob = zeros(D, M, D, D, N)
        # dc_hd_glob[dnm, :, dnm, :, :] += Xmn[nax, :, :, :] / ll / ll
        # dc_hd_glob[dnm, :, :, dnm, :] += Xmn[nax, :, :, :] / ll / ll
        for i in 1:D
            dc_hd_glob[i, :, i, :, :] += Xnm ./ ll / ll
            dc_hd_glob[i, :, :, i, :] += Xnm ./ ll / ll
        end

        # print(dc_hd_glob[0, 0, 0, :, 0])
        # print(dc_hd_glob.reshape(2, N, 2))
        Khd = (dc_hd_glob + dc_hd_diag) + dc_hd_back
        # MxDxDxN x Mx1x1xN -> MxDxDxN
        Khg *= K[:, nax, nax, :]
        # DxMxDxDxN * 1xMx1x1xN
        Khd *= K[nax, :, nax, nax, :]
        # MxDDxN -> M x DDN
        Khg = reshape(Khg, M, D * D * N)
        # DxMxDDxN -> DM x DDN
        Khd = reshape(Khd, D * M, D * D * N)
        # print(Khd.shape)
        return vcat(Khg, Khd)  # (D+1)M x DDN
    end
    println("Sizes K, Kgd, Kdg, Kdd", size(K), size(Kgd), size(Kdg), size(Kdd))
    Kext = [K Kdg;
            Kgd Kdd]  # (D+1)M x (D+1)N
    if noise
        noise_f = get(hyperparameters, "sigma_n^e", 0.)
        noise_df = get(hyperparameters, "sigma_n^f", 0.)
        noise = [noise_f .* ones(M)..., noise_df .* ones(D * M)...]
        return Kext .+ noise .* I((D + 1) * M)
    end
    return Kext
end

struct SE
    hyperparameters::Dict{String, Float64}
    function SE(hyperparameters)
        new(hyperparameters)
    end
end
"""
Xm : DxM
Xn : DxN
ll : D
"""
function (ker::SE)(Xm::Array{Float64, 2}, Xn::Array{Float64, 2};
                            orig=false,
                            noise=false,
                            hyperparameters=nothing,
                            gradient_only=false,
                            hessian_only=false,
                            potential_only=false)

    hyperparameters = hyperparameters == nothing ? ker.hyperparameters :
                                                    hyperparameters

    nax = [CartesianIndex()]
    # Xn = Xn == nothing ? copy(Xm) : Xn
    D, M = size(Xm) # Xm : DxM
    N = size(Xn)[2] # Xn : DxN
    ll = [get(hyperparameters, "l$i", 1) for i in 1:D]
    σ_f = get(hyperparameters, "sigma_f", 1.)
    Xmn = (Xm[:, :, nax] .- Xn[:, nax, :]) ./ ll        # DxMxN
    dists = dropdims(sum(Xmn.^2, dims=1), dims=1)      # MxN
    K = σ_f .* exp.(-0.5 * dists)      # M x N
    # for m=1:M
    #     for n=1:N
    #         ll = l[m] * l[n]
    #         K[m, n] = sigma_f * exp(-(Xm[m] - Xn[n])^2 / 2ll)
    #     end
    # end
    if orig
        noise_f = get(hyperparameters, "noise_f", 0)
        return K + (noise_f * I(M))
    end

    # Derivative coefficient D x M x N
    dc_gd = -Xmn
    # DxMxN x 1xMxN -> DxMxN -> DMxN
    Kgd = dc_gd .* K[nax, :, :]
    Kgd = vcat([Kgd[i, :, :] for i in 1:size(Kgd, 1)]...)
    if potential_only
        return vcat(K, Kgd)                  # (D+1)M x N
    end
    # DxMxN * 1xMxN -> MxDN
    Kdg = -dc_gd .* K[nax, :, :]
    Kdg = hcat([Kdg[i, :, :] for i in 1:size(Kdg, 1)]...)
    # DxMxN -> MxDxN
    Xnm = permutedims(Xmn, [2, 1, 3])
    # DxMx1xN  * 1xMxDxN  -> D x M x D x N
    dc_dd_glob = -Xmn[:, :, nax, :] .* Xnm[nax, :, :, :]
    # dc_dd_glob = np.einsum('inm, jnm -> injm', dc_gd, -dc_gd)
    # ∂_mn exp(Xn - Xm)^2
    dc_dd_diag = I(D)[:, nax, :, nax] ./ ll
    # DxMxDxN - DxMxDxN
    Kdd = @. (dc_dd_glob + dc_dd_diag) * K[nax, :, nax, :]
    # DM x DN
    Kdd = reshape(permutedims(Kdd, [2, 1, 4, 3]), D * M, D * N)
    if gradient_only
        # (D+1)M x DN
        return vcat(Kdg, Kdd)
    end
    if hessian_only
        # Delta _ dd
        dnm = [1:D;]
        # MxDxN
        dc_dg = -Xnm
        # MxDx1xN * Mx1xDxN -> MxDxDxN
        dc_hg_glob = -Xnm[:, :, nax, :] .* Xnm[:, nax, :, :]
        # dc_hg_glob = np.einsum('inm, jnm -> nijm', -dc_gd, -dc_gd)
        # DxMxN -> MxDxN -> MxDDxN
        dc_hg_diag = zeros(M, D, D, N)
        # ∂_mm K(Xm,Xn) MxDxDxN
        # dc_hg_diag[:, dnm, dnm, :] = -1 ./ ll
        for i in 1:D
            dc_hg_diag[:, i, i, :] .= -1 / ll[i]
        end
        # MxDxDxN + MxDxDxN -> MxDxDxN
        Khg = dc_hg_glob + dc_hg_diag
        # Bacground term: ∂_mmn K(Xm,Xn) DxMx1x1xN * 1xMxDxDxN -> DxMxDxDxN
        dc_hd_back = dc_gd[:, :, nax, nax, :] .* Khg[nax, :, :, :, :]

        # Diagonal term : DxMxDxDxN * 1xMxDxDxN -> DxMxDxDxN
        dc_hd_diag = zeros(D, M, D, D, N)
        # Global term :
        dc_hd_glob = zeros(D, M, D, D, N)

        # dc_hd_glob[dnm, :, dnm, :, :] += Xmn[nax, :, :, :] ./ ll
        # dc_hd_glob[dnm, :, :, dnm, :] += Xmn[nax, :, :, :] ./ ll
        for i in 1:D
            dc_hd_glob[i, :, i, :, :] += Xnm ./ ll[i]
            dc_hd_glob[i, :, :, i, :] += Xnm ./ ll[i]
        end
        # print(dc_hd_glob[0, 0, 0, :, 0])
        # print(dc_hd_glob.reshape(2, N, 2))
        Khd = (dc_hd_glob + dc_hd_diag) + dc_hd_back
        # MxDxDxN x Mx1x1xN -> MxDxDxN
        Khg .*= K[:, nax, nax, :]
        # DxMxDxDxN * 1xMx1x1xN
        Khd .*= K[nax, :, nax, nax, :]
        # MxDDxN -> MxNxDD -> M x NDD
        Khg = reshape(permutedims(Khg, [1, 4, 3, 2]), M, D * D * N)
        # DxMxDDxN -> MxDxNxDD -> MD x NDD
        Khd = reshape(permutedims(Khd, [2, 1, 5, 4, 3]), D * M, D * D * N)
        # print(Khd.shape)
        return vcat(Khg, Khd)  # (D+1)M x DDN
    end
    Kext = [K Kdg;
            Kgd Kdd]  # (D+1)M x (D+1)N
    if noise
        noise_f = get(hyperparameters, "sigma_n^e", 0.)
        noise_df = get(hyperparameters, "sigma_n^f", 0.)
        noise = [noise_f .* ones(M)..., noise_df .* ones(D * M)...]
        return Kext .+ noise .* I((D + 1) * M)
    end
    return Kext
end

"""
Gaussian only with limited potential
"""
mutable struct GaussianLite <: Kernel
    hyperparameters::Any
    data::Any
    function GaussianLite(;hyperparameters=nothing, data=nothing)
        new(hyperparameters, data)
    end
end

function (ker::GaussianLite)(Xm::Array{Float64, 2}, Xn::Array{Float64, 2}, noise::Bool)
    @assert noise
    D, M = size(Xm)
    noise_f =  ker.hyperparameters[end]
    K = ker(Xm, Xn)
    return K + (noise_f * I(M))
end

function (ker::GaussianLite)(Xm::Array{Float64, 2}, Xn::Array{Float64, 2};
                             gradients=false, hessian=false)

    D, M = size(Xm) # Xm : DxM
    A = D > 3 ? D ÷ 3 : 1
    N = size(Xn)[2] # Xn : DxN

    ll = ker.hyperparameters[1:end-2]
    σ_f = ker.hyperparameters[end-1]

    Xmn = (Xm[:, :, nax] .- Xn[:, nax, :]) ./ ll       # DxMxN
    dists = dropdims(sum(Xmn.^2, dims=1), dims=1)      # MxN
    K = σ_f .* exp.(-0.5 * dists)                      # M x N

    if gradients || hessian
        # DxMxN * 1xMxN -> MxDN
        # Derivative coefficient D x M x N
        dc_gd = -Xmn ./ ll
        Kdg = -dc_gd .* K[nax, :, :]
        # Kdg = hcat([Kdg[i, :, :] for i in 1:size(Kdg, 1)]...)
        # DxMxN -> MxDxN
        Kdg = permutedims(Kdg, [2, 1, 3])
        # Xnm = permutedims(Xmn, [2, 1, 3])
        # # DxMx1xN  * 1xMxDxN  -> D x M x D x N
        # dc_dd_glob = -Xmn[:, :, nax, :] .* Xnm[nax, :, :, :]
        # # dc_dd_glob = np.einsum('inm, jnm -> injm', dc_gd, -dc_gd)
        # # ∂_mn exp(Xn - Xm)^2
        # dc_dd_diag = I(D)[:, nax, :, nax] ./ ll
        # # DxMxDxN - DxMxDxN
        # Kdd = @. (dc_dd_glob + dc_dd_diag) * K[nax, :, nax, :]
        # # DM x DN
        # Kdd = reshape(permutedims(Kdd, [2, 1, 4, 3]), D * M, D * N)
        if gradients
            return reshape(Kdg, M, D*N)
        end
    end

    if hessian
        # Delta _ dd
        # dnm = [1:D;]
        # MxDxN
        dc_gd = permutedims(-Xmn ./ ll, [2, 1, 3])
        # dc_gd = -Xnm
        # Loop A :=> Mx3x1xN * Mx1x3xN -> Mx3x3AxN
        dc_hg_glob = zeros(Float64, M, 3, D, N)
        for a=1:A
            i, j = 3a-2, 3a
            dc_hg_glob[:, :, i:j, :] = dc_gd[:, i:j, nax, :] .* dc_gd[:, nax, i:j, :]
        end
        # dc_hg_glob = np.einsum('inm, jnm -> nijm', -dc_gd, -dc_gd)
        # DxMxN -> MxDxN -> MxDDxN
        # ∂_mm K(Xm,Xn) MxDxDxN
        # dc_hg_diag[:, dnm, dnm, :] = -1 ./ ll
        # graph = []
        dc_hg_diag = zeros(M, 3, D, N)
        for a in 1:A
            # nidx = graph[a, :]
            nidx = [1, 2, 3]
            xyz = [3a-2, 3a-1, 3a]
            dc_hg_diag[:, nidx, xyz, :] .= -1 / ll[xyz]
        end
        # Mx3xDxN + Mx3xDxN -> Mx3xDxN
        Khg = dc_hg_glob + dc_hg_diag
        Khg = Khg .* K[:, nax, nax, :]
        return reshape(Khg, M, 3D*N)
    end

    return K
end

"""
Gaussian only with limited potential
"""
mutable struct GaussianLite2 <: Kernel
    hyperparameters::Any
    data::Any
    function GaussianLite2(;hyperparameters=nothing, data=nothing)
        new(hyperparameters, data)
    end
end

function (ker::GaussianLite2)(Xm::Array{Float64, 2}, Xn::Array{Float64, 2}, noise::Bool)
    @assert noise
    D, M = size(Xm)
    M′ = M < 10 ? M : 10
    noise_f =  ker.hyperparameters[end]
    noise_e =  ker.hyperparameters[end-1]
    noise = zeros(M + D*M′, M + D*M′)
    for d=1:M
        noise[d, d] = noise_e
    end
    for d=M+1:M+D*M′
        noise[d, d] = noise_f
    end
    Kext = ker(Xm, Xn)
    return Kext + noise
end

function (ker::GaussianLite2)(Xm::Array{Float64, 2}, Xn::Array{Float64, 2};
                              potentials=false, gradients=false, hessian=false,
                              covariance=false)

    D, M = size(Xm) # Xm : DxM
    A = D > 3 ? D ÷ 3 : 1
    N = size(Xn)[2] # Xn : DxN

    ll = ker.hyperparameters[1:end-2]
    σ_f = ker.hyperparameters[end-1]

    Xmn = (Xm[:, :, nax] .- Xn[:, nax, :]) ./ ll       # DxMxN
    dists = dropdims(sum(Xmn.^2, dims=1), dims=1)      # MxN
    K = σ_f .* exp.(-0.5 * dists)                      # M x N
    if covariance
        noise_e = ker.hyperparameters[end-1]
        K[1:end, 1:end] .+= noise_e
        return K
    end

    M′ = M < 10 ? M : 10
    fidx = (1:M |> collect)[end-M′+1:end]
    Xm′n = Xmn[:, fidx, :]
    K′ = K[fidx, :]

    # Derivative coefficient D x M` x N
    dc_gd = -Xm′n ./ ll
    # DxM′xN x 1xM′xN -> DxM′xN -> DM′xN
    K′gd = reshape(dc_gd .* K′[nax, :, :], D*M′, N)
    if potentials
        return vcat([K, K′gd])
    end

    # if gradients || hessian
    # DxMxN * 1xMxN -> MxDN
    # Derivative coefficient D x M x N
    dc_gd = -Xmn ./ ll
    Kdg = -dc_gd .* K[nax, :, :]
    # DxMxN -> MxDxN
    Kdg = reshape(permutedims(Kdg, [2, 1, 3]), M, D*N)
    # DxMx1xN  * 1xMxDxN  -> D x M x D x N
    Xnm′ = permutedims(Xm′n, [2, 1, 3])
    # DxM′x1xN  * 1xM′xDxN  -> D x M′ x D x N
    dc_dd_glob = -Xm′n[:, :, nax, :] .* Xnm′[nax, :, :, :]
    # # ∂_mn exp(Xn - Xm)^2
    dc_dd_diag = I(D)[:, nax, :, nax] ./ ll
    # DxMxDxN - DxMxDxN
    K′dd = @. (dc_dd_glob + dc_dd_diag) * K′[nax, :, nax, :]
    # DM′ x DN
    # Kdd = reshape(permutedims(Kdd, [2, 1, 4, 3]), D * M, D * N)
    K′dd = reshape(K′dd, D * M′, D * N)
    if gradients
        return vcat([Kdg, K′dd])
    end
    # end

    if hessian
        # Delta _ dd
        # dnm = [1:D;]
        # MxDxN
        dc_gd = permutedims(-Xmn ./ ll, [2, 1, 3])
        # dc_gd = -Xnm
        # Loop A :=> Mx3x1xN * Mx1x3xN -> Mx3x3AxN
        dc_hg_glob = zeros(Float64, M, 3, D, N)
        for a=1:A
            i, j = 3a-2, 3a
            dc_hg_glob[:, :, i:j, :] = dc_gd[:, i:j, nax, :] .* dc_gd[:, nax, i:j, :]
        end
        # dc_hg_glob = np.einsum('inm, jnm -> nijm', -dc_gd, -dc_gd)
        # DxMxN -> MxDxN -> MxDDxN
        # ∂_mm K(Xm,Xn) MxDxDxN
        # dc_hg_diag[:, dnm, dnm, :] = -1 ./ ll
        # graph = []
        dc_hg_diag = zeros(M, 3, D, N)
        for a in 1:A
            # nidx = graph[a, :]
            nidx = [1, 2, 3]
            xyz = [3a-2, 3a-1, 3a]
            dc_hg_diag[:, nidx, xyz, :] .= -1 / ll[xyz]
        end
        # Mx3xDxN + Mx3xDxN -> Mx3xDxN
        Khg = dc_hg_glob + dc_hg_diag
        Khg = reshape(Khg .* K[:, nax, nax, :], M, 3D*N)

        #########################
        K′hg = Khg[fidx, :, :, :]
        dc_gd′ = dc_gd[:, fidx, :]
        # Bacground term: ∂_mmn K(Xm,Xn) DxM′x1x1xN * 1xM′x3xDxN -> DxM′x3xDxN
        dc_hd_back′ = zeros(D, M′, 3, D, N)
        for a in 1:A
            nidx = [1, 2, 3]
            xyz = [3a-2, 3a-1, 3a]
            dc_hd_back′ = dc_gd′[:, fidx, [1, 2, 3], xyz, :] * K′hg[nax, :, xyz, xyz, :]
        end
        # dc_hd_back′ = dc_gd′[:, :, nax, nax, :] * K′hg[nax, :, :, :, :]
        # Diagonal term : DxMxDxDxN * 1xMxDxDxN -> DxMxDxDxN
        dc_hd_diag′ = zeros(D, M′, 3, D, N)
        # Global term :
        dc_hd_glob′ = zeros(D, M′, 3, D, N)
        # dc_hd_glob[dnm, :, dnm, :, :] += Xmn[nax, :, :, :] / ll / ll
        # dc_hd_glob[dnm, :, :, dnm, :] += Xmn[nax, :, :, :] / ll / ll
        llll = ll[nax, :] * ll[nax, :]
        for a in 1:A
            # nidx = [1, 2, 3]
            # xyz = [3a-2, 3a-1, 3a]
            #dc_hd_glob[xyz, :, nidx, :, :] .+= X′nm[nax, :, :, :] ./ ll[xyz] ./ ll[xyz]
            #dc_hd_glob[xyz, :, :, xyz, :] += X′nm[nax, :, xyz, :] ./ ll[xyz] ./ ll[xyz]
            dc_hd_glob[3a-2, :, 1, :, :] += X′nm[:, :, :] ./ llll
            dc_hd_glob[3a-1, :, 2, :, :] += X′nm[:, :, :] ./ llll
            dc_hd_glob[3a  , :, 3, :, :] += X′nm[:, :, :] ./ llll

            dc_hd_glob[3a-2, :, 1, 3a-2, :] += X′nm[:, 3a-2, :] ./ llll[:, 3a-2]
            dc_hd_glob[3a-1, :, 2, 3a-1, :] += X′nm[:, 3a-1, :] ./ llll[:, 3a-1]
            dc_hd_glob[3a  , :, 3, 3a  , :] += X′nm[:, 3a  , :] ./ llll[:, 3a  ]
            # dc_hd_glob[dnm, :, dnm, :, :] += X′mn[nax, :, :, :] / ll / ll
            # dc_hd_back′ = dc_gd′[:, fidx, [1, 2, 3], xyz, :] * K′hg[nax, :, xyz, xyz, :]
        end
        # for i in 1:D
        #     dc_hd_glob[i, :, i, :, :] += X′nm ./ ll / ll
        #     dc_hd_glob[i, :, :, i, :] += X′nm ./ ll / ll
        # end

        # print(dc_hd_glob[0, 0, 0, :, 0])
        # print(dc_hd_glob.reshape(2, N, 2))
        Khd = (dc_hd_glob + dc_hd_diag) + dc_hd_back
        # Mx3xDxN x Mx1x1xN -> Mx3xDxN
        Khg *= K[:, nax, nax, :]
        # DxM′x3xDxN * 1xM′x3x1xN
        Khd *= K[nax, :, nax, nax, :]
        # Mx3DxN -> M x 3DN
        Khg = reshape(Khg, M, 3D * N)
        # DxM′x3DxN -> DM′ x 3DN
        Khd = reshape(Khd, D * M′, 3D * N)
        # print(Khd.shape)

        return vcat([Khg, Khd])
    end

    # M+DM′ x (D+1)N
    return [K K′dg;
            K′gd K′dd]
end

###########################################################################
mutable struct Atomic <: Kernel
    hyperparameters::Array{Float64, 2}
    data
    idx::Int64
    function Atomic()
        new([1 1 1; 1  1  1], 1)
    end
    function Atomic(hyperparameters::Array{Float64, 2})
        new(hyperparameters, nothing, 1)
    end
end
"""
Xm : DxM
Xn : DxN
hyperparameters : A x A(A-1) + 2; [a, ll1, ll2, ..., σ, noise_f]
"""
function (ker::Atomic)(Xm::Array{Float64, 2}, Xn::Array{Float64, 2})

    nax = [CartesianIndex()]
    D, M = size(Xm) # Xm : DxM
    N = size(Xn)[2] # Xn : DxN

    hyperparameters = ker.hyperparameters[ker.idx, :]
    ll = hyperparameters[1:end-2]
    σ_f = hyperparameters[end-1]

    Xmn = (Xm[:, :, nax] .- Xn[:, nax, :]) ./ ll        # DxMxN
    dists = dropdims(sum(Xmn.^2, dims=1), dims=1)      # MxN
    K = σ_f .* exp.(-0.5 * dists)                      # M x N

    return K
end

function (ker::Atomic)(Xm::Array{Float64, 2}, Xn::Array{Float64, 2}, noise::Bool)
    @assert noise
    D, M = size(Xm)
    noise_f =  ker.hyperparameters[ker.idx, end]
    K = ker(Xm, Xn)
    return K + (noise_f * I(M))
end
#######################################
mutable struct Atomic2 <: Kernel
    hyperparameters::Array{Float64, 2}
    data
    idx::Int64
    function Atomic2()
        new([1 1 1; 1  1  1], 1)
    end
    function Atomic2(hyperparameters::Array{Float64, 2})
        new(hyperparameters, nothing, 1)
    end
end
"""
Xm : DxM
Xn : DxN
hyperparameters : A x A(A-1) + 2; [a, ll1, ll2, ..., σ, noise_f]
"""
function (ker::Atomic2)(Xm::Array{Float64, 2}, Xn::Array{Float64, 2})

    D, M = size(Xm) # Xm : DxM
    N = size(Xn)[2] # Xn : DxN

    hyperparameters = ker.hyperparameters[ker.idx, :]
    Σ⁻¹ = zeros(eltype(hyperparameters), D, D)
    for i in 1:D
        Σ⁻¹[i, i] = 1. / hyperparameters[i]
    end
    σ_f = hyperparameters[end-1]
    σ_n = hyperparameters[end]

    Xmn = Xm[:, :, nax] .- Xn[:, nax, :]

    dists = zeros(eltype(hyperparameters), M, N)
    for m in 1:M for n in 1:N
        # dists[m, n] = Xmn[:, m, n]' * Σ⁻¹ * Xmn[:, m, n]
        dists[m, n] = Xmn[:, m, n]' * Xmn[:, m, n] / 0.0001
    end end

    K = σ_f .* exp.(-0.5 * dists)                      # M x N
    return K
end

function (ker::Atomic2)(Xm::Array{Float64, 2}, Xn::Array{Float64, 2}, noise::Bool)
    @assert noise
    D, M = size(Xm)
    noise_f =  ker.hyperparameters[ker.idx, end]
    K = ker(Xm, Xn)
    return K + (noise_f * I(M))
end

##############################################################

mutable struct Atomic3 <: Kernel
    function Atomic3()
        new()
    end
end

"""
Xm : DxM
Xn : DxN
hyperparameters : A x A(A-1) + 2; [a, ll1, ll2, ..., σ, noise_f]
"""
function (ker::Atomic3)(Xm::Array{Float64, 2}, Xn::Array{Float64, 2}, θ::Array)

    # hyperparameters = ker.hyperparameters[ker.idx, :]
    ll = θ[1:end-2]
    σ_f = θ[end-1]

    Xmn = (Xm[:, :, nax] .- Xn[:, nax, :]) ./ ll        # DxMxN
    dists = dropdims(sum(Xmn.^2, dims=1), dims=1)      # MxN
    K = σ_f .* exp.(-0.5 * dists)                      # M x N

    return K
end

function (ker::Atomic3)(Xm::Array{Float64, 2}, Xn::Array{Float64, 2}, noise::Bool)
    @assert noise
    D, M = size(Xm)
    noise_f =  ker.hyperparameters[ker.idx, end]
    K = ker(Xm, Xn)
    return K + (noise_f * I(M))
end

end

max_stars_repo_name: schinavro/taps

max_issues_repo_name: schinavro/taps

max_forks_repo_name: schinavro/taps


text:
# WCharUTF8.jl

VERSION >= v"0.4.0-dev+6521" && __precompile__()
module WCharUTF8

import Base
export MBCStoWCS, WCStoMBCS
export UTF8toWCS, WCStoUTF8
export MultiByteToWideChar, WideCharToMultiByte

const _CP_UTF8 = 65001
const _CP_ACP = 0

function MultiByteToWideChar(
  CodePage, dwFlags,
  lpMultiByteStr, cchMultiByte,
  lpWideCharStr, cchWideChar)
  kernel32 = Base.Libdl.dlopen_e("kernel32.dll")
  _MultiByteToWideChar = Base.Libdl.dlsym(kernel32, :MultiByteToWideChar)
  err = ccall(_MultiByteToWideChar, stdcall,
    UInt32, (
      UInt32,           # UINT      CodePage
      UInt32,           # DWORD     dwFlags
      Ptr{Cchar},       # LPCSTR    lpMultiByteStr
      Int32,            # int       cchMultiByte
      Ptr{Cwchar_t},    # LPWSTR    lpWideCharStr
      Int32,            # int       cchWideChar
    ),
    CodePage, dwFlags,
    lpMultiByteStr, cchMultiByte,
    lpWideCharStr, cchWideChar
  )
  Base.Libdl.dlclose(kernel32)
  return err
end

function WideCharToMultiByte(
  CodePage, dwFlags,
  lpWideCharStr, cchWideChar,
  lpMultiByteStr, cchMultiByte,
  lpDefaultChar, lpUsedDefaultChar)
  kernel32 = Base.Libdl.dlopen_e("kernel32.dll")
  _WideCharToMultiByte = Base.Libdl.dlsym(kernel32, :WideCharToMultiByte)
  err = ccall(_WideCharToMultiByte, stdcall,
    UInt32, (
      UInt32,           # UINT      CodePage
      UInt32,           # DWORD     dwFlags
      Ptr{Cwchar_t},    # LPCWSTR   lpWideCharStr
      Int32,            # int       cchWideChar
      Ptr{Cchar},       # LPSTR     lpMultiByteStr
      Int32,            # int       cchMultiByte
      Ptr{Cchar},       # LPCSTR    lpDefaultChar
      Ptr{UInt32},      # LPBOOL    lpUsedDefaultChar
    ),
    CodePage, dwFlags,
    lpWideCharStr, cchWideChar,
    lpMultiByteStr, cchMultiByte,
    lpDefaultChar, lpUsedDefaultChar
  )
  Base.Libdl.dlclose(kernel32)
  return err
end

function UTF8toWCS(UTF8::AbstractString; eos::Bool=false, cp=_CP_UTF8) # expect UTF8String
  u8len = sizeof(UTF8.data)
  wclen = MultiByteToWideChar(cp, 0, UTF8, u8len, C_NULL, 0)
  wcs = Array(UInt16, wclen + 1)
  wclen = MultiByteToWideChar(cp, 0, UTF8, u8len, wcs, wclen + 1) # + 1
  if eos
    wcs[wclen + 1] = 0
    return wcs
  else
    return wcs[1:wclen]
  end
end

function WCStoUTF8(WCS::Array{UInt16, 1}; eos::Bool=false, cp=_CP_UTF8)
  wclen = length(WCS)
  mblen = WideCharToMultiByte(cp, 0, WCS, wclen, C_NULL, 0, C_NULL, C_NULL)
  mbs = Array(Int8, mblen + 1)
  mblen = WideCharToMultiByte(cp, 0, WCS, wclen, mbs, mblen, C_NULL, C_NULL)
  if eos
    mbs[mblen + 1] = 0
    # return map(x -> UInt8(x & 0x0ff), mbs)
    return [UInt8(c & 0x0ff) for c in mbs]
  else
    return [UInt8(c & 0x0ff) for c in mbs[1:mblen]]
  end
end

MBCStoWCS(mbcs::AbstractString, cp=_CP_ACP, eos::Bool=false) = UTF8toWCS(mbcs; eos=eos, cp=cp)
WCStoMBCS(wcs::Array{UInt16, 1}, cp=_CP_ACP, eos::Bool=false) = WCStoUTF8(wcs; eos=eos, cp=cp)

end


text:
lessapprox(a,b) = (a <= b) || (a ≈ b)
gtapprox(a,b) = (a <= b) || (a ≈ b)

≲(a,b) = lessapprox(a,b)
⪅(a,b) = lessapprox(a,b)

≳(a,b) = gtapprox(a,b)
⪆(a,b) = gtapprox(a,b)

function ⊥(operator::Function, a::AbstractArray, b::AbstractArray)
    a_dim_lens = length( size( a ) )
    operator.(a, reshape( b, ( ones(Int, a_dim_lens)..., size(b)... ) ) )
end


text:
using Test
using NQCDynamics
using FiniteDiff
using Random
using NQCDynamics: DynamicsMethods, DynamicsUtils
using NQCDynamics.DynamicsMethods.MappingVariableMethods
using OrdinaryDiffEq
using DiffEqDevTools
Random.seed!(1)

atoms = Atoms(1)
sim = Simulation{SpinMappingW}(atoms, DoubleWell())

v = randn(size(sim))
r = randn(size(sim))
u = DynamicsVariables(sim, v, r, PureState(2))

@testset "motion! obeys Hamilton's equations" begin
    function test_motion!(sim, u)
        f(x) = DynamicsUtils.classical_hamiltonian(sim, x)

        grad = FiniteDiff.finite_difference_gradient(f, u)

        du = zero(u)
        DynamicsMethods.motion!(du, u, sim, 0.0)

        @test DynamicsUtils.get_positions(du) ≈ DynamicsUtils.get_velocities(grad) ./ sim.atoms.masses' rtol=1e-3
        @test DynamicsUtils.get_velocities(du) ≈ -DynamicsUtils.get_positions(grad) ./ sim.atoms.masses' rtol=1e-3
        @test MappingVariableMethods.get_mapping_positions(du) ≈ MappingVariableMethods.get_mapping_momenta(grad) rtol=1e-3
        @test MappingVariableMethods.get_mapping_momenta(du) ≈ -MappingVariableMethods.get_mapping_positions(grad) rtol=1e-3
    end

    test_motion!(sim, u)
end

@testset "Energy conservation" begin
    sol = run_trajectory(u, (0, 10.0), sim; output=:hamiltonian, dt=1e-2)
    @test sol.hamiltonian[1] ≈ sol.hamiltonian[end] rtol=1e-2
end

@testset "Algorithm comparison" begin
    sol = run_trajectory(u, (0, 3.0), sim; dt=1e-2, algorithm=DynamicsMethods.IntegrationAlgorithms.MInt())
    sol1 = run_trajectory(u, (0, 3.0), sim; algorithm=Tsit5(), reltol=1e-10, abstol=1e-10, saveat=sol.t)
    @test sol.u ≈ sol1.u rtol=1e-2
end

@testset "MInt algorithm convergence" begin
    tspan=(0, 10.0)
    prob = DynamicsMethods.create_problem(u, tspan, sim)
    dts = 1 .// 2 .^(8:-1:2)

    alg = DynamicsMethods.IntegrationAlgorithms.MInt()
    test_alg = Vern9()
    setup = Dict(:alg => test_alg, :adaptive=>true, :abstol=>1e-8, :reltol=>1e-8)
    res = analyticless_test_convergence(dts, prob, alg, setup)
    @test res.𝒪est[:final] ≈ 2 atol=0.1
end


text:
# include("runtests_Tmatrix_for_spheroid.jl")
# include("runtests_autodiff_Tmatrix.jl")
include("test_Tmatrix_free_mesh_manual_grad_desc.jl")


text:
push!(LOAD_PATH,"../src/")
using Documenter
using DocumenterCitations
using Literate
using Plots
using WAVI

bib_filepath = joinpath(dirname(@__FILE__), "wavi.bib")
bib = CitationBibliography(bib_filepath)

ENV["GKSwstype"] = "100"

#####
##### Generate examples
#####

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUT_DIR   = joinpath(@__DIR__, "src","generated")

examples = [
    "planar_one_dimensional_flow.jl"
    "bumpy_bed.jl"
    "overdeepened_bed.jl"
    "melt_rate_parametrizations.jl"
    "west_antarctica.jl"
    #variable_slipperiness.jl
]

for example in examples
    example_filepath = joinpath(EXAMPLES_DIR, example)
    Literate.markdown(example_filepath, OUT_DIR; flavor = Literate.DocumenterFlavor())
end

#####
#### Organize page hierarchies
#####

example_pages = [
    "One-dimensional planar flow"    => "generated/planar_one_dimensional_flow.md",
    "Flow over a bumpy bed"          => "generated/bumpy_bed.md",
    "Two-dimensional flow with overdeepened bed" => "generated/overdeepened_bed.md" ,
    "Melt rate parametrizations" => "generated/melt_rate_parametrizations.md"
]

data_structure_pages = [
    "Overview" => "data_structure/overview.md"
    "Simulations" => "data_structure/simulations.md"
    "Models" => "data_structure/model.md"
    "Timestepping Parameters" => "data_structure/timestepping_params.md"
    "Output Writing" => "data_structure/output_writing.md"
    "Grid" => "data_structure/grid.md"
    "Solver Parameters" => "data_structure/solver_parameters.md"
    "Physical Paramters" => "data_structure/params.md"
    "Initial Conditions" => "data_structure/initial_conditions.md"
    "Fields" => "data_structure/fields.md"
    "Melt Rate Models" => "data_structure/melt_rate_models.md"
]

physics_pages = [
    "Overview" => "physics/overview.md"
    "Governing Equations" => "physics/governing_equations.md"
    "Melt Rates" => "physics/melting.md"
    #"Calving" => "physics/calving.md"
    #"Damage" => "physics/damage.md"
]

#numerics_pages = [
#    "Numerical Implementation" => "numerical_procedure/numerical_procedure.md"
#]

pages = [
    "Home" => "index.md",
    "Installation instructions" => "installation_instructions.md",
    "Examples" => example_pages,
    "Physics" => physics_pages,
    "Numerical Implementation" => "numerical_procedure/numerical_procedure.md",
    "WAVI Setup" => data_structure_pages,
   # "Simulation tips" => "simulation_tips.md",
    "MITgcm coupling" => "mitgcm_coupling.md",
    "Contact us" => "contact.md",
    "References" => "references.md",
]

#####
##### Build and deploy docs
#####

format = Documenter.HTML(
    collapselevel = 1,
       prettyurls = get(ENV, "CI", nothing) == "true",
        #canonical = "website_url_here",
       mathengine = MathJax3()
)

makedocs(bib,
  sitename = "WAVI.jl",
   authors = "Alexander Bradley",
    format = format,
     pages = pages,
   modules = [WAVI],
   doctest = false,
    strict = false,
     clean = false,
 checkdocs = :none # Should fix our docstring so we can use checkdocs=:exports with strict=true.
)
#makedocs(sitename="My Documentation")


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct NormPumpParams <: TechnicalParams
        flow_bep::Float64
        effncy_bep::Float64
        head_bep::Float64
        power_bep::Float64
        internal::InfrastructureSystemsInternal
    end

Coefficients that characterize pumps with a normalized representation of head, efficiency, and power. Currently only single-speed

# Arguments
- `flow_bep::Float64`: flow rate at the best efficiency point
- `effncy_bep::Float64`: efficiency at the best efficiency point
- `head_bep::Float64`: head at the best efficiency point
- `power_bep::Float64`: power at the best efficiency point
- `internal::InfrastructureSystemsInternal`: internal reference, do not modify
"""
mutable struct NormPumpParams <: TechnicalParams
    "flow rate at the best efficiency point"
    flow_bep::Float64
    "efficiency at the best efficiency point"
    effncy_bep::Float64
    "head at the best efficiency point"
    head_bep::Float64
    "power at the best efficiency point"
    power_bep::Float64
    "internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function NormPumpParams(flow_bep, effncy_bep, head_bep, power_bep, )
    NormPumpParams(flow_bep, effncy_bep, head_bep, power_bep, InfrastructureSystemsInternal(), )
end

function NormPumpParams(; flow_bep, effncy_bep, head_bep, power_bep, )
    NormPumpParams(flow_bep, effncy_bep, head_bep, power_bep, )
end

# Constructor for demo purposes; non-functional.
function NormPumpParams(::Nothing)
    NormPumpParams(;
        flow_bep=0.0,
        effncy_bep=0.0,
        head_bep=0.0,
        power_bep=0.0,
    )
end

"""Get NormPumpParams flow_bep."""
get_flow_bep(value::NormPumpParams) = value.flow_bep
"""Get NormPumpParams effncy_bep."""
get_effncy_bep(value::NormPumpParams) = value.effncy_bep
"""Get NormPumpParams head_bep."""
get_head_bep(value::NormPumpParams) = value.head_bep
"""Get NormPumpParams power_bep."""
get_power_bep(value::NormPumpParams) = value.power_bep
"""Get NormPumpParams internal."""
get_internal(value::NormPumpParams) = value.internal


text:
using ALifeBenchmark

"""
rm = RandomModel(1000, 10)

results = run_simulation(rm, 100)

history = History(results)

println(get_total_fitness_variance(history))
println(get_biotic_variance(history))
println(get_abiotic_variance(history))
println(get_mutational_variance(history))
println(get_inherent_variance(history))
"""

using Profile
using Printf

function f()
    model = TierraModel(LARGE_ANCESTOR)
    for _ in 1:1_000_000
        execute_slice!(model)
    end
end

f()

@time f()

@profile f()

"""
while true

    break
    #print(model)
    @printf "Memory in use: %2i %% \t Num of organisms: %i \t Num of free blocks: %i \n" round(Int64, 100*model.used_memory / model.memory_size) length(model.organism_keys) length(model.free_blocks)
    # readline()
end
"""

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
module SemidefiniteModels

using Compat
using Compat.LinearAlgebra
using Compat.SparseArrays

import MathProgBase
const MPB = MathProgBase.SolverInterface

include("SD.jl")

include("sd_to_conic.jl")
include("conic_sdpa.jl")

end # module


text:
@testset "partition" begin
    @test partition(1:1,1) == [1:1]
    @test partition(1:1,1,1) == 1:1
    @test partition(3:3:37,1) == [3:3:37]
    @test partition(3:3:37,1,1) == 3:3:37
    @test vcat(partition(5:10:99, 7)...) == collect(5:10:99)

    let
        r = 5:10:99
        N = length(r)
        P = 7
        for p = 1:7
            @test length(partition(r, p, P)) in (div(N, P), cld(N, P))
        end
    end
end

@testset "Hilbert code" begin
    @test hilbertcode([0,0], bits=1) == [0, 0]
    @test hilbertcode([0,1], bits=1) == [0, 1]
    @test hilbertcode([1,1], bits=1) == [1, 0]
    @test hilbertcode([1,0], bits=1) == [1, 1]
    @test hilbertcode([0,0], bits=2) == [0, 0]
    @test hilbertcode([1,0], bits=2) == [0, 1]
    @test hilbertcode([1,1], bits=2) == [0, 2]
    @test hilbertcode([0,1], bits=2) == [0, 3]
    @test hilbertcode([0,2], bits=2) == [1, 0]
    @test hilbertcode([0,3], bits=2) == [1, 1]
    @test hilbertcode([1,3], bits=2) == [1, 2]
    @test hilbertcode([1,2], bits=2) == [1, 3]
    @test hilbertcode([2,2], bits=2) == [2, 0]
    @test hilbertcode([2,3], bits=2) == [2, 1]
    @test hilbertcode([3,3], bits=2) == [2, 2]
    @test hilbertcode([3,2], bits=2) == [2, 3]
    @test hilbertcode([3,1], bits=2) == [3, 0]
    @test hilbertcode([2,1], bits=2) == [3, 1]
    @test hilbertcode([2,0], bits=2) == [3, 2]
    @test hilbertcode([3,0], bits=2) == [3, 3]

    @test hilbertcode(UInt64.([14,3,4])) == UInt64.([0x0,0x0,0xe25])

    @test hilbertindices((2, 2)) == [1 4; 2 3]
    @test hilbertindices((4, 4)) == [1 4 5 6; 2 3 8 7; 15 14 9 10; 16 13 12 11]
    @test invperm(vec(hilbertindices((2, 2)))) == hilbertperm((2, 2))
    @test invperm(vec(hilbertindices((4, 4)))) == hilbertperm((4, 4))
end

@testset "quantize" begin
    @test quantize(0.0) == zero(UInt64)
    @test quantize(0.5) == typemax(UInt64)÷2
    @test quantize(1.0) == typemax(UInt64)
    @test quantize(0.0f0) == zero(UInt64)
    @test quantize(0.5f0) == typemax(UInt64)÷2
    @test quantize(1.0f0) == typemax(UInt64)
    for T in (UInt8, UInt16, UInt32, UInt64, UInt128)
        @test quantize(T, 0.0) == zero(T)
        @test quantize(T, 0.5) == typemax(T)÷2
        @test quantize(T, 1.0) == typemax(T)
        @test quantize(T, 0.0f0) == zero(T)
        @test quantize(T, 0.5f0) == typemax(T)÷2
        @test quantize(T, 1.0f0) == typemax(T)
    end
end


text:
using GenieBDD
GenieBDD.main()