text:
function evolve_system(	particle_type::String,
						R::Array{Float64, 2},
						Lx::Float64,
						Ly::Float64,
						Lz::Float64,
						X::Array{Float64, 1},
						Y::Array{Float64, 1},
						Z::Array{Float64, 1},
						Q0::Array{Float64, 1},
						Q1::Array{Float64, 1},
						Q2::Array{Float64, 1},
						Q3::Array{Float64, 1},
						A11::Array{Float64, 1},
						A12::Array{Float64, 1},
						A13::Array{Float64, 1},
						A21::Array{Float64, 1},
						A22::Array{Float64, 1},
						A23::Array{Float64, 1},
						A31::Array{Float64, 1},
						A32::Array{Float64, 1},
						A33::Array{Float64, 1},
						position_constraint_axis::String,
						position_constraint_lower::Float64,
						position_constraint_upper::Float64,
						orientation_axis::Array{Float64, 1},
						orientation_constraint_axis::Array{Float64, 1},
						orientation_constraint_lower::Float64,
						orientation_constraint_upper::Float64,
						sigma_translation::Float64,
						sigma_translation_max::Float64,
						sigma_rotation::Float64,
						sigma_rotation_max::Float64,
						sigma_ratio::Float64,
						number_of_relaxation_sweeps_max::Int64,
						number_of_equilibration_sweeps::Int64)

	is_relaxed::Bool = false

	number_of_particles::Int64 = size(R, 1)

	acceptance_probability_target::Float64 = 0.25

	# Pre-computed max radii (i.e. radius of bounding sphere).
	RMAX::Array{Float64, 1} = zeros(number_of_particles)
	if particle_type == "sphere"
		for current_particle = 1:number_of_particles
			RMAX[current_particle] = R[current_particle, 1]
		end
	elseif particle_type == "ellipse"
		for current_particle = 1:number_of_particles
			RMAX[current_particle] = maximum( R[current_particle, 1:2] )
		end
	elseif particle_type == "ellipsoid"
		for current_particle = 1:number_of_particles
			RMAX[current_particle] = maximum( R[current_particle, 1:3] )
		end
	elseif particle_type == "cuboid"
		for current_particle = 1:number_of_particles
			RMAX[current_particle] = sqrt(R[current_particle, 1]^2 + R[current_particle, 2]^2 + R[current_particle, 3]^2)
		end
	end

	# Preallocation.
	x_star::Float64 = 0.0
	y_star::Float64 = 0.0
	z_star::Float64 = 0.0
	q0_star::Float64 = 0.0
	q1_star::Float64 = 0.0
	q2_star::Float64 = 0.0
	q3_star::Float64 = 0.0
	a11_star::Float64 = 0.0
	a12_star::Float64 = 0.0
	a13_star::Float64 = 0.0
	a21_star::Float64 = 0.0
	a22_star::Float64 = 0.0
	a23_star::Float64 = 0.0
	a31_star::Float64 = 0.0
	a32_star::Float64 = 0.0
	a33_star::Float64 = 0.0

	energy_system::Float64 = 0.0
	energy_particle::Float64 = 0.0
	energy_particle_star::Float64 = 0.0
	current_relaxation_sweep::Int64 = 0
	current_equilibration_sweep::Int64 = 0
	acceptance_probability_translation::Float64 = 0.0
	acceptance_probability_rotation::Float64 = 0.0
	acceptance_probability_combined::Float64 = 0.0
	xAB::Float64 = 0.0
	yAB::Float64 = 0.0
	zAB::Float64 = 0.0
	overlapfun::Float64 = 0.0

	is_position_ok::Bool = false

	orientation_axis_rotated::Array{Float64, 1} = zeros(3)
	angle_to_orientation_constraint_axis::Float64 = 0.0
	is_orientation_ok::Bool = false

	energy_system = 1.0
	while (energy_system > 0.0 && current_relaxation_sweep < number_of_relaxation_sweeps_max) || current_equilibration_sweep < number_of_equilibration_sweeps
		if energy_system > 0.0
			current_relaxation_sweep += 1
			#println(join(("Relaxation sweep ", string(current_relaxation_sweep))))
		else
			current_equilibration_sweep += 1
			#println(join(("Equlibration sweep ", string(current_equilibration_sweep), " out of ", string(number_of_equilibration_sweeps))))
		end

		acceptance_probability_translation = 0.0
		acceptance_probability_rotation = 0.0

		energy_system = 0.0

		for currentA = 1:number_of_particles
			# Compute current local energy.
			energy_particle = 0.0
			for currentB = [1:currentA-1 ; currentA+1:number_of_particles]
				xAB = signed_distance_mod(X[currentA], X[currentB], Lx)
				yAB = signed_distance_mod(Y[currentA], Y[currentB], Ly)
				zAB = signed_distance_mod(Z[currentA], Z[currentB], Lz)

				if xAB^2 + yAB^2 + zAB^2 < (RMAX[currentA] + RMAX[currentB])^2
					if particle_type == "sphere"
						overlapfun = (RMAX[currentA] + RMAX[currentB])^2 - (xAB^2 + yAB^2 + zAB^2)
						energy_particle += overlapfun
					elseif particle_type == "ellipse"
						overlapfun = overlap_ellipse(xAB, yAB, zAB, A11[currentA], A12[currentA], A13[currentA], A21[currentA], A22[currentA], A23[currentA], A31[currentA], A32[currentA], A33[currentA], A11[currentB], A12[currentB], A13[currentB], A21[currentB], A22[currentB], A23[currentB], A31[currentB], A32[currentB], A33[currentB])

						if overlapfun < 1.0
							energy_particle += (1.0 - overlapfun)^2
						end
					elseif particle_type == "ellipsoid"
						overlapfun = overlap_ellipsoid(xAB, yAB, zAB, A11[currentA], A12[currentA], A13[currentA], A21[currentA], A22[currentA], A23[currentA], A31[currentA], A32[currentA], A33[currentA], A11[currentB], A12[currentB], A13[currentB], A21[currentB], A22[currentB], A23[currentB], A31[currentB], A32[currentB], A33[currentB], R[currentA, 1]^2 * R[currentA, 2]^2 * R[currentA, 3]^2)

						if overlapfun < 1.0
							energy_particle += (1.0 - overlapfun)^2
						end
					elseif particle_type == "cuboid"
						overlapfun = overlap_cuboid(xAB, yAB, zAB, A11[currentA], A12[currentA], A13[currentA], A21[currentA], A22[currentA], A23[currentA], A31[currentA], A32[currentA], A33[currentA], A11[currentB], A12[currentB], A13[currentB], A21[currentB], A22[currentB], A23[currentB], A31[currentB], A32[currentB], A33[currentB], R[currentA, 1], R[currentA, 2], R[currentA, 3], R[currentB, 1], R[currentB, 2], R[currentB, 3])

						energy_particle += overlapfun
					end

				end
			end

			# Generate random proposal position and compute new local energy with translation.
			is_position_ok = false
			while !is_position_ok
				(x_star, y_star, z_star) = generate_proposal_position(X[currentA], Y[currentA], Z[currentA], Lx, Ly, Lz, sigma_translation)
				if position_constraint_axis == "x"
					if position_constraint_lower * Lx <= x_star <= position_constraint_upper * Lx
						is_position_ok = true
					end
				elseif position_constraint_axis == "y"
					if position_constraint_lower * Ly <= y_star <= position_constraint_upper * Ly
						is_position_ok = true
					end
				elseif position_constraint_axis == "z"
					if position_constraint_lower * Lz <= z_star <= position_constraint_upper * Lz
						is_position_ok = true
					end
				end
			end

			energy_particle_star = 0.0
			for currentB = [1:currentA-1 ; currentA+1:number_of_particles]
				xAB = signed_distance_mod(x_star, X[currentB], Lx)
				yAB = signed_distance_mod(y_star, Y[currentB], Ly)
				zAB = signed_distance_mod(z_star, Z[currentB], Lz)

				if xAB^2 + yAB^2 + zAB^2 < (RMAX[currentA] + RMAX[currentB])^2
					if particle_type == "sphere"
						overlapfun = (RMAX[currentA] + RMAX[currentB])^2 - (xAB^2 + yAB^2 + zAB^2)
						energy_particle_star += overlapfun
					elseif particle_type == "ellipse"
						overlapfun = overlap_ellipse(xAB, yAB, zAB, A11[currentA], A12[currentA], A13[currentA], A21[currentA], A22[currentA], A23[currentA], A31[currentA], A32[currentA], A33[currentA], A11[currentB], A12[currentB], A13[currentB], A21[currentB], A22[currentB], A23[currentB], A31[currentB], A32[currentB], A33[currentB])

						if overlapfun < 1.0
							energy_particle_star += (1.0 - overlapfun)^2
						end
					elseif particle_type == "ellipsoid"
						overlapfun = overlap_ellipsoid(xAB, yAB, zAB, A11[currentA], A12[currentA], A13[currentA], A21[currentA], A22[currentA], A23[currentA], A31[currentA], A32[currentA], A33[currentA], A11[currentB], A12[currentB], A13[currentB], A21[currentB], A22[currentB], A23[currentB], A31[currentB], A32[currentB], A33[currentB], R[currentA, 1]^2 * R[currentA, 2]^2 * R[currentA, 3]^2)

						if overlapfun < 1.0
							energy_particle_star += (1.0 - overlapfun)^2
						end
					elseif particle_type == "cuboid"
						overlapfun = overlap_cuboid(xAB, yAB, zAB, A11[currentA], A12[currentA], A13[currentA], A21[currentA], A22[currentA], A23[currentA], A31[currentA], A32[currentA], A33[currentA], A11[currentB], A12[currentB], A13[currentB], A21[currentB], A22[currentB], A23[currentB], A31[currentB], A32[currentB], A33[currentB], R[currentA, 1], R[currentA, 2], R[currentA, 3], R[currentB, 1], R[currentB, 2], R[currentB, 3])

						energy_particle_star += overlapfun
					end

				end
			end

			if energy_particle_star <= energy_particle
				X[currentA] = x_star
				Y[currentA] = y_star
				Z[currentA] = z_star

				acceptance_probability_translation += 1.0
				energy_particle = energy_particle_star
			end

			# Generate random proposal orientation and compute new local energy with rotation.
			if particle_type != "sphere"
				is_orientation_ok = false
				while !is_orientation_ok
					(q0_star, q1_star, q2_star, q3_star) = generate_proposal_orientation(Q0[currentA], Q1[currentA], Q2[currentA], Q3[currentA], sigma_rotation)
					(a11_star, a12_star, a13_star, a21_star, a22_star, a23_star, a31_star, a32_star, a33_star) = rotation_matrix(q0_star, q1_star, q2_star, q3_star)
					orientation_axis_rotated[1] = a11_star * orientation_axis[1] + a12_star * orientation_axis[2] + a13_star * orientation_axis[3]
					orientation_axis_rotated[2] = a21_star * orientation_axis[1] + a22_star * orientation_axis[2] + a23_star * orientation_axis[3]
					orientation_axis_rotated[3] = a31_star * orientation_axis[1] + a32_star * orientation_axis[2] + a33_star * orientation_axis[3]
					angle_to_orientation_constraint_axis = acos(orientation_axis_rotated[1] * orientation_constraint_axis[1] + orientation_axis_rotated[2] * orientation_constraint_axis[2] + orientation_axis_rotated[3] * orientation_constraint_axis[3])
					if orientation_constraint_lower <= angle_to_orientation_constraint_axis <= orientation_constraint_upper
						is_orientation_ok = true
					end
				end

				if particle_type == "ellipse"
					(a11_star, a12_star, a13_star, a21_star, a22_star, a23_star, a31_star, a32_star, a33_star) = characteristic_matrix_ellipse(q0_star, q1_star, q2_star, q3_star, R[currentA, 1], R[currentA, 2])
				elseif particle_type == "ellipsoid"
					(a11_star, a12_star, a13_star, a21_star, a22_star, a23_star, a31_star, a32_star, a33_star) = characteristic_matrix_ellipsoid(q0_star, q1_star, q2_star, q3_star, R[currentA, 1], R[currentA, 2], R[currentA, 3])
				elseif particle_type == "cuboid"
					(a11_star, a12_star, a13_star, a21_star, a22_star, a23_star, a31_star, a32_star, a33_star) = rotation_matrix(q0_star, q1_star, q2_star, q3_star)
				end

				energy_particle_star = 0.0
				for currentB = [1:currentA-1 ; currentA+1:number_of_particles]
					xAB = signed_distance_mod(X[currentA], X[currentB], Lx)
					yAB = signed_distance_mod(Y[currentA], Y[currentB], Ly)
					zAB = signed_distance_mod(Z[currentA], Z[currentB], Lz)

					if xAB^2 + yAB^2 + zAB^2 < (RMAX[currentA] + RMAX[currentB])^2
						if particle_type == "ellipse"
							overlapfun = overlap_ellipse(xAB, yAB, zAB, a11_star, a12_star, a13_star, a21_star, a22_star, a23_star, a31_star, a32_star, a33_star, A11[currentB], A12[currentB], A13[currentB], A21[currentB], A22[currentB], A23[currentB], A31[currentB], A32[currentB], A33[currentB])

							if overlapfun < 1.0
								energy_particle_star += (1.0 - overlapfun)^2
							end
						elseif particle_type == "ellipsoid"
							overlapfun = overlap_ellipsoid(xAB, yAB, zAB, a11_star, a12_star, a13_star, a21_star, a22_star, a23_star, a31_star, a32_star, a33_star, A11[currentB], A12[currentB], A13[currentB], A21[currentB], A22[currentB], A23[currentB], A31[currentB], A32[currentB], A33[currentB], R[currentA, 1]^2 * R[currentA, 2]^2 * R[currentA, 3]^2)

							if overlapfun < 1.0
								energy_particle_star += (1.0 - overlapfun)^2
							end
						elseif particle_type == "cuboid"
							overlapfun = overlap_cuboid(xAB, yAB, zAB, a11_star, a12_star, a13_star, a21_star, a22_star, a23_star, a31_star, a32_star, a33_star, A11[currentB], A12[currentB], A13[currentB], A21[currentB], A22[currentB], A23[currentB], A31[currentB], A32[currentB], A33[currentB], R[currentA, 1], R[currentA, 2], R[currentA, 3], R[currentB, 1], R[currentB, 2], R[currentB, 3])

							energy_particle_star += overlapfun
						end
					end
				end

				if energy_particle_star <= energy_particle
					Q0[currentA] = q0_star
					Q1[currentA] = q1_star
					Q2[currentA] = q2_star
					Q3[currentA] = q3_star

					A11[currentA] = a11_star
					A12[currentA] = a12_star
					A13[currentA] = a13_star
					A21[currentA] = a21_star
					A22[currentA] = a22_star
					A23[currentA] = a23_star
					A31[currentA] = a31_star
					A32[currentA] = a32_star
					A33[currentA] = a33_star

					acceptance_probability_rotation += 1.0
					energy_particle = energy_particle_star
				end
			end

			energy_system += energy_particle
		end

		# Update sigma_translation and sigma_rotation based on acceptance probabilities.
		if particle_type == "sphere"
			acceptance_probability_translation /= number_of_particles
			if acceptance_probability_translation <= acceptance_probability_target
				sigma_translation *= 0.95
			else
				sigma_translation = min(1.05 * sigma_translation, sigma_translation_max)
			end
		else
			acceptance_probability_translation /= number_of_particles
			if acceptance_probability_translation <= acceptance_probability_target
				sigma_translation *= 0.95
			else
				sigma_translation = min(1.05 * sigma_translation, sigma_translation_max)
			end

			acceptance_probability_rotation /= number_of_particles
			if acceptance_probability_rotation <= acceptance_probability_target
				sigma_rotation *= 0.95
			else
				sigma_rotation = min(1.05 * sigma_rotation, sigma_rotation_max, sigma_ratio * sigma_translation)
			end

		#	if sigma_ratio >= 0.0
		#		acceptance_probability_combined = (acceptance_probability_translation + acceptance_probability_rotation) / (2.0 * number_of_particles)
		#		if acceptance_probability_combined <= acceptance_probability_target
		#			sigma_translation *= 0.95
		#		else
		#			sigma_translation = min(1.05 * sigma_translation, sigma_translation_max)
		#		end

		#		sigma_rotation = sigma_ratio * sigma_translation
		#	else
		#		acceptance_probability_translation /= number_of_particles
		#		if acceptance_probability_translation <= acceptance_probability_target
		#			sigma_translation *= 0.95
		#		else
		#			sigma_translation = min(1.05 * sigma_translation, sigma_translation_max)
		#		end

		#		acceptance_probability_rotation /= number_of_particles
		#		if acceptance_probability_rotation <= acceptance_probability_target
		#			sigma_rotation *= 0.95
		#		else
		#			sigma_rotation = min(1.05 * sigma_rotation, sigma_rotation_max)
		#		end
		#	end
		end

#		# Final computation of system energy after finishing the sweep. This needs to be performed because only local optimization is performed when the particles translate and rotate.
#		energy_system = 0.0
#		if particle_type == "sphere"
#			for currentA = 1:number_of_particles
#				for currentB = [1:currentA-1 ; currentA+1:number_of_particles]
#					xAB = signed_distance_mod(X[currentA], X[currentB], Lx)
#					yAB = signed_distance_mod(Y[currentA], Y[currentB], Ly)
#					zAB = signed_distance_mod(Z[currentA], Z[currentB], Lz)
#					overlapfun = (RMAX[currentA] + RMAX[currentB])^2 - (xAB^2 + yAB^2 + zAB^2)
#					if xAB^2 + yAB^2 + zAB^2 <= (RMAX[currentA] + RMAX[currentB])^2
#						energy_system += overlapfun
#					end
#				end
#			end
#		elseif particle_type == "ellipse"
#			for currentA = 1:number_of_particles
#				for currentB = [1:currentA-1 ; currentA+1:number_of_particles]
#					xAB = signed_distance_mod(X[currentA], X[currentB], Lx)
#					yAB = signed_distance_mod(Y[currentA], Y[currentB], Ly)
#					zAB = signed_distance_mod(Z[currentA], Z[currentB], Lz)
#					overlapfun = overlap_ellipse(xAB, yAB, zAB, A11[currentA], A12[currentA], A13[currentA], A21[currentA], A22[currentA], A23[currentA], A31[currentA], A32[currentA], A33[currentA], A11[currentB], A12[currentB], A13[currentB], A21[currentB], A22[currentB], A23[currentB], A31[currentB], A32[currentB], A33[currentB])
#					if overlapfun < 1.0
#						energy_system += (1.0 - overlapfun)^2
#					end
#				end
#			end
#		elseif particle_type == "ellipsoid"
#			for currentA = 1:number_of_particles
#				for currentB = [1:currentA-1 ; currentA+1:number_of_particles]
#					xAB = signed_distance_mod(X[currentA], X[currentB], Lx)
#					yAB = signed_distance_mod(Y[currentA], Y[currentB], Ly)
#					zAB = signed_distance_mod(Z[currentA], Z[currentB], Lz)
#					overlapfun = overlap_ellipsoid(xAB, yAB, zAB, A11[currentA], A12[currentA], A13[currentA], A21[currentA], A22[currentA], A23[currentA], A31[currentA], A32[currentA], A33[currentA], A11[currentB], A12[currentB], A13[currentB], A21[currentB], A22[currentB], A23[currentB], A31[currentB], A32[currentB], A33[currentB], R[currentA, 1]^2 * R[currentA, 2]^2 * R[currentA, 3]^2)
#					if overlapfun < 1.0
#						energy_system += (1.0 - overlapfun)^2
#					end
#				end
#			end
#		elseif particle_type == "cuboid"
#			for currentA = 1:number_of_particles
#				for currentB = [1:currentA-1 ; currentA+1:number_of_particles]
#					xAB = signed_distance_mod(X[currentA], X[currentB], Lx)
#					yAB = signed_distance_mod(Y[currentA], Y[currentB], Ly)
#					zAB = signed_distance_mod(Z[currentA], Z[currentB], Lz)
#					overlapfun = overlap_cuboid(xAB, yAB, zAB, A11[currentA], A12[currentA], A13[currentA], A21[currentA], A22[currentA], A23[currentA], A31[currentA], A32[currentA], A33[currentA], A11[currentB], A12[currentB], A13[currentB], A21[currentB], A22[currentB], A23[currentB], A31[currentB], A32[currentB], A33[currentB], R[currentA, 1], R[currentA, 2], R[currentA, 3], R[currentB, 1], R[currentB, 2], R[currentB, 3])
#					energy_system += overlapfun
#				end
#			end
#		end

	end

	#println((energy_system, sigma_translation, sigma_rotation))

	if energy_system == 0.0
		is_relaxed = true
	end

	return (X, Y, Z, Q0, Q1, Q2, Q3, A11, A12, A13, A21, A22, A23, A31, A32, A33, sigma_translation, sigma_rotation, is_relaxed)
end

max_stars_repo_name: roding/whitefish

max_issues_repo_name: roding/whitefish

max_forks_repo_name: roding/whitefish


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

using Markdown
using InteractiveUtils

# ╔═╡ d794dbb2-0254-4808-80ea-75f184586ad8
begin
	using PlutoUI
	PlutoUI.TableOfContents()
end

# ╔═╡ 32785fd3-0adb-4b61-a8dc-c920904d66ea
using OrderedCollections

# ╔═╡ 3a242ea0-a518-11eb-0c88-d541f2d52fb8
md"""
# Dictionaries
"""

# ╔═╡ ea932952-cb4c-4876-9dc4-3a2b56cf4761
md"""
A dictionary is an **unordered collection** of key-value pairs of elements.

## Declaration

The best way to understand what it means is by looking at one:
"""

# ╔═╡ 61cb5c19-2c9b-4484-b95b-c54c2bc4066e
dict = Dict("🍌" => 4, "🍓" => 15)

# ╔═╡ 94d9dbe3-4042-4bdf-9d14-7acc165514cc
md"""
As we can see a dictionary is declared using the `Dict` function, and then elements are declared using the `key => value` pair syntax. The keys and values can be of any type, and their type can differ between one key or value and another (although it is not recommanded for performance reasons). For example we can use symbols and strings for the keys, and integer and string for the values:
"""

# ╔═╡ 874a6524-8f3c-4bc0-bc6e-9e8fec402e56
Dict(:🍌 => 4, "🍓" => "this is not the season", 12 => 12.0)

# ╔═╡ 45835257-6034-4d64-be18-68c966d96603
md"""
And here's an example with integers:
"""

# ╔═╡ a0be3db2-a1a0-4405-b453-87b933d76ce2
Dict(1234 => 4, 999 => 15)

# ╔═╡ 9f81f1e8-f3c2-4fa4-8469-2d605a418481
md"""
Because it is used to retreive a value, each key in a dictionary must be unique. You couldn't retreive a value if it shared the same index as one another.

Please note also that a dictionary is unordered, meaning we can't index into a dictionary by using the order we used to declare our `key=>value` pairs.
"""

# ╔═╡ de10d9d7-3669-4999-8309-f6c0452f5f10
md"""
We can also declare a dictionary in two steps using the zip function:

- first we declare our keys and values separately:
"""

# ╔═╡ 349fbaf2-1d6e-41c8-a05e-49a1f2d12ac1
mykeys = (:🍌, :🍓)

# ╔═╡ 449eb41a-b5ac-422e-95c7-8baf36bb21f0
myvalues = (4, 15)

# ╔═╡ ae9376e7-3fc0-4668-a495-4338e139385d
md"""
And then we use the zip function to create a `key=>value` pair for each pair of element inside our dictionary:
"""

# ╔═╡ bb28ae8d-5bc2-4056-9c75-9983362c6716
Dict(zip(mykeys,myvalues))

# ╔═╡ 602bc23f-1acc-44a5-8986-cadea9c7f461
md"""
#### Empty dictionary
"""

# ╔═╡ 73dac266-f783-4406-91c7-27a08c4e4e60
md"""
We can declare an empty dictionary using `Dict()`, and even constrain the types of the keys and the values that will be added as follows:
"""

# ╔═╡ 14ac3d30-c0bc-4543-8e61-ccb1bfcd9287
constrained_dict = Dict{String, Float64}()

# ╔═╡ e1856197-b467-465c-ba4c-032282c989eb
constrained_dict["a"] = 2

# ╔═╡ 83c526bb-5db2-4a30-bcfa-7476faaf9f55
constrained_dict["a"]

# ╔═╡ 7e1dbf59-df69-417a-8578-b34bb403848d
constrained_dict["b"] = "test"

# ╔═╡ 73055842-a630-4403-8233-f0a774258323
md"""
## Indexing
"""

# ╔═╡ 70339c78-0e43-404b-a7f3-b5e837915ebb
md"""Indexing a dictionary is made using `[` along with the key of choice. For example we would type `dict[key]` to get the value in `dict` stored in `key`:"""

# ╔═╡ 9ea3f3db-2421-4e43-a706-7b7f7c0b1959
dict["🍓"]

# ╔═╡ 7a00c43e-d7a8-40c0-bcfa-594b3f4e7161
md"""
Remember that dictionaries are not ordered, but we can still reproduce the behavior of an array just for fun using integers as keys:
"""

# ╔═╡ 26e9f223-690d-4063-8185-ef3ce0b5e622
array_dict = Dict(1 => "first value", 2 => "second value")

# ╔═╡ cca40559-57dd-4827-a334-cd437b6a475d
array_dict[2]

# ╔═╡ 18fc5ca7-febe-40a4-bfa8-ef5a8c0b547c
md"""
If we want to index a value that does not exist yet in the dictionary, it returns an error:
"""

# ╔═╡ ae2b5a5e-9a3a-49e4-8fff-cc5ef26ed173
array_dict[3]

# ╔═╡ 37c1ceae-e993-4b36-96df-c088d3ab006b
md"""
We can check if a key exists using `haskey`:
"""

# ╔═╡ a609fc58-3ade-41ac-8325-56a7a66e2a3b
haskey(array_dict,2)

# ╔═╡ 6e36762d-e0d6-45b1-958d-f0602bdd11a5
md"""
And if we're not sure if a key exist and want to return a default value if not, there's `get`:
"""

# ╔═╡ bee06854-8bbd-4c92-953d-46c18520707c
get(array_dict,1,"default value") # The key exists, so it returns its value

# ╔═╡ ddb11f15-b000-4276-9c06-30b9d1e9936c
get(array_dict,3,"default value") # The key does not exist, it returns the default one

# ╔═╡ 693fac53-b2b6-489b-a7aa-93ea5eab8165
md"""
If you want to get all values or all keys at once, use `keys` and `values` like for arrays:
"""

# ╔═╡ de829c74-37f4-4fe8-8e81-da3c729e2019
keys(dict)

# ╔═╡ b08a7e04-6106-4ff6-baa3-b153eb627eba
values(dict)

# ╔═╡ 296be6bc-0ee2-4d2b-b8e5-8b730a974162
md"""
!!! note
	As we can see `keys` and `values` return a `KeySet` and a `ValueIterator` respectively. But we can use `collect` to get them as arrays.
"""

# ╔═╡ 9f0b1f62-f9f2-45b8-8fda-fe02a99fa602
collect(values(dict))

# ╔═╡ f106e0b7-7b17-4564-8ba8-a9ba0c4c413a
md"""
## Filter
"""

# ╔═╡ 6a22a2b0-8c8d-4f52-99ae-d37bc2f1ef87
md"""
Filtering values is rather simple, just use the filer function 🙂. It works by providing a filtering function as the first argument:
"""

# ╔═╡ 4751e962-d173-4047-aeee-7c562d15ca8d
filter(x -> iseven(x.second), dict)

# ╔═╡ fae3a4ab-222a-4b4e-a617-9f8199194dc7
md"""
!!! note
	`x -> iseven(x.second)` is an anonymous function, we'll get into that in the future. What it does here is testing if the value of a `key=>value` pair is even or not, and returns `true` if so, and `false` if not. `filter` then filters out the values that are not even, *i.e.* the ones that return `false`.
"""

# ╔═╡ 6bbe9f7c-17b1-4566-89b4-6247fd937256
md"""
## Mutate
"""

# ╔═╡ d0371dc1-5f08-4f1d-a79f-f4bd954eb974
md"""
Dictionaries are easy to work with. We can add new elements using `push!`, merge two dictionaries using `merge` and delete elements using either `pop!` or `delete!`.

Let's see an example using the dict we previously declared. Here's a reminder of which elements it holds:
"""

# ╔═╡ 28ff01f9-681c-400b-b62b-e3524dcfe2c1
dict

# ╔═╡ b19e99d0-3729-4255-808b-4261fd68145a
md"""
#### Add elements
"""

# ╔═╡ c895188a-8b24-45d7-a862-18110c8b7fd1
md"""
Now let's add a new element:
"""

# ╔═╡ f7dbbd52-a2da-478f-89b0-e712f945103c
begin
	dict_tmp = copy(dict)
	push!(dict_tmp, "🍒" => 30)
end

# ╔═╡ 15ac5ad1-8843-4be6-9337-5c41fe8589cf
md"""

!!! note

	If we were outside of Pluto, we would simply use `push!(dict, "🍒" => 30)`. Pluto is designed around reactivity, so it is generally bad practice to mutate a variable inside a Pluto notebook. We would normally prefer assigning the results to a new variable. I use a trick here to still be able to show you an example with a mutating function (i.e. one that change the value of its input argument directly). The trick is to put the code in a `begin ... end` statement to make it as a single statement, and apply `push!` on a copy of `dict` (here `dict_tmp`) so `dict` is never changed. If it is too complicated, no worries, just consider this big blob of code just as if it was only written `push!(dict, "🍒" => 30)`.
"""

# ╔═╡ 9c7974c3-6e50-44d5-b5ec-cb535280e610
md"""
#### Merge dictionaries
"""

# ╔═╡ e04a56ef-7d7a-47da-9380-28e3fe9e4720
md"""
We can merge two dictionaries into a single one using `merge`:
"""

# ╔═╡ 8e132c33-2d5a-48d2-a6d7-c2d7fdc5b2b5
dict2 = Dict("🍆" => 2, "🍅" => 4)

# ╔═╡ 85fac12d-02f4-4598-884f-82a93b29dae0
merge(dict,dict2)

# ╔═╡ 8e27d819-d58a-42b9-bf63-e9ab6c01576e
md"""
#### Delete elements
"""

# ╔═╡ 93732551-1730-41da-ab86-dbaa3732dfa0
md"""
To delete an element from a dictionary, you can either use `pop!`:
"""

# ╔═╡ fda3a3b9-7e06-422c-9b46-668c923c53e7
begin
	dict_tmp2 = copy(dict)
	pop!(dict_tmp2, "🍌")
end

# ╔═╡ e6e1515d-a893-44e9-a5e3-c6e191dfaa2d
md"""
You can see it returns 4 here. That is because `pop!` returns the value of the element it just removed, and it is very useful sometimes! We can check the newly created `dict_tmp2` and see if there still is the banana element:
"""

# ╔═╡ 1d68302e-0fbf-4e60-8bd4-97d90c5ee0b1
dict_tmp2

# ╔═╡ 2b0361ec-b65b-426c-80a5-b5bfa95a6730
md"""
Nope, it was removed as expected.
"""

# ╔═╡ bab7580c-228e-4bd1-91c5-e476f1c97ad5
md"""
If you prefer to get the remaining elements instead of the one removed, you can use `delete!`:
"""

# ╔═╡ f62a1379-8afe-48a6-9466-a9235199cfdb
begin
	dict_tmp3 = copy(dict)
	delete!(dict_tmp3, "🍌")
end

# ╔═╡ cdd99268-f499-40ff-b9ad-f29bef3c05ed
dict_tmp3

# ╔═╡ a6a378c0-f3bc-4ce9-8e3d-29b3caf43ab5
md"""
## Ordered dictionaries
"""

# ╔═╡ 77496c2f-4e8f-44b6-bddb-1763ac179398
md"""
We saw earlier that dictionaries are unordered. But if you think unordered collections are too chaotic and you want to put an end to this madness, you can use ordered dictionaries from the `OrderedCollections.jl` package instead.

So what is an ordered dictionary? The only difference with a standard dictionary is that elements are ordered, meaning we assign them a position in the dictionary, and they stick to it. And if you add an element, it is added as the last element, as things should be. Think of it like an array with names for each element.
"""

# ╔═╡ ee62f664-1e9d-4d0d-bd89-c46210d0bf5e
begin
	order_dict = OrderedDict(dict_tmp)
	push!(order_dict, "🌲" => 1)
end

# ╔═╡ a75a0cee-1a9e-4186-a39b-6f7e1e3a9f89
md"""
!!! note
	If you have an error on the code above, it is certainly because you didn't activate the project before opening the pluto notebook, so the `OrderedCollections` package is not available to you. There's two solutions to this issue:

	1/ activate the environment before opening this notebook. To do so, enter the following command in Julia `] activate .`. Then instantiate the environment by entering `] instantiate`. This will work assuming you are working from the whole [julia_course](https://github.com/VEZY/julia_course) project.

	2/ Or add a new cell below with the following code: `import Pkg; Pkg.add("OrderedCollections")`.
"""

# ╔═╡ b3a97c75-07a8-404a-ae42-01a3e865af16
md"""
## Useful functions
"""

# ╔═╡ a43f99bc-9ead-40ab-8743-31923537b953
md"""
Union of two dictionaries:
"""

# ╔═╡ 603cd2c3-1437-437a-a931-caa893b5c05a
union(dict2, dict_tmp)

# ╔═╡ ab10a88a-8542-479b-8452-9612d25ef26e
md"""
Intersection of two dictionaries:
"""

# ╔═╡ 95aba192-6027-46d2-8cf8-6222ddbed0d7
intersect(dict, dict_tmp)

# ╔═╡ b0d12465-c1ba-4520-a0d4-321b86219845
md"""
Difference between two dictionaries:
"""

# ╔═╡ 11c0414b-ef75-4211-ad69-9e402f6061c9
setdiff(dict2, dict_tmp)

# ╔═╡ c3a5e91d-d7d0-4abe-b7e6-e0c08780657c
md"""
These functions are related to mathematical sets.

![](https://upload.wikimedia.org/wikipedia/commons/8/86/A_union_B.svg)[^1]

The inclusion-exclusion principle is used to calculate the size of the union of sets: the size of the union is the size of the two sets, minus the size of their intersection.

[^1]: See markdown link for source.  """

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
OrderedCollections = "~1.4.1"
PlutoUI = "~0.7.18"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "0ec322186e078db08ea3e7da5b8b2885c099b393"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.0"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "5efcf53d798efede8fee5b2c8b09284be359bf24"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.2"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "ae4bbcadb2906ccc085cf52ac286dc1377dceccc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "57312c7ecad39566319ccf5aa717a20788eb8c1f"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.18"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─d794dbb2-0254-4808-80ea-75f184586ad8
# ╟─3a242ea0-a518-11eb-0c88-d541f2d52fb8
# ╟─ea932952-cb4c-4876-9dc4-3a2b56cf4761
# ╠═61cb5c19-2c9b-4484-b95b-c54c2bc4066e
# ╟─94d9dbe3-4042-4bdf-9d14-7acc165514cc
# ╠═874a6524-8f3c-4bc0-bc6e-9e8fec402e56
# ╟─45835257-6034-4d64-be18-68c966d96603
# ╠═a0be3db2-a1a0-4405-b453-87b933d76ce2
# ╟─9f81f1e8-f3c2-4fa4-8469-2d605a418481
# ╟─de10d9d7-3669-4999-8309-f6c0452f5f10
# ╠═349fbaf2-1d6e-41c8-a05e-49a1f2d12ac1
# ╠═449eb41a-b5ac-422e-95c7-8baf36bb21f0
# ╟─ae9376e7-3fc0-4668-a495-4338e139385d
# ╠═bb28ae8d-5bc2-4056-9c75-9983362c6716
# ╟─602bc23f-1acc-44a5-8986-cadea9c7f461
# ╟─73dac266-f783-4406-91c7-27a08c4e4e60
# ╠═14ac3d30-c0bc-4543-8e61-ccb1bfcd9287
# ╠═e1856197-b467-465c-ba4c-032282c989eb
# ╠═83c526bb-5db2-4a30-bcfa-7476faaf9f55
# ╠═7e1dbf59-df69-417a-8578-b34bb403848d
# ╟─73055842-a630-4403-8233-f0a774258323
# ╟─70339c78-0e43-404b-a7f3-b5e837915ebb
# ╠═9ea3f3db-2421-4e43-a706-7b7f7c0b1959
# ╟─7a00c43e-d7a8-40c0-bcfa-594b3f4e7161
# ╠═26e9f223-690d-4063-8185-ef3ce0b5e622
# ╠═cca40559-57dd-4827-a334-cd437b6a475d
# ╟─18fc5ca7-febe-40a4-bfa8-ef5a8c0b547c
# ╠═ae2b5a5e-9a3a-49e4-8fff-cc5ef26ed173
# ╟─37c1ceae-e993-4b36-96df-c088d3ab006b
# ╠═a609fc58-3ade-41ac-8325-56a7a66e2a3b
# ╟─6e36762d-e0d6-45b1-958d-f0602bdd11a5
# ╠═bee06854-8bbd-4c92-953d-46c18520707c
# ╠═ddb11f15-b000-4276-9c06-30b9d1e9936c
# ╟─693fac53-b2b6-489b-a7aa-93ea5eab8165
# ╠═de829c74-37f4-4fe8-8e81-da3c729e2019
# ╠═b08a7e04-6106-4ff6-baa3-b153eb627eba
# ╟─296be6bc-0ee2-4d2b-b8e5-8b730a974162
# ╠═9f0b1f62-f9f2-45b8-8fda-fe02a99fa602
# ╟─f106e0b7-7b17-4564-8ba8-a9ba0c4c413a
# ╟─6a22a2b0-8c8d-4f52-99ae-d37bc2f1ef87
# ╠═4751e962-d173-4047-aeee-7c562d15ca8d
# ╟─fae3a4ab-222a-4b4e-a617-9f8199194dc7
# ╟─6bbe9f7c-17b1-4566-89b4-6247fd937256
# ╟─d0371dc1-5f08-4f1d-a79f-f4bd954eb974
# ╠═28ff01f9-681c-400b-b62b-e3524dcfe2c1
# ╟─b19e99d0-3729-4255-808b-4261fd68145a
# ╟─c895188a-8b24-45d7-a862-18110c8b7fd1
# ╠═f7dbbd52-a2da-478f-89b0-e712f945103c
# ╟─15ac5ad1-8843-4be6-9337-5c41fe8589cf
# ╟─9c7974c3-6e50-44d5-b5ec-cb535280e610
# ╟─e04a56ef-7d7a-47da-9380-28e3fe9e4720
# ╠═8e132c33-2d5a-48d2-a6d7-c2d7fdc5b2b5
# ╠═85fac12d-02f4-4598-884f-82a93b29dae0
# ╟─8e27d819-d58a-42b9-bf63-e9ab6c01576e
# ╟─93732551-1730-41da-ab86-dbaa3732dfa0
# ╠═fda3a3b9-7e06-422c-9b46-668c923c53e7
# ╟─e6e1515d-a893-44e9-a5e3-c6e191dfaa2d
# ╠═1d68302e-0fbf-4e60-8bd4-97d90c5ee0b1
# ╟─2b0361ec-b65b-426c-80a5-b5bfa95a6730
# ╟─bab7580c-228e-4bd1-91c5-e476f1c97ad5
# ╠═f62a1379-8afe-48a6-9466-a9235199cfdb
# ╠═cdd99268-f499-40ff-b9ad-f29bef3c05ed
# ╟─a6a378c0-f3bc-4ce9-8e3d-29b3caf43ab5
# ╟─77496c2f-4e8f-44b6-bddb-1763ac179398
# ╠═32785fd3-0adb-4b61-a8dc-c920904d66ea
# ╠═ee62f664-1e9d-4d0d-bd89-c46210d0bf5e
# ╟─a75a0cee-1a9e-4186-a39b-6f7e1e3a9f89
# ╟─b3a97c75-07a8-404a-ae42-01a3e865af16
# ╟─a43f99bc-9ead-40ab-8743-31923537b953
# ╠═603cd2c3-1437-437a-a931-caa893b5c05a
# ╟─ab10a88a-8542-479b-8452-9612d25ef26e
# ╠═95aba192-6027-46d2-8cf8-6222ddbed0d7
# ╟─b0d12465-c1ba-4520-a0d4-321b86219845
# ╠═11c0414b-ef75-4211-ad69-9e402f6061c9
# ╟─c3a5e91d-d7d0-4abe-b7e6-e0c08780657c
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "libtiff"
version = v"4.0.9"

# Collection of sources required to build libtiff
sources = [
    "https://gitlab.com/libtiff/libtiff.git" =>
    "020bd2fd3b4235b2175dafa9633d2b191c50e7cd",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd libtiff/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain
make
make install
exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libtiff", :libtiff),
    LibraryProduct(prefix, "libtiffxx", :libtiffxx)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.2/build_Zlib.v1.2.11.jl",
    "https://github.com/SimonDanisch/LibpngBuilder/releases/download/1.0.0/build_libpng.v1.6.31.jl",
    "https://github.com/SimonDanisch/LibJPEGBuilder/releases/download/v9b/build_libjpeg.v9.0.0-b.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
"""
A LatticeBond is a struct which represents
a single bond in a geometrical lattice or
else on interaction graph defining a physical
model such as a quantum Hamiltonian.

LatticeBond has the following data fields:
* s1::Int -- number of site 1
* s2::Int -- number of site 2
* x1::Int -- x coordinate of site 1
* y1::Int -- y coordinate of site 1
* x2::Int -- x coordinate of site 2
* y2::Int -- y coordinate of site 2
* type::String -- optional description of bond type
"""
struct LatticeBond
  s1::Int
  s2::Int
  x1::Float64
  y1::Float64
  x2::Float64
  y2::Float64
  type::String
end

"""
    LatticeBond(s1::Int,s2::Int)

    LatticeBond(s1::Int,s2::Int,
                x1::Real,y1::Real,
                x2::Real,y2::Real,
                type::String="")

Construct a LatticeBond struct by
specifying just the numbers of sites
1 and 2, or additional details including
the (x,y) coordinates of the two sites and
an optional type string.
"""
function LatticeBond(s1::Int, s2::Int)
  return LatticeBond(s1, s2, 0.0, 0.0, 0.0, 0.0, "")
end

function LatticeBond(
  s1::Int, s2::Int, x1::Real, y1::Real, x2::Real, y2::Real, bondtype::String=""
)
  cf(x) = convert(Float64, x)
  return LatticeBond(s1, s2, cf(x1), cf(y1), cf(x2), cf(y2), bondtype)
end

"""
Lattice is an alias for Vector{LatticeBond}
"""
const Lattice = Vector{LatticeBond}

"""
    square_lattice(Nx::Int,
                   Ny::Int;
                   kwargs...)::Lattice

Return a Lattice (array of LatticeBond
objects) corresponding to the two-dimensional
square lattice of dimensions (Nx,Ny).
By default the lattice has open boundaries,
but can be made periodic in the y direction
by specifying the keyword argument
`yperiodic=true`.
"""
function square_lattice(Nx::Int, Ny::Int; kwargs...)::Lattice
  yperiodic = get(kwargs, :yperiodic, false)
  yperiodic = yperiodic && (Ny > 2)
  N = Nx * Ny
  Nbond = 2N - Ny + (yperiodic ? 0 : -Nx)
  latt = Lattice(undef, Nbond)
  b = 0
  for n in 1:N
    x = div(n - 1, Ny) + 1
    y = mod(n - 1, Ny) + 1
    if x < Nx
      latt[b += 1] = LatticeBond(n, n + Ny, x, y, x + 1, y)
    end
    if Ny > 1
      if y < Ny
        latt[b += 1] = LatticeBond(n, n + 1, x, y, x, y + 1)
      end
      if yperiodic && y == 1
        latt[b += 1] = LatticeBond(n, n + Ny - 1, x, y, x, y + Ny)
      end
    end
  end
  return latt
end

"""
    triangular_lattice(Nx::Int,
                       Ny::Int;
                       kwargs...)::Lattice

Return a Lattice (array of LatticeBond
objects) corresponding to the two-dimensional
triangular lattice of dimensions (Nx,Ny).
By default the lattice has open boundaries,
but can be made periodic in the y direction
by specifying the keyword argument
`yperiodic=true`.
"""
function triangular_lattice(Nx::Int, Ny::Int; kwargs...)::Lattice
  yperiodic = get(kwargs, :yperiodic, false)
  yperiodic = yperiodic && (Ny > 2)
  N = Nx * Ny
  Nbond = 3N - 2Ny + (yperiodic ? 0 : -2Nx + 1)
  latt = Lattice(undef, Nbond)
  b = 0
  for n in 1:N
    x = div(n - 1, Ny) + 1
    y = mod(n - 1, Ny) + 1

    # x-direction bonds
    if x < Nx
      latt[b += 1] = LatticeBond(n, n + Ny)
    end

    # 2d bonds
    if Ny > 1
      # vertical / y-periodic diagonal bond
      if (n + 1 <= N) && ((y < Ny) || yperiodic)
        latt[b += 1] = LatticeBond(n, n + 1)
      end
      # periodic vertical bond
      if yperiodic && y == 1
        latt[b += 1] = LatticeBond(n, n + Ny - 1)
      end
      # diagonal bonds
      if x < Nx && y < Ny
        latt[b += 1] = LatticeBond(n, n + Ny + 1)
      end
    end
  end
  return latt
end


text:
# ---
# title: 1357. Apply Discount Every n Orders
# id: problem1357
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Design
# link: <https://leetcode.com/problems/apply-discount-every-n-orders/description/>
# hidden: true
# ---
#
# There is a sale in a supermarket, there will be a `discount` every `n`
# customer.
# There are some products in the supermarket where the id of the `i-th` product
# is `products[i]` and the price per unit of this product is `prices[i]`.
# The system will count the number of customers and when the `n-th` customer
# arrive he/she will have a `discount` on the bill. (i.e if the cost is `x` the
# new cost is `x - (discount * x) / 100`). Then the system will start counting
# customers again.
# The customer orders a certain amount of each product where `product[i]` is the
# id of the `i-th` product the customer ordered and `amount[i]` is the number of
# units the customer ordered of that product.
#
# Implement the `Cashier` class:
#
#   * `Cashier(int n, int discount, int[] products, int[] prices)` Initializes the object with `n`, the `discount`, the `products` and their `prices`.
#   * `double getBill(int[] product, int[] amount)` returns the value of the bill and apply the discount if needed. Answers within `10^-5` of the actual value will be accepted as correct.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["Cashier","getBill","getBill","getBill","getBill","getBill","getBill","getBill"]
#     [[3,50,[1,2,3,4,5,6,7],[100,200,300,400,300,200,100]],[[1,2],[1,2]],[[3,7],[10,10]],[[1,2,3,4,5,6,7],[1,1,1,1,1,1,1]],[[4],[10]],[[7,3],[10,10]],[[7,5,3,1,6,4,2],[10,10,10,9,9,9,7]],[[2,3,5],[5,3,2]]]
#     **Output**
#     [null,500.0,4000.0,800.0,4000.0,4000.0,7350.0,2500.0]
#     **Explanation**
#     Cashier cashier = new Cashier(3,50,[1,2,3,4,5,6,7],[100,200,300,400,300,200,100]);
#     cashier.getBill([1,2],[1,2]);                        // return 500.0, bill = 1 * 100 + 2 * 200 = 500.
#     cashier.getBill([3,7],[10,10]);                      // return 4000.0
#     cashier.getBill([1,2,3,4,5,6,7],[1,1,1,1,1,1,1]);    // return 800.0, The bill was 1600.0 but as this is the third customer, he has a discount of 50% which means his bill is only 1600 - 1600 * (50 / 100) = 800.
#     cashier.getBill([4],[10]);                           // return 4000.0
#     cashier.getBill([7,3],[10,10]);                      // return 4000.0
#     cashier.getBill([7,5,3,1,6,4,2],[10,10,10,9,9,9,7]); // return 7350.0, Bill was 14700.0 but as the system counted three more customers, he will have a 50% discount and the bill becomes 7350.0
#     cashier.getBill([2,3,5],[5,3,2]);                    // return 2500.0
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 10^4`
#   * `0 <= discount <= 100`
#   * `1 <= products.length <= 200`
#   * `1 <= products[i] <= 200`
#   * There are **not** repeated elements in the array `products`.
#   * `prices.length == products.length`
#   * `1 <= prices[i] <= 1000`
#   * `1 <= product.length <= products.length`
#   * `product[i]` exists in `products`.
#   * `amount.length == product.length`
#   * `1 <= amount[i] <= 1000`
#   * At most `1000` calls will be made to `getBill`.
#   * Answers within `10^-5` of the actual value will be accepted as correct.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using EnsembleInference

@testset "Utility functions" begin
    @testset "kron_delta" begin
        @test EnsembleInference.kron_delta(0, 0)
        @test EnsembleInference.kron_delta(1, 1)
        @test EnsembleInference.kron_delta(2, 2)
        @test !EnsembleInference.kron_delta(0, 1)
        @test !EnsembleInference.kron_delta(1, 0)
        @test !EnsembleInference.kron_delta(0, 2)
        @test !EnsembleInference.kron_delta(2, 0)
        @test !EnsembleInference.kron_delta(1, 2)
        @test !EnsembleInference.kron_delta(2, 1)
    end
end


text:
@doc raw"""
    radial_dipole_moment(a, b)

Compute the radial dipole moment between the orbitals `a` and `b` according to

```math
\begin{align}
\tag{BS63.2}
R_{n \ell}^{n',\ell-1} &=
ABCD,
&
n \neq n'\\
\tag{BS63.5}
R_{n,\ell-1}^{n\ell} &=
R_{n\ell}^{n,\ell-1} =
\frac{3}{2}
n\sqrt{n^2-\ell^2},
& n = n',
\end{align}
```
where
```math
\begin{aligned}
A &= \frac{(-)^{n'-\ell}}{4(2\ell-1)!}, \\
B &= \sqrt{\frac{(n+\ell)!(n'+\ell-1)!}{(n-\ell-1)!(n'-\ell)!}}, \\
C &= \frac{(4nn')^{\ell+1}(n-n')^{n+n'-2\ell-2}}{(n+n')^{n+n'}}, \\
D &= {_2}F_1(-n_r, -n_r', 2\ell, x) -
\left(\frac{n-n'}{n+n'}\right)^2
{_2}F_1(-n_r-2, -n_r', 2\ell, x), \\
n_r &= n-\ell-1, \quad
n_r' = n' - \ell, \quad
x = -\frac{4nn'}{(n-n')^2},
\end{aligned}
```
and ``{_2}F_1(a,b,c;x)`` is a hypergeometric function that can be
computed used
[HyperGeometricFunctions.jl](https://github.com/JuliaMath/HypergeometricFunctions.jl).

"""
function radial_dipole_moment(a::Orbital, b::Orbital)
    a == b && return 0.0

    n,ℓ = a.n,a.ℓ
    n′,ℓ′ = b.n,b.ℓ

    abs(ℓ-ℓ′) ≠ 1 && return 0.0
    ℓ < ℓ′ && return radial_dipole_moment(b, a)

    # Eq. (BS63.5)
    n == n′ && return 3/2*n*√(n^2-ℓ^2)

    A = powneg1(n′-ℓ)/(4*factorial(2ℓ-1))
    B = √(factorial_ratio(n+ℓ,n-ℓ-1)*factorial_ratio(n′+ℓ-1,n′-ℓ))
    C = (4*n*n′)^(ℓ+1)*(n-n′)^(n+n′-2ℓ-2)/((n+n′)^(n+n′))

    nᵣ = n-ℓ-1
    nᵣ′ = n′-ℓ
    x = -(4n*n′)/((n-n′)^2)
    D = _₂F₁(-nᵣ, -nᵣ′, 2ℓ, x) - ((n-n′)/(n+n′))^2*_₂F₁(-nᵣ-2, -nᵣ′, 2ℓ, x)

    A*B*C*D
end

radial_dipole_moment(a::SpinOrbital, b::SpinOrbital) =
    radial_dipole_moment(a.orb, b.orb)

"""
    dipole_moment(a, b[, component=:z])

Compute the dipole moment between orbitals `a` and `b` along the
Cartesian direction indicated by `component`. The radial integral is
computed using [`radial_dipole_moment`](@ref) and the angular integral
using
[AngularMomentumAlgebra.jl](https://github.com/JuliaAtoms/AngularMomentumAlgebra.jl).

```jldoctest
julia> a,b,c,d = SpinOrbital(o"1s", 0, half(1)), SpinOrbital(o"2s", 0, half(1)), SpinOrbital(o"2p", 0, half(1)), SpinOrbital(o"2p", 1, half(1))
(1s₀α, 2s₀α, 2p₀α, 2p₁α)

julia> dipole_moment(a, c)
0.7449355390278029

julia> dipole_moment(b, c)
2.999999999999999

julia> dipole_moment(a, d)
0.0

julia> dipole_moment(a, d, :x)
-0.5267489711934153
```

"""
function dipole_moment(a::SpinOrbital, b::SpinOrbital, component=:z)
    d = AngularMomentumAlgebra.Dipoles.𝐫[(x = 1, y = 2, z = 3)[component]]
    v = dot(a,d,b)
    iszero(v) && return 0.0
    r = radial_dipole_moment(a, b)
    r*v.operators[1][2]
end

"""
    dipole_matrix(orbitals[, component=:z])

Return the sparse matrix whose elements correspond to the dipole
moments between all of the `orbitals` as measured along the Cartesian
direction indicated by `component`.

# Examples

```jldoctest
julia> dipole_matrix(sos"1[s] 2[s-p]")
10×10 SparseArrays.SparseMatrixCSC{Float64,Int64} with 8 stored entries:
  [7, 1]  =  0.744936
  [8, 2]  =  0.744936
  [7, 3]  =  3.0
  [8, 4]  =  3.0
  [1, 7]  =  0.744936
  [3, 7]  =  3.0
  [2, 8]  =  0.744936
  [4, 8]  =  3.0

julia> dipole_matrix(sos"1[s] 2[s-p]", :y)
10×10 SparseArrays.SparseMatrixCSC{Complex{Float64},Int64} with 16 stored entries:
  [5 ,  1]  =  0.0+0.526749im
  [9 ,  1]  =  0.0+0.526749im
  [6 ,  2]  =  0.0+0.526749im
  [10,  2]  =  0.0+0.526749im
  [5 ,  3]  =  0.0+2.12132im
  [9 ,  3]  =  0.0+2.12132im
  [6 ,  4]  =  0.0+2.12132im
  [10,  4]  =  0.0+2.12132im
  [1 ,  5]  =  -0.0-0.526749im
  [3 ,  5]  =  -0.0-2.12132im
  [2 ,  6]  =  -0.0-0.526749im
  [4 ,  6]  =  -0.0-2.12132im
  [1 ,  9]  =  -0.0-0.526749im
  [3 ,  9]  =  -0.0-2.12132im
  [2 , 10]  =  -0.0-0.526749im
  [4 , 10]  =  -0.0-2.12132im
```
"""
function dipole_matrix(orbitals::AbstractVector{<:SpinOrbital}, component=:z)
    d = AngularMomentumAlgebra.Dipoles.𝐫[(x = 1, y = 2, z = 3)[component]]

    T = component == :y ? ComplexF64 : Float64
    norb = length(orbitals)

    sparse_builder(T, norb, norb) do set!
        for (i,a) in enumerate(orbitals)
            for (j,b) in enumerate(orbitals)
                v = dot(a,d,b)
                iszero(v) && continue
                r = radial_dipole_moment(a.orb, b.orb)
                set!(i, j, r*v.operators[1][2])
            end
        end
    end
end

export radial_dipole_moment, dipole_moment, dipole_matrix


text:
module Indices
using ..AbstractLattice

export index, gridcoords, neighbor_index
import Base: length, size

# Base type for all indexing kernels
abstract Indexing
typealias GridCoords{I <: Integer} Vector{I}

include("cartesian.jl")
include("periodic.jl")

length(d::Indexing) = prod(d.dimensions)
size(d::Indexing) = d.dimensions
function size(d::Indexing, i::Integer)
    @assert i ≥ 1 && i ≤ length(d.dimensions)
    d.dimensions[i]
end

function index{N, I <: Integer}(kernel::Indexing, indices::NTuple{N, I})
    index(kernel::Indexing, collect(indices))
end

function gridcoords(kernel::Indexing, index::Integer)
    result = ones(typeof(index), length(kernel.dimensions))
    index -= 1
    for i in length(kernel.dimensions):-1:1
       result[i] += floor(typeof(index), index // kernel.strides[i])
       index %= kernel.strides[i]
    end
    result
end

function neighbor_index(indexing::Indexing, site::GridCoords, direction::GridCoords)
    @assert length(site) == length(direction)
    const neighbor = site + direction
    @assert all(neighbor .> 0) && all(neighbor .<= size(indexing))
    index(indexing, neighbor)
end
function neighbor_index(indexing::Indexing, site::Integer, direction::GridCoords)
    @assert site > 0 && site < length(indexing)
    neighbor_index(indexing, gridcoords(indexing, site), direction)
end

include("cached.jl")

end


text:
using Gadfly, DataArrays, RDatasets

plot(dataset("Zelig", "macro"), x="Year", y="Country", color="GDP", Geom.rectbin)


text:
function getUserInput(T = String)
    if T == String
        return readline()
    else
    try
        return parse(T,readline())
    catch
        println("Sorry, I could not interpret your answer. Please try again")
        getUserInput(T,msg)
        end
    end
end

function brack(s, n, f)
    x = [["(", "[", "{", "<"], [")", "]", "}", ">"]];
    for i = 1:4
        if s == x[n][i]
            if f
                return s
            else
                return i
            end
        end
    end
    return " "
end

function brackets(s)
    x = String[];
    y = 0;
    for i = 1:length(s)
        a = brack(s[i], 1, true);
        if a != " "
            resize!(x, y + 1);
            x[y + 1] = a;
            y += 1;
        else
            b = brack(s[i], 2, true);
            if b != " "
                if y == 0
                    return "Brackets are uneven!"
                elseif brack(x[y], 1, false) == brack(s[i], 2, false)
                    deleteat!(x, length(x));
                    y -= 1;
                else
                    return "Brackets are uneven!"
                end
            end
        end
    end
    if y != 0
        return "Brackets are uneven!"
    else
        return "Brackets are even!"
    end
end

println(brackets(split(getUserInput(String), "")));

max_stars_repo_name: MediocreBoris/gitgud

max_issues_repo_name: MediocreBoris/gitgud

max_forks_repo_name: MediocreBoris/gitgud


text:
using Revise
using Shearlab
using ScatteringTransform
using Plots

using cartoonLike
using LinearAlgebra

# compare the two transforms on a random matrix
m1=2; n1=50; p1=55
X1 = randn(n1,p1)
layer1 = stParallel(m1, size(X1), typeBecomes=eltype(X1))
scatteringSheared = ScatteringTransform.sheardec2D(X1, layer1.shears[1])

shears = Shearlab.getshearletsystem2D(size(X1)..., 2)
tinyShears = Shearlab.getshearletsystem2D(50,55,2)
tinyOutput = Shearlab.SLsheardec2D(X1,tinyShears)
norm((scatteringSheared)/norm(scatteringSheared,Inf) - tinyOutput/norm(tinyOutput,Inf))
plot(heatmap(scatteringSheared[:,:,1]), heatmap(real.(tinyOutput)[:,:,1]), heatmap(real.(scatteringSheared - tinyOutput)[:,:,1]), layout=(3,1))

n = 400
points = [-.25 0;
          -.26 -.3;
           .1 -.25;
           .3  .3;
          -.3  .25]
X1,xLocs,yLocs = given_points(points, cart = cartoonLikeImage{Float64}(xSize=n, ySize=n))
heatmap(xLocs, yLocs, X1)
layer1 = stParallel(m1, size(X1), typeBecomes=eltype(X1))
scatteringSheared = ScatteringTransform.sheardec2D(X1, layer1.shears[1])

shears = Shearlab.getshearletsystem2D(size(X1)..., 2)
tinyShears = Shearlab.getshearletsystem2D(400, 400, 2)
tinyOutput = Shearlab.SLsheardec2D(X1, tinyShears)

norm((scatteringSheared)/norm(scatteringSheared,Inf) - tinyOutput/norm(tinyOutput,Inf))
# it agrees very closely on a more realistic image that isn't just gaussian noise everywhere.
k = 7; plot(heatmap(scatteringSheared[:,:,k]), heatmap(real.(tinyOutput)[:,:,k]), heatmap(real.(scatteringSheared - tinyOutput)[:,:,k]), layout=(3,1))

subsamp = 1.5; nScales = 3
m2=2; n2=201; p2=325
X2 = randn(n2,p2)
layer2 = stParallel(m2, size(X2); subsample = subsamp, nScale = nScales,
                          typeBecomes = eltype(X2))
# A more carefully constructed test that also tests a different element type
m3=2
X3 = zeros(Float32, 1, 50, 50)
X3[1, 13:37, 13:37] = ones(Float32, 25, 25)
layer3 = stParallel(m3, size(X3)[end-1:end], typeBecomes=eltype(X3))


text:
module Graphs

abstract type AbstractGraph{T} end

mutable struct DirectedGraph{T} <: AbstractGraph{T}
    vertices::Vector{T}
    edges::Matrix{Bool}
    vertices_index::Dict{T, Int}
end

function Base.show(io::IO, g::DirectedGraph)
    println(io, "DirectedGraph{$(eltype(g))}")

    isempty(g) && return

    println(io, "Vertices")
    for (i, v) in enumerate(g.vertices)
        println(io, "[$i] $v")
    end

    println(io, "Edges")
    println(io, g.edges)
end

function Base.eltype(g::DirectedGraph{T}) where {T}
    T
end

function Base.isempty(g::DirectedGraph) :: Bool
    return nv(g) == 0
end

function DirectedGraph(::Type{T}) where {T}
    DirectedGraph(Vector{T}(), Matrix{Bool}(undef, 0, 0), Dict{T, Int}())
end

"Number of Edges"
ne(g::DirectedGraph) = sum(g.edges)

"Number of Vertices"
nv(g::DirectedGraph) = length(g.vertices)

function has_vertex(g::DirectedGraph{T}, v::T) where {T}
    return haskey(g.vertices_index, v)
end

function get_vertex_index(g::DirectedGraph{T}, v::T) where {T}
    return g.vertices_index[v]
end

function get_vertex_from_index(g::DirectedGraph, ind::Integer)
    return g.vertices[ind]
end

function add_vertex!(g::DirectedGraph{T}, v::T) where {T}
    @assert !has_vertex(g, v) "Vertex $v is already in graph"
    push!(g.vertices, v)
    vertex_index = length(g.vertices)
    g.vertices_index[v] = vertex_index
    augment_edges_matrix!(g, vertex_index)

    nothing
end

function augment_edges_matrix!(g::DirectedGraph, vertex_index::Integer)
    @assert size(g.edges) == (vertex_index - 1, vertex_index - 1)
    new_edges = similar(g.edges, vertex_index, vertex_index)

    # copies old content
    for j in 1:(vertex_index - 1), i in 1:(vertex_index - 1)
        new_edges[i, j] = g.edges[i, j]
    end

    # reset new content
    for i in 1:vertex_index
        new_edges[i, vertex_index] = false
        new_edges[vertex_index, i] = false
    end

    # point to the new matrix
    g.edges = new_edges

    nothing
end

function add_edge!(g::DirectedGraph, source::Integer, dest::Integer)
    g.edges[source, dest] = true
    nothing
end

function add_edge!(g::DirectedGraph{T}, source::T, dest::T) where {T}
    add_edge!(g, get_vertex_index(g, source), get_vertex_index(g, dest))
end

has_edge(g::DirectedGraph, source::Integer, dest::Integer) :: Bool = g.edges[source, dest]

function has_edge(g::DirectedGraph{T}, source::T, dest::T) :: Bool where {T}
    return has_edge(g, get_vertex_index(g, source), get_vertex_index(g, dest))
end

function remove_edge!(g::DirectedGraph, source::Integer, dest::Integer)
    g.edges[source, dest] = false
    nothing
end

function remove_edge!(g::DirectedGraph{T}, source::T, dest::T) where {T}
    remove_edge!(g, get_vertex_index(g, source), get_vertex_index(g, dest))
end

abstract type GraphIteratorResult end

struct ReturnIndex <: GraphIteratorResult
end

struct ReturnVertex <: GraphIteratorResult
end

struct GraphIterator{R<:GraphIteratorResult, G<:AbstractGraph, T}
    graph::G
    edges_view::T
    result::R
end

Base.length(iter::GraphIterator) = sum(iter.edges_view)
Base.eltype(iter::GraphIterator{ReturnIndex}) = Int
Base.eltype(iter::GraphIterator{ReturnVertex, G}) where {T, G<:AbstractGraph{T}} = T

function view_on_edges(g::DirectedGraph, vertex_index::Integer, fwd::Bool)
    if fwd
        return view(g.edges, vertex_index, :)
    else
        return view(g.edges, :, vertex_index)
    end
end

function GraphIterator(g::AbstractGraph, vertex_index::Integer, fwd::Bool, result::R) where {R<:GraphIteratorResult}
    GraphIterator(
            g,
            view_on_edges(g, vertex_index, fwd),
            result
        )
end

each_adjacent_vertex_index(g::AbstractGraph, vertex_index::Integer, fwd::Bool) = GraphIterator(g, vertex_index, fwd, ReturnIndex())
each_adjacent_vertex(g::AbstractGraph, vertex_index::Integer, fwd::Bool) = GraphIterator(g, vertex_index, fwd, ReturnVertex())
adjacent_vertex_count(g::AbstractGraph, vertex_index::Integer, fwd::Bool) = length(each_adjacent_vertex_index(g, vertex_index, fwd))

function adjacent_vertex_count(g::AbstractGraph{T}, vertex::T, fwd::Bool) where {T}
    return adjacent_vertex_count(g, get_vertex_index(g, vertex), fwd)
end

predecessors_count(g::AbstractGraph, vertex) = adjacent_vertex_count(g, vertex, false)
successors_count(g::AbstractGraph, vertex) = adjacent_vertex_count(g, vertex, true)

function Base.iterate(iter::GraphIterator{R, G, T}, state=1) where {R, G, T}

    if state > length(iter.edges_view)
        return nothing
    end

    for vertex_index in state:length(iter.edges_view)
        if iter.edges_view[vertex_index]
            # edge exists
            next_vertex_index = vertex_index + 1
            if R == ReturnIndex
                return (vertex_index, next_vertex_index)
            else
                return (get_vertex_from_index(iter.graph, vertex_index), next_vertex_index)
            end
        end
    end

    return nothing
end

function each_successor(g::AbstractGraph{T}, from_vertex::T) where {T}
    return each_successor(g, get_vertex_index(g, from_vertex))
end

function each_successor(g::AbstractGraph, from_vertex_index::Integer)
    return each_adjacent_vertex(g, from_vertex_index, true)
end

function each_predecessor(g::AbstractGraph{T}, from_vertex::T) where {T}
    return each_predecessor(g, get_vertex_index(g, from_vertex))
end

function each_predecessor(g::AbstractGraph, from_vertex_index::Integer)
    return each_adjacent_vertex(g, from_vertex_index, false)
end

successors(g::AbstractGraph, from_vertex) = Set(each_successor(g, from_vertex))
predecessors(g::AbstractGraph, from_vertex) = Set(each_predecessor(g, from_vertex))

#
# BFS
#

mutable struct BFSVertexInfo
    color::Symbol
    distance::Int
    pred::Union{Nothing, Int} # predecessor vertex index
end

"Breadth-First Search"
function bfs(g::AbstractGraph, source_vertex_index::Integer, fwd::Bool) :: Vector{BFSVertexInfo}
    vertices_count = nv(g)

    # setup
    bfs_vertices = Vector{BFSVertexInfo}(undef, vertices_count)
    for i in 1:vertices_count
        if i == source_vertex_index
            bfs_vertices[i] = BFSVertexInfo(:gray, 0, nothing)
        else
            bfs_vertices[i] = BFSVertexInfo(:white, typemax(Int), nothing)
        end
    end
    queue = Vector{Int}()
    push!(queue, source_vertex_index)

    # loop
    while !isempty(queue)
        vertex_index = pop!(queue)

        for adjacent_vertex_index in each_adjacent_vertex_index(g, vertex_index, fwd)
            bfs_vertex = bfs_vertices[adjacent_vertex_index]
            if bfs_vertex.color == :white
                bfs_vertex.color = :gray
                bfs_vertex.distance = bfs_vertices[vertex_index].distance + 1
                bfs_vertex.pred = vertex_index
                push!(queue, adjacent_vertex_index)
            end
            bfs_vertices[vertex_index].color = :black
        end
    end

    return bfs_vertices
end

function bfs(g::AbstractGraph{T}, source_vertex::T, fwd::Bool) where {T}
    return bfs(g, get_vertex_index(g, source_vertex), fwd)
end

"Returns the Set of Vertices that are reachable from `source_vertex`."
function reachable_vertices(g::AbstractGraph{T}, source_vertex, fwd::Bool) where {T}
    result = Set{T}()

    for (vertex_index, bfs_vertex) in enumerate(bfs(g, source_vertex, fwd))
        if is_reachable(bfs_vertex)
            push!(result, get_vertex_from_index(g, vertex_index))
        end
    end

    return result
end

is_reachable(v::BFSVertexInfo) = v.color != :white

#
# Dominance
#

"""
In a flow graph with entry node V1,
node Vi dominates node Vj if and only if
Vi lies on every path from V1 to Vj.
"""
function dominator_sets(g::AbstractGraph)

    local dom::Vector{Set{Int}}
    vertices_count = nv(g)

    if vertices_count == 0
        error("Can't execute dominance analysis on an empty graph")
    elseif vertices_count == 1
        return Set([1])
    end

    let
        # initial state: Dom(n0) = {n0}, Dom(n) = N, n != n0, N = all vertices
        all_vertices = Set([ i for i in 1:vertices_count ])
        dom = [ (i == 1 ? Set([1]) : copy(all_vertices)) for i in 1:vertices_count ]
    end

    @assert isempty(predecessors(g, 1)) "First vertex $(get_vertex_from_index(1)) can't have a predecessor."

    local changed::Bool = true

    while changed
        changed = false

        for i in 2:vertices_count
            temp = ∪(Set([i]), ∩([ dom[j] for j in each_adjacent_vertex_index(g, i, false) ]...))
            if temp != dom[i]
                dom[i] = temp
                changed = true
            end
        end
    end

    return dom
end

function strict_dominator_sets(g::AbstractGraph)
    sets = dominator_sets(g)
    for (i, set) in enumerate(sets)
        pop!(set, i)
    end
    return sets
end

"""
Returns the index of the nearest reachable vertex from a given `set`
Throws error if there is a vertex in `set` that is not reachable from `vertex_index`.
"""
function nearest_reachable_vertex_from_set(g::AbstractGraph, vertex_index::Integer, set::Set{Int}, fwd::Bool) :: Int
    bfs_vertices = bfs(g, vertex_index, fwd)

    distance_to_vertex = Dict{Int, Int}()
    for vertex in set
        bfs_vertex = bfs_vertices[vertex]
        @assert is_reachable(bfs_vertex) "vertex $vertex is not reachable."
        distance_to_vertex[bfs_vertex.distance] = vertex
    end

    nearest_distance = min(keys(distance_to_vertex)...)
    return distance_to_vertex[nearest_distance]
end

function idominators(g::AbstractGraph)
    result = Vector{Union{Nothing, Int}}(undef, nv(g))

    for (i, set) in enumerate(strict_dominator_sets(g))
        if isempty(set)
            result[i] = nothing
        elseif length(set) == 1
            result[i] = pop!(set)
        else
            result[i] = nearest_reachable_vertex_from_set(g, i, set, false)
        end
    end

    return result
end

function dominance_frontiers(g::AbstractGraph)
    idoms = idominators(g)
    vertices_count = nv(g)

    local runner::Int

    # initial state -> empty sets for all vertices
    result = [ Set{Int}() for i in 1:vertices_count ]

    for vertex_index in 1:vertices_count
        if predecessors_count(g, vertex_index) > 1
            for pred_vertex_index in each_adjacent_vertex_index(g, vertex_index, false)
                runner = pred_vertex_index
                while runner != idoms[vertex_index]
                    push!(result[runner], vertex_index)
                    runner = idoms[runner]
                end
            end
        end
    end

    return result
end

end


text:
function VerticesToCentroids!( kc::Matrix{Float64}, kv::Matrix{Float64} )
    @tturbo for j=1:size(kc,2)
        for i=1:size(kc,1)
            kc[i,j] = 0.25*(kv[i,j] + kv[i+1,j] + kv[i,j+1] + kv[i+1,j+1])
        end
    end
end
export VerticesToCentroids!

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

function CentroidsToVertices!( kv::Matrix{Float64}, kc::Matrix{Float64}, ncx::Int64, ncy::Int64, BC::BoundaryConditions )

    @tturbo kv .= 0.0

    # First make extended arry using BC's
    kc_ex   = zeros( ncx+2, ncy+2)
    @tturbo kc_ex[2:end-1,2:end-1] .= kc
    if BC.periodix==1 # Periodic case
        kc_ex[  1,2:end-1] .= kc[end,:]
        kc_ex[end,2:end-1] .= kc[  1,:]
    else # other cases
        kc_ex[  1,2:end-1] .= kc[  1,:]
        kc_ex[end,2:end-1] .= kc[end,:]
    end
    kc_ex[  :,  1] .= kc_ex[  :,    2]
    kc_ex[  :,end] .= kc_ex[  :,end-1]

    # Conductivity averaging
    VerticesToCentroids!( kv, kc_ex )
end
export CentroidsToVertices!

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

function AverageConductivity( kc::Matrix{Float64}, ncx::Int64, ncy::Int64, Th_BC::Thermal_BC )

    # Security
    if Th_BC.type_W==0 && Th_BC.type_E!=0
        error("Both sides should be periodic!")
    end
    kx        = zeros( ncx+1, ncy  )
    ky        = zeros( ncx  , ncy+1)
    kc_ex     = zeros( ncx+2, ncy  )
    kc_ex[2:end-1,:] .= kc
    if Th_BC.type_W==0 && Th_BC.type_E==0 # Periodic case
        kc_ex[  1,:] .= kc[end,:]
        kc_ex[end,:] .= kc[  1,:]
    else # other cases
        kc_ex[  1,:] .= kc[  1,:]
        kc_ex[end,:] .= kc[end,:]
    end

    # Conductivity averaging x
    @tturbo for j=1:ncy
        for i=1:ncx+1
            kx[i,j] += 0.5*(kc_ex[i,j] + kc_ex[i+1,j])
        end
    end

    kc_ex     = zeros( ncx, ncy+2)
    kc_ex[:,2:end-1] .= kc
    kc_ex[:,  1] .= kc[:,  1]
    kc_ex[:,end] .= kc[:,end]

    # Conductivity averaging y
    @tturbo for i=1:ncx
        for j=1:ncy+1
            ky[i,j] += 0.5*(kc_ex[i,j] + kc_ex[i,j+1])
        end
    end

    return kx, ky
end
export AverageConductivity

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

@views function GenerateMesh( xmin::Float64, xmax::Float64, ymin::Float64, ymax::Float64, ncx::Int64, ncy::Int64 )
    dx, dy    = (xmax-xmin)/ncx, (ymax-ymin)/ncy
    xc        = LinRange( xmin+dx/2, xmax-dx/2, ncx  )
    yc        = LinRange( ymin+dy/2, ymax-dy/2, ncy  )
    xce       = LinRange( xmin-dx/2, xmax+dx/2, ncx+2)
    yce       = LinRange( ymin-dy/2, ymax+dy/2, ncy+2)
    xv        = LinRange( xmin     , xmax     , ncx+1)
    yv        = LinRange( ymin     , ymax     , ncy+1)
    return dx, dy, xc, yc, xce, yce, xv, yv
end
export GenerateMesh

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

@views function GenerateMesh!( dom::ModelDomain )
    """
    This functions generate x and y coordinates of mesh centroids and vertices, it also computes the grid spacing
    """
    dom.dx, dom.dy    = (dom.xmax-dom.xmin)/dom.ncx, (dom.ymax-dom.ymin)/dom.ncy
    dom.xc        = LinRange( dom.xmin+dom.dx/2, dom.xmax-dom.dx/2, dom.ncx  )
    dom.yc        = LinRange( dom.ymin+dom.dy/2, dom.ymax-dom.dy/2, dom.ncy  )
    dom.xce       = LinRange( dom.xmin-dom.dx/2, dom.xmax+dom.dx/2, dom.ncx+2)
    dom.yce       = LinRange( dom.ymin-dom.dy/2, dom.ymax+dom.dy/2, dom.ncy+2)
    dom.xv        = LinRange( dom.xmin     , dom.xmax     , dom.ncx+1)
    dom.yv        = LinRange( dom.ymin     , dom.ymax     , dom.ncy+1)
    return
end
export GenerateMesh!

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

@views function PureShearBoxUpdate!( dom::ModelDomain, BC::BoundaryConditions, params::ModelParameters )
    dom.xmin    += mean(BC.Vx.Dir_W) * params.dt
    dom.xmax    += mean(BC.Vx.Dir_E) * params.dt
    dom.ymin    += mean(BC.Vy.Dir_S) * params.dt
    dom.ymax    += mean(BC.Vy.Dir_N) * params.dt
    BC.Vx.Dir_W .=-params.Ebg*dom.xmin*ones(dom.ncy+2)
    BC.Vx.Dir_E .=-params.Ebg*dom.xmax*ones(dom.ncy+2)
    BC.Vy.Dir_S .= params.Ebg*dom.ymin*ones(dom.ncx+2)
    BC.Vy.Dir_N .= params.Ebg*dom.ymax*ones(dom.ncx+2)
    dom.L = dom.xmax - dom.xmin
    println("Box pure shear deformation: ", (dom.L-dom.L0)/dom.L0*100, "%" )
    GenerateMesh!( dom )
    return dom.dx, dom.dy
end
export PureShearBoxUpdate!


text:
@testset "utils.jl" begin
    @testset "submat" begin

        @testset "PDMat" begin
            M = PDMat(TEST_MATRICES["Positive definite"])
            M_dense = Matrix(M)

            subM = submat(M, [1, 3])
            @test subM isa PDMat
            @test Matrix(subM) == M_dense[[1, 3], [1, 3]]
        end

        @testset "PDDiagMat" begin
            M = PDiagMat(diag(TEST_MATRICES["Positive definite"]))
            M_dense = Matrix(M)

            subM = submat(M, [1, 3])
            @test subM isa PDiagMat
            @test Matrix(subM) == M_dense[[1, 3], [1, 3]]
        end

        @testset "PDSparseMat" begin
            M = PDSparseMat(sparse(TEST_MATRICES["Positive definite"]))
            M_dense = Matrix(M)

            subM = submat(M, [1, 3])
            @test subM isa PDSparseMat
            @test Matrix(subM) == M_dense[[1, 3], [1, 3]]
        end

        @testset "ScalMat" begin
            M = ScalMat(3, 0.1)
            M_dense = Matrix(M)

            subM = submat(M, [1, 3])
            @test subM isa ScalMat
            @test Matrix(subM) == M_dense[[1, 3], [1, 3]]
            @test_throws BoundsError submat(M, [1, 10])
        end

        @testset "PSDMat" begin
            SM = TEST_MATRICES["Positive semi-definite"]
            pivoted = cholesky(SM, Val(true); check=false)
            M = PSDMat(SM, pivoted)
            M_dense = Matrix(M)

            subM = submat(M, [1, 3])
            @test subM isa PSDMat
            @test Matrix(subM) == M_dense[[1, 3], [1, 3]]
        end

        @testset "WoodburyPDMat" begin
            A = randn(4, 2)
            D = Diagonal(randn(2).^2 .+ 1)
            S = Diagonal(randn(4).^2 .+ 1)
            x = randn(size(A, 1))

            W = WoodburyPDMat(A, D, S)
            W_dense = PDMat(Symmetric(Matrix(W)))

            subW = submat(W, [1, 3])
            @test subW isa WoodburyPDMat
            @test subW.A == W.A[[1, 3], :]
            @test subW.D == W.D
            @test subW.S == Diagonal(W.S.diag[[1, 3]])
            @test Matrix(subW) == W_dense[[1, 3], [1, 3]]
        end
    end
end


text:
#Context contains all necessary structures for program building -> solving

mutable struct Context
    ctl
end


text:
# # "Unstable API" Analysis
#
# ```@meta
# CurrentModule = JET
# ```

# ## Motivation
#
# Julia doesn't have any facilities to truly hide module internals.
# This means, we can always access to whatever defined within a module and use it freely,
# but some of them may be considered as the module's "internal"s and subject to changes.
# When possible, we want to avoid their usages for better maintainability in the future.
# But the problem is, how can we automatically find them already used in an existing code ?
#
# This analysis is motivated by [this discussion](https://github.com/JuliaLang/julia/pull/40745#issuecomment-850876150).

# ## Implementation

# Let's define "unstable API" s such that, they're
# - undefined binding, or
# - not `export`ed nor documented, if defined
# and now we can implement such analyzer that detects code that matches the definition
# above using JET.jl's pluggable-analysis framework.
#
# The implementation below is _almost sound_, under the assumption that the bindings are
# resolved statically.
# One thing to note is that, the analysis implements an heuristic to avoid false positives
# from "language intrinsics", for example, `Base.indexed_iterate` and `Base.Broadcast.broadcasted`.
# They're _usually_ introduced into your code implicitly by Julia's iteration protocols and such,
# and we're not responsible for their details (thus not interested in their usages).
# But the problem is that the analyzer below doesn't distinguish those introduced by the
# language and those written by ourselves, and in the latter case we're certainly uses
# "unstable API" under the definition above.

using JET.JETInterface   # to load APIs of the pluggable analysis framework
using JET                # to use analysis entry points
const CC = Core.Compiler # to inject a customized report pass

# First off, we define `UnstableAPIAnalyzer`, which is a new [`AbstractAnalyzer`](@ref) and will
# implement the customized report pass

struct UnstableAPIAnalyzer{T} <: AbstractAnalyzer
    state::AnalyzerState
    is_target_module::T
end
function UnstableAPIAnalyzer(;
    is_target_module = ==(@__MODULE__),
    jetconfigs...)
    return UnstableAPIAnalyzer(AnalyzerState(; jetconfigs...), is_target_module)
end
JETInterface.AnalyzerState(analyzer::UnstableAPIAnalyzer) = analyzer.state
JETInterface.AbstractAnalyzer(analyzer::UnstableAPIAnalyzer, state::AnalyzerState) =
    UnstableAPIAnalyzer(state, analyzer.is_target_module)
JETInterface.ReportPass(analyzer::UnstableAPIAnalyzer) = UnstableAPIAnalysisPass()
JETInterface.get_cache_key(analyzer::UnstableAPIAnalyzer) = AnalyzerState(analyzer).param_key

# Next, we overload some of `Core.Compiler`'s [abstract interpretation](@ref abstractinterpret) methods,
# and inject a customized analysis pass (here we gonna name it `UnstableAPIAnalysisPass`).
# In this analysis, we are interested in whether a binding that appears in a target code is
# an "unstable API" or not, and we can simply check if each abstract element appeared during
# abstract interpretation meets our criteria of "unstable API".
# For that purpose, it's suffice to overload `Core.Compiler.abstract_eval_special_value`
# and `Core.Compiler.builtin_tfunction`.
# To inject a report pass, we use [`ReportPass(::AbstractAnalyzer)`](@ref) interface.

struct UnstableAPIAnalysisPass <: ReportPass end

function CC.abstract_eval_special_value(analyzer::UnstableAPIAnalyzer, @nospecialize(e), vtypes::CC.VarTable, sv::CC.InferenceState)
    if analyzer.is_target_module(sv.mod) # we care only about what we wrote
        ReportPass(analyzer)(UnstableAPI, analyzer, sv, e)
    end

    ## recurse into JET's default abstract interpretation routine
    return Base.@invoke CC.abstract_eval_special_value(analyzer::AbstractAnalyzer, e, vtypes::CC.VarTable, sv::CC.InferenceState)
end

function CC.builtin_tfunction(analyzer::UnstableAPIAnalyzer, @nospecialize(f), argtypes::Vector{Any}, sv::CC.InferenceState)
    if f === getfield
        if length(argtypes) ≥ 2
            a1, a2 = argtypes[1:2]
            if isa(a1, Core.Const) && (v1 = a1.val; isa(v1, Module))
                if isa(a2, Core.Const) && (v2 = a2.val; isa(v2, Symbol))
                    if analyzer.is_target_module(sv.mod) || # we care only about what we wrote, but with relaxed filter
                       (parent = sv.parent; isa(parent, CC.InferenceState) && analyzer.is_target_module(parent.mod))
                        ReportPass(analyzer)(UnstableAPI, analyzer, sv, GlobalRef(v1, v2))
                    end
                end
            end
        end
    end

    ## recurse into JET's default abstract interpretation routine
    return Base.@invoke CC.builtin_tfunction(analyzer::AbstractAnalyzer, f, argtypes::Vector{Any}, sv::CC.InferenceState)
end

# Additionally, we can cut off the performance cost involved with Julia's native compiler's optimizations passes:
CC.may_optimize(analyzer::UnstableAPIAnalyzer) = return false

# Now we implement the body of our analysis.
# We define "unstable API"s such that they're:
# 1. undefined binding, or
# 2. not `export`ed nor documented, if defined
# and we're not interested in any other program properties other than whether our code contains "unstable API"s or not.

# So in our report pass, we would like to ignore all the reports implemented by JET.jl by default
(::UnstableAPIAnalysisPass)(T::Type{<:InferenceErrorReport}, analyzer, state, @nospecialize(spec_args...)) = return

# but except the report of undefined global references (i.e. `GlobalUndefVarErrorReport`).
# This overload allow us to find code that falls into the category 1.
function (::UnstableAPIAnalysisPass)(T::Type{JET.GlobalUndefVarErrorReport}, analyzer, state, @nospecialize(spec_args...))
    JET.BasicPass()(T, analyzer, state, spec_args...) # forward to JET's default report pass
end

# And now we will define new [`InferenceErrorReport`](@ref) report type `UnstableAPI`,
# which represents the category 2, and implement a report pass to detect it.

@reportdef struct UnstableAPI <: InferenceErrorReport
    g::GlobalRef
end
function JETInterface.get_msg(::Type{UnstableAPI}, sv, g::GlobalRef)
    (; mod, name) = Base.resolve(g) # resolve to original name
    return "$mod.$name is unstable !"
end

function (::UnstableAPIAnalysisPass)(::Type{UnstableAPI}, analyzer::UnstableAPIAnalyzer, sv, @nospecialize(e))
    if isa(e, GlobalRef)
        isdefined(e.mod, e.name) || return false # this global reference falls into the category 1, should be caught by `GlobalUndefVarErrorReport` instead

        (; mod, name) = Base.resolve(e) # this reference will be safely resolved
        analyzer.is_target_module(mod) && return # we don't care about what we defined ourselves

        if isunstable(mod, name)
            add_new_report!(sv.result, UnstableAPI(sv, e))
        end
    end
end

# In the report pass above, `isunstable` will take the heavy lifting to find "unstable API"s.
# Here we will implement `isunstable` according to the definition above but with some heuristics
# to exclude language intrinsics, which can automatically be included into our code and aren't
# usually of our interest.

function isunstable(mod, name)
    ## exclude language intrinsics
    mod === Core && return false
    x = getfield(mod, name)
    x isa Core.Builtin && return false
    (x === Base.indexed_iterate || x === Base.SizeUnknown) && return false # iteration protocol
    (x === Base.Iterators.Filter || x === Base.Iterators.Flatten) && return false # iterator protocol
    x === Base.Broadcast.broadcasted && return false # broadcase protocol
    x === Base.kwerr && return false # ignore keyword lowering

    return !isexported(mod, name) && !hasdoc(mod, name)
end

function isexported(mod, name)
    (; mod, name) = Base.resolve(GlobalRef(mod, name))
    return Base.isexported(mod, name)
end

## adapted from https://github.com/JunoLab/CodeTools.jl/blob/56e7f0b514a7476864c27523bcf9d4bc04699ce1/src/summaries.jl#L24-L34

using Base.Docs
function hasdoc(mod, name)
    binding = Docs.Binding(mod, name)
    for m in Docs.modules
        meta = Docs.meta(m)
        haskey(meta, binding) && return true
        (; mod, var) = binding
        isdefined(mod, var) && haskey(meta, getfield(mod, var)) && return true
    end
    return false
end

# ## Usages

# Now we find "unstable API"s in your code using [JET's analysis entry points](@ref jetanalysis-entry)
# with passing `UnstableAPIAnalyzer` as the `analyzer` configuration.

# ### Simple cases

# Let's first use the [interactive analysis entries](@ref jetanalysis-interactive-entry) and
# try simple test cases.

# `UnstableAPIAnalyzer` can find an "unstable" function:
function some_reflection_code(@nospecialize(f))
    return any(Base.hasgenerator, methods(f)) # Base.hasgenerator is unstable
end
@report_call analyzer=UnstableAPIAnalyzer some_reflection_code(sin)

# `UnstableAPIAnalyzer` can find an "unstable" global variable:
module foo; bar = 1 end
report_call((Any,); analyzer=UnstableAPIAnalyzer) do a
    foo.bar + a # foo.bar is unstable
end

# `UnstableAPIAnalyzer` can detect "unstable API"s even if they're imported binding or
# nested reference (, which will be resolve to `getproperty`)
import Base: hasgenerator
report_call((Any,); analyzer=UnstableAPIAnalyzer) do mi
    ## NOTE every function call appearing here is unstable
    ci = hasgenerator(mi) ? Core.Compiler.get_staged(mi) : Base.uncompressed_ast(mi)
end

# ### Analyze a real-world package

# Finally we can use [JET's top-level analysis entry points](@ref jetanalysis-toplevel-entry) to analyze
# a whole script or package.
#
# Here we will run `UnstableAPIAnalyzer` on [IRTools.jl](https://github.com/FluxML/IRTools.jl),
# which uses `Base.isgenerated`, which is renamed to `Base.hasgenerator` in Julia v1.7 and
# invoked the discussion at <https://github.com/JuliaLang/julia/pull/40745#issuecomment-850876150>.
# Especially, it uses `Base.isgenerator` [here](https://github.com/FluxML/IRTools.jl/blob/1f3f43be654a41d0db154fd16b31fdf40f30748c/src/reflection/reflection.jl#L49),
# and you can see the analyzer correctly detects it if you run the following code with IRTools@v0.4.2 installed.

using Pkg #src
if "IRTools" in keys(Pkg.project().dependencies) #src
is_irtools(mod) = occursin("IRTools", string(Symbol(mod))) # module context will be virtualized by `report_package`, thus use string match
report_package("IRTools"; analyzer=UnstableAPIAnalyzer, is_target_module=is_irtools)
else #src
@warn "IRTools isn't installed in the current environment at $(Pkg.project().path)" #src
end #src

# ```
# ═════ 59 possible errors found ═════
# ┌ @ /Users/aviatesk/.julia/packages/IRTools/aSVI5/src/reflection/reflection.jl:39 Core.kwfunc(IRTools.Inner.invoke_meta)(Core.apply_type(Core.NamedTuple, (:world,))(Core.tuple(world)), IRTools.Inner.invoke_meta, T)
# │┌ @ /Users/aviatesk/.julia/packages/IRTools/aSVI5/src/reflection/reflection.jl:69 IRTools.Inner.#invoke_meta#6(world, _3, T)
# ││┌ @ /Users/aviatesk/.julia/packages/IRTools/aSVI5/src/reflection/reflection.jl:74 Core.kwfunc(IRTools.Inner.meta)(Core.apply_type(Core.NamedTuple, (:types, :world))(Core.tuple(S, world)), IRTools.Inner.meta, T)
# │││┌ @ /Users/aviatesk/.julia/packages/IRTools/aSVI5/src/reflection/reflection.jl:38 IRTools.Inner.#meta#1(types, world, _3, T)
# ││││┌ @ /Users/aviatesk/.julia/packages/IRTools/aSVI5/src/reflection/reflection.jl:43 Base._methods_by_ftype
# │││││ Base._methods_by_ftype is unstable !: Base._methods_by_ftype
# ││││└─────────────────────────────────────────────────────────────────────────────────
# ││││┌ @ /Users/aviatesk/.julia/packages/IRTools/aSVI5/src/reflection/reflection.jl:49 Base.isgenerated
# │││││ variable Base.isgenerated is not defined: Base.isgenerated
# ││││└─────────────────────────────────────────────────────────────────────────────────
# ││││┌ @ /Users/aviatesk/.julia/packages/IRTools/aSVI5/src/reflection/reflection.jl:49 Base.uncompressed_ast
# │││││ Base.uncompressed_ast is unstable !: Base.uncompressed_ast
# ││││└─────────────────────────────────────────────────────────────────────────────────
# ││││┌ @ /Users/aviatesk/.julia/packages/IRTools/aSVI5/src/reflection/reflection.jl:54
# ... # many other "unstable API"s detected
# ```


text:
#=
Illustrates vectors in the plane.

@author : Spencer Lyon <spencer.lyon@nyu.edu>
          Victoria Gregory <victoria.gregory@nyu.edu>

@date: 07/09/2014
=#

using Plots
pyplot()
using LaTeXStrings

function plane_fig()

  vecs = ([2, 4], [-3, 3], [-4, -3.5])
  x_vals = zeros(2, length(vecs))
  y_vals = zeros(2, length(vecs))
  labels = []

  # Create matrices of x and y values, labels for plotting
  for i = 1:length(vecs)
    v = vecs[i]
    x_vals[2, i] = v[1]
    y_vals[2, i] = v[2]
    labels = [labels; (1.1 * v[1], 1.1 * v[2], "$v")]
  end

  plot(x_vals, y_vals, arrow=true, color=:blue,
       legend=:none, xlims=(-5, 5), ylims=(-5, 5),
       annotations = labels, xticks=-5:1:5, yticks=-5:1:5)
  vline!([0], color=:black)
  hline!([0], color=:black)
  plot!(foreground_color_axis=:white, foreground_color_text=:white,
        foreground_color_border=:white)

end

function scalar_multiply()  # illustrate scalar multiplication

  x = [2, 2]
  scalars = [-2, 2]

  # Create matrices of x and y values, labels for plotting
  x_vals = zeros(2, 1 + length(scalars))
  y_vals = zeros(2, 1 + length(scalars))
  labels = []
  x_vals[2, 3] = x[1]
  y_vals[2, 3] = x[2]
  labels = [labels; (x[1] + 0.4, x[2] - 0.2, L"$x$")]

  # Perform scalar multiplication, store results in plotting matrices
  for i = 1:length(scalars)
    s = scalars[i]
    v = s .* x
    x_vals[2, i] = v[1]
    y_vals[2, i] = v[2]
    labels = [labels; (v[1] + 0.4, v[2] - 0.2, LaTeXString("\$$s x\$"))]
  end

  plot(x_vals, y_vals, arrow=true, color=[:red :red :blue],
       legend=:none, xlims=(-5, 5), ylims=(-5, 5),
       annotations = labels, xticks=-5:1:5, yticks=-5:1:5)
  vline!([0], color=:black)
  hline!([0], color=:black)
  plot!(foreground_color_axis=:white, foreground_color_text=:white,
        foreground_color_border=:white)

end

# Plot the first figure --- three vectors in the plane

plane_fig()

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
using WeaponSystems
using Documenter

#DocMeta.setdocmeta!(WeaponSystems, :DocTestSetup, :(using WeaponSystems); recursive=true)
push!(LOAD_PATH,"../src/")
makedocs(;
    modules=[WeaponSystems],
    #authors="Irene Ndindabahizi, Royal Millitary Academy, Belgium",
    #repo="https://github.com/YayeIrene/WeaponSystems.jl/blob/{commit}{path}#{line}",
    sitename="WeaponSystems.jl",
  #=  format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://YayeIrene.github.io/WeaponSystems.jl",
        assets=String[],
    ),=#
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/YayeIrene/WeaponSystems.jl.git",
)


text:
using JLD
using Random
import Distributions.Bernoulli

let
  accuracy_exact = zeros(6)
  timing_exact = zeros(6)
  nnz_exact = zeros(6)

  accuracy_predict = zeros(6,3)
  timing_predict = zeros(6,3)
  nnz_predict = zeros(6,3)

  ########################################
  # Setting up the problem
  ########################################
  Random.seed!(123)
  include("BEMUtils.jl")
  λ = 1.0
  @show Threads.nthreads()

  n_test = 100
  test_points = points([1.0, 1.0, 1.0] .- [1.0, 0.1, 0.1] .* rand(3, n_test))

  n_charges = 1000
  charges_weights = [sign.(rand(Bernoulli(0.8), n_charges) .- 0.5);
                     sign.(rand(Bernoulli(0.3), n_charges) .- 0.5)]
  charges_points = points([(rand(3, n_charges) .* [1.0, 1.0, 0.1] .+ [0.0, 0.0, 1.0])';
                           (rand(3, n_charges) .* [1.0, 0.1, 1.0] .+ [0.0, 1.0, 0.0])']')

  u_true = assembleMatrix(test_points, charges_points) * charges_weights
  ########################################

  for q ∈ 3:8
    for ρ ∈ [1.0, 2.0, 3.0]
    hwv = haarCube(q)

    ########################################
    # Loading the "approximate approximation"
    ########################################

      println("Loading ordering from file")
      ld = load("./out/jld/BEM_sort_q_$(q).jld")
      P = ld["P"]
      revP = ld["revP"]
      distances = ld["distances"]
      hwv = hwv[P]

      println("loading boundary data from file")
      ld = load("./out/jld/BEM_boundaryData_q_$(q).jld")
      boundaryData = ld["boundaryData"]

      println("loading prediction from file")
      ld = load("./out/jld/BEM_predict_q_$(q)_rho_$(Int(ρ * 10)).jld")
      u_predict = ld["u_predict"]
      time_predict = ld["time_predict"]
      length_rowval = ld["length_rowval"]

      accuracy_predict[q - 2, Int(ρ)] = mean((u_true - u_predict).^2)
      timing_predict[q - 2, Int(ρ)] = time_predict
      nnz_predict[q - 2, Int(ρ)] = length_rowval
    ########################################

    ########################################
    # Loading the "exact approximation" if existant
    ########################################
    if isfile("./out/jld/BEM_exact_q_$(q).jld")
      ld = load("./out/jld/BEM_exact_q_$(q).jld")
      Θ = ld["Θ"]
      time_assemble_dense = ld["timing_assemble_dense"]

      time_compute_dense = @elapsed u_BEM_exact = assembleMatrix(test_points, hwv) * (Θ \ boundaryData)
      time_dense = time_compute_dense + time_assemble_dense

      @show q
      @show accuracy_exact[q - 2] = mean((u_true - u_BEM_exact).^2)
      timing_exact[q - 2] = time_compute_dense + time_dense
    else
      accuracy_exact[q - 2] = NaN
      timing_exact[q - 2] = NaN
    end
    nnz_exact[q - 2] = (length(hwv) + 1) * length(hwv) / 2
    ########################################
    end
  end
save("./out/jld/BEM_Plot.jld",
     "accuracy_exact", accuracy_exact,
     "timing_exact", timing_exact,
     "nnz_exact", nnz_exact,
     "accuracy_predict", accuracy_predict,
     "timing_predict", timing_predict,
     "nnz_predict", nnz_predict)
end


text:
function time_adaptive_grid(mparams, capital, co2, prior_loc, prior_scale,
    damage, loc_quad, dam_quad, states)
"""
Adapt grid to simulated trajectories: learning
"""

    smax = Array{Float64}(undef, mparams[:time_horizon] + 1, 5)
    smin = Array{Float64}(undef, mparams[:time_horizon] + 1, 5)

    smin[1:mparams[:sim_length_adapt] + 1, 1], smax[1:mparams[:sim_length_adapt] + 1, 1] = get_bounds(capital, 1e-6)
    smin[1:mparams[:sim_length_adapt] + 1, 2], smax[1:mparams[:sim_length_adapt] + 1, 2] = get_bounds(co2, 1e-6)
    smin[1:mparams[:sim_length_adapt] + 1, 3], smax[1:mparams[:sim_length_adapt] + 1, 3] = get_bounds(prior_loc, 1e-6)
    smin[1:mparams[:sim_length_adapt] + 1, 4], smax[1:mparams[:sim_length_adapt] + 1, 4] = get_bounds(prior_scale, 1e-6)
    smin[1:mparams[:sim_length_adapt] + 1, 5], smax[1:mparams[:sim_length_adapt] + 1, 5] = get_bounds(damage, 1e-6)

    # Add buffer
    smax[:, [1 2]] *= 1.05
    smax[:, 4] *= 1.1
    smin *= .95

    smax[:, 5] .= mparams[:initial_smax][5]

    for t = mparams[:sim_length_adapt]+1:mparams[:time_horizon]+1
        smax[t,:] = mparams[:initial_smax]
        smin[t,:] = mparams[:initial_smin]
    end

    return smax, smin

end

function time_adaptive_grid(mparams, capital, co2, damage, dam_quad, states)
"""
Adapt grid to simulated trajectories: no learning
"""

    smax = Array{Float64}(mparams[:time_horizon] + 1, 3)
    smin = Array{Float64}(mparams[:time_horizon] + 1, 3)

    smin[1:mparams[:sim_length_adapt] + 1, 1], smax[1:mparams[:sim_length_adapt] + 1, 1] = get_bounds(capital, 1e-6)
    smin[1:mparams[:sim_length_adapt] + 1, 2], smax[1:mparams[:sim_length_adapt] + 1, 2] = get_bounds(co2, 1e-6)
    smin[1:mparams[:sim_length_adapt] + 1, 3], smax[1:mparams[:sim_length_adapt] + 1, 3] = get_bounds(damage, 1e-6)

    # Add buffer
    smax[:, [1 2]] *= 1.05
    smin *= .95

    smax[:, 3] = mparams[:initial_smax][3]

    for t = mparams[:sim_length_adapt]+1:mparams[:time_horizon]+1
        smax[t, :] = mparams[:initial_smax]
        smin[t, :] = mparams[:initial_smin]
        smax[t, 3] = mparams[:initial_smax][3]
    end

    return smax, smin

end

max_stars_repo_name: irudik/optimal-climate-policy-aej

max_issues_repo_name: irudik/optimal-climate-policy-aej

max_forks_repo_name: irudik/optimal-climate-policy-aej


text:
# Create example data
s = LinRange(-4, 4, 50)
v = cos.(s) .+ cos.(s)' .+ cos.(reshape(s,1,1,:))
# Draw the isosurface
isosurface(v, 0.5, color=0x99ffcc)


text:
using YaoQX
using Test
using TestSetExtensions

@testset ExtendedTestSet "YaoQX.jl" begin
    @includetests ARGS
end


text:
export ForwardBatcher, next!, current, release!

type ForwardBatcher{C, A} <: Batcher
    ctx::C
    inputs::A
    targets::A
    size::Int
    index::Int
    inputSeg::A
    targetSeg::A
end

function ForwardBatcher(ctx, inputs, targets, size)
    @assert dims(inputs, 0) == dims(targets, 0)
    @assert size < dims(inputs, 0)
    @assert dims(inputs, 0) % size == 0
    A = typeof(inputs)
    ForwardBatcher(ctx, inputs, targets, size, 0, array(ctx), array(ctx))
end

function release!(br::ForwardBatcher)
    release(br.inputSeg)
    release(br.targetSeg)
end

function next!(br::ForwardBatcher)
    scope!(br.ctx) do this
        startPos = br.index * br.size
        endPos = startPos + br.size - 1

        br.inputSeg[] = br.inputs[seq(br.ctx, startPos, endPos), :]
        br.targetSeg[] = br.targets[seq(br.ctx, startPos, endPos), :]

        if endPos == dims(br.inputs, 0) - 1
            br.index = 0
        else
            br.index += 1
        end
    end
end

current(br::ForwardBatcher) = (isEmpty(br.inputSeg) && error("Bacher is not initialized."); (br.inputSeg, br.targetSeg))

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
using Images

import Plots; using Plots
pyplot()

import Statistics; using Statistics
import TiffImages; using TiffImages

#stack = TiffImages.load("multispec_newsetup2.tif", mmap = true)
filename = "C:/Users/Alexander Collins/Desktop/Code/Data00005.tif"
stack = TiffImages.load(filename, mmap = true)

function plotimagemean(stack)
	N = length(stack[1,1,:])

	[Float64(sum(stack[:,:,i]) / length(stack[:,:,1])) for i in 1:N]
end

meanvals = plotimagemean(stack)

plot(1:length(stack[1,1,:]), meanvals ./ maximum(meanvals), xlabel = "Frame", ylabel = "Normalised Average Intensity", legend = false, dpi = 300, grid = false)

#savefig("C:\\Users\\Alexander Collins\\Desktop\\Code\\counts.png")


text:
struct RealizedMeta
    start_time::Dates.DateTime
    resolution::Dates.TimePeriod
    len::Int
    start_offset::Int
    end_offset::Int
    interval_len::Int
    realized_timestamps::AbstractVector{Dates.DateTime}
end

function RealizedMeta(
    res::SimulationProblemResults;
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    existing_timestamps = get_timestamps(res)
    interval = existing_timestamps.step
    resolution = get_resolution(res)
    interval_len = Int(interval / resolution)
    realized_timestamps = get_realized_timestamps(res, start_time=start_time, len=len)

    result_start_time = existing_timestamps[findlast(
        x -> x .<= first(realized_timestamps),
        existing_timestamps,
    )]
    result_end_time = existing_timestamps[findlast(
        x -> x .<= last(realized_timestamps),
        existing_timestamps,
    )]

    len = length(result_start_time:interval:result_end_time)

    start_offset = length(result_start_time:resolution:first(realized_timestamps))
    end_offset = length(
        (last(realized_timestamps) + resolution):resolution:(result_end_time + interval - resolution),
    )

    return RealizedMeta(
        result_start_time,
        resolution,
        len,
        start_offset,
        end_offset,
        interval_len,
        realized_timestamps,
    )
end

function get_realization(
    result_values::Dict{
        OptimizationContainerKey,
        SortedDict{Dates.DateTime, DataFrames.DataFrame},
    },
    meta::RealizedMeta,
)
    realized_values = Dict{OptimizationContainerKey, DataFrames.DataFrame}()
    for (key, result_value) in result_values
        results_concat = Dict{Symbol, Vector{Float64}}()
        for (step, (t, df)) in enumerate(result_value)
            first_id = step > 1 ? 1 : meta.start_offset
            last_id =
                step == meta.len ? meta.interval_len - meta.end_offset : meta.interval_len
            for colname in propertynames(df)
                colname == :DateTime && continue
                col = df[!, colname][first_id:last_id]
                if !haskey(results_concat, colname)
                    results_concat[colname] = col
                else
                    results_concat[colname] = vcat(results_concat[colname], col)
                end
            end
        end
        realized_values[key] = DataFrames.DataFrame(results_concat, copycols=false)
        DataFrames.insertcols!(
            realized_values[key],
            1,
            :DateTime => meta.realized_timestamps,
        )
    end
    return realized_values
end


text:
# A `Relation` should support `row ∈ relation` (and `r1 ⨝ r2`)
abstract type Relation{names}; end

@inline names(::Relation{_names}) where {_names} = _names
@inline labels(::Relation{names}) where {names} = labels(names...)
@pure labels(names::Symbol...) = _labels(names...)
@inline _labels(name::Symbol) = (Label{name}(),)
@inline _labels(name::Symbol, other_names...) = (Label{name}(), _labels(other_names...)...)


text:
if VERSION < v"0.4-"
   using Dates, TimeSeries, DataFrames
else
   using Base.Dates, TimeSeries, DataFrames
end

module Quandl

if VERSION < v"0.4-"
   using Dates, TimeSeries, DataFrames, JSON, Requests
else
   using Base.Dates, TimeSeries, DataFrames, JSON, Requests
end

export quandlget,
       quandl,
       quandlsearch,
       set_auth_token,
       interactivequandl

include("api.jl")
include("timearray.jl")
include("dataframe.jl")

# Create empty auth token if none exists
if !isfile(Pkg.dir("Quandl/token/auth_token"))
    set_auth_token("")
end

end