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

using Markdown
using InteractiveUtils

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

# ╔═╡ 9bc349e8-efb0-11ea-1392-cd735828af33
# packaged loads, etc
begin
	using Colors
	using DelimitedFiles
	using NIfTI
	using Images
	using ImageView
	using StatsBase
	using Plots
	using PlutoUI
	md"""loading packages..."""
end

# ╔═╡ 47a60d48-058e-11eb-3140-03a5c22a27ee
html"<button onclick=present()>Present</button>"

# ╔═╡ ab3afe3a-0684-11eb-090d-ffeea5fcadec
html"""

"""

# ╔═╡ 730d087c-067b-11eb-0427-e9425b7c1ed6
md"""# Images and colourmaps"""

# ╔═╡ 8029070e-efb0-11ea-1d0a-7f402f2e9989
md"""## What are images and colourmaps?

- a brief explanations of RGB colours
- what are **intensity images**
- how can be superimpose two images by using different colour maps?

"""

# ╔═╡ 06433d72-058b-11eb-0565-c3e0f393ef5d
md"""## How do colours work on a computer?"""

# ╔═╡ 38bf1198-049f-11eb-0b21-334cd2fe39d3
md"""Colours used on computer displays are often defined in terms of the proportion of **red**, **green**, and **blue** the contain """

# ╔═╡ 542b0706-058b-11eb-0ae3-d57aa9ead256
[ RGB(1, 0, 0), RGB(0, 1, 0), RGB(0, 0, 1)]

# ╔═╡ bfb16300-058e-11eb-38d2-3572d8db913f
md"""## Controlling strength of each colour"""

# ╔═╡ ec3816c4-058b-11eb-2d8b-a1ed496f153d
md"""we can change the amount of each by multiplying by a number between 0 and 1

multiplying each color by a number varies the **amount**: $(@bind colorMult Slider(0:0.1:1, default = 1.0, show_value=true)) """

# ╔═╡ 0e5ddcca-058c-11eb-3068-73441b317478
[ RGB(colorMult, 0, 0), RGB(0, colorMult, 0), RGB(0, 0, colorMult)]

# ╔═╡ cb5f17f4-067b-11eb-3f8b-3ff5d0aa088c
md"""
- sometimes the colour values are displays as fractions $[0..1]$
- on the computer, those values are often stored between $[0..255]$ (**8 bits**)
"""

# ╔═╡ 62d614fc-058f-11eb-1fc0-09e86807181e
md"""## Mixing RGB creates unique colours """

# ╔═╡ 8baff8a4-058f-11eb-0882-9f305cf65b7c
md"""We can create many different colors by independently changing values and then adding:

The following proportions of

`red` $(@bind redM Slider(0:0.1:1, default = 0.5, show_value=true))

`green` $(@bind greenM Slider(0:0.1:1, default = 0, show_value=true))

`blue` $(@bind blueM Slider(0:0.1:1, default = 0, show_value=true))

lead the following colour:
"""

# ╔═╡ cc6e657e-058f-11eb-2314-5b6df2f908e9
begin
	cPatch =  [RGB(redM,greenM, blueM) for i=range(0,1, length=2), j= range(0,1, length=2)]
	heatmap(cPatch,
		xticks=:none,
		yticks=:none, title="RGB = $([redM, greenM, blueM])")

end

# ╔═╡ 70e80930-04a1-11eb-055b-152d2006d067
md"""## Example: mixing red with...

number of steps in color map $(@bind nSteps Slider(5:5:15, show_value=true)) """ # number of steps in ramp"""

# ╔═╡ 1b79dd34-041d-11eb-037a-93335f4b4116
r = [RGB(i, 0, 0) for i=range(0,1, length=nSteps), j= range(0,1, length=nSteps)];

# ╔═╡ 50856a72-05ba-11eb-0f1f-7dc909270ad5
heatmap(r, xticks=:none, yticks=:none)

# ╔═╡ 82d57e14-04a1-11eb-2190-9daace12efa8
md"""## ... more green/blue

results in a grid of colours that have more red $\downarrow$ and

1. more green towards $\rightarrow$ right
2. more blue towards $\rightarrow$ right"""

# ╔═╡ 6f3591f2-049f-11eb-2a97-ad43e34fb000
begin
	g = [RGB(0,j, 0) for i=range(0,1, length=nSteps), j= range(0,1, length=nSteps)]
	b = [RGB(0,0, j) for i=range(0,1, length=nSteps), j= range(0,1, length=nSteps)]
	rg = [RGB(i,j, 0) for i=range(0,1, length=nSteps), j= range(0,1, length=nSteps)]
	rb = [RGB(i,0, j) for i=range(0,1, length=nSteps), j= range(0,1, length=nSteps)]
	heatmap(g, xticks=:none, yticks=:none, layout=(2,2), subplot=1,
		title="green")
	heatmap!(b, xticks=:none, yticks=:none, layout=(2,2), subplot=2,
		title="blue")
	heatmap!(rg, xticks=:none, yticks=:none, layout=(2,2), subplot=3,
		title="red + green")
	heatmap!(rb, xticks=:none, yticks=:none, layout=(2,2), subplot=4,
		title="red + blue")

end

# ╔═╡ 79e91b4e-0654-11eb-2fac-777bdb7c41f5
rgbTable = [RGB(i,j, k) for i=range(0,1, length=nSteps),
		j=range(0,1, length=nSteps),
		k= range(0,1, length=nSteps)];

# ╔═╡ 0d3ab51c-0590-11eb-1d21-1f6fde64b56f
md""" ## So $\rightarrow$ 3 numbers to specify a colour

### Question

What if we only want have 1 number for each pixel (or voxel in three dimensions)?

For example, the signal **intensity** from magnetic resonance imaging (MRI) is a single number at each point.

### Answer

Change `r`, `g`, and `b` values by same amount... this results in

[black ... grey ... white]
"""

# ╔═╡ 3c37796e-058c-11eb-330f-c9f07956d74e
md""" ## grey """

# ╔═╡ 2a2a7d84-058f-11eb-27bb-db743f9910d5
md"""$(@bind greyMult Slider(0:0.1:1, default = 0.7, show_value=true)) """

# ╔═╡ 69bc09c0-058c-11eb-36a6-bd6e942b5e08
md"""... if we mix `R` , `G` and `B` together in equal amounts, we get GREY"""

# ╔═╡ 81adabb0-058c-11eb-3ee6-959bbf5f0c5d
 begin
	rgb = [RGB(greyMult,greyMult, greyMult) for i=range(0,1, length=2), j= range(0,1, length=2)]
	heatmap(rgb, xticks=:none, yticks=:none, title="RGB = $(repeat([greyMult],1,3))")
end

# ╔═╡ 266fa81e-067e-11eb-19f6-37bcc6c623e3
md"""## Colourmaps?!

The rule for which colour goes with which value...

$\begin{matrix}
0 & 1 & \dots & 254 & 255 \\
\text{black} & \text{dark grey} & \dots & \text{light grey} & \text{white}
\end{matrix}$

"""

# ╔═╡ 1c814864-0686-11eb-1c6f-2f0fa3d874cf
md"""or **reds**"""

# ╔═╡ c8c8b590-0685-11eb-037a-d17f5e4c3ba4
reverse(cgrad(:reds))

# ╔═╡ 26100fbe-0686-11eb-10eb-6752e71e2481
md"""or **jet** """

# ╔═╡ 3a80a800-0686-11eb-2da5-49c8cb23254d
cgrad(:jet1)

# ╔═╡ 8b8056b0-0686-11eb-15ae-c32b9e0034a2
md"""or **hsv**, **viridis**, ..."""

# ╔═╡ c00be5ca-0686-11eb-1a62-c75bb04147df
cgrad(:hsv)

# ╔═╡ c4cd3778-0686-11eb-169e-57ca0d655f5b
cgrad(:viridis)

# ╔═╡ 46874f50-0686-11eb-0e2b-8f54c14e5e8a
md"""[a whole list of them](http://docs.juliaplots.org/latest/generated/colorschemes/#scientific) online
"""

# ╔═╡ 7f61ac58-0684-11eb-182a-ad7903933dd5
md"""## Displaying a slice w/ different colourmaps """

# ╔═╡ 49141f5e-06fa-11eb-1cb0-7ffc29a11b48
md"""## Superimposing anatomy + statistical result

in practice, we often want to merge two images with different color maps, say *gray* (for anatomical imformation) and *inferno* (for statistical results)

"""

# ╔═╡ 1386a3c6-06fb-11eb-2abe-29cfbf505e0e
stat_image = niread("thresh_zstat3.nii");

# ╔═╡ f30741d0-06fa-11eb-290b-375c807e13c6
bg_image = niread("tmean.nii");

# ╔═╡ 3ebec912-06fc-11eb-23de-cdb28471c46d
md"""slice: $(@bind someSlice Slider(1:size(bg_image,3), default = 12, show_value=true))
"""

# ╔═╡ 9c1b27d4-0684-11eb-235b-eb019285d84d
html"""

"""

# ╔═╡ 27115482-0682-11eb-19d8-c344581ec8af
colormap_names = ["greys","reds", "jet", "viridis", "hsv", "inferno"];

# ╔═╡ 7d531b12-0682-11eb-3417-4f641024d076
@bind sliceColor Select(colormap_names)

# ╔═╡ 4aa28d70-067c-11eb-1bfc-179c255fcdd0
md"""## Some fun notes...

somet things you can do with `Colors.jl` in a Pluto notebook.

"""

# ╔═╡ 6381fc42-0654-11eb-3d50-9fd2b4d241b7
md"""## the whole table of r, g, b...

mixing $(nSteps) red by $(nSteps) green by $(nSteps) blue results in $(nSteps.^3) colour tiles """

# ╔═╡ 9d530386-0654-11eb-13ea-4dbc0c7d9f03
begin
	allColors2d = reshape(rgbTable, nSteps, nSteps*nSteps )
	heatmap(allColors2d, xticks=:none, yticks=:none)
end

# ╔═╡ d12a8c42-0680-11eb-225d-f9d14982eef4
md"""## Showing colors and fMRI response

two different views of the same timecourse data:

- time (in volumes) runs from top $\downarrow$
- left view: fMRI response in colors (dark: low, bright: high)
- right view: a plot (just rotated by 90$^\circ$)

"""

# ╔═╡ 8fa5f8da-04a1-11eb-2941-65dd79fdba84
md"""### Colour picking in notebook """

# ╔═╡ 9d485144-0588-11eb-3996-0fc8abf6ed9d
 @bind someColor html"""<input type="color" id="favcolor" name="favcolor" value="#ff0000">"""

# ╔═╡ bd6e4a08-0588-11eb-0a97-d1d98b30d894
md"""the colour is stored a variable in **hexadecimal format: $(someColor)"""

# ╔═╡ c35d782a-04a0-11eb-2d29-296d9f33ecd9
md"""### Named colors in the `Colors.jl` package """

# ╔═╡ 409aa96c-04a0-11eb-377a-2732f7fecbaa
[keys(Colors.color_names)...]

# ╔═╡ 0b4e44d4-04a1-11eb-3e57-53f8cdb029a7
# get the color names out of the dictionary
color_names = sort([keys(Colors.color_names)...]);

# ╔═╡ 2c6ad37c-04a0-11eb-28b7-8b55c2f60cb3
@bind the_color Select(color_names);

# ╔═╡ f3661c44-049f-11eb-1ff7-6f0d6d96d8fe
RGB(Colors.color_names[the_color]./255...)

# ╔═╡ 9f76bc18-0421-11eb-2d94-f7aa0f2e0998
md"""## math typsetting etc

an equation

\begin{equation}
	\sum_{i=0}^{\infty}\frac{1}{x^i}
\end{equation}
"""

# ╔═╡ 44de79f0-f450-11ea-0f93-15cb49d19311
md"""## An interactive 3d viewer"""

# ╔═╡ 2b2f6be4-f07f-11ea-1459-775749ca949f
md"""## Some helper functions """

# ╔═╡ 141dde62-efb1-11ea-34df-cd5e838d8129
begin
	# read data in from file
	data = NIfTI.niread("filtered_func_data.nii");
	md"""loading in `data` """
end

# ╔═╡ a8ab8ee6-f083-11ea-1316-7fbe4e125914
begin
md"""
x  $(@bind xLoc Slider(1:size(data,1), default=42, show_value=true));
y  $(@bind yLoc Slider(1:size(data,2), default=10, show_value=true)); z  $(@bind zLoc Slider(1:size(data,3), default=8, show_value=true))
"""
end

# ╔═╡ b247bd76-f07e-11ea-1e07-577643ff20af
begin
	heatmap(repeat(data[xLoc,yLoc,zLoc,:],1,1),
		c=:greys,
		layout=(1,2), subplot=1,
	title = "timecourse at [$(xLoc), $(yLoc), $(zLoc)]")
	xaxis!(showaxis=:y, subplot=1)
	yaxis!(flip=true, ylabel="time (volumes)")
	plot!(data[xLoc,yLoc,zLoc,:], 1:size(data,4) ,subplot=2)

	xaxis!(xlabel="fMRI response", xflip=false, subplot=2)
	yaxis!(showaxis=:xy, yflip=true, subplot=2)
end

# ╔═╡ f0ef494e-f80c-11ea-1163-b954c1739162
view(data, :, :, 12, 1)

# ╔═╡ 13a8fe64-efb4-11ea-1ac6-dbd49975dddf
minMax = percentile(data[:], [1; 99])

# ╔═╡ 8df7fcb6-f602-11ea-3f15-ff7c241ce9d4
plotSlice = function(zl, cmap="greys")
	heatmap(view(data, :,:, zl, 1),
		aspect_ratio=:equal,
		axis=:off,
		clim=tuple(minMax...),
		c=Symbol(cmap),
	title = "Slice # $(zl) with colors: $(cmap)")
end

# ╔═╡ 9dc04240-0681-11eb-3db5-ad9e7e1ee810
plotSlice(zLoc,sliceColor)

# ╔═╡ 71606d26-efb4-11ea-12e2-75f85b9ea7d8
scaleFunc = scaleminmax(minMax[1], minMax[2])

# ╔═╡ 3ed9060c-06fb-11eb-1440-a168e684e874
begin
	bgSlice = bg_image[:,:,someSlice]
	statSlice = stat_image[:,:, someSlice]
	statSlice[statSlice .< 2.3] .= NaN
	# minAndMax = percentile(bgSlice, (5, 95))
	#Gray.(scaleFunc.(displayData)
	gImage = Gray.(scaleFunc.(bgSlice))
	statScale = scaleminmax(0, 8.0)
	sImage =  statSlice; #statScale.(statSlice)
	cImage = gImage .+ sImage
	heatmap(gImage, title="slice $(someSlice), anatomy + statistical map")
	heatmap!(sImage, aspect_ratio=1,c=:hot, clims=(2.3,8))
	xaxis!(showaxis=:no)

end

# ╔═╡ 072ce484-efb2-11ea-36ee-29a0147cfc04
"""
displaySlice

	return and display a slice in a particular orientation.

"""
displaySlice = function(e, s = 1, orientation = 1, scaling = scaleFunc)
	# - can index in +1st dimension of dataset
	# so even if dims = 3, indexing into 1 of dim=4 will work.
	#if (ndims(e) == 4)
	#	e = mean(e; dims = 4)
	#	println("console message - averaging across time")
	#end
	if  (orientation == 1)
		displayData = view(e, s, :,:,1)
	elseif (orientation == 2)
		displayData = view(e, :,s,:,1)
	elseif( orientation == 3)
		displayData = view(e, :, :, s, 1)
	else
		error("uhoh!")
	end
	rotl90(Gray.(scaling.(displayData)))
end

# ╔═╡ bb15657a-f5f4-11ea-2555-65b86fec3444
[displaySlice(data, i, 2, scaleFunc) for i in 1:10:size(data, 2)]

# ╔═╡ c6f20f76-f5f4-11ea-05d2-aba2612911f4
[displaySlice(data, i, 1, scaleFunc) for i in 1:10:size(data, 1)]

# ╔═╡ f289d412-05be-11eb-172b-316cebd67e80
gr()
#plotly()

# ╔═╡ 25d9fd54-f450-11ea-3b2a-2f872e88a372
md"""z coord $(@bind zSlider Slider(1:size(data,3)))"""

# ╔═╡ 1fed409a-f450-11ea-21be-8d44768dea80
displaySlice(data, zSlider, 3, scaleFunc)

# ╔═╡ Cell order:
# ╟─47a60d48-058e-11eb-3140-03a5c22a27ee
# ╟─ab3afe3a-0684-11eb-090d-ffeea5fcadec
# ╟─730d087c-067b-11eb-0427-e9425b7c1ed6
# ╟─8029070e-efb0-11ea-1d0a-7f402f2e9989
# ╟─06433d72-058b-11eb-0565-c3e0f393ef5d
# ╟─38bf1198-049f-11eb-0b21-334cd2fe39d3
# ╟─542b0706-058b-11eb-0ae3-d57aa9ead256
# ╟─bfb16300-058e-11eb-38d2-3572d8db913f
# ╟─ec3816c4-058b-11eb-2d8b-a1ed496f153d
# ╟─0e5ddcca-058c-11eb-3068-73441b317478
# ╟─cb5f17f4-067b-11eb-3f8b-3ff5d0aa088c
# ╟─62d614fc-058f-11eb-1fc0-09e86807181e
# ╟─8baff8a4-058f-11eb-0882-9f305cf65b7c
# ╟─cc6e657e-058f-11eb-2314-5b6df2f908e9
# ╟─70e80930-04a1-11eb-055b-152d2006d067
# ╟─1b79dd34-041d-11eb-037a-93335f4b4116
# ╟─50856a72-05ba-11eb-0f1f-7dc909270ad5
# ╟─82d57e14-04a1-11eb-2190-9daace12efa8
# ╟─6f3591f2-049f-11eb-2a97-ad43e34fb000
# ╟─79e91b4e-0654-11eb-2fac-777bdb7c41f5
# ╟─0d3ab51c-0590-11eb-1d21-1f6fde64b56f
# ╟─3c37796e-058c-11eb-330f-c9f07956d74e
# ╟─2a2a7d84-058f-11eb-27bb-db743f9910d5
# ╟─69bc09c0-058c-11eb-36a6-bd6e942b5e08
# ╟─81adabb0-058c-11eb-3ee6-959bbf5f0c5d
# ╟─266fa81e-067e-11eb-19f6-37bcc6c623e3
# ╟─1c814864-0686-11eb-1c6f-2f0fa3d874cf
# ╟─c8c8b590-0685-11eb-037a-d17f5e4c3ba4
# ╟─26100fbe-0686-11eb-10eb-6752e71e2481
# ╟─3a80a800-0686-11eb-2da5-49c8cb23254d
# ╟─8b8056b0-0686-11eb-15ae-c32b9e0034a2
# ╟─c00be5ca-0686-11eb-1a62-c75bb04147df
# ╟─c4cd3778-0686-11eb-169e-57ca0d655f5b
# ╟─46874f50-0686-11eb-0e2b-8f54c14e5e8a
# ╟─7f61ac58-0684-11eb-182a-ad7903933dd5
# ╟─9dc04240-0681-11eb-3db5-ad9e7e1ee810
# ╟─7d531b12-0682-11eb-3417-4f641024d076
# ╟─49141f5e-06fa-11eb-1cb0-7ffc29a11b48
# ╟─3ebec912-06fc-11eb-23de-cdb28471c46d
# ╟─3ed9060c-06fb-11eb-1440-a168e684e874
# ╟─1386a3c6-06fb-11eb-2abe-29cfbf505e0e
# ╟─f30741d0-06fa-11eb-290b-375c807e13c6
# ╟─9c1b27d4-0684-11eb-235b-eb019285d84d
# ╟─27115482-0682-11eb-19d8-c344581ec8af
# ╟─4aa28d70-067c-11eb-1bfc-179c255fcdd0
# ╟─6381fc42-0654-11eb-3d50-9fd2b4d241b7
# ╟─9d530386-0654-11eb-13ea-4dbc0c7d9f03
# ╟─d12a8c42-0680-11eb-225d-f9d14982eef4
# ╟─b247bd76-f07e-11ea-1e07-577643ff20af
# ╟─a8ab8ee6-f083-11ea-1316-7fbe4e125914
# ╟─8fa5f8da-04a1-11eb-2941-65dd79fdba84
# ╟─bd6e4a08-0588-11eb-0a97-d1d98b30d894
# ╟─9d485144-0588-11eb-3996-0fc8abf6ed9d
# ╟─c35d782a-04a0-11eb-2d29-296d9f33ecd9
# ╟─409aa96c-04a0-11eb-377a-2732f7fecbaa
# ╟─f3661c44-049f-11eb-1ff7-6f0d6d96d8fe
# ╟─2c6ad37c-04a0-11eb-28b7-8b55c2f60cb3
# ╟─0b4e44d4-04a1-11eb-3e57-53f8cdb029a7
# ╟─9f76bc18-0421-11eb-2d94-f7aa0f2e0998
# ╠═f0ef494e-f80c-11ea-1163-b954c1739162
# ╠═bb15657a-f5f4-11ea-2555-65b86fec3444
# ╠═c6f20f76-f5f4-11ea-05d2-aba2612911f4
# ╟─44de79f0-f450-11ea-0f93-15cb49d19311
# ╟─1fed409a-f450-11ea-21be-8d44768dea80
# ╟─2b2f6be4-f07f-11ea-1459-775749ca949f
# ╠═8df7fcb6-f602-11ea-3f15-ff7c241ce9d4
# ╟─9bc349e8-efb0-11ea-1392-cd735828af33
# ╟─141dde62-efb1-11ea-34df-cd5e838d8129
# ╟─13a8fe64-efb4-11ea-1ac6-dbd49975dddf
# ╠═71606d26-efb4-11ea-12e2-75f85b9ea7d8
# ╟─072ce484-efb2-11ea-36ee-29a0147cfc04
# ╟─f289d412-05be-11eb-172b-316cebd67e80
# ╟─25d9fd54-f450-11ea-3b2a-2f872e88a372

max_stars_repo_name: schluppeck/hands-on-brain-data

max_issues_repo_name: schluppeck/hands-on-brain-data

max_forks_repo_name: schluppeck/hands-on-brain-data


text:
module ResultTypes

export Result, ErrorResult, unwrap, unwrap_error, iserror

type Result{T, E<:Exception}
    result::Nullable{T}
    error::Nullable{E}
end

Result{T}(x::T) = Result{T, ErrorException}(Nullable{T}(x), Nullable{ErrorException}())
Result{T, E}(x::T, ::Type{E}) = Result{T, E}(Nullable{T}(x), Nullable{E}())

function ErrorResult{T, E<:Exception}(::Type{T}, e::E)
    Result{T, E}(Nullable{T}(), Nullable{E}(e))
end

function ErrorResult{T}(::Type{T}, e::AbstractString)
    Result{T, ErrorException}(Nullable{T}(), Nullable{ErrorException}(ErrorException(e)))
end

function ErrorResult{T}(::Type{T})
    Result{T, ErrorException}(Nullable{T}(), Nullable{ErrorException}(ErrorException("")))
end

function unwrap{T, E}(r::Result{T, E})::T
    if !isnull(r.result)
        return get(r.result)
    elseif !isnull(r.error)
        throw(get(r.error))
    else
        error("Empty Result{$T, $E} type")
    end
end

function unwrap_error{T, E}(r::Result{T, E})::E
    if !isnull(r.error)
        return get(r.error)
    else
        error("$r is not an ErrorResult")
    end
end

Base.convert{T, S, E}(::Type{T}, r::Result{S, E})::T = unwrap(r)

function Base.convert{T, S, E}(::Type{Result{S, E}}, x::T)
    Result{S, E}(Nullable{S}(convert(S, x)), Nullable{E}())
end

function Base.convert{T, E}(::Type{Result{T, E}}, e::E)
    Result{T, E}(Nullable{T}(), Nullable{E}(e))
end

function Base.show{T, E}(io::IO, r::Result{T, E})
    if iserror(r)
        print(io, "ErrorResult(", T, ", ", unwrap_error(r), ")")
    else
        print(io, "Result(", unwrap(r), ")")
    end
end

iserror(r::Result) = !isnull(r.error)

end # module

ext: jl


text:
"""
# Divide-and-conquer algorithms

- Multiplication
- Recurrence relations
    - Master theorem
    - Binary search
- Mergesort
    - An n log n lower bound for sorting
- Medians
    - A randomized divide-and-conquer algorithm for selection
    - Efficiency analysis
    - The Unix sort command
- Matrix multiplication
- The fast Fourier transform
    - Why multiply polynomials?
    - An alternative representation of polynomials
    - Evaluation by divide-and-conquer
    - Interpolation
    - A matrix reformulation
    - Interpolation resolved
    - Inversion formula
    - A closer look at the fast Fourier transform
    - The definitive FFT algorithm
    - The fast Fourier transform unraveled
    - The slow spread of a fast algorithm
"""
module Chapter2

export
    karatsuba,
    karatsuba2,
    mergesort,
    iterativemergesort,
    polynomialmultiplication

"""
    karatsuba(x, y)

Using the Karatsuba algorithm, compute the product of two integers `x` and `y`.

# Examples

```jldoctest
julia> karatsuba(2412, 3231)
7793172
```
"""
function karatsuba(x::Int, y::Int)
    if ndigits(x) == 1 || ndigits(x) == 1
        return x * y
    end

    n = max(ndigits(x), ndigits(y))
    m = n ÷ 2

    xL, xR = x ÷ 10^m, x % 10^m
    yL, yR = y ÷ 10^m, y % 10^m

    P1 = karatsuba(xL, yL)
    P2 = karatsuba(xR, yR)
    P3 = karatsuba((xL + xR), (yL + yR))

    return (P1 * 10^(2m)) + ((P3 - P1 - P2) * 10^m) + P2
end

"""
    karatsuba2(x, y)

Using the Karatsuba algorithm, compute the product of two integers `x` and `y`.

# Examples

```jldoctest
julia> karatsuba(2412, 3231)
7793172
```
"""
function karatsuba2(x::Int, y::Int)
    if ndigits(x, base=2) == 1 || ndigits(y, base=2) == 1
        return x * y
    end

    n = max(ndigits(x, base=2), ndigits(y, base=2))
    m = n ÷ 2

    xL, xR = x ÷ 2^m, x % 2^m
    yL, yR = y ÷ 2^m, y % 2^m

    P1 = karatsuba(xL, yL)
    P2 = karatsuba(xR, yR)
    P3 = karatsuba((xL + xR), (yL + yR))

    return (P1 * 2^n) + ((P3 - P1 - P2) * 2^m) + P2
end

"""
    mergesort(a)

Sort an array of numbers `a` using the merge sort algorithm.

# Examples

```jldoctest
julia> mergesort([10, 2, 5, 3, 7, 13, 1, 6])
8-element Vector{Int64}:
  1
  2
  3
  5
  6
  7
 10
 13

julia> mergesort([1, 4.54312, 3.4321, π])
4-element Vector{Float64}:
 1.0
 3.141592653589793
 3.4321
 4.54312
```
"""
function mergesort(a::AbstractVector)
    if length(a) > 1
        m = length(a) ÷ 2
        return merge(mergesort(a[1:m]), mergesort(a[(m + 1):end]))
    else
        return a
    end
end

function merge(x::AbstractVector, y::AbstractVector)
    if length(x) == 0
        return y
    end
    if length(y) == 0
        return x
    end
    if x[1] <= y[1]
        return pushfirst!(merge(x[2:end], y), x[1])
    else
        return pushfirst!(merge(x, y[2:end]), y[1])
    end
end

function iterativemergesort(a::AbstractVector)
    Q = []
    for i ∈ 1:length(a)
        push!(Q, [a[i]])
    end
    while length(Q) > 1
        push!(Q, merge(eject(Q), eject(Q)))
    end
    return eject(Q)
end

function polynomialmultiplication(A::AbstractVector{T}, B::AbstractVector{T}) where {T <: Real}
    product = zeros(length(A) + length(B))

    for i ∈ 1:length(A)
        for j ∈ 1:length(B)
            product[i + j] += A[i] * B[j]
        end
    end

    return product
end

# TODO: FFT
# function FFT(A, ω)
#
# end

end # module


text:
using GraphMatFun, LsqFit, LinearAlgebra,Random
include("exp_reset_all.jl");
#include("simulationtools.jl");
include("simulationtools.jl");

f=exp
m=6
target_n=200 # Measure error wrt this discr
n=50;        # Optimize wrt this discr
rho=2.7;  # Increased SID rho m=5

its=5;
droptol0=1e-10;
opt_kwargs=Dict(:logger=> 0,:γ0 => 0.5,:droptol => droptol0,
                :linlsqr => :real_svd,:maxit => its)
base = init_state_mult(f,rho,n,m; eltype=Complex{BigFloat})
base.params[:n]=n;
base.params[:target_n]=target_n;
base.params[:opt_kwargs]=opt_kwargs;
base.params[:kickit_mode]=4;

oldrho="2_22";
datadir=joinpath("..","..","data","exp");
# First the optimization free
filename=joinpath(datadir,"exp_sid_m$(m).cgr");
sid_org=init_state_file!(deepcopy(base),:sid,filename,showmeta=true)

#filename=joinpath("exp_bbc_m$(m)_rho$oldrho.cgr");
#bbc_org=init_state_file!(deepcopy(base),:bbc,filename,showmeta=true)

filename=joinpath(datadir,"exp_sastre_m$(m).cgr");
sastre_org=init_state_file!(deepcopy(base),:sastre,filename,showmeta=true)

#sastre_org.cref=get_degopt_crefs(sastre_org.graph);

filename=joinpath(datadir,"exp_ps_m$(m)_rho$oldrho.cgr");
ps_org=init_state_file!(deepcopy(base),:ps,filename,showmeta=true)

filename=joinpath(datadir,"exp_mono_m$(m)_rho$oldrho.cgr");
mono_org=init_state_file!(deepcopy(base),:mono,filename,showmeta=true)

filename=joinpath(datadir,"exp_mono_m6_opt_rho2_7_warmstart.cgr");
#filename="/home/jarl/jobb/src/matfun/simulations/newgraphs/exp_m$(m)_mono+GN_$oldrho.cgr";
mono0=init_state_file!(deepcopy(base),:mono,filename,showmeta=true)
(mono,simlist,commandlist)=
     run_sequence(mono0,"dddddddddddddddddddddsssGGGssdddssssdddddsssssssdddsdddsssszq");
#"dddddddddddddddddddddsssGGGSSSddddSkSSdSKKKKKKsSdddgggggggggSSSSSSSSsssssssssssSSDsSSSSSSSSSSDSSSSGGGGGGGGGSSSdSdkDSSztkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkSSq"
filename=joinpath(datadir,"exp_ps_m$(m)_opt_rho$oldrho.cgr");
ps0=init_state_file!(deepcopy(base),:ps,filename,showmeta=true)
ps0.params[:kickit_mode]=2;
(ps,simlist,commandlist)=
     run_sequence(ps0,
     "sdddddddddSdddSdSddddSSNSkSSSGGGSddSddSSddSSzddSddkSkkkkkkkkkkkSdSgggggggSSSSSSSSssGGGGGGGSSSzSddssssddsddszk");

filename=joinpath(datadir,"exp_sastre_m$(m)_opt_rho$oldrho.cgr");
sastre0=init_state_file!(deepcopy(base),:sastre,filename,showmeta=true)
(sastre,simlist,commandlist)=
     run_sequence(sastre0,
                  "ddddddddddsdddddddddddssGGGssddddsdddddsssdddssszq");

filename=joinpath(datadir,"exp_sid_m$(m)_opt_rho$oldrho.cgr");
sid0=init_state_file!(deepcopy(base),:sid,filename,showmeta=true)
(sid,simlist,commandlist)=
     run_sequence(sid0,
     "ddddddddddsdddddddddddssGGGssddddsdddddsssdddsssdddsssddddsssssdsddszq")

#filename=joinpath("exp_bbc+GN_m$(m)_rho$oldrho.cgr");
#bbc0=init_state_file!(deepcopy(base),:bbc,filename,showmeta=true)
#(bbc,simlist,commandlist)=
#     run_sequence(bbc0,
#     "sdddddsdddddssssssddddsssssssddddsssssssssssssdddssssssssssddddssssssssdsggggggsdssssssssssssssssGGGGGsssGGGsssq")
#

include("exp_print_all.jl");
println("Run include(exp_save_all.jl) if you want to save");

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using Images, FileIO
using StatsBase
using AverageShiftedHistograms
using DSP
using FFTW

#using CSV
#using DataFrames

#Pre processes an image
#takes a 2D array
#returns a 2D array

"""
    pre_proc(im)

convolve image with a disk-shaped structure element
"""
function pre_proc(im)
    img = Gray.(im) .> 0.05
    bw = Float64.(img)
    strel = [0 1 1 1 0;
             1 1 1 1 1;
             1 1 1 1 1;
             1 1 1 1 1;
             0 1 1 1 0]
    R = findall(x->x==true, strel)
    LocalFilters::erode!(bw, R)
    LocalFilters::dilate!(bw, R)
    #localfilter!(bw, :, min, 3)
    #localfilter!(bw, :, max, 3)
    im2 = nonzero(bw)
    return im2
end

"""
    nonzero(image)

evaluate where an image is
nonzero
"""
function nonzero(image)
    dims = size(image)
     im = Array{Bool}(undef, dims[1], dims[2])
    for i = 1:dims[1]
        for j = 1:dims[2]
            im[i,j] = !(image[i,j]>0)
        end
    end
    return im
end

#Sloppy euclidian distance computation
#takes a 1D array of CartesianIndex tuples
#returns a 1D array of Float32
function pdist(indices)

    dims = size(indices)[1]

    lens = Array{Float32}(undef, Nlen(size(indices)[1]+1))
    idx = 1

    for i = 1:dims[1]
            for j = i:dims[1]
                lens[idx] = cart(indices[i], indices[j])
                idx += 1
            end
    end

    return lens
end

#Computes the euclidian distance between a pair of CartesianIndex tuples
#takes a tuple of two CartesianIndex tuples
#returns a float
function cart(a, b)
    return sqrt((a[1]-b[1])^2 + (a[2]-b[2])^2)
end

#Computes distances between a list of indices, and histograms these
#values across a 1:100 interval with bins of width 1

#takes a list of CartesianIndices
#returns an AverageShiftedHistograms ash object (see methods on github)
#TODO make the range adjustable, add optimisation to shortcut for cluster functions
#(computing this many histograms is wasteful vs one big one)
function im_stat(indx)
    dims = size(indx)[1]
    lens = zeros(Float32, dims)
    primer = Array{Float32}([1])
    s1 = ash(primer, rng=0:1:100.)
    idx = 1
    @inbounds @fastmath  Threads.@threads for i = 1:dims[1]
        #a = indx[i]
            @inbounds @fastmath for j = i:((dims)[1])
                lens[j] = cart(indx[i], indx[j])
               #print("\n i = ", i, " j = ", j)
            end
        #idx += 1

        merge!(s1, ash(lens[idx:end], rng=0:1:100.))
    end

    return s1
end

function im_stat1(indx)
    dims = size(indx)[1]
    lens = zeros(Float32, dims)
    primer = Array{Float32}([1])
    s1 = ash(primer, rng=1:1:100.)
    idx = 1
    @inbounds @fastmath for i = 1:dims[1]
        #a = indx[i]
            @simd for j = idx:((dims)[1])
                lens[j] .= cart(indx[i], indx[j])
               #print("\n i = ", i, " j = ", j)
            end
        idx += 1

        merge!(s1, ash(lens[idx:end], rng=1:1:100.))
    end

    return s1
end

#Computes distances between a list of indices, and histograms these
#values across a 1:100 interval with bins of width 1
#does this after computing the 8-connectivity set of clusters
#does this for all clusters. Also computes the SN of the image if not provided

#takes an image
#returns an AverageShiftedHistograms ash object (see methods on github)
#TODO tuple this with SN as output
function tpcof(im, SN=0...)

    if SN == 0
        full_arr = trues((1:size(im)[1], 1:size(im)[2]))

        ind = findall(x->x==true, full_arr) #CartesianIndices((1:size(im)[1], 1:size(im)[2]))
        #ind = CartesianIndices((size(im)[1], size(im)[2]))
        Nind = length(ind)
        SN = (2 .* im_stat(ind).counts) ./ length(im)
        print("\n")

    end

    indices = findall(x->x==true, im)
    BN = (2 .* im_stat(indices).counts) ./ length(indices)
    S2 = (length(indices)/length(im)) .* (BN ./ SN)
    #print(length(indices), "    ", length(im) , "\n")

    #S2 =  (length(indices)/length(im)).*((BN.counts)./ (2 .* SN.counts ./))
    S2[1] = length(indices)/length(im)

    return S2, length(indices), SN, 0:1:100.
end

#Computes distances between a list of indices, and histograms these
#values across a 1:100 interval with bins of width 1
#Also computes the SN of the image if not provided

#takes an image
#returns an AverageShiftedHistograms ash object (see methods on github)
#TODO tuple this with SN as output
function tpclf(im, SN=0...)

    label = label_components(im)
    comps = component_subscripts(label)[2:end]
    bbox = component_boxes(label)
    primer = Array{Float32}([1])
    s1 = ash(primer, rng=1:1:100.)

    if SN == 0
        ind = CartesianIndices((1:size(im)[1], 1:size(im)[2]))
        Nind = length(ind)
        SN = im_stat(ind)
    end

    for cluster in comps
        merge!(s1, im_stat(cluster))
    end
    #print(size(xy(s1))
    C2 = (Nind./length(im).*(xy(s1)[2]./ xy(SN)[2]))

    return C2
end

"""
Compute a binary disk-shaped structure element
given radius r on the basis of (areal) interpolation

"""
function disk_strel(r)
    dim = (2*r)-1
    arr = zeros(Bool, dim, dim)
    coords = CartesianIndices(arr)
    #dump(coords)
    middle = (length(coords) ÷ 2) +1
    mid = coords[middle]
    for idx in coords
        if cart(idx, mid) < 5
            arr[idx] = true
        end

    end
    return arr
    #centre = CartesianIndices(arr[r,r])

end

function polydisp_iso_circs(dims, r, N)
    x = dims[1]
    y = dims[2]
    off = (r-1)
    im = zeros(Bool, x+(off*2), y+(off*2))
    strel = disk_strel(r)
    #carts = findall(x->x==true, view(im, off:x, off:y))
    for i= 1:N
        x1 = rand(off+1:x)
        y1 = rand(off+1:y)
        im[(x1-off):(x1+off), (y1-off):(y1+off)] .= im[(x1-off):(x1+off), (y1-off):(y1+off)] .| strel
        #print(x1,"\t", y1, "\n")
    end
    #dump(carts)
    #strel = disk_strel(r)
    #centres = sample(carts, 80)
    #for coord in centres
        #x = coord[1]
        #y = coord[2]
        #print(coord)
        #im[(x-off):(x+off), (y-off):(y+off)] = strel

return im
end

"""
Generate isotropically distributed circles that
may overlap. Keep the implied phase fraction
ϕ below 0.3 if possible

"""
function monodisp_circ(dims, R, N)
    x = dims[1]
    y = dims[2]
    im = zeros(Bool, dims)
    xrand = rand(1:x, N)
    yrand = rand(1:y, N)
    coords = hcat(xrand, yrand)
    for coord in zip(xrand, yrand)
        #print(coord, "\n")
        for i in 1:x, j in 1:y

            if cart((i,j), coord) < R
                im[i,j] = true
            end
        end
    end
    return im

end

function S2_analyt(r, R, phi)
    #v2dv1=2*(r>2*R)+2/pi*(pi+r/2/R.*sqrt(1-r.^2/4/R^2)-acos(r/2/R)).*(2*R>=r);
    squarebit = real(sqrt.(Complex.((1 .- r .^2 ./ 4 ./ (R^2)))))
    f = (r ./ 2 ./ R)
    cosbit = acos.(f .* (f .<= 1 )) .* (f .<= 1 )
    leftbit = 2 .* (r .> (2*R))
    rightbit = 2 .* R .>=r
    midbit = r ./ 2 ./ R
    total = leftbit .+ 2 ./ pi .* (pi .+ (midbit .* squarebit) .-  cosbit) .* rightbit

    eta=-log(phi);

    S2=exp.(-total .* eta);

    return S2
end

function S2_fft(I)
    n,m = size(I)
    I_fft = fft(I)
    abs.(fftshift(ifft(I_fft .* conj(I_fft)))) ./ (n*m)
end

function S2_fft_test1(I,P)
    n,m = size(I)
    I_fft = fft(I, P)
    abs.(fftshift(ifft(I_fft .* conj(I_fft)))) ./ (n*m)
end

function S2_fft_test(I)

    I_fft = rfft(I)

    n,m = size(I_fft)
    abs.(fftshift(ifft(I_fft .* conj(I_fft)))) ./ (n*m)
end

function C2_pdist(im, SN=0...)

    label = label_components(im)
    comps = component_subscripts(label)[2:end]
    bbox = component_boxes(label)
    primer = Array{Float32}([1])
    s1 = ash(primer, rng=0:1:100.)

    if SN == 0
        ind = CartesianIndices((size(im)))
        Nind = length(ind)
        SN = im_stat(ind)
    end

    for cluster in comps
        merge!(s1, blas_stat(cluster))
    end
    #print(size(xy(s1))
    C2 = (Nind./length(im).*(xy(s1)[2]./ xy(SN)[2]))
    for i= 1:length(C2)
    if isnan(C2[i])
        C2[i] = 0
    end
end

    return C2
end

function C2_fft_MT(I)

    l = label_components(I)
    bbox = component_boxes(l)[2:end]
    C2= zeros(size(I)[1], size(I)[2], Threads.nthreads())
    #A = zeros(Threads.nthreads())
    alloc = falses(size(I)[1], size(I)[2], Threads.nthreads())
    Threads.@threads for cluster in bbox
        #print(S2_fft(cluster))
        xbound = cluster[1][1]:cluster[2][1]
        ybound = cluster[1][2]:cluster[2][2]
        alloc[xbound, ybound, Threads.threadid()] = I[xbound, ybound]
        #print(size(alloc[:,:,Threads.threadid()]))
        #print("\n")
        #print(size(C2[:,:, Threads.threadid()]))
        C2[:,:, Threads.threadid()] = C2[:,:, Threads.threadid()] .+ S2_fft(alloc[:,:,Threads.threadid()])
        #if mod(i , 100) == 0
        #    imshow(C2)
        #end
        alloc[xbound, ybound, Threads.threadid()] .= false
        #C2 = C2 .+ S2_fft(cluster)
    end

    return C2
end

function C2_fft_MT_test(I)

    P = plan_fft(I)
    l = label_components(I)
    bbox = component_boxes(l)[2:end]
    C2= zeros(size(I)[1], size(I)[2], Threads.nthreads())
    alloc = falses(size(I)[1], size(I)[2], Threads.nthreads())

    Threads.@threads for cluster in bbox
        xbound = cluster[1][1]:cluster[2][1] #find bounds
        ybound = cluster[1][2]:cluster[2][2]
        alloc[xbound, ybound, Threads.threadid()] = I[xbound, ybound] #bind values to preallocated array
        C2[:,:, Threads.threadid()] = C2[:,:, Threads.threadid()] .+ S2_fft(alloc[:,:,Threads.threadid()])
        alloc[xbound, ybound, Threads.threadid()] .= false

    end

    return C2
end

function blas_stat(indx)
    F = Tuple.(indx)
    x, y= first.(F), last.(F)
    x1, y1 = copy(x), copy(y)
    len = size(x)[1]
    #print(len)
    xout = zeros(len)
    yout = zeros(len)
    dist = zeros(len)
    primer = Array{Float32}([1])
    s1 = ash(primer, rng=0:1:100.)
    @inbounds @fastmath Threads.@threads for i = 1:len
        dist .= sqrt.((x .- x1[1]) .^2 .+ (y .- y1[1]) .^2)
        ash!(s1, view(dist, i:len))
    end
    return s1
end

function C2_fft(I)
    l = label_components(I)
    bbox = component_boxes(l)[2:end]
    C2= zeros(size(I)[1], size(I)[2])
    #A = zeros(Threads.nthreads())
    alloc = falses(size(I)[1], size(I)[2])
    for cluster in bbox
        #print(S2_fft(cluster))
        xbound = cluster[1][1]:cluster[2][1]
        ybound = cluster[1][2]:cluster[2][2]
        alloc[xbound, ybound] = I[xbound, ybound]
        #print(size(alloc[:,:,Threads.threadid()]))
        #print("\n")
        #print(size(C2[:,:, Threads.threadid()]))
        C2 = C2 .+ S2_fft(alloc[:,:])
        #if mod(i , 100) == 0
        #    imshow(C2)
        #end
        alloc[xbound, ybound] .= false
        #C2 = C2 .+ S2_fft(cluster)
    end
    return C2
end

"""
    make_rand_im(philen, dims)

create a random binary array of size dims with
philen white pixels set within it. Also return
a coordinate list for white and black pixels
"""
function make_rand_im(philen, dims)
    idx = CartesianIndices((dims))
    #println(length(idx))
    wpix = sample(idx, philen, replace=false)
    #println(length(wpix))
    #println(wpix)
    bpix = setdiff(idx, wpix)
    im = zeros(Bool, (dims))

    for pix in wpix
        im[pix] = true
    end
    return im, (wpix, bpix)
end

"""
    get_rand_pix((pix))::Tuple{Array{CartesianIndex{2},1},Array{CartesianIndex{2},1}}

Select a random white and black pixel from a tupled
pair of lists of cartesian coordinates
"""
function get_rand_pix((pix))::Tuple{Array{CartesianIndex{2},1},Array{CartesianIndex{2},1}}
    wpick = [sample(pix[1])]
    bpick = [sample(pix[2])]
    #print(wpick, "\t", bpick, "\n")

    return wpick, bpick
end

function ediff_S2_fft(S2_guess,  S2)
    return sum((S2 .- S2_guess) .^2)
end

function ediff_C2(C2_guess, C2)
    return sum((((C2 .- C2_guess) .* 2) .^2))
end

function ediff_S2(S2_guess, S2)
    return sum((((S2 .- S2_guess) .* 2) .^2))
end

function ediff_test(C2_guess, C2, S2_guess,  S2)
    return (ediff_C2(C2_guess, C2), ediff_S2_fft(S2_guess, S2))
end

"""
    cluster_stat(clusters, step::Float64, maxrng)

Fast compution of S2 count for an array of arrays of cartesian indices.
Multithreaded in the sense that
"""
function cluster_stat(clusters, step::Float64, maxrng)
    s1_cache_threads = [zeros(Int64,maxrng) for i in 1:Threads.nthreads()];
    @inbounds @fastmath Threads.@threads for indx in clusters
        s1_cache_threads[Threads.threadid()] .+= blas_stat_st2(indx, step, maxrng)
    end
    f = zeros(Int64, length(s1_cache_threads[1]))
    for i in 1:Threads.nthreads(); f = f .+ s1_cache_threads[i]; end
    return f
end

function findlens(list, lower, upper)::Tuple{Int64,Int64}
    mymap = map( x->(lower<=length(x)<=upper), list)
    left = findfirst(isequal(1), mymap)
    right = findlast(isequal(1), mymap)
    #println("left = ", left, "\tright = ", right)
    if left === nothing
        return (0,0)
    else
    return (left,right)
    end
end

"""
    C2_pdist4(im, maxrng::Int64, SN::Array{Float64,1}=Array{Float64}(undef, 0))

Fast total computation of the C2 count for an array of arrays of cartesian indices
This version splits the problem by the type of cluster, selecting an appropriate
algorithm to suit the scale of each component.
"""
function C2_pdist4(im, maxrng::Int64, SN::Array{Float64,1}=Array{Float64}(undef, 0))

    comps = component_subscripts(label_components(im, trues(3,3)))[2:end]
    sort!(comps, by = x -> (length(x)))
    #println(length.(comps))
    singles = findlens(comps, 1, 1)
    #dump(singles)
    smalls = findlens(comps, 2, 20)
    bigs = findlens(comps, 21, length(im))
    #println("singles = ", singles, " smalls = ", smalls, " bigs = ", bigs)
    #bbox = component_boxes(label)
    #primer = Array{Float32}([])
    #maxrng = 100#min(100, size(im)[1], size(im)[2])
    BN = fit(Histogram, Array{Float32}([]), 0:1.:maxrng, closed=:left).weights

    if isempty(SN)#SN == 0.0
        SN = SN_comp(size(im),  maxrng)
        #ind = vec(CartesianIndices((size(im))))

        #SN = 2 .* blas_stat4(ind, 1.0, maxrng) ./ length(im)
    end
    #println("there are ", length(comps), " clusters")
    #println("SN = ", SN)
    #plot(SN)
    if singles !==(0,0) #Accumulate BN for all single-pixel clusters (dist=0 implicitly)
        BN[1] += 2*singles[2]
        #println("trig'ed\t", length(singles))
    end

    if smalls !==(0,0)
        BN .= (BN .+ (2 .* cluster_stat((comps[smalls[1]:smalls[2]]), 1.0, maxrng)))
    end

    if bigs!==(0,0)
        for cluster in (@view comps[bigs[1]:bigs[2]])
            BN .= (BN .+ (2 .* blas_stat4(cluster, 1.0, maxrng)))
            #break
        end
    end
    BN =  (BN ./ sum(im))
    C2 = ((sum(im)/length(im)) .* (BN ./ SN))
    for i= 1:length(C2)
        if isnan(C2[i])
            C2[i] = 0
        end
    end
    return C2, BN
end

"""
    SN_comp(dims::Tuple{Int64,Int64},  maxrng::Int64)

Wrapper for ```blas_stat5``` acting as a function barrier
"""
function SN_comp(dims::Tuple{Int64,Int64},  maxrng::Int64)::Array{Int64,1}
    2 .* blas_stat5(vec(CartesianIndices(dims)), maxrng)
end

"""
    inner_blas2(x, y, dist, len)

compute the lower triangular self-interaction matrix of a list of N items.
gives a flat vector of F32, of length (N(N-1)/2). Interaction here is computing
the L2 norm, but could be otherwise.
"""
function inner_blas2(x::Array{Float32,1}, y::Array{Float32,1}, dist::Array{Float32,1}, len::UInt64)
    @fastmath @inbounds for i = 1:len
        #println("left = ", (triang(i-1)+1), "\tright = ", triang(i))
        fragment!(view(dist, (triang(i-1)+1):triang(i)), view(x, 1:i), view(y, 1:i), x[i], y[i])
    end
end

"""
    blas_stat_st2(indx, step, maxrng)

For a set of N points, compute the unique distance between all possible pairs of
points. Bins distance measurements with a histogram of bin width step and length maxrng
Single threaded implementation, greedy with memory. Good for smaller computations
"""
function blas_stat_st2(indx, step::Float64, maxrng)
    F = Tuple.(indx)
    x, y= Float32.(first.(F)), Float32.(last.(F))
    dist = Vector{Float32}(undef,triang(length(x)))
    len::UInt64 = length(x)
    inner_blas2(x, y, dist, len)
    #print(length(fit(Histogram, dist, 0:step:(maxrng)).weights))
    return fit(Histogram, dist, 0:step:(maxrng)).weights
end

# function blas_stat_st2(x::Array{Float64,1}, y::Array{Float64,1}, step::Float64, maxrng)
#     F = Tuple.(indx)
#     x, y= Float32.(first.(F)), Float32.(last.(F))
#     dist = Vector{Float32}(undef,triang(length(x)))
#     len::UInt64 = length(x)
#     inner_blas2(x, y, dist, len)
#     #print(length(fit(Histogram, dist, 0:step:(maxrng)).weights))
#     return fit(Histogram, dist, 0:step:(maxrng)).weights
# end

function triang(n::UInt64)::UInt64
    return (n*(n+1))÷2
end

function triang(n::Int64)::Int64
    return (n*(n+1))÷2
end

"""
Compute the L2 norm between points (xi, yi) and the arrays (arrx, arry) and write
to preallocated array dist
"""
function fragment!(dist, arrx, arry, xi, yi)
    dist .= sqrt.((arrx .- xi) .^2 .+ (arry .- yi) .^2)
end

"""
    blas_stat4(indx, step::Float64, maxrng)::Array{Int64,1}

Computes a very fast multithreaded probabilistic S2 count. Employs kernel density
estimation in place of a true histogram. Gives hilarious results when used
for reconstruction.
"""
function blas_stat4(indx
    , step::Float64, maxrng)::Array{Int64,1}
    F = Tuple.(indx)
    x, y= Float32.(first.(F)), Float32.(last.(F))
    len = size(x)[1]
    primer = Array{Float32}([0])
    _s1_cache_threads = [ash(primer, rng=0:step:(maxrng-1), m =1) for i in 1:Threads.nthreads()];
    _dist_cache_threads = [Vector{Float32}(undef,len) for i in 1:Threads.nthreads()]

    @inbounds @fastmath Threads.@threads for i = 1:len
        fragment!(view(_dist_cache_threads, Threads.threadid())[1], x, y, x[i], y[i] )
        println(view(_dist_cache_threads[Threads.threadid()], i:len))
        ash!(view(_s1_cache_threads, Threads.threadid())[1], view(_dist_cache_threads[Threads.threadid()], i:len))
    end
    f = zeros(Int64, length(_s1_cache_threads[1].counts))
    for i in 1:Threads.nthreads(); f .= f .+ _s1_cache_threads[i].counts; end
    f[1] = f[1] - Threads.nthreads()
    return f
end

"""
Compute the S2 count (multithreaded, O(N) memory)

For a list of N CartesianIndex tuples, compute and bin the L2 norm between all pixel
pairs in the list.
"""
function blas_stat5(indx, maxrng)
    F = Tuple.(indx)
    x, y= Float32.(first.(F)), Float32.(last.(F))
    len = size(x)[1]
    primer = Array{Float32}([])
    cnt_thrd = [zeros(Int64,maxrng) for i in 1:Threads.nthreads()];
    dist_thrd = [zeros(Float32, len) for i in 1:Threads.nthreads()]

    @inbounds @fastmath Threads.@threads for i = 1:len
        fragment!(view(dist_thrd, Threads.threadid())[1], x, y, x[i], y[i] )
        view(cnt_thrd, Threads.threadid())[1] .+=  fit(Histogram, view(dist_thrd[Threads.threadid()], i:len), 0:1:maxrng).weights
    end
    f = zeros(Int64, length(cnt_thrd[1]))
    for i in 1:Threads.nthreads(); f = f .+ cnt_thrd[i]; end
    #f[1] = f[1] - Threads.nthreads()
    return f
end

"""
Get list of all surface pixels from a 2d binary image using convolution to discriminate
"""
function surf_opt2(im::Array{Bool,2})::Tuple{Array{CartesianIndex{2},1},Array{CartesianIndex{2},1}}
    edge_kern::Array{Int64} = [1 1 1; 1 10 1; 1 1 1]
    arrim::Array{Int64,2} = conv(im, edge_kern)
    return (findall(in(10:17), arrim), findall(in(1:8), arrim))
end

"""
Get list of all surface pixels from a 2d binary image using convolution to discriminate
Sanitised version to prevent the function returning nothing
"""
function surf_opt(im::Array{Bool,2})::Tuple{Array{CartesianIndex{2},1},Array{CartesianIndex{2},1}}
    edge_kern::Array{Int64} = [1 1 1; 1 10 1; 1 1 1]
    arrim::Array{Int64,2} = conv(im, edge_kern)[2:end-1, 2:end-1]
    #display(heatmap(arrim, color=:grays, aspect_ratio=1))
    #imshow(im)
    #imshow(arrim)
    bb::Array{CartesianIndex{2},1} = findall(in(1:8), arrim)
    wb::Array{CartesianIndex{2},1} = findall(in(10:17), arrim)

    if !isnothing(bb) & ! isnothing(wb)
        return (wb, bb)
    else
        return  (CartesianIndices(collect(1:1)), CartesianIndices(collect(1:1)))
    end
end

"""
Update tuple of 2 cartesian indices by swapping elements from each list at linear
index wpick (for 1) and bpick (for 0)
"""
function swap_pix(guess, pix, wpick, bpick)
    #println("swap_pix : \t\twpix length = ", length(pix[1])," bpix length = ", length(pix[2]) )
    guess[pix[1][wpick]] = false ##THIS IS THE DANGER ZONE
    guess[pix[2][bpick]] = true ########THIS IS NOT SYMMETRIC

    pix[1][wpick], pix[2][bpick] = pix[2][bpick], pix[1][wpick]

    return (pix)
end

"""
    find_equivalent(shortlist, pix, pick_s)

Computes equivalent position of unique entry in a shortlist to an entry in a master
list. Uses a much more sensible tupled pick - should use that everywhere
"""
function find_equivalent(shortlist, pix, pick_s)
    wcart = shortlist[1][pick_s[1]]
    bcart = shortlist[2][pick_s[2]]
    #println("find_equivalent : \twpick = ", pick_s[1], "\tbpick = ", pick_s[2])
    #println("find_equivalent : \twpix len = ", length(pix[1]), "\tbpix len = ", length(pix[2]))
    wrand = findfirst(x -> x == wcart, pix[1])
    brand = findfirst(x -> x == bcart, pix[2])
    #println("find_equivalent : \twrand = ", wrand, "\tbrand = ", brand)
    return (wrand, brand)
end

"""
Find the intersection between two tupled lists of cartesian indices
We use this to create a segment of a list that maintains ordering
"""
function surf_rand(surf_res, pix)
    return (intersect(surf_res[1], pix[1]), (intersect(surf_res[2], pix[2])))
end

"""
Selects random elements in a pair of lists
"""
function pick_rand_coord(list::Tuple{Array{CartesianIndex{2},1},Array{CartesianIndex{2},1}})::Tuple{Int64, Int64}
    wpick = rand(1:length(list[1]))
    #println("pick_rand_coord : \twlen = ", length(list[1]), "\tblen = ", length(list[2]))
    bpick = rand(1:length(list[2]))
    #println("pick_rand_coord : \twpick = ", wpick, "\tbpick = ", bpick)
    return wpick, bpick
end

"""
Fetch the region around a pixel selection; give a reduced window if near an edge
"""
function fetch_locale(im, pix, kern, wpick, bpick)
    re, be = size(im)[1], size(im)[2]
    wp = pix[1][wpick]
    bp = pix[2][bpick]

    #println(kern)
    #println("wp = ",wp)
    kern = check_edge(wp, re, be)
    #println(kern)
    #println(wp .+ kern)
    a = view(im, wp .+ kern)
    kern = check_edge(bp, re, be)
    #println("bp = ", bp)
    #println(kern)
    #println(bp .+ kern)
    b = view(im, bp .+ kern)
    return (a, b)
end

"""
Prevent fetch_locale from colliding with image edges by providing a reduced
window if an edge is collided with
"""
function check_edge(p, re, be)
    kern = CartesianIndices((-1:1, -1:1))
    l1 = -1
    r1 = 1
    l2 = -1
    r2 = 1
    if p[1] ==1
        l1 = 0
        #kern = kern .= CartesianIndices((0:1, -1:1))
    end
    if p[1] ==re
        r1 = 0
       # kern = CartesianIndices((-1:0, -1:1))
    end
    if p[2] ==1
        l2 = 0
        #kern = CartesianIndices((-1:1, 0:1))
    end
    if p[2] ==be
        r2 = 0
        #kern = CartesianIndices((-1:1, -1:0))
    end
    kern = CartesianIndices((l1:r1, l2:r2))

    return kern
end

"""
Find all unique values present in a small window into the region adjacent to a
pixel selection. Used to interrogate a cluster labelled image
"""
function get_region_names(regions)
    w1 = unique(regions[1][1])
    w2 = unique(regions[2][1])
    b1 = unique(regions[1][2])
    b2 = unique(regions[2][2])
    return ((w1, b1),(w2, b2))
end

"""
Deduct the contribution of a given set of pixels from a previously calculated
S2 count
"""
function subtract_cluster(hist, cluster)
    hist .= hist .- blas_stat_st2(cluster, 1.0, 100)
end

"""
Examine which clusters are present in multiple windows simultaneously to determine
which clusters are due to change from a pixel swap
"""
function find_unique_regions(names)
    olds = filter(x->x!=0,union(names[1][1], names[1][2]))
    news = filter(x->x!=0,union(names[2][1], names[2][2]))

    return (olds, news)
end

"""
Compute the C2 count in an image. Serves as a starting point for a reconstruction
"""
function C2_initialise(im, maxrng::Int64)
    comps = component_subscripts(label_components(im, trues(3,3)))[2:end]
    BN = zeros(Int64, maxrng)
    #for comp in comps
    #println(length(comp))
    #end
    BN .= (BN .+ 2 .* (cluster_stat((comps), 1.0, maxrng)))
    return BN
end

"""
Sloppy computation to find the contribution to C2 count from the old state and the
new, and hence the difference between them.
"""
function compute_change!(subtract, add, clusts_old, clusts_new, maxrng)
    for clust in clusts_old; subtract .= subtract .+ blas_stat5(clust, maxrng);  end
    for clust in clusts_new; add      .= add .+ blas_stat5(clust, maxrng); end
    return  add, subtract
end

"""
Eliminates NaN and Inf in an array
"""
function naninf(arr)
    for i=1:length(arr)
        if isnan(arr[i]) | isinf(arr[i])
            arr[i] = 0
        end
    end
    return arr
end

"""
Computes the C2 count change between a pair of images with one pixel swapped between them
"""
function update_C2_BN(BN, pix, original::Array{Int64, 2}, modified::Array{Int64, 2}, wpick, bpick, maxrng)
    kern = CartesianIndices((-1:1, -1:1))
    names = get_region_names((fetch_locale(original, pix, kern, wpick, bpick), fetch_locale(modified, pix, kern, wpick, bpick)))
    old, new = find_unique_regions(names)
    #println("old unique regions = ", old)
    #println("new unique regions = ", new)
    subtract = zeros(Int64, length(BN))
    add = zeros(Int64, length(BN))
    #println("BN is of length ",  length(BN))
    #clusts_old = component_subscripts(original) #we can replace this with findall for a particular index or indices of interest
    #clusts_new = component_subscripts(modified)
    clusts_old, clusts_new = get_clusters(old, new, original, modified)
    #println(size(clusts_new), "\t", size(clusts_old) ,"\t", clusts_new, "\t",clusts_old)
    add, subtract = compute_change!(subtract, add,  clusts_old, clusts_new, maxrng)

    return  (add .- subtract) .* 2
end

"""
Fuses clusters to reduce the scale of the problem
"""
function get_clusters(old, new, original, modified)
    clusts_old = [CartesianIndex{2}[]]
    clusts_new=  [CartesianIndex{2}[]]
    for clust in old; push!(clusts_old, findall(x->x==clust, original)); end
    for clust in new; push!(clusts_new, findall(x->x==clust, modified)); end
    popfirst!(clusts_old)
    popfirst!(clusts_new)
    #println(clusts_old, "\n", clusts_new)
    return clusts_old, clusts_new
end

"""
Compute the update to S2 count. Find the contribution from
all white pixel's distance from wpick and bpick, and the difference between them
"""
function update_S2_BN(BN, pix, wpick, bpick, philen, maxrng)
    wp = pix[1][wpick]
    bp = pix[2][bpick]
    F = Tuple.(pix[1])
    x, y= Float32.(first.(F)), Float32.(last.(F))
    dist1 = zeros(Float32, philen) #this can depend on philen
    dist2 = zeros(Float32, philen)

    #pix[2][bpick] = pix[1][wpick]# we need to actually modify x and y, not pix
    x[wpick],y[wpick] = bp[1], bp[2]
    add = fit(Histogram, fragment!(dist1, x, y, wp[1], wp[2]), 0:1.0:maxrng).weights
    x[wpick],y[wpick] = wp[1], wp[2]
    #pix[1][wpick] = pix[2][bpick]
    subtract = fit(Histogram, fragment!(dist2, x, y, bp[1], bp[2]), 0:1.0:maxrng).weights
    #println("ash is of length ", length(subtract))
    return (add .- subtract) .* 2
end

"""
Find the S2 count for an image
"""
function S2_initialise(pix, maxrng)
     (2 .* (blas_stat5(pix[1], maxrng)))
end

"""
Normalise S2 count into S2 proper by dividing by BN and ϕ
"""
function S2_finalise(SN, BN, im)
    dims = size(im)
    BN =  (BN ./ sum(im))
    S2 = ((sum(im)/length(im)) .* (BN ./ (SN ./ length(im))))
end

"""
Loads and thresholds an image
"""
function loadim(fname::String)::Array{Bool, 2}
    img::Array{Bool, 2} = (Gray.(load(fname))) .< 0.50
end


text:
"""
     EmbeddingPositions

If the embedding was constructed from multivariate data, which of the
original data variables was mapped to which variable in the embedding?

## Fields
- **`positions`**: position[i] returns the index of the variable that goes
    into the i-th coordinate axis of the embedding.
"""
struct EmbeddingPositions
    positions::AbstractVector{Int}
end

export EmbeddingsPositions


text:
abstract type Potential end

apply_potential(p::Potential, r) = p.f.(r)

struct PseudoHS <: Potential
    f::Function
end

PseudoHS() = PseudoHS(pseudohs)
PseudoHS(λ::Real) = PseudoHS((x) -> pseudohs(x; λ=λ))

function pseudohs(r; λ=50.0)
    # FIXME: This formula only applies for a 3D system
    dT = 1.4543 + 1.199 / (λ^1.0545)
    A = λ * (λ / (λ - 1.0))^(λ - 1.0)
    B = λ / (λ - 1.0)
    eps_energy = 1.0 / dT

    ur = 0.0

    if r < B
        xr = 1.0 / (r^λ)
        xa = 1.0 / (r^(λ - 1.0))
        ur = A * eps_energy * (xr - xa) + eps_energy
    else
        ur = 0.0
    end

    return ur
end

struct HardSphere <: Potential
    f::Function

    HardSphere() = new(hard_sphere)
end

function hard_sphere(r)
    ur = 0.0

    if r < 1.0
        ur = Inf
    else
        ur = 0.0
    end

    return ur
end

struct SquareWell <: Potential
    f::Function
end

SquareWell(λ::Real) = SquareWell((x) -> square_well(x; λ=λ))

function square_well(r; λ=1.5)
    uij = 0.0

    if r < 1.0
        uij = Inf
    elseif 1.0 <= r < λ
        uij = -1.0
    else
        uij = 0.0
    end

    return uij
end

struct LennardJones <: Potential
    f::Function

    LennardJones() = new(lj)
end

function lj(r)
    r6 = 1.0 / r^6
    r12 = r6 * r6

    return 4.0 * (r12 - r6)
end


text:
const OnewayDict = Dict{Symbol,Union{Real,Missing,Array{Union{Real,Missing},1}}}

""" Interpretes the oneway tags of a OSM way object. Returns an array of strings.
"""
function oneway(way::Dict)::Dict
    if way["type"] != "way"
        throw(ArgumentError("`way` must be an OSM way object.") )
    end

    tags = get(way, "tags", "_")
    if tags != "_"
        oneway = get(tags, "oneway", "_")
         onewayconditional = get(tags, "oneway:conditional", "_")
        if onewayconditional != "_"
            # interpret conditional oneway tag
            out = parseonewayconditional(onewayconditional)
            # ckeck if oneway and oneway:conditional match
            if !(oneway == "reversible" || (oneway == "no" || oneway == "_") )
                @warn "The `oneway` and `oneway:conditional` tag are conflicting." oneway = oneway condition = out
            end

        # common streets
        elseif oneway == "no" || oneway == "alternating"
            out = Dict{Symbol, Dict}(:uv => OnewayDict(:weight => missing),
                                     :vu => OnewayDict(:weight => missing) )

        # oneway streets
        elseif oneway == "yes"
            out = Dict{Symbol, Dict}(:uv => OnewayDict(:weight => missing))
        elseif oneway == "-1"
            out = Dict{Symbol, Dict}(:vu => OnewayDict(:weight => missing))
        elseif oneway == "reversible"
            out = Dict{Symbol, Dict}()
            @warn "A `way` with a `reversible` tag but without a `oneway:conditional` tag was discovered, it is ignored." ID=way["id"]

        # if not marked as oneway but one of the following tags is given, it is oneway
        elseif oneway == "_" && any(get(tags, "junction", "_") .== ["roundabout", "circular"])
            out = Dict{Symbol, Dict}(:uv => OnewayDict(:weight => missing))
        elseif oneway == "_" && get(tags, "highway", "_") == "motorway"
            out = Dict{Symbol, Dict}(:uv => OnewayDict(:weight => missing))
        elseif oneway == "_"
            out = Dict{Symbol, Dict}(:uv => OnewayDict(:weight => missing),
                                     :vu => OnewayDict(:weight => missing))
        else
            out = Dict{Symbol, Dict}()
            @warn "An unhandled oneway tag was discovered therefore the according way is ignored." ID=way["id"]
        end
    else
        out = Dict{Symbol, Dict}()
        @warn "A `way` object without `tags` is discovered it is ignored." ID=way["id"]
    end

    return out
end

""" Interpretes the oneway:conditional part of the oneway information. Returns a dictinary.
"""
function parseonewayconditional(onewayconditional::AbstractString)::Dict
    condition = parseconditionaltag(onewayconditional)
    condikeys = collect(keys(condition))
    for ck = condikeys
        if length(collect(keys(condition[ck][:rules]))) > 1
            @error "Handling `oneway:conditional` with more the one rule per condition is not implemented." condition = onewayconditional
        end
        value = condition[ck][:value]
        rule = condition[ck][:rules][:rule1]

        uvweight = repeat(Union{Real,Missing}[missing], 168)
        vuweight = repeat(Union{Real,Missing}[missing], 168)

        if value != "yes"
            vuweight[rule[:hoursofweek] .+ 1] .= Inf
        elseif value != "-1"
            uvweight[rule[:hoursofweek] .+ 1] .= Inf
        else
            @error "`one way:conditional` tags with following value, are not supported." value = value
        end

    end
    if all(uvweight .=== missing)
        uvweight = missing
    elseif all(uvweight .=== Inf)
        uweight = Inf
    end
    if all(vuweight .=== missing)
        vuweight = missing
    elseif all(vuweight .=== Inf)
        vuweight = Inf
    end
    out = Dict{Symbol, Dict}(:uv => OnewayDict(:weight => uvweight),
                             :vu => OnewayDict(:weight => vuweight) )

    return out
end


text:
module SceneTrees

export SceneNode, Path

struct Path <: AbstractVector{String}
    entries::Vector{String}
end

Base.size(p::Path) = size(p.entries)
Base.IndexStyle(::Type{Path}) = IndexLinear()
Base.convert(::Type{Path}, x::AbstractVector{<:AbstractString}) = Path(x)

Base.joinpath(p::Path, s...) = foldl(joinpath, p, s)
Base.joinpath(p::Path, s::Symbol) = joinpath(p, String(s))
Base.joinpath(p::Path, s::AbstractString) = _joinpath(p, split(s, '/'))

_joinpath(p::Path, s::AbstractVector{<:AbstractString}) = foldl(_joinpath, p, s)
_joinpath(p::Path, s::AbstractString) = isempty(s) ? Path(String[]) : Path(vcat(p.entries, s))

Base.show(io::IO, p::Path) = print(io, string('/', join(p.entries, '/')))
Base.getindex(p::Path, i::Int) = p.entries[i]
Base.getindex(p::Path, v::AbstractVector) = Path(p.entries[v])
Base.setindex(p::Path, x, i::Int) = p.entries[i] = x

mutable struct SceneNode
    object::Nullable{Vector{UInt8}}
    transform::Nullable{Vector{UInt8}}
    children::Dict{String, SceneNode}
end

SceneNode() = SceneNode(nothing, nothing, Dict{String, SceneNode}())

Base.getindex(s::SceneNode, name::AbstractString) = get!(SceneNode, s.children, name)
function Base.delete!(s::SceneNode, name::AbstractString)
    if haskey(s.children, name)
        delete!(s.children, name)
    end
end

function Base.getindex(s::SceneNode, path::Path)
    if length(path.entries) == 0
        return s
    else
        return s[path[1]][path[2:end]]
    end
end

function Base.delete!(s::SceneNode, path::Path)
    parent = s[path[1:end-1]]
    child = path[end]
    delete!(parent, child)
end

function Base.foreach(f::Function, s::SceneNode)
    f(s)
    for (k, v) in s.children
        foreach(f, v)
    end
end

end


text:
# src/devCtrlAux.jl
# Created by Alejandro Villoria and Anabel Ovide, June 3, 2021
# MIT license
# Auxiliary functions for sub-module QCCDDevControl

using .QCCDevControl_Types

"""
Helper function for `linear_transport()`.
Removes the ion from the origin chain, adds it to the destination chain,
and sets `destination` to `nothing` if it has arrived to its final destination.
# Arguments
* `ion` - Ion to be moved.
* `origin` - Current zone the ion is in.
* `destination` - Zone the ion is going to.
"""
function _move_ion(ion ::Qubit,
    origin ::Union{GateZone, AuxZone, LoadingZone},
    destination ::Union{GateZone, AuxZone, LoadingZone})

  # Remove ion from origin
  if origin.zoneType === :loadingZone
    origin.hole = nothing
  else
    index = origin.end0 === destination.id ? 1 : length(origin.chain)
    deleteat!(origin.chain, index)
  end

  # Add ion to destination and change its position
  if destination.zoneType === :loadingZone
    destination.hole = ion.id
  else
    destination.end0 === origin.id ? pushfirst!(destination.chain, [ion.id]) :
                                    push!(destination.chain, [ion.id])
  end
  ion.position = destination.id

  # Remove destination to ion if it has arrived to its destination
  if ion.destination === destination.id
    ion.destination = nothing
  end
end

"""
Helper function for `swap()`.
Swaps the ions in the chain.
# Arguments
* `qdc` - Control device struct.
* `ion1_idx` - First ion to be swaped.
* `ion2_idx` - Second ion to be swaped.
"""
function _swap_ions(qdc::QCCDevControl, ion1_idx:: Int, ion2_idx:: Int)
  zone = giveZone(qdc, qdc.qubits[ion1_idx].position)
  pos1 = nothing
  pos = 0
  for (num,i) in enumerate(zone.chain)
    pos1 = findall(x->x==ion1_idx, i)
    if !isempty(pos1)
      pos = num
      pos1 = pos1[1]
      break
    end
  end
  pos2 = collect(Iterators.flatten(map( y -> findall(x->x==ion2_idx, y), zone.chain)))[1]
  zone.chain[pos][pos1] = ion2_idx
  zone.chain[pos][pos2] = ion1_idx
end

"""
Helper function to compute the time and actualize time in the control device.
# Arguments
* `qdc` - Control device struct.
* `t::Time_t` — time at which the operation commences.
* `op_t::Time_t` — time the opration last.
Returns the actual time after the operation.
"""
function compute_time(qdc::QCCDevControl, t::Time_t, op_t::Time_t)
  t₀ = t + op_t
  t₀ > t  || throw(Error("Error while computing time"))
  qdc.t_now = t₀
  return t₀
end

max_stars_repo_name: TartuQC/trapped-ion-qccd-simulator

max_issues_repo_name: TartuQC/trapped-ion-qccd-simulator

max_forks_repo_name: TartuQC/trapped-ion-qccd-simulator


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

using Markdown
using InteractiveUtils

# ╔═╡ 7950392c-f279-11ea-16cd-117125c04757
using Plots

# ╔═╡ 6be0aea8-f273-11ea-20a8-b9e192d2df7c
md"""# PW - N
Langton
* werk lokaal
* ken uw computer (waar staat mijn data/waar staat mijn notebook/waar ben ik aan het werk?)
"""

# ╔═╡ 76aac224-f273-11ea-38d8-81a3efd97984
# waar ben ik an het werken
pwd()

# ╔═╡ f1476d90-f274-11ea-20a1-1dd4670be6d9
# wat bestaat er in mijn pwd?
readdir(pwd())

# ╔═╡ 930cbd5a-f273-11ea-1ef5-1324cc695386
# relatief path
readdir("./data/")

# ╔═╡ 13990aea-f275-11ea-363d-8b94278d165e
# absoluut path
readdir("/Users/bart/Documents/Stack/ES313/ES313.jl/Exercises/data/")

# ╔═╡ b6f1587a-f273-11ea-0cd3-5b035ae9d0e4
r = readlines("./data/Langtonsrules.txt")

# ╔═╡ 1886dc3e-f276-11ea-234d-6d4368d5e6eb
split(r[1],"")

# ╔═╡ c1ff5738-f275-11ea-3dd2-b9e753043688
begin
	res = parse.(Int,split(r[1],""))
	current = res[1]
	neighbors = res[2:5]
	next = res[6]
	println("current $(current) with neighbords $(neighbors) becomes $(next)")
end

# ╔═╡ 8097b168-f276-11ea-389b-358431906bc8
function getrules(path::String)
	r = readlines(path)
	# initiate rule dict
	d = Dict()
	for rule in r
		# parse rule
		res = parse.(Int,split(rule,""))
		current = res[1]
		neighbors = res[2:5]
		next = res[6]
		# add rule to dict
		for i in 0:3
			println("current $(current) with neighbords $(neighbors) becomes $(next)")
			get!(get!(d,current,Dict(circshift(neighbors,i)=>next)), circshift(neighbors,i), next)
		end
	end
	d
end

# ╔═╡ be1dcdb8-f278-11ea-0794-f7497b0283e0
d =getrules("./data/Langtonsrules.txt")

# ╔═╡ 13d8d032-f278-11ea-21a1-7327e1e5f8ed
d[1][[1,0,2,2]]

# ╔═╡ ec80dad8-f278-11ea-0f69-7fc4b77a4853
readlines("./data/Langtonstart.txt")

# ╔═╡ 064d8dee-f279-11ea-393a-154c73a1bc9b
# inlezen startsitutie
begin
	s = readlines("./data/Langtonstart.txt")
	start = permutedims( hcat( [parse.(Int,split(line,"")) for line in s]...) )
end

# ╔═╡ 7d196e40-f279-11ea-09be-1b2902d25369
begin
	mycmap = cgrad([RGBA(0/255,0/255,0/255),
    RGBA(0/255,0/255,255/255),
    RGBA(255/255,0/255,0/255),
    RGBA(0/255,255/255,0/255),
    RGBA(255/255,255/255,0/255),
    RGBA(255/255,0/255,255/255),
    RGBA(255/255,255/255,255/255),
    RGBA(0/255,255/255,255/255)]);
	heatmap(start, yflip=true, color=mycmap)
end

# ╔═╡ 6de438f0-f7eb-11ea-0932-cfa0f78d7cd4
md"""resterende taken van vandaag:
* implementeren van startsituatie in groter veld
* updates realiseren
* animeren
"""

# ╔═╡ f039e0e0-f7f3-11ea-3385-bf1224bda003
size(start)

# ╔═╡ 3f646880-f7f3-11ea-1b48-9386ca64e247
# creëren vd wereld - start in midden
function genworld(start, dims)
	# input check
	ssize = size(start)
	if any(ssize .>= dims)
		error("world dimensions not OK")
	end
	# initialisatie
	world = zeros(Int,dims...)
	# vullen
	i = round(Int,(dims[1] - ssize[1])/2) + 1
	j = round(Int,(dims[2] - ssize[2])/2) + 1
	world[i:i+ssize[1]-1, j:j+ssize[2]-1] = start

	return world
end

# ╔═╡ 42662f70-f7f6-11ea-27d5-7d9c6a3d289a
size(start)

# ╔═╡ af920f40-f7f3-11ea-0487-1352825d5983
world = genworld(start,(14,20))

# ╔═╡ 66be144e-f7f6-11ea-31a0-f9f5407a8e79

# ╔═╡ 4a098770-f7f8-11ea-2335-0d281a74ae3d
"""
	updateworld!(world, rules)

update the current world (modifying it)
"""
function updateworld!(world, rules)
	# nieuwe staat initialiseren
	newworld = zeros(Int, size(world))
	# voor iedere locatie nieuwe toestand bepalen
	for i in 2:size(world,1)-1
		for j in 2:size(world,2)-1
			newworld[i,j] = rules[world[i,j]][[world[i-1,j],world[i,j+1],world[i+1,j], world[i,j-1]]]
		end
	end
	@show world .== newworld
	world = copy(newworld)
end

# ╔═╡ 9a2138c0-f7f8-11ea-3a7d-87737a676985
begin
	updateworld!(world, d)
	world
end

# ╔═╡ e2ec52b0-f7f8-11ea-0f27-237ed18a5dc6
typeof(updateworld!(world,d))

# ╔═╡ bb636410-f7f6-11ea-177d-8d31331fe198
function updateworld(world, rules)
	# nieuwe staat initialiseren
	newworld = zeros(Int, size(world))
	# voor iedere locatie nieuwe toestand bepalen
	for i in 2:size(world,1)-1
		for j in 2:size(world,2)-1
			newworld[i,j] = rules[world[i,j]][[world[i-1,j],world[i,j+1],world[i+1,j], world[i,j-1]]]
		end
	end
	return newworld
end

# ╔═╡ 7d265b70-f7f7-11ea-3893-a3262c84ccbe
heatmap(updateworld(world, d), yflip=true, color=mycmap, size=(200,250))

# ╔═╡ b80b6a00-f7f7-11ea-3be5-9b94c64f465a
heatmap(world, yflip=true, color=mycmap, size=(200,250))

# ╔═╡ a5799ef0-f7f9-11ea-01fa-5b10e3cd3d24
function changeme!(A)
	B = similar(A)
	for i in eachindex(A)
		B[i] = A[i] + sum(i)
	end
	A = deepcopy(B)
end

# ╔═╡ a551cb9e-f7f9-11ea-345f-0f23a31af889
A = ones(Int, 4,4)

# ╔═╡ a543e8f0-f7f9-11ea-23ad-c195f4ba1172
for i in 1:10
	changeme!(A)
	println(A)
end

# ╔═╡ 0b0d09f0-f7fa-11ea-2091-5b2b8608911e
changeme!(A)

# ╔═╡ Cell order:
# ╠═6be0aea8-f273-11ea-20a8-b9e192d2df7c
# ╠═76aac224-f273-11ea-38d8-81a3efd97984
# ╠═f1476d90-f274-11ea-20a1-1dd4670be6d9
# ╠═930cbd5a-f273-11ea-1ef5-1324cc695386
# ╠═13990aea-f275-11ea-363d-8b94278d165e
# ╠═b6f1587a-f273-11ea-0cd3-5b035ae9d0e4
# ╠═1886dc3e-f276-11ea-234d-6d4368d5e6eb
# ╠═c1ff5738-f275-11ea-3dd2-b9e753043688
# ╠═8097b168-f276-11ea-389b-358431906bc8
# ╠═be1dcdb8-f278-11ea-0794-f7497b0283e0
# ╠═13d8d032-f278-11ea-21a1-7327e1e5f8ed
# ╠═ec80dad8-f278-11ea-0f69-7fc4b77a4853
# ╠═064d8dee-f279-11ea-393a-154c73a1bc9b
# ╠═7950392c-f279-11ea-16cd-117125c04757
# ╠═7d196e40-f279-11ea-09be-1b2902d25369
# ╟─6de438f0-f7eb-11ea-0932-cfa0f78d7cd4
# ╠═f039e0e0-f7f3-11ea-3385-bf1224bda003
# ╠═3f646880-f7f3-11ea-1b48-9386ca64e247
# ╠═42662f70-f7f6-11ea-27d5-7d9c6a3d289a
# ╠═af920f40-f7f3-11ea-0487-1352825d5983
# ╠═66be144e-f7f6-11ea-31a0-f9f5407a8e79
# ╠═9a2138c0-f7f8-11ea-3a7d-87737a676985
# ╠═e2ec52b0-f7f8-11ea-0f27-237ed18a5dc6
# ╠═4a098770-f7f8-11ea-2335-0d281a74ae3d
# ╠═bb636410-f7f6-11ea-177d-8d31331fe198
# ╠═7d265b70-f7f7-11ea-3893-a3262c84ccbe
# ╠═b80b6a00-f7f7-11ea-3be5-9b94c64f465a
# ╠═a5799ef0-f7f9-11ea-01fa-5b10e3cd3d24
# ╠═a551cb9e-f7f9-11ea-345f-0f23a31af889
# ╠═a543e8f0-f7f9-11ea-23ad-c195f4ba1172
# ╠═0b0d09f0-f7fa-11ea-2091-5b2b8608911e

max_stars_repo_path: Exercises/PWN (endstate).jl

max_issues_repo_path: Exercises/PWN (endstate).jl

max_forks_repo_path: Exercises/PWN (endstate).jl


text:
#--------------------------------------------------------------
# MultilevelTT package
# example code for https://arxiv.org/abs/1802.09062
#--------------------------------------------------------------
# Code for 1D Poisson equation example,
# using AMEn solver via Python interface (see TTPy.jl)
module PoissonAMEn

using TT, TTFEM, TTAlg
using PyCall, TTPy

export fullSystem, amen1D, amen2D

# explicit assembly of preconditioned system
# (as required by amen solver)
function fullSystem(L::Integer, D::Integer)
  C = ttBPXdnC(L, D)

  f1 = ttrhsdn(L)
  f = f1
  for k = 2:D
      f = ttkron(f, f1)
  end
  g = matvec(C, f)
  svdtrunc!(g, svd!(g), 1e-14)

  Q = ttBPXdnQLaplace(L, D, 1)

  B = matmat(transpose(Q), Q)
  for k = 2:D
      Q = ttBPXdnQLaplace(L, D, k)
      B = add!(B, 1., matmat(transpose(Q), Q))
  end

  p = size(B[L],1)
  q = size(B[L],4)
  T = reshape(B[L],(p*2^D*2^D,q))
  S = reshape(B[L+1],(q,1))
  B[L] = reshape(T*S,(p,2^D,2^D,1))
  pop!(B)

  return B, g
end

function amen1D(L::Integer, ε::AbstractFloat)
    B, g = fullSystem(L, 1)
    C = ttBPXdnC(L, 1)
    u = amen(B, g, g, ε)
    Cu = matvec(C, u)
    return u, Cu
end

function amen2D(L::Integer, ε::AbstractFloat,
            kickrank::Integer = 4, nswp::Integer = 20)
    B, g = fullSystem(L, 2)
    u = amen(B, g, g, ε, kickrank, nswp)
    B = 0
    gc()

    C = ttBPXdnC(L, 2)
    Cu = matvec(C, u)

    return u, Cu
end

end


text:
using Pkg
Pkg.add(PackageSpec(url="https://github.com/sbos/AdaGram.jl.git"))

max_stars_repo_name: benandow/UiRef

max_issues_repo_name: benandow/UiRef

max_forks_repo_name: benandow/UiRef


text:
# # Heat Equation
#
# This example demonstrates how to combine `OrdinaryDiffEq` with `DiffEqOperators` to solve a time dependent PDE.
# We consider the heat equation on the unit interval, with dirichlet boundary conditions:
# ∂ₜu = Δu
# u(x=0,t)  = a
# u(x=1,t)  = b
# u(x, t=0) = u₀(x)
#
# For `a = b = 0` and `u₀(x) = sin(2πx)` a solution is given by:
u_analytic(x, t) = sin(2*π*x) * exp(-t*(2*π)^2)
#
# We want to reproduce it numerically
#
using DiffEqOperators, OrdinaryDiffEq

nknots = 100
h = 1.0/(nknots+1)
knots = range(h, step=h, length=nknots)
ord_deriv = 2
ord_approx = 2

const Δ = CenteredDifference(ord_deriv, ord_approx, h, nknots)
const bc = Dirichlet0BC(Float64)

t0 = 0.0
t1 = 0.03
u0 = u_analytic.(knots, t0)

step(u,p,t) = Δ*bc*u
prob = ODEProblem(step, u0, (t0, t1))
alg = KenCarp4()
sol = solve(prob, alg)

using Test
@test u_analytic.(knots, t1) ≈ sol[end] rtol=1e-3


text:
module projection

"These are scaling terms to convert array subscripts to obj space"
XREF = cld(length(x1objspace),2)
YREF = cld(length(x1objspace),2)
Xidx = repeat(1:length(x1objspace), outer=(length(x2objspace)*length(x3objspace)))
Yidx = repeat(1:length(x2objspace), inner=(length(x1objspace)*length(x3objspace)))
Zidx = repeat(1:length(x3objspace), inner=(length(x1objspace)*length(x2objspace)))

"For shifting the images relative to object space"
SHIFTX = Float32.((Xidx.-XREF).*OSR)
SHIFTY = Float32.((Yidx.-YREF).*OSR)

"Resampling points -- just to ensure the center point remains consistent"
half1 = centerPT:-OSR:1
half2 = (centerPT + OSR):OSR:x1_imgspacelength
samples = vcat(half1,half2)
sort!(samples)

"Make lowpass 2D sinc filter for preprocess step prior to resampling"
sinc2d(x,y) = (sinc(x)*sinc(y))
x1 = range(-2, stop=2, length=13) #NOTE THIS ASSUMES OSR OF 3!!!
x2 = x1'
sincfilter = sinc2d.(x1,x2).*0.330

"Preallocation for image stacks"
imgsperlayer = Nnum*Nnum
H = makeHmatrix(psfWAVE_STACK[:,:,1],subpixelpitch,d,k0,lambda)
Hlayer = complex(zeros(size(MLARRAY,1),size(MLARRAY,2),imgsperlayer))
Himgs = zeros(length(samples),length(samples), Nnum*Nnum*length(x3objspace))

# First draft. You might be able to make it more concise...
function shiftimg()
    if x > 0
        if y > 0
            return @view
        elseif y < 0
            return @view
        else
            return @view
        end
    elseif x < 0
        if y > 0
            return @view
        elseif y < 0
            return @view
        else
            return @view
        end
    else
        if y > 0
            return @view
        elseif y < 0
            return @view
        else
            return @view
        end
    end
end

function ishiftimg()
    if x > 0
        if y > 0
            return @view
        elseif y < 0
            return @view
        else
            return @view
        end
    elseif x < 0
        if y > 0
            return @view
        elseif y < 0
            return @view
        else
            return @view
        end
    else
        if y > 0
            return @view
        elseif y < 0
            return @view
        else
            return @view
        end
    end
end

println("Computing LF PSFs (2/3)")
#NOTE: You can probably do this a lot simpler with careful/clever use of @view
#       I started to hash out the idea below. Verify. Probably need if-else for
# each shift case (positive, negative and zero). Could be a separate function that returns a view or range.
for layer in 1:objspace.zlen
    stackstart = 1 + (layer-1)*imgsperlayer
    stackend   = layer*imgsperlayer
    baseimg = psfWAVE_STACK[:,:,layer]
    for img in 1:imgsperlayer

        @views stackbin .= shift(baseimg, SHIFTX, SHIFTY) .* ishift(MLARRAY, SHIFTX, SHIFTY)

        #Fresnel 2d one step size of b0
        fft2!(stackbin)
        stackbin = stackbin .* H
        ifft2!(stackbin)

        @views finalbin .= abs2.(ishift(stackbin, SHIFTX, SHIFTY))

        "Non original code: addition of low pass filter before downsampling"
        stackbin = conv(stackbin, sincfilter)[croppedrng, croppedrng]

        Hlayer[:,:,img] = stackbin
    end
    Himgs[:,:,stackstart:stackend] .= Hlayer[samples,samples,:]
    println("Layer: " * string(layer) * " of " * string(objspace.zlen) * " complete.")
end

end # for loop timer

"Normalize the whole data set"
Himgs = Himgs./maximum(Himgs)

"""Threshold all the images (grouped by Z-plane) by the tolerance
value. Flooring all values below the cutoff to zero."""
tol = 0.005

for layer in 1:length(x3objspace)
   H4Dslice = Himgs[:,:,Zidx .== layer]
   H4Dslice[H4Dslice .< (tol*maximum(H4Dslice))] .= 0
   Himgs[:,:,Zidx .== layer] = H4Dslice
end

end # end projection module

#= unlikely needed if using phase space deconv
function calcht(Himgs, Nnum, x2objspace,Xidx, Yidx, Zidx)
    Himgsize = size(Himgs,0)
    x2length = length(x2objspace)
    lenslets = cld(Himgsize,Nnum) # How many lenslets wide is the image
    if mod(lenslets,1) == 0 # Pad the image area
        imgsize = (lenslets+1)*Nnum
    else
        imgsize = (lenslets+2)*Nnum
    end
    Nnumpts = Nnum^1
    imcenter = cld(imgsize,1)
    imcenterinit = imcenter - cld(Nnum,1)
    x0 = Xidx[Zidx.==0]
    y0 = Yidx[Zidx.==0]
    xpoints = imcenterinit .+ x0
    ypoints = imcenterinit .+ y0
    shiftx = cld(Nnum,1) .- xpoints
    shifty = cld(Nnum,1) .- ypoints
    cropstart = Int(imcenter - (Himgsize-2)/1)
    cropstop = Int(imcenter + (Himgsize-2)/1)

    projection = zeros(imgsize, imgsize, Nnumpts)
    for i in 0:Nnumpts
        projection[xpoints[i],ypoints[i],i] = 0.-1
    end
    Ht = zeros(size(Himgs))
    Hkernel = reverse(Himgs, dims=0)
    Hkernel .= reverse(Hkernel, dims=1)

    ### similar to reconstruction
    for i in 0:Nnumpts
        projection_GPU = (projection[:,:,i])
        Backprojection = (zeros(size(projection, 0), size(projection, 0), x2length ))

        for j in 0:Nnum
            for k in 0:Nnum
                for p in 0:x2length

                indexes = (Xidx .== j) .& (Yidx .== k) .& (Zidx .== p)
                Hkernel_GPU = (Himgs[:,:,indexes])

                Hkernel_GPU = convolve1(projection_GPU,Hkernel_GPU, AF_CONV_DEFAULT, AF_CONV_AUTO)

                 Backprojection[j:Nnum:end, k:Nnum:end, p] = Hkernel_GPU[j:Nnum:end, k:Nnum:end]

                end
            end
        end
        println("innerloop done")

        Backprojection = Backprojection[cropstart:cropstop, cropstart:cropstop, :]
        Backprojection = translate(Backprojection, shiftx[i], shifty[i], size(Backprojection,0), size(Backprojection,1))

        #back to host here
        indexes1 = (Xidx .== x0[i]) .& (Yidx .== y0[i])
        Ht[:,:,indexes1] = Array(Backprojection)

        println("full loop")
    end
    Backprojection = finalize(Backprojection)
    convolution = finalize(convolution)
    Hkernel_GPU = finalize(Hkernel_GPU)
    afgc(999)

    ### end similar to reconstruction

    return Ht
end
=#