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

# using Markdown
# using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ c0a30957-4c7b-4d7b-bfa9-c2fb691a077b
#=╠═╡ notebook_exclusive
begin
	using Revise
	using PlutoUtils
	using PlotlyBase
	using BenchmarkTools
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ 74975885-9a4e-4857-8135-9e4f69061caf
begin
	using DocStringExtensions
	using StaticArrays
	using LinearAlgebra
	using Dictionaries
	using SplitApplyCombine
end

# ╔═╡ 379613ec-0973-4000-ae8c-d7c33ddca18e
#=╠═╡ notebook_exclusive
md"""
# Packages
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 736b0cf6-bec2-4226-8ef4-70f6a865d34a
#=╠═╡ notebook_exclusive
ToC()
  ╠═╡ notebook_exclusive =#

# ╔═╡ f8243a65-9f5e-464e-bb06-0bb4f5131b8b
#=╠═╡ notebook_exclusive
md"""
# Exports
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 7f645e69-3334-44db-9ba1-9f2d3e0127a2
const ColorOrderDict{I} = Dictionary{SVector{2, I}, I}

# ╔═╡ 8660a7c4-eb78-4e7c-966b-d759df7f3dfa
#=╠═╡ notebook_exclusive
md"""
# Lattice Functions
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 71163795-9695-4f11-acc2-6e3838c8a158
#=╠═╡ notebook_exclusive
md"""
## generate\_regular_lattice
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ da97848f-a7ff-4f2d-b98d-e8bf1ccc3038
function lattice_generator(dx::T, dy::T, ds::T;x0::T = T(0), y0::T = T(0), M::Int = 70,N::Int = M) where T<:Real
	# Function to generate x position as function of row,column number m,n
	x(m, n) = m * dx + n * ds + x0
	# Function to generate y position as function of row,column number m,n
	y(n) = n * dy + y0
	# Generate the elements. For each row, shift the columns to always have the search domain around x=0
	gen = (SVector(x(m - round(Int,n * ds / dx), n), y(n)) for n in -N:N,m in -M:M)
	return gen
end

# ╔═╡ f8a53711-e07f-4b6b-84ea-803679496571
"""
Generate a regular lattice of points
$(TYPEDSIGNATURES)

# Arguments
- `dx` → element spacing on the x axis
- `dy` → element spacing on the y axis
- `ds` → displacement along x between rows of elements
- `f_cond::Function` → function of two arguments (`f(x,y) = ...`) returning `true` if element at position `x,y` must be kept and `false` otherwise

# Keyord Arguments
- `x0 = 0` → x coordinate of the origin of the lattice
- `y0 = 0` → y coordinate of the origin of the lattice
- `M::Int = 70` → Number of elements to generate per row of points before appliying the filtering function `f_cond`
- `N::Int = M` → Number of rows of points to generate before appliying the filtering function `f_cond`
"""
function generate_regular_lattice(dx::T, dy::T, ds::T, f_cond::Function = (x, y) -> true;x0::T = T(0), y0::T = T(0), M::Int = 70,N::Int = M) where T<:Real
	gen = lattice_generator(dx, dy, ds; x0, y0, M, N)
	return [x for x ∈ gen if f_cond(x...)]
end

# ╔═╡ ce22b91e-6bba-4312-a89b-1a78f84034d3
function regular_lattice_nelements(dx::T, dy::T, ds::T, f_cond::Function = (x, y) -> true;x0::T = T(0), y0::T = T(0), M::Int = 70,N::Int = M) where T<:Real
	gen = lattice_generator(dx, dy, ds; x0, y0, M, N)
	sum(x -> f_cond(x...), gen)
end

# ╔═╡ 0eb5d19d-b535-4cda-a89b-26ba295e2711
generate_regular_lattice(dx::Real,dy::Real,ds::Real,args...;kwargs...) = generate_regular_lattice(promote(dx,dy,ds)...,args...;kwargs...)

# ╔═╡ 2ba73b51-ecb9-4632-9f39-bdaeb8c5bd34
#=╠═╡ notebook_exclusive
@benchmark generate_square_lattice(1, (x,y) -> x^2 + y^2 < 100)
  ╠═╡ notebook_exclusive =#

# ╔═╡ f0834b38-8efe-4e77-b0f9-47e5b7595191
#=╠═╡ notebook_exclusive
md"""
## generate\_rect_lattice
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 9e5db472-f96a-4acb-96ae-024b5c73a93d
"""
    generate_rect_lattice(spacing_x::Real,spacing_y::Real[,f_cond];kwargs...)
# Summary
Generate a rectangular lattice of points (with different spacing among x and y directions)
# Arguments
- `spacing_x` → spacing between points on the x axis
- `spacing_y` → spacing between points on the y axis

See [`generate_regular_lattice`](@ref) for a description of `f_cond` and of  the keyword arguments
"""
generate_rect_lattice(spacing_x::Real,spacing_y::Real,args...;kwargs...) = generate_regular_lattice(spacing_x,spacing_y,0,args...;kwargs...)

# ╔═╡ c41f9f41-d8bd-4001-98cb-2ab788404b1b
#=╠═╡ notebook_exclusive
md"""
## generate\_square\_lattice
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 8787134f-9d14-4329-8dda-72557e3175b8
"""
`generate_square_lattice(spacing::Real[,f_cond];kwargs...)`
# Summary
Generate a square lattice of points (with equal spacing among x and y directions)
# Arguments
- `spacing` → spacing between points on both x and y axis

See [`generate_regular_lattice`](@ref) for a description of `f_cond` and of  the keyword arguments
"""
generate_square_lattice(spacing::Real,args...;kwargs...) = generate_regular_lattice(spacing,spacing,0,args...;kwargs...)

# ╔═╡ f589306c-919d-468a-a0fd-9367acc36a7b
#=╠═╡ notebook_exclusive
md"""
## generate\_hex\_lattice
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 6ca05079-4c0d-4c45-8486-a4291310189d
"""
`generate_hex_lattice(spacing::Real[,f_cond];kwargs...)`
# Summary
Generate a hexagonal lattice of points (with equal distance between neighboring points).
The hexagonal lattice generated by this function has distance between points on the same
column √3 times greater than the distance between points on the same row.
# Arguments
- `spacing` → spacing between points

See [`generate_regular_lattice`](@ref) for a description of `f_cond` and of  the keyword arguments
"""
generate_hex_lattice(spacing::Real,args...;kwargs...) = generate_regular_lattice(spacing .* (1,√3/2,.5)...,args...;kwargs...)

# ╔═╡ 243621c4-245a-4267-9bb4-568e673450fa
#=╠═╡ notebook_exclusive
generate_hex_lattice(1;M=20, x0 = .5) |> x -> scatter(x; mode="markers") |> Plot
  ╠═╡ notebook_exclusive =#

# ╔═╡ 059045e0-9acc-438d-b3f5-602f8d5892f7
#=╠═╡ notebook_exclusive
md"""
# Misc Functions
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 2f1f02c5-3ea5-40c1-8fae-704d150036e6
# Get the conversion from linear to db and viceversa
"""
$(TYPEDSIGNATURES)
Convert a number from linear to dB
"""
lin2db(x::Real) = 10log10(x)

# ╔═╡ 5be397fe-a531-423c-8be0-5d31df79dd2f
"""
$(TYPEDSIGNATURES)
Convert a number from dB to linear
"""
db2lin(x::Real) = 10^(x/10)

# ╔═╡ b2e80c33-bbfe-43ca-8795-c9d8d6fa52a9
# Convert between frequency and wavelength
"""
$(TYPEDSIGNATURES)
Get the wavelength (in m) starting from the frequency (in Hz)
"""
f2λ(f::Real) = c₀/f

# ╔═╡ 9165c4d4-69b5-456c-813c-4725feeb5b52
"""
$(TYPEDSIGNATURES)
Get the frequency (in Hz) starting from the wavelength (in m)
"""
λ2f(λ::Real) = c₀/λ

# ╔═╡ a5ca5a8a-8497-41e2-9af0-92db5db9ce73
#=╠═╡ notebook_exclusive
md"""
# Generate Colors
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 0ddf072d-009d-42f2-9a8f-f69fbab750c6
#=╠═╡ notebook_exclusive
md"""
The idea here is to find the optimal color ordering that starting from the color associated with the lattice point in (0,0) and populate the available colors in order to always selected the *un-picked* color that has the highest distance from the currently active colors (*picked* lattice points).

To do so, the idea is to verify the distane not just with the unique lattice points, but also with the neighbor lattice points on the direction of the center of the lattice parallelpiped.
To find the direction of the center of the parallelepiped, we subtract half of the lattice generating vectors from each point and compute the resulting point projection on the base of the parallelepipeid. The sign of the coefficients gives 4 possible quadrant which identify where to extend the lattice points for the computation of the closest active beam.
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 5db224d4-7379-4c8f-bcee-9cf00011d286
#=╠═╡ notebook_exclusive
md"""
## plot\_color\_basis
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 5f62c7d2-ebfb-43e6-b9e3-06ca78c99390
"""
$TYPEDSIGNATURES
Plot the two arrows of the two basis vectors identified by matrix F (one vector per column).
Plotly anotations are used for creating the arrows
"""
function plot_basis_vector(F)
	common_params = (
		ax = 0,
		ay = 0,
		axref = "x",
		ayref = "y",
		arrowhead = 1,
		arrowwidth = 2,
		xanchor = "right", # Needed to have the tail match the origin
		yanchor = "top", # Needed to have the tail match the origin
	)
	# Plot the first vector arrow
	a1 = attr(;
		x = F[1,1],
		y = F[2,1],
		common_params...
	)
	# Plot the first vector arrow
	a2 = attr(;
		x = F[1,2],
		y = F[2,2],
		common_params...
	)
	Plot(Layout(annotations = [a1, a2]))
end

# ╔═╡ 6107a3dc-26dd-4a0d-aeff-5eca2cd1dcd4
#=╠═╡ notebook_exclusive
@bind NN Slider(4:20)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 0d33b162-cc06-4c3b-8ade-5b40106dec0e
#=╠═╡ notebook_exclusive
NN
  ╠═╡ notebook_exclusive =#

# ╔═╡ 063db114-1b95-47d8-8f8b-26eaff8f9574
#=╠═╡ notebook_exclusive
md"""
## adjugate
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 8d4eb806-8af1-4127-8f72-6dd68f810eb5
"""
	adjugate(A::SMatrix{2,2,<:Number,4})
Compute the [adjugate matrix](https://en.wikipedia.org/wiki/Adjugate_matrix) for a 2x2 `SMatrix`.

Used to compute the integer vector modulo in 2D space.

See: [`mod`](@ref)
"""
adjugate(A::T) where T <: SMatrix{2, 2, <:Number, 4} = T(A[4], -A[2], -A[3], A[1])

# ╔═╡ d175246c-552a-4f0f-8415-2339e9af833c
#=╠═╡ notebook_exclusive
md"""
## Base.mod
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 087730c6-cf72-4133-8064-d5619ea4b188
"""
# Integer Vector Modulo
	mod(m::SVector{2,<:Integer}, M::SMatrix{2,2,<:Integer,4})
Compute the integer modulo vector operation of an integer 2D vector with respect to an integer 2x2 matrix. Taken from `[1]`, explots the adjugate matrix instead of the inverse to avoid rounding problems

See: [`adjugate`](@ref)

References:
- [[1]](https://doi.org/10.1109/TSP.2020.3023584) L. Xiao, X. -G. Xia and Y. -P. Wang, "Exact and Robust Reconstructions of Integer Vectors Based on Multidimensional Chinese Remainder Theorem (MD-CRT)," in IEEE Transactions on Signal Processing, vol. 68, pp. 5349-5364, 2020, doi: 10.1109/TSP.2020.3023584.
"""
function Base.mod(m::SVector{2,<:Integer}, M::SMatrix{2,2,<:Integer,4})
	M * mod.(adjugate(M)m, det(M))/det(M) |> SVector{2,Int}
end;

# ╔═╡ 745bb1c1-a312-4dbd-a29f-0836b7dbe8a7
#=╠═╡ notebook_exclusive
md"""
## get\_deterministic\_color\_order
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 76275b93-5668-4ac9-a6a2-2f4b30ca8ab3
#=╠═╡ notebook_exclusive
md"""
## get\_color\_illumination\_order
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 9ccf08d4-2ffb-4f36-a632-3b0ed4017d92
#=╠═╡ notebook_exclusive
let
	a = rand(SVector{2,Float64},100)
	@benchmark map(x -> Int.(sign.(x)), $a)
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ ca09052d-1fce-422d-9205-ae4e87dc4db4
#=╠═╡ notebook_exclusive
@benchmark get_color_illumination_order(50; lattice_type = :square)
  ╠═╡ notebook_exclusive =#

# ╔═╡ f700edf6-0c44-4507-bbf5-4c5dc02fa74c
#=╠═╡ notebook_exclusive
NNN = 16
  ╠═╡ notebook_exclusive =#

# ╔═╡ 019c9ab2-ef2b-4677-a5aa-d0363fffef72
#=╠═╡ notebook_exclusive
n5 = 17
  ╠═╡ notebook_exclusive =#

# ╔═╡ cffd2b81-66c1-4f50-948c-e38ec011105d
#=╠═╡ notebook_exclusive
N4 = 4
  ╠═╡ notebook_exclusive =#

# ╔═╡ 7785e6b1-881e-4088-82fe-3dad106b07be
#=╠═╡ notebook_exclusive
n4 = 4
  ╠═╡ notebook_exclusive =#

# ╔═╡ 9cd5609b-1e58-4d39-87ab-b3d7542de691
#=╠═╡ notebook_exclusive
n6, nc = 8, 13
  ╠═╡ notebook_exclusive =#

# ╔═╡ 89888da7-4912-4557-a375-0150f80ee703
#=╠═╡ notebook_exclusive
md"""
## new\_color\_order
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ f942e39c-0d93-4bc8-8ff2-112fed566014
"""
$TYPEDSIGNATURES
Find the new optimal color ordering when starting from color `n` rather than from 1.
# Arguments
- `n` → New starting color.
- `color_order_dict` → original optimal color order for the current coloring scheme obtained via [`get_color_illumination_order`](@ref).
- `F` → lattice generating matrix for the current coloring scheme, obtained using [`generate_F_reuse_matrix`](@ref).
"""
function new_color_order(n::Int, color_order_dict::ColorOrderDict, F)
	@assert n <= length(color_order_dict) "The provided color is higher than the length of the `color_order_dict`"
	# Find the point lattice point associated to the current color
	p_current = filter(x -> x == n, color_order_dict).indices.values[1]
	idict = sort(color_order_dict) |> x->Dictionary(x.values, x.indices.values)
	# Find the new color
	new_order = map(idict) do p
		v = mod(p + p_current, F)
		color_order_dict[v]
	end
	return new_order.values
end

# ╔═╡ df0602a2-d278-4fae-9957-30c85b55598a
#=╠═╡ notebook_exclusive
md"""
## get\_projection
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ a86b2ff6-1967-424f-86eb-40fb4288c8b5
"""
$TYPEDSIGNATURES
Compute the generic projection coefficients (coordinates) of `v1` on the base identified by matrix `B`, composed of the two generating vectors `v2` and `v3`

Return the projection coefficients as a StaticVector of 2 elements.

References:
- [`https://math.stackexchange.com/questions/148199/equation-for-non-orthogonal-projection-of-a-point-onto-two-vectors-representing`](https://math.stackexchange.com/questions/148199/equation-for-non-orthogonal-projection-of-a-point-onto-two-vectors-representing)
"""
function get_projection(v1, B)
	# extract the basis generating vectors
	# Compute the matrix of the linear system of equations
	A = B'B
	# Compute the target vector
	b = B'v1
	# solve the system
	x = A\b
end

# ╔═╡ 41817b1f-4c5e-4a55-93b7-520d6b71ea9c
function get_color_illumination_order(F, unique_points)
	N = length(unique_points) # Number of colors
	# We preallocate the array of the distances
	check_vec = Vector{eltype(unique_points)}(undef, N) # +3 because the origin is replicated 3 times
	# Preallocate the output to be modified
	out = Dictionary(unique_points, ones(Int,N))
	pts = deepcopy(unique_points)
	# Put the expanded origin points
	check_vec[1] = popat!(pts,1)
	@views @inbounds for n = 2:N
		# Find distance between all remaining points and currently populated points
		maxdist = 0.0
		id = 0
		for (i,p) ∈ enumerate(pts)
			mindist = Inf
			# Find the projection of on the generator axis
			h = get_projection(p, F)
			α, β = -Int(sign(h[1])), -Int(sign(h[2]))
			# Here we don't only check p but also the relevant closest lattice points
			for (nn,v) ∈ enumerate((SA[0,0], α*F[:,1], β*F[:,2], F * SA[α, β]))
				pp = v + p
				for c ∈ check_vec[1:n-1]
					dist = norm(pp-c)
					mindist = dist < mindist ? dist : mindist
				end
			end
			# Find the index of the point that has the maximum minimum distance
			maxdist, id = mindist > maxdist ? (mindist, i) : (maxdist, id)
		end
		# Put that vector in the check_vec and change it's order in the dict
		pt = popat!(pts, id)
		check_vec[n] = pt
		set!(out, pt, n)
	end
	out
end

# ╔═╡ 3f2a31d4-0fa8-40fa-9dc4-bd6a26d2ddc9
# Initialize the vector that contains the matrix to compute the beam coloring. We limit ourselves at 500 colors to start
const F_reuse_matrix = (square = SMatrix{2,2,Int,4}[], triangular = SMatrix{2,2,Int,4}[])

# ╔═╡ dbef7000-7c39-49f7-b24e-f0fb436eb54e
#=╠═╡ notebook_exclusive
md"""
## compute\_F\_cell
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 9183aa40-9dc1-4237-8bf5-de42c93b149f
#=╠═╡ notebook_exclusive
md"""
## \_coloring\_inner\_bruteforce!
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 6b9beb62-dc7e-4b8b-9b7c-8fee5b1da98f
function _coloring_inner_bruteforce!(T_mat, rot_mat, grid_max)
    max_colours = length(T_mat)
    check_vec = fill((typemax(Int),typemax(Int)),max_colours)
    @inline norm2(x) = sum(abs2.(x))
    @inbounds for x1 = 0:grid_max, y1 = grid_max:-1:-grid_max, x2 = 0:grid_max, y2 = 0:grid_max
        mat = @SMatrix [x1 y1;x2 y2]
        # Compute the determinant
        t_det = Int(det(mat))
        # Skip points which are not likely to give useful results
        if (t_det < 1)  || (t_det > max_colours) || (t_det > grid_max^2) || (maximum(abs.(mat)) > ceil(sqrt(t_det) + 3))
            continue
        end
        # Compute the angle between the basis vectors identified by [x1,y1] and [x2,y2]
        angle = abs(atan(y1,x1) - atan(y2,x2))
        # Skip cases where the angle between vectors is either too acute or too obtuse
        if abs(π/2-angle) > π/4
            continue
        end
        # Create temp variables for computation of the minimum distance
        dmat = mat*rot_mat
        # Compute frobenius norm and minimum squared distance for the candidate lattice generating matrix
        frobe = round(Int,norm2(dmat))
        # display(frobe)
        # Minimum squared distance is either the modulo of one of the vectors or the modulo of sum or difference of them
        dmin = round(Int,minimum((norm2(dmat[1,:]), norm2(dmat[2,:]), norm2(sum(dmat;dims=1)), norm2(diff(dmat;dims=1)))))
        # Check if the current realization is better than the saved one
        if isless((-dmin,frobe),check_vec[t_det])
            # Update the check_vec
            check_vec[t_det] = (-dmin,frobe)
            # Update the vector containing the generating matrices
            T_mat[t_det] = round.(Int, mat')
        end
    end
end

# ╔═╡ 14cb2a0b-2ea8-471b-987f-1647f1516992
## Here we have the functions for the coloring computation
function compute_F_cell(max_colours::Int;grid_max::Int=ceil(Int,sqrt(max_colours))+5)
    #=
    This function is used to compute all the possible 2x2 lattice generating matrices for possible coloring schemes up to 'max_colours' colors
    Computation is done with a brute-force approach, generating all possible 2x1 vectors with maximum elements up to grid_max.

    The very heuristic automatic number of the grid_max is valid up to 5000 colors
    =#

    @assert max_colours <= 5000 "A number of colors greater than 5000 is currently not supported"

    # Find the current length of the pre-computed vector
    current_length = length(F_reuse_matrix.square)
    if current_length >= max_colours
        # We already computed the function for the required number of colors
        return
    end
    n_missing = max_colours - current_length
    append!(F_reuse_matrix.square,Vector{SMatrix{2,2,Int,4}}(undef,n_missing))
    append!(F_reuse_matrix.triangular,Vector{SMatrix{2,2,Int,4}}(undef,n_missing))
    # Compute the matrix for the square lattice
    _coloring_inner_bruteforce!(F_reuse_matrix.square,I,grid_max)
    # Compute the matrix for the triangular lattice
    _coloring_inner_bruteforce!(F_reuse_matrix.triangular,@SMatrix([1 0;cosd(60) sind(60)]),grid_max)
end

# ╔═╡ 1155e836-99d0-4cc8-83d0-355a6ab6fcc0
#=╠═╡ notebook_exclusive
md"""
## generate\_F\_reuse\_matrix
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 1d44cf1c-11a5-4366-94f3-85b695c6ca12
function generate_F_reuse_matrix(N_Colours::Int=4; lattice_type::Symbol=:triangular, max_colours::Int=max(10,N_Colours))
	@assert lattice_type ∈ (:triangular, :square) "The lattice type has to be either :triangular or :square"
    compute_F_cell(max_colours)
    return getproperty(F_reuse_matrix,lattice_type)[N_Colours]
end

# ╔═╡ a7037a49-c2cb-48b6-bbf9-ca8150afcfbe
#=╠═╡ notebook_exclusive
FF = generate_F_reuse_matrix(:triangular, NN)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 259ea307-2862-42f3-866a-be8f4eb83cf3
#=╠═╡ notebook_exclusive
plot_basis_vector(FF)
  ╠═╡ notebook_exclusive =#

# ╔═╡ e6033b3e-51fa-4eeb-87ff-e5b5359aebc5
"""
$TYPEDSIGNATURES
Assign a given set of unique of N integer 2D vectors, assign a deterministic and consistent order for associating these vectors to indices (1:N) used for coloring.

The order is obtained by sorting the vectors first by the norm and then by the angle it makes with respect to the X-axis (first dimension)
"""
function get_deterministic_color_order(N::Int; lattice_type = :triangular)
	T = SVector{2,Int}
	# Get the F_matrix
	F = generate_F_reuse_matrix(N; lattice_type)
	# Find the range over which to generate the potential lattice points
	M = [F*SA[1,1] F]
	xstart, xend = extrema(M[1,:])
	ystart, yend = extrema(M[2,:])
	# pts = ((x,y) for x ∈ xstart:xend, y ∈ ystart:yend)
	# Create the vector of possible values
	V = zeros(T, N)
	# Start from n = 1 as 0,0 is for sure in the set
	n = 1
	@inbounds for x ∈ xstart:xend, y ∈ ystart:yend
		p = T(x,y)
		# Get the vector modulo the color generating matrix
		v = mod(p,F)
		if v ∉ V
			n += 1
			V[n] = v
			n == N && break
		end
	end
	# Check if we found all the vector
	@assert n == N "Not all unique vectors were found!"
	# We now sort this vector
	sort!(V; by = v -> (norm(v), mod(atan(reverse(v)...), 2π)))
end

# ╔═╡ 272234bf-ea90-4d24-b5f4-1cd6cdaea20b
#=╠═╡ notebook_exclusive
get_deterministic_color_order(150)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 5da650b5-9309-418b-92d9-3388b5385e9b
#=╠═╡ notebook_exclusive
fffvec = get_deterministic_color_order(NNN; lattice_type = :triangular)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 6579518a-2130-4e25-b4eb-966e60203f17
"""
$TYPEDSIGNATURES
Get the illumination order to achieve maximum possible distance between already illuminated colors and next illuminated color.

Returns a `Dictionary` associating each unique color vector obtained from [`get_deterministic_color_order`](@ref) to the ordering color.
"""
function get_color_illumination_order(N; lattice_type = :triangular)
	F = generate_F_reuse_matrix(N;lattice_type)
	vec = get_deterministic_color_order(N; lattice_type)
	get_color_illumination_order(F, vec), F
end

# ╔═╡ 0f605f2c-029b-4af8-8e52-2c7bc4c7626a
#=╠═╡ notebook_exclusive
get_color_illumination_order(4; lattice_type = :square)
  ╠═╡ notebook_exclusive =#

# ╔═╡ f97bb119-0248-4782-8cf9-cca61a666dbf
#=╠═╡ notebook_exclusive
let
	lattice_type = :triangular
	cdict, F = get_color_illumination_order(64; lattice_type)
	order = sort(cdict) |> x -> Dictionary(x.values, x.indices.values)
	pts = map(i -> order[i],1:n5)
	outline = scatter(
		[SA[0,0],F[:,1], F*SA[1,1], F[:,2], SA[0,0]];
		mode = "lines",
		line_color = "black",
	)
	data = [
		outline,
		scatter(pts; mode = "markers"),
	]
	Plot(data)
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ de86f516-b65a-445e-8988-4cfc9dafd000
#=╠═╡ notebook_exclusive
FFF = generate_F_reuse_matrix(NNN; lattice_type = :triangular)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 7f4b9a79-ce48-487f-be37-d884ee9db0e9
#=╠═╡ notebook_exclusive
dio = get_color_illumination_order(FFF, fffvec)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 39820d4a-d05f-4dbc-8d6d-7e20b579688c
#=╠═╡ notebook_exclusive
dio.indices.values
  ╠═╡ notebook_exclusive =#

# ╔═╡ 7bbafd36-6130-4f00-ac86-63b2bc467b77
#=╠═╡ notebook_exclusive
idio = sort(dio) |> x->Dictionary(x.values, x.indices.values)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 2a45ea3b-ea3e-4d5e-912d-36b0055cf307
#=╠═╡ notebook_exclusive
new_color_order(1, dio, FFF)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 759d6fd6-9649-4bea-b088-1cae1647ad3d
#=╠═╡ notebook_exclusive
map(x -> mod(x + SA[2, -2], FFF), idio)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 39948c5f-d472-48ce-8385-1c67b8f1580a
#=╠═╡ notebook_exclusive
md"""
## get\_L
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ f5f5a045-1c9b-40b8-b326-7d0e1336412c
"""
$TYPEDSIGNATURES
Get the lattice generating matrix as a function of the `lattice_type` and of the lattice `spacing`.

`lattice_type` has to be either `:triangular` or `:square`
"""
function get_L(spacing; lattice_type = :triangular)
	@assert lattice_type ∈ (:triangular, :square) "`lattice_type` has to be either `:triangular` or `:square`"
	if lattice_type === :triangular
		return SA[1 1/2;0 √3/2] * spacing
	else
		return SA[1 0;0 1] * spacing
	end
end

# ╔═╡ cd9ee96f-91a9-43cf-980c-6253a0c6018f
#=╠═╡ notebook_exclusive
md"""
## generate\_colors
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 16614e6a-2ef0-4b36-913d-6fd19440b60b
asdf = let
	N = 4
	CD, F = get_color_illumination_order(N; lattice_type = :triangular)
	spacing = .5
	BC = generate_hex_lattice(spacing; M = 100)
	L = spacing * SA[1 1/2; 0 √3/2]
	(BC, CD, L, F)
end

# ╔═╡ d3e49592-cee0-4414-9b92-ce2fb66529df
"""
$TYPEDSIGNATURES
Provide a the coloring breakdown for a given set of lattice points.
# Arguments
- `BeamCenters` → Array of SVector{2,<:Number} points identifying the beam lattice centers
- `colordict` → Dictionary of the optimal color order for the given coloring scheme, obtained using [`get_color_illumination_order`](@ref)
- `L` → Beam Lattice generating matrix, function of the beam spacing and lattice type and obtained using [`get_L`](@ref)
- `F` → Color Lattice generating matrix, function of the number of colors and beam lattice type and obtained using [`generate_F_reuse_matrix`](@ref)
# Keyword Arguments
- `first_color_idx` → Provide the index of the point in BeamCenters that should be associated with color 1.
- `first_color_coord` → Provide the coordinates of the point that should be associated with color 1, if both coord and idx are provided, the idx is ignored. Defaults to the giving color 1 to the point in (0,0).

# Note
The colors are ordered in such a way that if illuminating the colors in sequential order, the next color is always the one that has the maximum minimum distance with points of the previoulsy illuminated colors (so the one most likely to create lower interference).
"""
function generate_colors(BeamCenters, colordict, L, F; first_color_idx = nothing, first_color_coord = first_color_idx === nothing ? zero(eltype(BeamCenters)) : BeamCenters[first_color_idx])
	map(BeamCenters) do p
		T = SVector{2, Int}
		# Find the integer lattice representation
		l = inv(L)*(p - first_color_coord)
		# For some reason doing broadcast round would allocate
		v = SA[round(Int,l[1]), round(Int,l[2])]
		# Get the resulting color based on the modulo value of l w.r.t F
		color = colordict[mod(v,F)]
	end
end;

# ╔═╡ 7e68054e-4268-424e-b413-ef18baf832ac
"""
    generate_colors(BeamCenters::AbstractVector,N_Colours::Int=4;lattice_type::Symbol=:triangular)

Provide a the coloring breakdown for a given set of lattice points.
# Arguments
- `BeamCenters` → Vector of Tuple or StaticVectors expressing the U-V coordinates of each point in the lattice for which the coloring is being computed
- `N_Colours` → Number of colors to divide the lattice in. Defaults to `4`

# keyword Arguments
- `lattice_type` → Symbol that can either be `:triangular` or `:square`, idenifying the type of lattice. Defaults to `:triangular`
"""
function generate_colors(BeamCenters::AbstractVector,N_Colours::Int=4;first_color_coord=nothing,first_color_idx=nothing,precision_digits::Int=7,lattice_type::Symbol=:triangular)
    #=
    **************************************************************************
       Generate frequency colouring file
    **************************************************************************

     References:
     [1]   "On the frequency allocation for mobile radio telephone systems", C. de
           Almeida; R. Palazzo, Proceedings of 6th International Symposium on
           Personal, Indoor and Mobile Radio Communications, Year: 1995, Volume: 1, Pages: 96 - 99 vol.1
     [2]   P. Angeletti, "Simple implementation of vectorial modulo operation based
           on fundamental parallelepiped," in Electronics Letters, vol. 48, no. 3, pp. 159-160,
           February 2 2012. doi: 10.1049/el.2011.3667
     [3]   L. Xiao, X.  -G. Xia and Y. -P. Wang, "Exact and Robust
           Reconstructions of Integer Vectors Based on Multidimensional Chinese
           Remainder Theorem (MD-CRT)," in IEEE Transactions on Signal Processing,
           vol. 68, pp. 5349-5364, 2020, doi: 10.1109/TSP.2020.3023584

     Authors: Alberto Mengali, 2020, European Space Agency

    Input Arguments:
    BeamCenters:         A vector containing the U,V coordinates of the beam centers as tuples or staticvectors
    N_Colours:           An integer number specifying the number of colors to generate in the association
    first_color_coord:   A tuple or static vector containing the U,V coordinates of the beam that will have the first color
    first_color_idx:     The beam index of the beam containing the first color, either this variable of first_order_coord can be specified, not together
    precision_digits:    The number of digits to be used in the flooring operation
    =#

    # Check if either the first color coordinates or first color idx are given
    if first_color_coord !== nothing
        if first_color_idx !== nothing
            @warn "Both first_color_idx and first_color_coord were given, disregarding the idx variable"
        end
        first_color_idx = findfirst(x -> x == first_color_coord, BeamCenters)
        if first_color_idx === nothing
          error("The provided `first_color_coord` $(first_color_coord) is not part of the beam centers vector")
        end
    end
    if lattice_type ∉ (:triangular, :square)
        @error "The specified lattice type ($lattice_type) is not recognized, it should be either :triangular or :square"
    end

    n_beams = length(BeamCenters)

    # If only 1 color is requested, simply return the trivial result
    if N_Colours == 1
        Colours = ones(n_beams)
        nbeams_per_color = n_beams
        idxs = fill(fill(true,n_beams))
    end

    # Find the minimum distance in U and V
    minU_dist = minimum(diff(sort(first.(BeamCenters)) |> x -> unique!(y -> round(y;digits=precision_digits),x)))
    minV_dist = minimum(diff(sort(last.(BeamCenters)) |> x -> unique!(y -> round(y;digits=precision_digits),x)))
    if lattice_type === :triangular
        beamU_spacing = 2minU_dist
        beamV_spacing = 2minV_dist
        # Matrix to normalize the grid points in u-v into a integer grid with y(v) axis not pointing north but north-east with 60° inclination
        D = @SMatrix [1 -1/2;0 1]
    elseif lattice_type === :square
        beamU_spacing = minU_dist
        beamV_spacing = minV_dist
        D = @SMatrix [1 0;0 1]
    end
    # Get the coloring generating matrix
    F_reuse_matrix = generate_F_reuse_matrix(N_Colours;lattice_type)
    # Create the set that will contain the unique results
    unique_colors_vecs = SVector{2,Int}[]
    # Initialize the colors vector
    Colors = similar(BeamCenters, Int)
    @inbounds for (n,p₀) ∈ enumerate(BeamCenters)
        # Compute the integer beam centers
        p₀_normalized = p₀ ./ SA_F64[beamU_spacing, minV_dist]
        # Transform the beam centers from u-v coordinates in radians into integer indexes
        beam_index_vector = round.(Int,D*p₀_normalized);
        # Find the values of the beam indexes modulo F_reuse_matrix
        unique_beam_index = round.(Int,beam_index_vector .- (F_reuse_matrix*floor.(round.(inv(F_reuse_matrix)*beam_index_vector,digits=precision_digits))))
        # Check if this color has already been assigned/registered
        idx = findfirst(x -> x == unique_beam_index,unique_colors_vecs)
        if idx isa Nothing
            push!(unique_colors_vecs, unique_beam_index)
            cidx = length(unique_colors_vecs)
        else
            cidx = idx
        end
        Colors[n] = cidx
    end

    if first_color_idx !== nothing
      # Reshuffle the colors in order to have first_color_idx as color number 1
      ref_col = Colors[first_color_idx]
      @inbounds @simd for n ∈ eachindex(Colors)
        c = Colors[n]
        Colors[n] = mod(c - ref_col,N_Colours) + 1
      end
    end

    return Colors
end;

# ╔═╡ d8584d03-8eb4-4864-b646-a6a0656a2e12
begin
	export generate_regular_lattice, generate_square_lattice, generate_hex_lattice, generate_rect_lattice
	export generate_colors
	export f2λ, λ2f, db2lin, lin2db
	export get_color_illumination_order, get_L, new_color_order
end

# ╔═╡ 66417309-8d6d-47e9-b5ef-7fcc0cde9194
function plot_illumination(n, N, spacing; lattice_type = :triangular)
	cdict, F = get_color_illumination_order(N; lattice_type)
	# Get the max colors
	plotly_colors = PlutoUtils.color_order_64
	# Generate the lattice of the beams
	p = generate_hex_lattice(spacing; M = round(Int, 5/spacing))
	# Find the lattice generating matrix
	L = get_L(spacing; lattice_type)
	# Find the colors
	colors = generate_colors(p, cdict, L, F)
	order = group(i -> colors[i], 1:length(p))
	data = @views map(1:n) do i
		pts = p[order[i]]
		scatter(pts;mode = "markers", marker_color = plotly_colors[i], showlegend = false)
	end
	Plot(data)
end

# ╔═╡ ba3660ba-e274-4e11-9fd1-75aafbb0a776
#=╠═╡ notebook_exclusive
let
	n = n4
	Nmax = N4 .^ [1 2 3]
	spacing = [1 .5 .25]
	P = map(Nmax, spacing) do N,sp
		plot_illumination(min(n, N), N, sp)
	end
	vcat(P...)
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ 47ffd006-29b2-4774-a791-11b3f88a7aba
#=╠═╡ notebook_exclusive
let n = n6, spacing = .5
	N = 16
	lattice_type = :triangular
	cdict, F = get_color_illumination_order(N; lattice_type)
	# Get the max colors
	plotly_colors = PlutoUtils.color_order_64
	# Generate the lattice of the beams
	p = generate_hex_lattice(spacing; M = round(Int, 5/spacing))
	# Find the lattice generating matrix
	L = get_L(spacing; lattice_type)
	# Find the colors
	colors = generate_colors(p, cdict, L, F)
	order = group(i -> colors[i], 1:length(p))
	corder = new_color_order(nc, cdict, F)
	data = @views map(corder[1:n]) do i
		pts = p[order[i]]
		scatter(pts;mode = "markers", marker_color = plotly_colors[i], showlegend = false)
	end
	Plot(data)
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ c2a099fb-446c-4a39-b505-973223c38a27
#=╠═╡ notebook_exclusive
let
	N = 4
	lattice_type = :triangular
	lat = lattice_type === :square ? generate_square_lattice : generate_hex_lattice
	CD, F = get_color_illumination_order(N; lattice_type)
	spacing = 1
	BC = lat(spacing; M = 10)
	L = get_L(spacing; lattice_type)
	colors = generate_colors(BC, CD, L, F; first_color_coord = SA[-1/2,√3/2])
	# colors = generate_colors(BC, N; lattice_type)
	data = scatter(BC; mode = "markers", marker = attr(
		color = colors,
		colorscale = predefined_colorscale(N),
	))
	Plot(data)
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ 8160086a-6349-447c-87ae-880b02fa97f5
#=╠═╡ notebook_exclusive
md"""
# Tests
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 3ea0415c-af14-430c-bf7c-2c7d71b7a333
#=╠═╡ notebook_exclusive
let
	N_colors = 6
	lat = generate_hex_lattice(1; M = 10)
	colors = generate_colors(lat,N_colors)
	data = scatter(lat;mode="markers", marker_color = colors)
	Plot(data)
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ 65b74f4a-e2af-4caf-8719-5a59c6349bb9
#=╠═╡ notebook_exclusive
let
	N_colors = 4
	lat = generate_square_lattice(1; M = 10)
	colors = generate_colors(lat,N_colors; lattice_type = :square)
	data = scatter(lat;mode="markers", marker_color = colors)
	Plot(data)
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
PlutoUtils = "ed5d0301-4775-4676-b788-cf71e66ff8ed"
Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

[compat]
BenchmarkTools = "~1.2.0"
Dictionaries = "~0.3.15"
DocStringExtensions = "~0.8.6"
PlotlyBase = "~0.8.18"
PlutoUtils = "~0.5.5"
Revise = "~3.1.20"
SplitApplyCombine = "~1.2.0"
StaticArrays = "~1.2.13"
"""

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

julia_version = "1.7.0-rc2"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "abb72771fd8895a7ebd83d5632dc4b989b022b5b"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.2"

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

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

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

[[deps.BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
git-tree-sha1 = "61adeb0823084487000600ef8b1c00cc2474cd47"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.2.0"

[[deps.Chain]]
git-tree-sha1 = "cac464e71767e8a04ceee82a889ca56502795705"
uuid = "8be319e6-bccf-4806-a6f7-6fae938471bc"
version = "0.4.8"

[[deps.CodeTracking]]
deps = ["InteractiveUtils", "UUIDs"]
git-tree-sha1 = "9aa8a5ebb6b5bf469a7e0e2b5202cf6f8c291104"
uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
version = "1.0.6"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "a851fec56cb73cfdf43762999ec72eff5b86882a"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.15.0"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

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

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Dictionaries]]
deps = ["Indexing", "Random"]
git-tree-sha1 = "43ae37eac34e76ac97d1a7db28561243e7242461"
uuid = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
version = "0.3.15"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

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

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Glob]]
git-tree-sha1 = "4df9f7e06108728ebf00a0a11edee4b29a482bb2"
uuid = "c27321d9-0574-5035-807b-f59d2c89b15c"
version = "1.3.0"

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

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

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

[[deps.Indexing]]
git-tree-sha1 = "ce1566720fd6b19ff3411404d4b977acd4814f9f"
uuid = "313cdc1a-70c2-5d6a-ae34-0150d3930a38"
version = "1.1.1"

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

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

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

[[deps.JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
git-tree-sha1 = "e273807f38074f033d94207a201e6e827d8417db"
uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
version = "0.8.21"

[[deps.LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"

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

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

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

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

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

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

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

[[deps.LoweredCodeUtils]]
deps = ["JuliaInterpreter"]
git-tree-sha1 = "491a883c4fef1103077a7f648961adbf9c8dd933"
uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
version = "2.1.2"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

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

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

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

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

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

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

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

[[deps.Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

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

[[deps.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"

[[deps.PlotlyBase]]
deps = ["ColorSchemes", "Dates", "DelimitedFiles", "DocStringExtensions", "JSON", "LaTeXStrings", "Logging", "Parameters", "Pkg", "REPL", "Requires", "Statistics", "UUIDs"]
git-tree-sha1 = "180d744848ba316a3d0fdf4dbd34b77c7242963a"
uuid = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
version = "0.8.18"

[[deps.PlutoDevMacros]]
deps = ["MacroTools", "Requires"]
git-tree-sha1 = "f5cbbcaa26fe68fd63bc25f5a9e626b7f0e515b0"
uuid = "a0499f29-c39b-4c5c-807c-88074221b949"
version = "0.4.4"

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

[[deps.PlutoUtils]]
deps = ["Chain", "Colors", "DocStringExtensions", "Glob", "HypertextLiteral", "OrderedCollections", "PlutoDevMacros", "PlutoUI", "PrettyTables", "Reexport", "Requires", "StaticArrays", "UUIDs"]
git-tree-sha1 = "27589f8f79b24c096e92520b03cdab764d2b6a88"
uuid = "ed5d0301-4775-4676-b788-cf71e66ff8ed"
version = "0.5.5"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "d940010be611ee9d67064fe559edbb305f8cc0eb"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.2.3"

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

[[deps.Profile]]
deps = ["Printf"]
uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"

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

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

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

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[deps.Revise]]
deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"]
git-tree-sha1 = "41deb3df28ecf75307b6e492a738821b031f8425"
uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
version = "3.1.20"

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

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

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

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SplitApplyCombine]]
deps = ["Dictionaries", "Indexing"]
git-tree-sha1 = "dec0812af1547a54105b4a6615f341377da92de6"
uuid = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
version = "1.2.0"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

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

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

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

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

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

[[deps.UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

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

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

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

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

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

# ╔═╡ Cell order:
# ╟─379613ec-0973-4000-ae8c-d7c33ddca18e
# ╠═c0a30957-4c7b-4d7b-bfa9-c2fb691a077b
# ╠═74975885-9a4e-4857-8135-9e4f69061caf
# ╠═736b0cf6-bec2-4226-8ef4-70f6a865d34a
# ╟─f8243a65-9f5e-464e-bb06-0bb4f5131b8b
# ╠═d8584d03-8eb4-4864-b646-a6a0656a2e12
# ╠═7f645e69-3334-44db-9ba1-9f2d3e0127a2
# ╟─8660a7c4-eb78-4e7c-966b-d759df7f3dfa
# ╟─71163795-9695-4f11-acc2-6e3838c8a158
# ╠═f8a53711-e07f-4b6b-84ea-803679496571
# ╠═ce22b91e-6bba-4312-a89b-1a78f84034d3
# ╠═da97848f-a7ff-4f2d-b98d-e8bf1ccc3038
# ╠═0eb5d19d-b535-4cda-a89b-26ba295e2711
# ╠═2ba73b51-ecb9-4632-9f39-bdaeb8c5bd34
# ╟─f0834b38-8efe-4e77-b0f9-47e5b7595191
# ╠═9e5db472-f96a-4acb-96ae-024b5c73a93d
# ╟─c41f9f41-d8bd-4001-98cb-2ab788404b1b
# ╠═8787134f-9d14-4329-8dda-72557e3175b8
# ╟─f589306c-919d-468a-a0fd-9367acc36a7b
# ╠═6ca05079-4c0d-4c45-8486-a4291310189d
# ╠═243621c4-245a-4267-9bb4-568e673450fa
# ╟─059045e0-9acc-438d-b3f5-602f8d5892f7
# ╠═2f1f02c5-3ea5-40c1-8fae-704d150036e6
# ╠═5be397fe-a531-423c-8be0-5d31df79dd2f
# ╠═b2e80c33-bbfe-43ca-8795-c9d8d6fa52a9
# ╠═9165c4d4-69b5-456c-813c-4725feeb5b52
# ╟─a5ca5a8a-8497-41e2-9af0-92db5db9ce73
# ╟─0ddf072d-009d-42f2-9a8f-f69fbab750c6
# ╟─5db224d4-7379-4c8f-bcee-9cf00011d286
# ╠═5f62c7d2-ebfb-43e6-b9e3-06ca78c99390
# ╠═6107a3dc-26dd-4a0d-aeff-5eca2cd1dcd4
# ╠═0d33b162-cc06-4c3b-8ade-5b40106dec0e
# ╠═a7037a49-c2cb-48b6-bbf9-ca8150afcfbe
# ╠═259ea307-2862-42f3-866a-be8f4eb83cf3
# ╟─063db114-1b95-47d8-8f8b-26eaff8f9574
# ╠═8d4eb806-8af1-4127-8f72-6dd68f810eb5
# ╟─d175246c-552a-4f0f-8415-2339e9af833c
# ╠═087730c6-cf72-4133-8064-d5619ea4b188
# ╟─745bb1c1-a312-4dbd-a29f-0836b7dbe8a7
# ╠═e6033b3e-51fa-4eeb-87ff-e5b5359aebc5
# ╠═272234bf-ea90-4d24-b5f4-1cd6cdaea20b
# ╟─76275b93-5668-4ac9-a6a2-2f4b30ca8ab3
# ╠═6579518a-2130-4e25-b4eb-966e60203f17
# ╠═9ccf08d4-2ffb-4f36-a632-3b0ed4017d92
# ╠═41817b1f-4c5e-4a55-93b7-520d6b71ea9c
# ╠═0f605f2c-029b-4af8-8e52-2c7bc4c7626a
# ╠═ca09052d-1fce-422d-9205-ae4e87dc4db4
# ╠═f700edf6-0c44-4507-bbf5-4c5dc02fa74c
# ╠═de86f516-b65a-445e-8988-4cfc9dafd000
# ╠═5da650b5-9309-418b-92d9-3388b5385e9b
# ╠═7f4b9a79-ce48-487f-be37-d884ee9db0e9
# ╠═39820d4a-d05f-4dbc-8d6d-7e20b579688c
# ╠═2a45ea3b-ea3e-4d5e-912d-36b0055cf307
# ╠═7bbafd36-6130-4f00-ac86-63b2bc467b77
# ╠═759d6fd6-9649-4bea-b088-1cae1647ad3d
# ╠═019c9ab2-ef2b-4677-a5aa-d0363fffef72
# ╠═f97bb119-0248-4782-8cf9-cca61a666dbf
# ╠═cffd2b81-66c1-4f50-948c-e38ec011105d
# ╠═7785e6b1-881e-4088-82fe-3dad106b07be
# ╠═ba3660ba-e274-4e11-9fd1-75aafbb0a776
# ╠═66417309-8d6d-47e9-b5ef-7fcc0cde9194
# ╠═9cd5609b-1e58-4d39-87ab-b3d7542de691
# ╠═47ffd006-29b2-4774-a791-11b3f88a7aba
# ╟─89888da7-4912-4557-a375-0150f80ee703
# ╠═f942e39c-0d93-4bc8-8ff2-112fed566014
# ╟─df0602a2-d278-4fae-9957-30c85b55598a
# ╠═a86b2ff6-1967-424f-86eb-40fb4288c8b5
# ╠═3f2a31d4-0fa8-40fa-9dc4-bd6a26d2ddc9
# ╟─dbef7000-7c39-49f7-b24e-f0fb436eb54e
# ╠═14cb2a0b-2ea8-471b-987f-1647f1516992
# ╟─9183aa40-9dc1-4237-8bf5-de42c93b149f
# ╠═6b9beb62-dc7e-4b8b-9b7c-8fee5b1da98f
# ╟─1155e836-99d0-4cc8-83d0-355a6ab6fcc0
# ╠═1d44cf1c-11a5-4366-94f3-85b695c6ca12
# ╟─39948c5f-d472-48ce-8385-1c67b8f1580a
# ╠═f5f5a045-1c9b-40b8-b326-7d0e1336412c
# ╟─cd9ee96f-91a9-43cf-980c-6253a0c6018f
# ╠═16614e6a-2ef0-4b36-913d-6fd19440b60b
# ╠═d3e49592-cee0-4414-9b92-ce2fb66529df
# ╠═c2a099fb-446c-4a39-b505-973223c38a27
# ╠═7e68054e-4268-424e-b413-ef18baf832ac
# ╟─8160086a-6349-447c-87ae-880b02fa97f5
# ╠═3ea0415c-af14-430c-bf7c-2c7d71b7a333
# ╠═65b74f4a-e2af-4caf-8719-5a59c6349bb9
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: disberd/TelecomUtils

max_issues_repo_name: disberd/TelecomUtils

max_forks_repo_name: disberd/TelecomUtils


text:
struct Akima
    n::Int
    x::Vector{Int}
    y::Vector{Float64}
    b::Vector{Float64}
    c::Vector{Float64}
    d::Vector{Float64}
end

function Akima(x, y)
    n = length(x)
    @assert n == length(y)
    dx = diff(x)
    m = Array{Float64}(undef, n + 3)
    m[3:end-2] = diff(y) ./ dx
    m[2] = 2m[3] - m[4]
    m[1] = 2m[2] - m[3]
    m[end-1] = 2m[end-2] - m[end-3]
    m[end] = 2m[end-1] - m[end-2]

    # if m1 == m2 != m3 == m4, the slope at the breakpoint is not defined.
    # This is the fill value:
    b = 0.5 .* (m[4:end] .+ m[1:end-3])
    # get the denominator of the slope t
    dm = abs.(diff(m))
    f1 = dm[3:n+2]
    f2 = dm[1:n]
    f12 = f1 + f2
    # These are the mask of where the the slope at breakpoint is defined:
    ind = findall(f12 .> 1e-9 * maximum(f12))
    # Set the slope at breakpoint
    b[ind] = (f1[ind] .* m[ind .+ 1] .+
              f2[ind] .* m[ind .+ 2]) ./ f12[ind]
    # calculate the higher order coefficients
    c = (3.0 .* m[3:end-2] .- 2.0 .* b[1:end-1] .- b[2:end]) ./ dx
    d = (b[1:end-1] .+ b[2:end] .- 2.0 .* m[3:end-2]) ./ dx.^2

    Akima(n, x, y, b, c, d)
end

@inline function interpolate(a::Akima, xi)
    idx = searchsortedlast(a.x, trunc(Int, xi))
    idx == 0 && return a.y[1]
    idx == a.n && return a.y[end]

    wj = xi - a.x[idx]
    @evalpoly wj a.y[idx] a.b[idx] a.c[idx] a.d[idx]
end

ext: jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libtree_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libtree")
JLLWrappers.@generate_main_file("libtree", UUID("210123e1-2ab9-5e3e-92f0-8698273d9eb6"))
end  # module libtree_jll


text:
# This test follows the tutorial at:
# https://dslpitt.org/genie/wiki/SMILE_Tutorial_Code_2:_Performing_Inference_with_a_Bayesian_Network

using Smile

# create a net
net = Network()

# add a success node
success = add_node(net, DSL_CPT, "Success")
somenames = IdArray()
add(somenames, "Success")
add(somenames, "Failure")
set_number_of_outcomes(definition(get_node(net, success)), somenames)

# add a forecast node
forecast = add_node(net, DSL_CPT, "Forecast")
flush(somenames)
add(somenames, "Good")
add(somenames, "Moderate")
add(somenames, "Poor")
set_number_of_outcomes(definition(get_node(net, forecast)), somenames)

# add an arc from "Success" to "Forecast" to represent conditional independence
add_arc(net, success, forecast)

# add the distributions
theprobs = DoubleArray()
set_size(theprobs, 2)
set_at(theprobs, 0, 0.2)
set_at(theprobs, 1, 0.8)
set_definition(definition(get_node(net, success)), theprobs)

thecoordinates = SysCoordinates(definition(get_node(net, forecast)))
set_unchecked_value(thecoordinates, 0.4); next(thecoordinates)
set_unchecked_value(thecoordinates, 0.4); next(thecoordinates)
set_unchecked_value(thecoordinates, 0.2); next(thecoordinates)
set_unchecked_value(thecoordinates, 0.1); next(thecoordinates)
set_unchecked_value(thecoordinates, 0.3); next(thecoordinates)
set_unchecked_value(thecoordinates, 0.6)

# use clustering algorith
set_default_BN_algorithm(net, DSL_ALG_BN_LAURITZEN)

# suppose we want to compute P("Forecast" = Moderate)
# update the network
update_beliefs(net)

# get the result value

# DSL_sysCoordinates theCoordinates(*theNet.GetNode(forecast)->Value());
# DSL_idArray *theNames;
# theNames = theNet.GetNode(forecast)->Definition()->GetOutcomesNames();
# int moderateIndex = theNames->FindPosition("Moderate"); // should be 1
# theCoordinates[0] = moderateIndex;
# theCoordinates.GoToCurrentPosition();

# // get P("Forecast" = Moderate)
# double P_ForecastIsModerate = theCoordinates.UncheckedValue();
# printf("P(\"Forecast\" = Moderate) = %f\n",P_ForecastIsModerate);

thecoordinates = SysCoordinates(value(get_node(net, forecast)))
thenames = get_outcome_names(definition(get_node(net, forecast)))
moderateIndex = find_position(thenames, "Moderate")
@assert(moderateIndex == 1)

thecoordinates[0] = moderateIndex
go_to_current_position(thecoordinates)

# get P("Forecast" = Moderate)
P_ForecastIsModerate = unchecked_value(thecoordinates)
@printf("P(\"Forecast\" = Moderate) = %f\n", P_ForecastIsModerate)
@assert(isapprox(P_ForecastIsModerate, 0.32))

# now we want to compute P("Success" = Failure | "Forecast" = Good)
# first, introduce the evidence
# 0 is the index of state [Good]

set_evidence(value(get_node(net, forecast)), 0)

# update the network again
update_beliefs(net)

# get the result value
link_to(thecoordinates, value(get_node(net, success)))
thecoordinates[0] = 1 # 1 is the index of state [Failure]
go_to_current_position(thecoordinates)
P_SuccIsFailGivenForeIsGood = unchecked_value(thecoordinates)
@printf("P(\"Success\" = Failure | \"Forecast\" = Good) = ")
@printf("%f\n",P_SuccIsFailGivenForeIsGood)
@assert(isapprox(P_SuccIsFailGivenForeIsGood, 0.5))

# now we want to compute P("Success" = Success | "Forecast" = Poor)
# first, clear the evidence in node "Forecast"
clear_evidence(value(get_node(net, forecast)))

# introduce the evidence in node "Success"
# 2 is the index of state [Poor]
set_evidence(value(get_node(net, forecast)), 2)

# update the network again
update_beliefs(net)

# get the result value
link_to(thecoordinates, value(get_node(net, success)))
thecoordinates[0] = 0 # 0 is the index of state [Success]
go_to_current_position(thecoordinates)
P_SuccIsSuccGivenForeIsPoor = unchecked_value(thecoordinates)
@printf("P(\"Success\" = Success | \"Forecast\" = Poor) = ")
@printf("%f\n", P_SuccIsSuccGivenForeIsPoor)
@assert(abs(P_SuccIsSuccGivenForeIsPoor - 0.0769) < 1e-2)


text:
using Colors

include("helpers/page.jl")
include("helpers/doc.jl")

intro = md"""
$(title(2, "Layout"))

$(vskip(1em))

Escher provides primitives like `hbox`, `vbox`, `hskip`, `vskip`, and `flex` for laying out tiles into grids. Complex layouts can be composed from smaller parts. For higher-order layouts such as tabs, pages, menus and collapsibles see [here](layout2-api.html).

$(vskip(1em))
# Functions
$(vskip(1em))

$(
    showdocs([width, height, size, container, minwidth, minheight, maxwidth,
              maxheight, clip, pad, flow, hbox, vbox, shrink, grow, flex,
              packitems, packacross, packlines])
)
"""

function main(window)
    docpage(intro)
end


text:
# Autogenerated wrapper script for ODEInterface_jll for x86_64-apple-darwin-libgfortran4
export libbvp_m_proxy, libbvpsol, libbvpsol_i32, libcolnew, libcolnew_i32, libddeabm, libddeabm_i32, libddebdf, libddebdf_i32, libdop853, libdop853_i32, libdopri5, libdopri5_i32, libodex, libodex_i32, libradau, libradau5, libradau5_i32, libradau_i32, librodas, librodas_i32, libseulex, libseulex_i32

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("ODEInterface")
JLLWrappers.@declare_library_product(libbvp_m_proxy, "@rpath/bvp_m_proxy.dylib")
JLLWrappers.@declare_library_product(libbvpsol, "@rpath/bvpsol.dylib")
JLLWrappers.@declare_library_product(libbvpsol_i32, "@rpath/bvpsol_i32.dylib")
JLLWrappers.@declare_library_product(libcolnew, "@rpath/colnew.dylib")
JLLWrappers.@declare_library_product(libcolnew_i32, "@rpath/colnew_i32.dylib")
JLLWrappers.@declare_library_product(libddeabm, "@rpath/ddeabm.dylib")
JLLWrappers.@declare_library_product(libddeabm_i32, "@rpath/ddeabm_i32.dylib")
JLLWrappers.@declare_library_product(libddebdf, "@rpath/ddebdf.dylib")
JLLWrappers.@declare_library_product(libddebdf_i32, "@rpath/ddebdf_i32.dylib")
JLLWrappers.@declare_library_product(libdop853, "@rpath/dop853.dylib")
JLLWrappers.@declare_library_product(libdop853_i32, "@rpath/dop853_i32.dylib")
JLLWrappers.@declare_library_product(libdopri5, "@rpath/dopri5.dylib")
JLLWrappers.@declare_library_product(libdopri5_i32, "@rpath/dopri5_i32.dylib")
JLLWrappers.@declare_library_product(libodex, "@rpath/odex.dylib")
JLLWrappers.@declare_library_product(libodex_i32, "@rpath/odex_i32.dylib")
JLLWrappers.@declare_library_product(libradau, "@rpath/radau.dylib")
JLLWrappers.@declare_library_product(libradau5, "@rpath/radau5.dylib")
JLLWrappers.@declare_library_product(libradau5_i32, "@rpath/radau5_i32.dylib")
JLLWrappers.@declare_library_product(libradau_i32, "@rpath/radau_i32.dylib")
JLLWrappers.@declare_library_product(librodas, "@rpath/rodas.dylib")
JLLWrappers.@declare_library_product(librodas_i32, "@rpath/rodas_i32.dylib")
JLLWrappers.@declare_library_product(libseulex, "@rpath/seulex.dylib")
JLLWrappers.@declare_library_product(libseulex_i32, "@rpath/seulex_i32.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libbvp_m_proxy,
        "lib/bvp_m_proxy.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbvpsol,
        "lib/bvpsol.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbvpsol_i32,
        "lib/bvpsol_i32.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcolnew,
        "lib/colnew.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcolnew_i32,
        "lib/colnew_i32.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libddeabm,
        "lib/ddeabm.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libddeabm_i32,
        "lib/ddeabm_i32.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libddebdf,
        "lib/ddebdf.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libddebdf_i32,
        "lib/ddebdf_i32.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdop853,
        "lib/dop853.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdop853_i32,
        "lib/dop853_i32.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdopri5,
        "lib/dopri5.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdopri5_i32,
        "lib/dopri5_i32.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libodex,
        "lib/odex.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libodex_i32,
        "lib/odex_i32.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libradau,
        "lib/radau.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libradau5,
        "lib/radau5.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libradau5_i32,
        "lib/radau5_i32.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libradau_i32,
        "lib/radau_i32.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        librodas,
        "lib/rodas.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        librodas_i32,
        "lib/rodas_i32.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libseulex,
        "lib/seulex.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libseulex_i32,
        "lib/seulex_i32.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Autogenerated wrapper script for blaspp_jll for i686-linux-musl-cxx11
export libblaspp

using CompilerSupportLibraries_jll
using OpenBLAS32_jll
JLLWrappers.@generate_wrapper_header("blaspp")
JLLWrappers.@declare_library_product(libblaspp, "libblaspp.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, OpenBLAS32_jll)
    JLLWrappers.@init_library_product(
        libblaspp,
        "lib/libblaspp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module LDJFranklin

import Franklin
using Franklin: globvar, locvar, pagevar
using Dates
using FranklinUtils
using LDJ
using JSON
using IterTools: chain
using DataStructures: OrderedDict
using Base.Unicode: titlecase
using Memoization
using Dates: year, now
using URIs: URI

using HTTP: get
using URIs: URI
import Base.convert

const BOOKS = []

export ldj_trans, ldjfranklin_hfuncs, hfun_ldj_website, hfun_ldj_book, hfun_ldj_author, hfun_ldj_publisher, hfun_ldj_crumbs,
    hfun_ldj_place, hfun_ldj_search, hfun_ldj_webpage, hfun_ldj_webpage_tag, hfun_ldj_library, hfun_insert_library

@memoize function ldj_publisher()
    orgschema(;name=globvar(:website_title),
              url=globvar(:website_url),
              logo=globvar(:logo; default=globvar(:author_image)),
              contact="info",
              email=globvar(:email),
              sameas=[globvar(:github), globvar(:twitter)])
end
function hfun_ldj_website(k="")
    website(globvar(:website_url),
            locvar(:author), year(now());
            image=globvar(:author_image)) |>
                wrap_ldj
end

function hfun_ldj_search(k="")
    search(locvar(:website_url);
           parts=(path="/search", query=(q="{input}",))) |> wrap_ldj
end

function hfun_ldj_place(k="")
    place(;country=locvar(:country), region=locvar(:region))|> wrap_ldj
end

function hfun_ldj_author(k=""; wrap=true)
    author(;name=locvar(:author),
           email=locvar(:email),
           image=locvar(:author_image),
           description=(isdefined( @__MODULE__, :author_bio) ? author_bio() : ""),
           sameAs=[locvar(:github), locvar(:twitter)]) |>
               x -> wrap_ldj(x, wrap)
end

function hfun_ldj_publisher(k=""; wrap=true)
    hfun_ldj_author(k; wrap)
end

@memoize function ldj_get_languages()
    languages((lang for (lang, _) in locvar(:languages)))
end

function hfun_ldj_webpage()
    webpage(id=locvar(:fd_full_url),
            title=locvar(:title),
            url=locvar(:fd_full_url),
            mtime=locvar(:fd_mtime_raw),
            selector= ".franklin-content",
            description=locvar(:rss_description),
            keywords=locvar(:tags),
            access_mode=locvar(:accessMode),
            access_sufficient=locvar(:accessModeSufficient),
            access_summary="Visual elements are tentatively described.",
            image=locvar(:images; default=globvar(:author_image)),
            lang=globvar(:lang),
            created=locvar(:fd_ctime),
            props=["availableLanguage" => ldj_get_languages(),
                   "author" => hfun_ldj_author(;wrap=false),
                   "audience" => "cool people",
                   "publisher" => ldj_publisher(),
                   "mentions" => locvar(:mentions)]
                    ) |> wrap_ldj
end

function hfun_ldj_webpage_tag()
    tag = locvar(:fd_tag; default="")
    tag_path = joinpath(globvar(:tag_page_path), tag)
    mdate = joinpath(Franklin.path(:site), tag_path, "index.html") |> mtime |>  Dates.unix2datetime
    name = splitext(tag_path)[1]
    url = joinpath(globvar(:website_url), replace(name, r"(index|404)$" => ""))
    webpage(id=url,
            title=locvar(:title),
            url=url,
            mtime=mdate,
            selector= isempty(tag) ? "#tag_cloud" : ".tag-content",
            description=locvar(:rss_description),
            keywords=locvar(:fd_tag),
            access_mode=locvar(:accessMode),
            access_sufficient=locvar(:accessModeSufficient),
            access_summary="Visual elements are tentatively described.",
            image=locvar(:images; default=globvar(:author_image)),
            lang=globvar(:lang),
            created=locvar(:fd_ctime),
            props=["availableLanguage" => ldj_get_languages(),
                   "author" => hfun_ldj_author(;wrap=false),
                   "audience" => "cool people",
                   "publisher" => ldj_publisher(),
                   "mentions" => locvar(:mentions)]
                    ) |> x -> wrap_ldj(x, true, "ldj-webpage")
end

@inline function nostring(str::Union{Nothing, AbstractString})
	isnothing(str) ? "" : str
end

@doc "file path must be relative to the project directory, assumes the published website is under '__site/'"
function ldj_trans(file_path, src_url, trg_url, lang)
    translation(;src_url, trg_url, lang,
                title=nostring(pagevar(file_path, :title)),
                mtime=nostring(pagevar(file_path, :fd_mtime_raw)),
                selector=".franklin-content",
                description=pagevar(file_path, :rss_description),
                keywords=pagevar(file_path, :tags),
                image=pagevar(file_path, :images),
                headline=pagevar(file_path, :title),
                translator_name=globvar(:translator_name),
                translator_url=globvar(:translator_url),
                props=["mentions" => locvar(:mentions)]) |> wrap_ldj
end

@doc "create breadcrumbs schema for posts, requires a function to generate breadcrumbs"
function hfun_ldj_crumbs(args)
    func = getfield(Main, Symbol(args[1]))
    func() |> breadcrumbs |> x -> wrap_ldj(x, true, "ldj-breadcrumbs")
end

@doc "create a book structure"
function hfun_ldj_book(args...; kwargs...)
	book(args...; kwargs...)
end

@inline function get_date()
	let d = fr.locvar(:date)
        if d === EMPTY_DATE
            DateTime(year(now), 1, 1)
        else
            DateTime(d, DATE_FORMAT)
        end
    end
end

# function lx_book(name, author, url, sameas="", _)
function lx_book(com, _)
    args = lxproc(com)
    name, author, url, tags, comments  = strip.(args) |>
        x -> split(x, r"; ?")
    # name, author, url = split(args, "\" \"")
    push!(BOOKS, (name=name, author=author, url=url, tags=tags, comments=comments, sameas=""))
    "[" * name * "](" * url * ")"
end

convert(::Type{Symbol}, s::String) = Symbol(s)
@doc "Fetch books list from a calibre content server"
function calibredb_books_server(server="http://localhost:8099";
                                query=Dict("library_id" => "",
                                           "num" => typemax(Int32)))
    get(server * "/interface-data/books-init";  query) |>
        res -> String(res.body) |>
        JSON.parse |>
        r -> r["metadata"]
end

@doc "generate LDJ data from a list of books"
function hfun_ldj_library()
    empty!(BOOKS)
    if locvar(:fd_rpath) === "reads/index.md"
        # add calibre books to library
        calibre_books = calibredb_books_server(locvar(:calibre_server);
                                               query=Dict("library_id" => locvar(:calibre_library),
                                                          "num" => 10000))
        for (_, book) in calibre_books
            push!(BOOKS, (name=book["title"], author=book["author_sort"],
                          url="", tags=book["tags"], sameas=""))
        end
        bookfeed(BOOKS) |> wrap_ldj
    else
        ""
    end
end

@doc "create html lists by grouping books based on their tags"
function hfun_insert_library(groups=[])
    c = IOBuffer()
    if isempty(groups)
        lists = Dict(["all" => Vector{String}()])
        setgroup! = book -> begin
            "read" ∈ book.tags && push!(lists["all"], book.name)
        end
    else
        lists = OrderedDict([g => [] for g in groups])
        setgroup! = book -> begin
            for g in groups
                if "read" ∈ book.tags && g ∈ book.tags
                    push!(lists[g], book)
                end
            end
        end
    end
    for book in BOOKS
        setgroup!(book)
    end
    write(c, "
")
    for (group_name, group_list) in lists
        if length(group_list) > 0
            write(c, "
$(titlecase(group_name))
")
            write(c, "
")
            for book in group_list
                write(c, "
")
                println(c, book.name, "
 - ", book.author,"
")
                write(c, "
")
            end
            write(c, "
")
        end
    end
    write(c, "
")
    ret = String(take!(c))
    close(c)
    ret
end

# function calibre_books_cli(library="http://localhost:8099", fields="authors,isbn,publisher,tags,title")
#     @assert !isnothing(Sys.which("calibredb"))
#     # since calibre is a python package, and is installed system-wide,
#     # make sure we are not using julia python env
#     let pythonpath = ENV["PYTHONPATH"]
#         delete!(ENV, "PYTHONPATH")
#         l = read(`calibredb --with-library $library list --for-machine -f $fields`, String)
#         ENV["PYTHONPATH"] = pythonpath
#         l |> JSON.parse
#     end
# end

end


text:
using RCall

__jhitandrun__ = RCall.rimport("hitandrun")
using .__jhitandrun__

rgenerate_controls = reval("generate_controls = function(XY, XY_obs, threshold){

  XY_unknown = 1 - sign(XY_obs)

  row_constraints = matrix(0, nrow = nrow(XY_unknown), ncol = prod(dim(XY_obs)))
  for(i in 1:nrow(XY)){
    for(j in 1:ncol(XY)){
      row_constraints[i,((i-1)*ncol(XY)+j)] = XY_unknown[i,j]
    }
  }

  col_constraints = matrix(0, nrow = ncol(XY_unknown), ncol = prod(dim(XY_obs)))
  for(i in 1:nrow(XY)){
    for(j in 1:ncol(XY)){
      col_constraints[j,((i-1)*ncol(XY)+j)] = XY_unknown[i,j]
    }
  }

  A1 = rbind(row_constraints, col_constraints)
   A1 = A1[,apply(A1,2,function(x) sum(x) > 0)]
   A2 = diag(-1, ncol(A1))
   A3 = diag(1, ncol(A1))

   b1 = c(rowSums(XY) - rowSums(XY_obs), colSums(XY) - colSums(XY_obs))
   b2 = rep(-1, nrow(A2))
   b3 = rep(threshold, nrow(A3))

   dirs1 = rep('=', length(b1))
   dirs2 = rep('<=', nrow(A2))
   dirs3 = rep('<=', nrow(A3))

   A = rbind(A1, A2, A3)
   b = c(b1, b2, b3)
   dirs = c(dirs1, dirs2, dirs3)

  to_keep = apply(A,1,sum) != 0

  A = A[to_keep,]
  b = b[to_keep]
  dirs = dirs[to_keep]

  to_keep

  controls = list(constr = A, rhs = as.numeric(b), dir = dirs)

  controls

}")

rmapping = reval("mapping = function(sampling, XY, XY_obs, threshold){
  x = apply(sampling,2,mean)
  XY_unknown = 1 - sign(XY_obs)
  k=1
  for(i in 1:nrow(XY)){
    for(j in 1:ncol(XY)){
      if(XY_unknown[i,j] == 1){
        XY_obs[i,j] = x[k]
        k = k+1
      }
    }
  }

  XY_obs

}")

rminus1 = reval("minus1 = function(ctrls){
  list(constr = ctrls[['constr']][-1,], rhs = ctrls[['rhs']][-1], dir = ctrls[['dir']][-1])
}")

function method_sampling_mcmc(kx, ky, n, percentile_cutoff, rho; verbose=true)
  verbose && @info "generating data"
  data = gendata(kx, ky, n, percentile_cutoff, rho)
  verbose && @info "data has size $(size(data.gen)) with $(count(ismissing, data.obs)) missings"
  A = data.gen
  A_obs = @. ifelse(ismissing(data.obs), 0, data.obs)
  A_obs = complete_trival(A, A_obs, data.threshold)
  if any(A_obs .== 0)
      rXY = robject(A)
      rXY_obs = robject(A_obs)
      rthreshold = robject(data.threshold)
      rctrls = rcall(rgenerate_controls, XY=rXY, XY_obs = rXY_obs, threshold = rthreshold)
      rctrls1 = rcall(rminus1, ctrls = rctrls)
      try
          rsampling = rcall(shakeandbake, constr = rctrls)
          rA_obs  = rcall(rmapping, sampling = rsampling, XY = rXY, XY_obs = rXY_obs, threshold = rthreshold)
          A_obs = rcopy(rA_obs)
          mae = mean_absolute_error(data, A_obs)
          return (data=data, blanks=count(ismissing, data.obs), mae=mae, fit = A_obs, vc = cramer_v(data.gen), full=true)
      catch
          try
              rsampling = rcall(shakeandbake, constr = rctrls, eliminate = false)
              rA_obs  = rcall(rmapping, sampling = rsampling, XY = rXY, XY_obs = rXY_obs, vc = cramer_v(data.gen), threshold = rthreshold)
              A_obs = rcopy(rA_obs)
              mae = mean_absolute_error(data, A_obs)
              return (data=data, blanks=count(ismissing, data.obs), mae=mae, fit = A_obs, vc = cramer_v(data.gen), full=true)
          catch
              rsampling = rcall(shakeandbake, constr = rctrls1, eliminate = false)
              rA_obs  = rcall(rmapping, sampling = rsampling, XY = rXY, XY_obs = rXY_obs, threshold = rthreshold)
              A_obs = rcopy(rA_obs)
              mae = mean_absolute_error(data, A_obs)
              return (data=data, blanks=count(ismissing, data.obs), mae=mae, fit = A_obs, vc = cramer_v(data.gen), full = false)
          end
      end
  else
      mae = mean_absolute_error(data, A_obs)
      return (data=data, blanks=count(ismissing, data.obs), mae=mae, fit = A_obs, vc = cramer_v(data.gen), full=true)
  end
end

max_stars_repo_name: Michal-Kot/FillMissingEntries

max_issues_repo_name: Michal-Kot/FillMissingEntries

max_forks_repo_name: Michal-Kot/FillMissingEntries


text:
#=
Normal SWAP-Test algorithm related Functions that help verify the validity of MSQR training results.
=#
export SWAPtest, SWAPtrain!, SWAPtestRes

"""
    SWAPtestRes{overlap::Float64, witnessOp::PutBlock, reg::ArrayReg, circuit::ChainBlock}
Fields:
\n`overlap::Float64`:     Overlap bewteen generated register and target register.
\n`witnessOp::PutBlock`:  The witness(measure) operator of the `SWAPtest` circuit.
\n`reg::ArrayReg`: The register goes through `SWAPtest`. nqubits(SWAPtest.reg) = 1 + 2*nBitT.
\n`circuit::ChainBlock`:  Circuit of `SWAPtest`.
"""
struct SWAPtestRes
    overlap::Float64     # Overlap bewteen generated register and target register.
    witnessOp::PutBlock  # The witness(measure) operator of the `SWAPtest` circuit.
    reg::ArrayReg # The register goes through `SWAPtest`. nqubits(SWAPtest.reg) = 1 + 2*nBitT.
    circuit::ChainBlock  # Circuit of `SWAPtest`.
end

"""
    SWAPtest(regG::ArrayReg, regT::ArrayReg; nMeasure::Int64=1, ϕ::Real=0, useCuYao::Bool=CUDA_ON)
    ->
    SWAPtestRes{overlap::Float64, witnessOp::PutBlock, reg::ArrayReg, circuit::ChainBlock}
SWAP test function which get the overlap between target register(regT) and generated register(reg0).
\n`regG::ArrayReg`: Generated quantum state to compare with target quantum state.
\n`regT::ArrayReg`: Target quantum state.
"""
function SWAPtest(regG::ArrayReg, regT::ArrayReg; nMeasure::Int64=1, ϕ::Real=0, useCuYao::Bool=CUDA_ON)
    nMcheck(nMeasure, regG, regT)
    reg1 = repeat(copy(regG), nMeasure)
    reg2 = repeat(copy(regT), nMeasure)
    reg0 = zero_state(1, nbatch=nMeasure)
    regA = join(reg0, reg1, reg2)
    useCuYao && (regA = regA |> CuYao.cu)
    resS = SWAPtest(regA, ϕ=ϕ)
    regA = resS.reg
    useCuYao && (regA = regA |> CuYao.cpu)
    SWAPtestRes(resS.overlap, resS.witnessOp, regA, resS.circuit)
end
"""
    SWAPtest(regA::ArrayReg; nMeasure::Int64=1, ϕ::Real=0)
    ->
    SWAPtestRes{overlap::Float64, witnessOp::PutBlock, reg::ArrayReg, circuit::ChainBlock}
Method 2 of `SWAPtest`:
\n`regA::ArrayReg`: Input register(state) for SWAP Test circuit. `regA = join(zero_state(1), reg1, reg2)`
"""
function SWAPtest(regA::ArrayReg; nMeasure::Int64=1, ϕ::Real=0)
    regA = copy(repeat(regA, nMeasure))
    nBitA = nqubits(regA)
    nBitT = Int((nBitA - 1) / 2)
    witnessOp = put(nBitA, nBitA=>Z)
    SWAPblock = chain(nBitA, [control(nBitA, nBitA, (i,i-nBitT)=>SWAP) for i=nBitA-1:-1:nBitT+1])
    SWAPcircuit = chain(nBitA, put(nBitA, nBitA=>H), put(nBitA, nBitA=>shift(ϕ)), SWAPblock, put(nBitA, nBitA=>H))
    regA |> SWAPcircuit
    overlaps = expect(witnessOp, regA)
    overlap = mean(overlaps) |> real
    SWAPtestRes(overlap, witnessOp, regA, SWAPcircuit)
end
"""
    SWAPtest(regTar::ArrayReg, circuit::ChainBlock; nMeasure::Int64=1, ϕ::Real=0, useCuYao::Bool=CUDA_ON)
    ->
    SWAPtestRes{overlap::Float64, witnessOp::PutBlock, reg::ArrayReg, circuit::ChainBlock}
Method 3 of `SWAPtest`:
\n`regTar::ArrayReg`: Target quantum state.
\n`circuit::ChainBlock`: circuit to generate the state(regG) to compare with `regTar`.
"""
function SWAPtest(regTar::ArrayReg, circuit::ChainBlock; nMeasure::Int64=1, ϕ::Real=0, useCuYao::Bool=CUDA_ON)
    n = nqubits(regTar)
    regT = copy(regTar)
    regG = zero_state(nqubits(regT)) |> circuit
    SWAPtest(regG, regT, nMeasure=nMeasure, ϕ=ϕ, useCuYao=useCuYao)
end
function SWAPtest(circuit::ChainBlock; regAll::ArrayReg, nMeasure::Int64=1, ϕ::Real=0)
    nMcheck(nMeasure, regAll)
    regA = copy(regAll)
    nBitT = Int((nqubits(regA) - 1) / 2)
    regA |> focus!(nBitT+1:2nBitT...) |> circuit |> relax!(nBitT+1:2nBitT...)
    SWAPtest(regA, nMeasure=nMeasure, ϕ=ϕ)
end

"""
    SWAPtrain!(regTar::ArrayReg, circuit::ChainBlock, nTrain::Union{Int64, :auto}; nMeasure::Int64=1, Gmethod::Union{String, Tuple{String,Float64}}="Qdiff", GDmethod=("default",0.01), show::Bool=false, useCuYao::Bool=CUDA_ON, ConvTh::Tuple{Float64, Float64}=(5e-4, 1e-3))
    ->
    overlaps::Array{Float64,1}
SWAP-Test training function. This function will change the parameters of differentiable gates in circuit. When set `nTrain = :auto`, trigger the automaic training ieration.
"""
function SWAPtrain!(regTar::ArrayReg, circuit::ChainBlock, nTrain::Union{Int64, Symbol}; nMeasure::Int64=1,
                    Gmethod::Union{String, Tuple{String,Float64}}="Qdiff", GDmethod=("default",0.01),
                    show::Bool=false, useCuYao::Bool=CUDA_ON, ConvTh::Tuple{Float64, Float64}=(5e-4, 1e-3))
    if show
        cPar = MPSDCpar(circuit)
        nBitT = cPar.nBitA
        vBit = cPar.vBit
        rBit = cPar.rBit
        depth = cPar.depth
        println("\nSWAP Training Parameters:")
        println("nBitT=$(nBitT) vBit=$(vBit) rBit=$(rBit) depth=$(depth)")
        println("nMeasure=$(nMeasure) nTrain=$(nTrain) GDmethod=$(GDmethod)\n")
        println("Initial overlap = $(SWAPtest(regTar, circuit, nMeasure=nMeasure, useCuYao=useCuYao).overlap)")
    end
    reg = zero_state(nqubits(regTar)+1, nbatch = nMeasure)
    regT = repeat(regTar, nMeasure)
    regA = copy(join(reg, regT))
    useCuYao == true && (regA = regA |> cu)
    if typeof(nTrain) == Int64
        res = train!(nTrain, circuit, Tmethod = circuit->SWAPtest(circuit, regAll = regA),
                     Gmethod=Gmethod, GDmethod=GDmethod, show=show)
    elseif nTrain == :auto
        res = train!(circuit, Tmethod = circuit->SWAPtest(circuit, regAll = regA),
                     Gmethod=Gmethod, GDmethod=GDmethod, show=show, ConvTh=ConvTh)
    end
    res
end


text:
using Thebes, Luxor, Random

include(dirname(pathof(Thebes)) * "/../data/moreobjects.jl")

function mygfunction(vertices, faces, labels, action=:stroke)
    if !isempty(faces)
        @layer begin
            for (n, p) in enumerate(faces)
                sethue("white")
                poly(p, action)
            end
        end
    end
end

function drawgeodesic(object, cpos, cscale, rotx, roty, rotz, eased)
    eyepoint(Point3D(1200, 1200, 200))

    c = scaleby!(object, cscale.x, cscale.y, cscale.z)
    moveby!(c, cpos)
    theta = rescale(eased, 0, 1, 0, 2pi)
    rotateby!(c, Point3D(0, 0, 0), theta, theta, theta)
    pin(c, gfunction= mygfunction)
end

function backdrop(scene, framenumber)
    pl = box(O, scene.movie.width, scene.movie.height, vertices=true)
    # start at bottom left
    mesh1 = mesh(pl, [
    "midnightblue",
        "grey80",
        "grey70",
        "midnightblue",
    ])
    setmesh(mesh1)
    poly(pl, :fill)
end

function frame1(scene, framenumber)
    sethue("black")
    setline(0.5)
    setlinejoin("bevel")
    setopacity(0.6)
    eased_n = scene.easingfunction(framenumber, 0, 1, scene.framerange.stop)
    # object, position, scale, rotation
    drawgeodesic(deepcopy(object), Point3D(0, 0, 0), Point3D(350, 350, 350), 0, 0, 0, eased_n)
end

geodesicmovie = Movie(500, 500, "geodesic")
object = sortfaces!(make(geodesic, "geodesic"))

animate(geodesicmovie, [
    Scene(geodesicmovie, backdrop, 1:50),
    Scene(geodesicmovie, frame1, 1:50, easingfunction=easeinoutsine)], creategif=true,
    pathname="/tmp/geodesic.gif")


text:
using RoutingNetworks
using Base.Test

# write your own tests here
n = squareNetwork(3)

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


text:
# # Conway's game of life

# ![](game_of_life.gif)

# [Game of life on wikipedia](https://en.wikipedia.org/wiki/Conway%27s_Game_of_Life).

# It is also available from the `Models` module as [`Models.game_of_life`](@ref).

using Agents, Plots
using Random # hide
gr(); # hide
nothing # hide

# ## 1. Define the rules

# Rules of Conway's game of life: DSRO (Death, Survival, Reproduction, Overpopulation).
# Cells die if the number of their living neighbors is <D or >O,
# survive if the number of their living neighbors is ≤S,
# come to life if their living neighbors are  ≥R and ≤O.
rules = (2, 3, 3, 3)
nothing # hide

# ## 2. Build the model

# First, define an agent type. It needs to have the compulsary `id` and `pos` fields, as well as an `status` field that is `true` for cells that are alive and `false` otherwise.

mutable struct Cell <: AbstractAgent
    id::Int
    pos::Dims{2}
    status::Bool
end

# The following function builds a 2D cellular automaton. `rules` is of type `Tuple{Int,Int,Int,Int}` representing DSRO.

# `dims` is a tuple of integers determining the width and height of the grid environment.
# `metric` specifies whether cells connect to their diagonal neighbors.

# This function creates a model where all cells are "off".

function build_model(; rules::Tuple, dims = (100, 100), metric = :chebyshev, seed = 120)
    space = GridSpace(dims; metric)
    properties = Dict(:rules => rules)
    model = ABM(Cell, space; properties, rng = MersenneTwister(seed))
    idx = 1
    for x in 1:dims[1]
        for y in 1:dims[2]
            add_agent_pos!(Cell(idx, (x, y), false), model)
            idx += 1
        end
    end
    return model
end
nothing # hide

# Now we define a stepping function for the model to apply the rules to agents.

function ca_step!(model)
    new_status = fill(false, nagents(model))
    for agent in allagents(model)
        nlive = nlive_neighbors(agent, model)
        if agent.status == true && (nlive ≤ model.rules[4] && nlive ≥ model.rules[1])
            new_status[agent.id] = true
        elseif agent.status == false && (nlive ≥ model.rules[3] && nlive ≤ model.rules[4])
            new_status[agent.id] = true
        end
    end

    for k in keys(model.agents)
        model.agents[k].status = new_status[k]
    end
end

function nlive_neighbors(agent, model)
    neighbor_positions = nearby_positions(agent, model)
    all_neighbors =
        Iterators.flatten(ids_in_position(np, model) for np in neighbor_positions)
    sum(model[i].status == true for i in all_neighbors)
end
nothing # hide

# now we can instantiate the model:
model = build_model(rules = rules, dims = (50, 50))

# Let's make some random cells on
for i in 1:nagents(model)
    if rand(model.rng) < 0.2
        model.agents[i].status = true
    end
end

# ## 3. Animate the model

# We use the [`plotabm`](@ref) function to create an animation.

cd(@__DIR__) #src
ac(x) = x.status == true ? :black : :white
anim = @animate for i in 0:100
    i > 0 && step!(model, dummystep, ca_step!, 1)
    p1 = plotabm(model; ac = ac, as = 3, am = :square, showaxis = false)
end
nothing # hide

# We can now save the animation to a gif.

gif(anim, "game_of_life.gif", fps = 5)

max_stars_repo_name: seantfish/chatbot

max_issues_repo_name: seantfish/chatbot

max_forks_repo_name: seantfish/chatbot


text:
using CameraModels
using Test
using StaticArrays

struct SomeTestModel <: CameraModel end

@testset "Test sensorsize using rows and columns." begin
    CameraModels.rows(m::SomeTestModel) = 11
    CameraModels.columns(m::SomeTestModel) = 22

    model = SomeTestModel()
    @test sensorsize(model) == SVector{2}(22,11)
end

include("CameraModelTestBench.jl")

include("Pinhole.jl")


text:
@testset "ConvertEltype" begin
    @test (ConvertEltype <: Augmentor.AffineOperation) == false
    @test (ConvertEltype <: Augmentor.ImageOperation) == false
    @test (ConvertEltype <: Augmentor.Operation) == true

    @testset "constructor" begin
        @test_throws MethodError ConvertEltype()
        @test typeof(@inferred(ConvertEltype(Float64))) <: ConvertEltype <: Augmentor.Operation
        @test typeof(@inferred(ConvertEltype(RGB))) <: ConvertEltype <: Augmentor.Operation
        @test typeof(@inferred(ConvertEltype(RGB{N0f8}))) <: ConvertEltype <: Augmentor.Operation
        @test str_show(ConvertEltype(Float64)) == "Augmentor.ConvertEltype(Float64)"
        @test str_show(ConvertEltype(RGB)) == "Augmentor.ConvertEltype(RGB{Any})"
        @test str_show(ConvertEltype(Gray{N0f8})) == "Augmentor.ConvertEltype(Gray{N0f8})"
        @test str_showconst(ConvertEltype(Float64)) == "ConvertEltype(Float64)"
        @test str_showconst(ConvertEltype(RGB{N0f8})) == "ConvertEltype(RGB{N0f8})"
        @test str_showcompact(ConvertEltype(Float64)) == "Convert eltype to Float64"
        @test str_showcompact(ConvertEltype(Gray)) == "Convert eltype to Gray{Any}"
    end
    @testset "eager" begin
        @test Augmentor.supports_eager(ConvertEltype) === true
        @test Augmentor.supports_eager(ConvertEltype{Float64}) === true
        img_out = convert(Array{Gray{Float32}}, rect)
        imgs = [
            (Float32, rect),
            (Float32, OffsetArray(rect, -2, -1)),
            (Gray{Float32}, rect),
            (Gray{Float32}, Float64.(rect)),
            (Gray{Float32}, reshape(view(rect,:,:), 2,3)),
            (Gray{Float32}, RGB{N0f8}.(rect)),
            (Gray{Float32}, Augmentor.prepareaffine(rect)),
            (Gray{Float32}, OffsetArray(rect, -2, -1)),
            (Gray{Float32}, view(rect, IdentityRange(1:2), IdentityRange(1:3))),
            (RGB{Float32}, rect),
            (RGB{Float32}, OffsetArray(rect, -2, -1))
        ]
        @testset "single image" begin
            for (T, img_in) in imgs
                res = @inferred(Augmentor.applyeager(ConvertEltype(T), img_in))
                out = T.(img_out)
                @test res ≈ out
                @test typeof(res) <: typeof(out)
            end
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(ConvertEltype) === false
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(ConvertEltype) === false
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(ConvertEltype) === true
        @test @inferred(Augmentor.supports_lazy(ConvertEltype{Float64})) === true
        @test @inferred(Augmentor.supports_lazy(typeof(ConvertEltype(Gray)))) === true
        @test @inferred(Augmentor.supports_lazy(typeof(ConvertEltype(Gray{N0f8})))) === true
        let img = @inferred(Augmentor.applylazy(ConvertEltype(Gray{Float32}), OffsetArray(rect,-2,-1)))
            @test parent(parent(img)) === rect
            @test typeof(img) <: ReadonlyMappedArray{Gray{Float32},2}
            @test axes(img) == (OffsetRange(-1:0), OffsetRange(0:2))
            @test img[0,0] isa Gray{Float32}
            @test collect(img) == convert.(Gray{Float32}, rect)
        end
        let img = @inferred(Augmentor.applylazy(ConvertEltype(Gray{Float32}), Augmentor.prepareaffine(rect)))
            @test parent(parent(parent(parent(img)))) === rect
            @test typeof(img) <: ReadonlyMappedArray{Gray{Float32},2}
            @test axes(img) === (1:2, 1:3)
            @test img[1,1] isa Gray{Float32}
            @test collect(img) == convert.(Gray{Float32}, rect)
        end
        let img = @inferred(Augmentor.applylazy(ConvertEltype(Gray{Float32}), view(rect, IdentityRange(1:2), IdentityRange(1:3))))
            @test parent(parent(img)) === rect
            @test typeof(img) <: ReadonlyMappedArray{Gray{Float32},2}
            @test axes(img) === (1:2, 1:3)
            @test img[1,1] isa Gray{Float32}
            @test collect(img) == convert.(Gray{Float32}, rect)
        end
        let img = @inferred(Augmentor.applylazy(ConvertEltype(Gray{Float32}), rgb_rect))
            @test parent(img) === rgb_rect
            @test axes(img) === (Base.OneTo(2), Base.OneTo(3))
            @test typeof(img) <: ReadonlyMappedArray{Gray{Float32},2}
            @test img == convert.(Gray{Float32}, rgb_rect)
        end
        let img = @inferred(Augmentor.applylazy(ConvertEltype(Float32), checkers))
            @test parent(img) === checkers
            @test axes(img) === (Base.OneTo(3), Base.OneTo(5))
            @test typeof(img) <: ReadonlyMappedArray{Float32,2}
            @test img == convert(Array{Float32}, checkers)
        end
        let img = @inferred(Augmentor.applylazy(ConvertEltype(RGB{N0f8}), checkers))
            @test parent(img) === checkers
            @test typeof(img) <: ReadonlyMappedArray{RGB{N0f8},2}
            @test img == convert.(RGB, checkers)
        end
        let img = @inferred(Augmentor.applylazy(ConvertEltype(RGB{Float64}), checkers))
            @test parent(img) === checkers
            @test typeof(img) <: ReadonlyMappedArray{RGB{Float64},2}
            @test img == convert.(RGB{Float64}, checkers)
        end
    end
    @testset "view" begin
        @test Augmentor.supports_view(ConvertEltype) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(ConvertEltype) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(ConvertEltype) === false
    end
end


text:
using Documenter, SparseRegression

makedocs(
    format = Documenter.HTML(),
    sitename = "SparseRegression.jl",
    authors = "Josh Day",
    clean = true,
    pages = [
        "index.md",
        "usage.md",
        "algorithms.md"
    ]
)

deploydocs(
    repo   = "github.com/joshday/SparseRegression.jl.git",
)


text:
# This file was generated, do not modify it. # hide
model = SimplePipe()
lambdas = range(model, :(clf.lambda), lower=1e-3, upper=100, scale=:log10)
tm = TunedModel(model=SimplePipe(), ranges=lambdas, measure=cross_entropy)
mtm = machine(tm, Xtrain, ytrain)
fit!(mtm)
best_pipe = fitted_params(mtm).best_model


text:
using Chatter
using Base.Test

struct Foo <: Chatter.Source end
struct BarMessage <: Chatter.Message
    msg::AbstractString
    BarMessage(msg) = new(msg)
end

client = ChatterClient(Foo())
for message in ["A", "B", "C", "D", "E", "F", "G"]
    push!(client.messages, BarMessage(message))
end

@test map(m -> m.msg, get!(client)) == ["A"]
@test map(m -> m.msg, get!(client; n=2)) == ["B", "C"]
@test map(m -> m.msg, flush!(client)) == ["D", "E", "F", "G"]
@test isempty(client.messages)


text:
#=
NOTICE:
This file cannot just be run as a script as there
are required manual steps listed in the comments.
=#

# Define a method for a function
# Notice that this is a *very* simple function
# so we are not dealing with linking in this sample

addone(x) = x + 1

# Trigger a specialization to be created

precompile(addone, Tuple{Int64})

# Generate and output LLVM textual IR to a file

open("sample.ll", "w") do io
    code_llvm(io, addone, (Int64,); dump_module=true)
end

# Take a look at the file.
# You will notice two functions getting defined.
# The one named jfptr_addone_somenumber is the wrapper.
# The one named julia_addone_someothernumber is the real one.
# Both are needed.
# You may want to fiddle with the literal number added
# by the julia_addone_somenumber function so you can tell
# later that you are running the shared library code
# rather than the jit code.

# In your shell convert the LLVM ll code to bc format
# You will have to figure out what your julia build directory is

julia_build_dir/usr/tools/llvm-as sample.ll -o sample.bc
julia_build_dir/usr/tools/clang -shared -fpic sample.bc -o sample.so

# Dig down in the metadata of the function
# to get to the method specializations

s = methods(addone).ms[1].specializations

# The cache is a linked list of code instances
# Since we only caused one instance to be created
# we know we are looking at the right one ;)

c = s.func.cache

# Here are the two function pointers needed
# to call the method with normal julia code

c.invoke
c.specptr

# If we had a native compiled shared library we could load it
using Libdl

# Get pointers to the functions
# (Look in the sample.ll file created aboe to get the numbers)

so = dlopen("somepath", someflags)
invoke_fptr = dlsym(so, "jfptr_addone_somenumber")
specptr_fptr = dlsym(so, "julia_addone_someothernumber")

# Save the pointers into the code instance from above

c.invoke = invoke_fptr
c.specptr = specptr_fptr

# Call the version of the method from the shared library

addone(Int64(1))

max_stars_repo_name: TimoLarson/NativeCache

max_issues_repo_name: TimoLarson/NativeCache

max_forks_repo_name: TimoLarson/NativeCache


text:
using CairoMakie
x(t; a = 5, b = 7, c = 1.5) = (a - b)*cos(t) + c*cos(t*(a/b - 1))
y(t; a = 5, b = 7, c = 1.5) = (a - b)*sin(t) - c*sin(t*(a/b - 1))
t = 0:0.001:20π
with_theme(theme_dark()) do
    fig, = lines(x.(t), y.(t); color = ("#08F7FE",0.8), linewidth = 0.8,
        axis = (; aspect = DataAspect(), xlabel = "x", ylabel = "y"),
        figure = (; resolution = (600,600)))
    [lines!(x.(t), y.(t); color = ("#08F7FE", 0.08), linewidth = 1 + 0.8*i)
        for i in 1:13]
    save("lines_hypotrochoid.png", fig, px_per_unit = 2)
    display(fig)
end