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

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

# ╔═╡ ccbc7c67-e16a-48f1-a927-c002257fcc54
begin
	import Pkg
	Pkg.add(["Images","TestImages","FileIO",
			"TiledIteration",
			"PlutoUI", "Plots",
			"MosaicViews", "PerceptualColourMaps", "VideoIO"])
	using Images,TestImages,LinearAlgebra,TiledIteration,PlutoUI,Plots
	using MosaicViews, PerceptualColourMaps, FFTW, FileIO, Random
end

# ╔═╡ cbd2ea36-7ef6-4ea8-82ac-34b3e96ed8d5
md"""
### Origin and brief references
> Paul Constantine and I looked at PCA of large collections of images in the context of computing Tall and Skinny QR factorizations where we first saw discrete fourier transform-like structure as we show below for a single image. Paul had originally noticed smoothness in terms of parameterized solutions to linear systems and parameterized computational science problems, which is closely related.
> In a different vein, Lek-Heng Lim often sharply critiqued work that blindly deploys data-as-model frameworks without sufficient thought. In this case, the matrix-of-tiles is at least a vector-space of samples, so it has better properties _as a linear operator_ or _linear space_. (Which is why we are able to find structure inside it.)
> Finally, Vaneet Agrawal mentioned _reshaping_ video data (originally a "3d" matrix) into a 5d array to improve low-rank modeling. This was in the context of work on low-tubal rank matrix completition, which uses a tensor computing framework based on ideas by Misha Kilmer and subsequent work by Kilmer, Bramen, and Hao. (See also our own take on these ideas in a paper with Gleich, Greif, and Varah.) Here we are really doing the same thing.
> And of course, there are influences from all the neural network patch or tile convolution models in all of the recent deep learning on images setups.
"""

# ╔═╡ 59c37536-fff2-480a-80bc-f97125e117af
md"## Initializing packages
This adds all the packages you need to run these examples.
"

# ╔═╡ 2c8a2738-77b1-4803-91d3-9b5172e4a94a
md"""# We need to start with a picture.

$(@bind im_choice Select(["webcam","upload","url","coffee","cat","earth","mandrill","house","blobs"] |> x -> x .=> x,default="coffee"))

The box above allows you to choose one!
- `webcam`: see the webcam image area below. You can take your own picture!
- `upload`: see the upload area below. You can send us a picture from your computer. (I have no idea what happens to it, but we don't intentially do anything nefarious.)
- `url`: you have the picture on the internet somewhere we can download it from. Great, tell us the URL and we'll grab it for you.
- or use one of the named images from `TestImages.jl`.

My favorite is `coffee`. Nothing like a fantastic coffee!
"""

# ╔═╡ 3112233a-543d-48f5-823f-5b1e1cafcb61
begin
	function process_raw_camera_data(raw_camera_data)
	# the raw image data is a long byte array, we need to transform it into something
	# more "Julian" - something with more _structure_.

	# The encoding of the raw byte stream is:
	# every 4 bytes is a single pixel
	# every pixel has 4 values: Red, Green, Blue, Alpha
	# (we ignore alpha for this notebook)

	# So to get the red values for each pixel, we take every 4th value, starting at
	# the 1st:
	reds_flat = UInt8.(raw_camera_data["data"][1:4:end])
	greens_flat = UInt8.(raw_camera_data["data"][2:4:end])
	blues_flat = UInt8.(raw_camera_data["data"][3:4:end])

	# but these are still 1-dimensional arrays, nicknamed 'flat' arrays
	# We will 'reshape' this into 2D arrays:

	width = raw_camera_data["width"]
	height = raw_camera_data["height"]

	# shuffle and flip to get it in the right shape
	reds = reshape(reds_flat, (width, height))' / 255.0
	greens = reshape(greens_flat, (width, height))' / 255.0
	blues = reshape(blues_flat, (width, height))' / 255.0

	# we have our 2D array for each color
	# Let's create a single 2D array, where each value contains the R, G and B value of
	# that pixel

	RGB.(reds, greens, blues)
end

	function camera_input(;max_size=200, default_url="https://i.imgur.com/SUmi94P.png")
"""

			<button id="stop" title="Stop video">✖</button>
			<button id="pop-out" title="Pop out/pop in">⏏</button>

		<video playsinline autoplay></video>

		<button id="shutter" title="Click to take a picture">📷</button>

		Enable webcam

""" |> HTML
end

	md"""The code to handle the camera input is hidden. It's still a bit raw. It saves the images in `im_webcam`"""
end

# ╔═╡ 2ac8cb5d-ffff-4561-bc08-2ffc368338e2
md""" ### Get your own webcam image.

From <https://github.com/mitmath/18S191/blob/Fall20/lecture_notebooks/week1/01-images-intro.jl>

$(@bind raw_camera_data camera_input(;max_size=1000))
"""

# ╔═╡ 03604f42-66ce-4a39-b1b1-f5fd49bfccdb
im_webcam = process_raw_camera_data(raw_camera_data)

# ╔═╡ e6e10ee8-2bcb-4743-86be-ffac6b591b96
md""" ### Upload a picture.
This uses the `PlutoUI.jl` upload code.

$(@bind image_upload_data FilePicker([MIME("image/jpg"),MIME("image/jpeg")]))
"""

# ╔═╡ 2d4125ac-5c10-4150-bc00-e2cfcbfaf5ad
im_upload = load(IOBuffer(image_upload_data["data"]))

# ╔═╡ 52d8a212-0b34-4306-b410-e3d3dce7de3f
md""" ### Give a URL

$(@bind image_url TextField((120,1);default="https://user-images.githubusercontent.com/6933510/110924885-d7f1b200-8322-11eb-9df7-7abf29c8db7d.png"))

This defaults to the the MIT Computational Thinking tree image.
"""

# ╔═╡ fd4b3fef-9df2-4253-81b4-98a9c59a0176
im_url = load(download(image_url))

# ╔═╡ 13b28401-d752-4dae-9ddf-bfd59f3ab995
md""" ## The image we'll use for our SVD analysis.

We crop this slightly to make it a multiple of 64 to make our tile-based analysis more straightforward.

"""

# ╔═╡ 08c9f59e-4408-409a-b324-36dedc005c65
begin
	imagename = im_choice
	if im_choice == "webcam"
		im = im_webcam
	elseif im_choice == "upload"
		im = im_upload
	elseif im_choice == "url"
		im = im_url
	elseif im_choice == "coffee"
		im = testimage("coffee")
	elseif im_choice == "cat"
		im = testimage("chelsea")
	elseif im_choice == "earth"
		im = testimage("earth_apollo17") # 3000x3000-ish
		im = imresize(im, ratio=0.25)
	elseif im_choice == "house"
		im = testimage("house")
	elseif im_choice == "mandrill"
		im = testimage("mandrill")
	elseif im_choice == "blobs"
		im = testimage("blobs")
	else
		@error("Please make up your mind!")
	end

	function crop_to_multiple(im, factor; center::Bool = true)
		newsize = div.(size(im), factor).*factor
		offsets = div.(mod.(size(im), factor), center ? 2 : 1)
		return im[range.(offsets.+1, newsize.+offsets, step=1)...]
	end
	im = crop_to_multiple(im, 64)
end

# ╔═╡ 4632bc52-9106-4d05-8f6f-a189e577388d
md"""
But we aren't going to handle color approximations here. (There is another workbook on this that gets deeper into those weeds.) So we are really approximating the following image."""

# ╔═╡ 6a13f71e-983b-4d48-8787-8d58abbebbf5
Gray.(im)

# ╔═╡ 17dae8a4-010a-4816-b5a4-c9448d3cf8b1
md"""
In the standard explaination, we pick a rank for the SVD-based approximation, and then display that approximation.

**Rank for SVD approximation** (Usually you see something around 15-20.)
- `rsvd` $(@bind rsvd Slider(1:50,default=15))
"""

# ╔═╡ fddb68e0-adc2-4ff5-908d-e0204bf074e3
md""" `rsvd` = $(rsvd)"""

# ╔═╡ c372842c-a529-4ad1-8edb-47d755018dbc
begin
	#= This is the standard Image-as-Matrix SVD example. =#
  	A = Float64.(Gray.(im))
  	U, S, V = svd(A);
end; # output hidden because it isn't that interesting

# ╔═╡ dbace9f6-6b02-4e91-88a4-92fe27ec063d
begin
  	function standard_svd_approx(r::Int)
    	B = U[:,1:r]*Diagonal(S[1:r])*V[:,1:r]'
    	nparam = length(U[:,1:r])+ length(V[:,1:r])
		return B, nparam
  	end
  	svd_approx_gray, nparam_gray = standard_svd_approx(rsvd)
  	Gray.(svd_approx_gray) # show this
end

# ╔═╡ 037bc839-97e0-44d9-a000-c3cb4a5dcf99
with_terminal() do
  println("Rank $rsvd approximation")
  println("  Error % = $(100*norm(A-svd_approx_gray)/norm(A))")
  println("  parameters = $nparam_gray")
end

# ╔═╡ 36dccd6c-0024-4f49-88b0-f6b5b2901f5d
md"""
## Number of parameters vs. Error

It's super handy to look at the amount of error compared with the number of parameters. In this case, we look at error as a percent of total image. This is just the 2-norm difference of the image - approximation as a vector (also called Frobenius 2-norm), divided by the norm of the image.
"""

# ╔═╡ 822c2c32-a3e1-40a3-a98f-6fa11a660035
begin
  function get_nparam_approx_data(reference_im, r, approx_fun)
	approx_im, nparam = approx_fun(r)
	return nparam, 100*norm(approx_im-reference_im)/norm(reference_im)
  end
  map(r->get_nparam_approx_data(A, r, standard_svd_approx), 1:20)
end

# ╔═╡ 74cc32bb-dc19-4a1e-9b6a-79fb782e00b3
begin
	scatter(map(r->get_nparam_approx_data(A, r, standard_svd_approx), 1:20), label="Image-as-Matrix SVD")
	ylabel!("Error %")
	xlabel!("Parameters")
end

# ╔═╡ 11a12eb5-2928-47c1-87c5-2faed7e1ba02
md""" # Image as Matrix-of-Tiles SVD Example

The entire point of this workbook is to convince you there is a better approach based on thinking of the matrix as a collection of image tiles or image patches. This leads to _better_ approximation results but also leads to _deeper_ questions and _more interesting analysis_. In particular, we will see the discrete cosine transform (DCT), which is used in the JPEG format, arise naturally via the SVD!

The idea is that we divide the image into little tiles, and then represent the images as a matrix where each column is a _vector_ representation of the little tile.

Here is what this looks like for 16x16 tiles. So each column of the matrix will be `16*16 = 256` sized.
"""

# ╔═╡ eae7193c-f287-4f62-bb50-cfc998638ba9
begin
""" This is a really crappy implementation to show tiles. This is done just to make it possible. """
function show_tiles(im, tilesz::Integer)
  img = Gray.(im)
  imc = RGB.(img) # convert back to RGB
  for i=tilesz:tilesz:size(imc,1) # show lines
    imc[i,:] .= colorant"yellow"
  end
  for j=tilesz:tilesz:size(imc,2) # show lines
    imc[:,j] .= colorant"yellow"
  end
  imc
end
show_tiles(im, 16)
end

# ╔═╡ c3fd4712-7802-4685-aa5b-4b9d8f81aeab
md""" #### Pick the tile size below

> We only support a few tile sizes because we set this up to be simple where the image has to evenly divide into tiles.

- Tile size $(@bind tilesz_str Select(string.([4,8,16,32,64]) .=> string.([4,8,16,32,64]), default="8"))
"""

# ╔═╡ cf8ba3f6-3347-49e5-ad5a-92518a4755b7
tilesz = parse(Int, tilesz_str)

# ╔═╡ e953db96-42a3-4506-ad68-3e84cabc47bd
md"""
The following code builds the `pixels-by-tiles` matrix `X`. It uses the `TiledIteration` package to iterate over each tile and then lays it into a column of the matrix by just `hcat`ing all the output from the tile iterator.

This is done on the color image to show the different color regions.
"""

# ╔═╡ 606d5856-7051-4633-893a-c33f691dd8ac
begin
	function tile_reshape(im,tilesz)
	  hcat(
			map(tile -> vec(im[tile...]),
				TileIterator(axes(im), (tilesz,tilesz))
			)...
		)

	end
	X = tile_reshape(im, tilesz)
end

# ╔═╡ 46d589b2-87a1-48a3-86a3-f784a17d6724
md"""We actually compute on the greyscale matrix as before."""

# ╔═╡ 90233a0d-c04b-4216-8b0b-9d7c1912150f
Gray.(X)

# ╔═╡ 02811d5e-7948-4cf1-9330-dd95e7bcd55a
md"""
On important note is that we aren't making any more data when we do this. The values of `X` and `A` are the same, and there are exactly the same number of them. The key difference is that those in `X` are structured in terms of tiles. This allows us to exploit different structure in the image.
"""

# ╔═╡ 3632bf49-d6bf-4de0-9450-5233517db4ef
length(X), length(A)

# ╔═╡ 342cc50a-81e9-4688-829a-e7018cff30d6
md"""
As before, we compute the SVD of the matrix X. The output `tU` means `tiled U` as a short-hand. But we need to remember that this needs to be interpreted differently.
"""

# ╔═╡ 25215211-7c33-41db-a327-7c2f45500b3c
tU,tS,tV = svd(Float64.(Gray.(X))); # output hidden because it's boring

# ╔═╡ 2ae1b465-2f67-483c-bd56-af0b859940f7
md"""
Here, we pick a rank for the Tiled SVD approximation.

**Rank for Tiled SVD approximation** (Usually you see something around 3-5.)
- `rtile` $(@bind rtile Slider(1:25,default=3))

The reconstruction code given a tile size is very different. We first build an rank `rtile` approximation of the matrix `X` in the matrix `Y`. Then we read each column of `Y` back into a matrix `B` that approximates the matrix `A` by using the TileIterator again.
"""

# ╔═╡ c2233159-59f4-4516-859c-5c25170d53aa
md" `rtile` = $(rtile)"

# ╔═╡ 3e7902ab-4f75-4d14-901b-d2e380996c8f
begin
	function tiled_svd_approx(r::Int, tilesz::Int; tU, tS, tV)
		Y = tU[:,1:r]*Diagonal(tS[1:r])*tV[:,1:r]'
		B = copy(A)
		for (i,tile) in enumerate(TileIterator(axes(A), (tilesz,tilesz)))
		  B[tile...] .= reshape(Y[:,i], (tilesz,tilesz))
		end
		return B, length(tU[:,1:r])+length(tV[:,1:r])
	end
	tile_svd_approx_gray, nparam_tile_gray = tiled_svd_approx(rtile,tilesz; tU, tS, tV)
	Gray.(tile_svd_approx_gray) # show this
end

# ╔═╡ 49bd0281-151c-4306-bcf8-2eb8c618e04d
begin
punchline_str = begin
	iam_approx,nparam_iam = standard_svd_approx(rsvd)
	iat_approx,nparam_iat = tiled_svd_approx(rtile,tilesz; tU, tS, tV)
	iam_error = 100*norm(iam_approx-A)/norm(A)
	iat_error =  100*norm(iat_approx-A)/norm(A)
	md"""
| Image-as-Matrix SVD Approximation | Images-as-$tilesz x $tilesz Tiles Approximation |
|:---------:|:-----------:|
| $(Gray.(iam_approx)) | $(Gray.(iat_approx)) |
| $(nparam_iam) parameters | $(nparam_iat) parameters |
| Rank $(rsvd) approximation | Rank $(rtile) approximation |
| $(round(iam_error;digits=2)) % Error | $(round(iat_error;digits=2)) % Error |
"""
end
md"> Expand to see code to generate punchline figures"
end

# ╔═╡ b9473b23-957b-4d23-8ad5-893c6c4063cf
md"""# Pedagogical image-as-matrix compared with image-as-tiles SVD examples

A common example in linear algebra classes or matrix computations or numerical computing classes is to take the SVD of a matrix that represents an image. Typically, the matrix is just taken _to be_ the image, which we call _image-as-matrix_ here. Then we create a low-rank approximation (from the SVD) of that matrix, and show the result as an image.

This example is curious for a few reasons.
- First, the _image-as-matrix_ view creates a matrix that is a rather odd linear operator. It's synatically valid, of course. This leads into a long discussion of _what is a matrix_ (but not _what is **the** matrix!) that we'll avoid.
- Second, the goal is to approximate the image! We'll see there is a better strategy to approximate the underlying images (and that makes much more sense as a linear operator).

The punchline is below.

$(punchline_str)
"""

# ╔═╡ 8c059291-0564-4ec7-80a0-13432962254a
with_terminal() do
  println("Tiled rank $rtile approximation")
  println("  Error % = $(100*norm(A-tile_svd_approx_gray)/norm(A))")
  println("  parameters = $nparam_tile_gray")
end

# ╔═╡ b688c31c-5394-40e8-9e6c-2f3fb64fffaf
md""" ## The number of parameters when `rtile = rsvd` are not the same!
Note that for the `rtile=3` approximation (the default set in this notebook) there are 10,560 parameters. The `rsvd=15` approximation (the default set here) there are 14,400 parameters for the coffee image. (This will vary depending on what image you have.

So what we want to look at is Error % as a function of tile-size, rank, and number of parameters. That's what we show below.

The standard Image-as-Matrix approximation is shown in red. We get better approximations for the tiled examples (on the coffee image) if the curve is below the red curve. This happens for many different tile sizes. But it does not for `64x64` tiles. This is because there are just too few such tiles for a small image. If you want to see those work, try a bigger image like `earth`.
"""

# ╔═╡ 3d9f7a23-b5d3-49ef-ae6f-d3644ba0bcf8
begin
	plot(map(r->get_nparam_approx_data(A, r, standard_svd_approx), 1:20), label="Image-as-Matrix SVD",marker=:dot,linewidth=3,color=:red)

	for ptilesz = [4,8,16,32,64]
		pX = tile_reshape(A, ptilesz)
		ptU,ptS,ptV = svd(pX)
		plot!(map(r->get_nparam_approx_data(A, r, r->tiled_svd_approx(r,ptilesz;tU=ptU,tV=ptV,tS=ptS)), 1:min(20,min(size(pX)...))), label="Image-as-$(ptilesz)x$(ptilesz)-Tiles SVD " ,marker=:dot)
	end

	ylabel!("Error %")
	xlabel!("Parameters")
	xlims!(1,length(U[:,1:20]) + length(V[:,1:20]))
end

# ╔═╡ 860b01bd-47e0-41f2-bd05-f5f3957b8370
md""" Of course, it's often useful to look at the picture the other way. For a fixed level of error?, is the tiled approximation better than the Image-as-Matrix approximation?

What we want to see for a given error level (x-axis) which approximation has fewer parameters. So we take their ratio. In this case, let's look at the ratio of so that ratios above 1 show that the tiled approximation is better.
"""

# ╔═╡ 9fed5b5e-0106-484f-898c-1d664033b339
begin
	# Here we want to make the opposite plot, and show the number of parameters used for a given approximation level.
	# This uses the fixed-tile size above (tilesz)
	svd_results = map(r->get_nparam_approx_data(A, r, standard_svd_approx), 1:20)
	tiled_svd_results = map(r->get_nparam_approx_data(A, r,
								r->tiled_svd_approx(r,tilesz; tU,tS,tV)),
							1:min(20,min(size(X)...)))
	relapprox = Tuple{Float64,Float64}[]
	for (nparam,approx_level) in svd_results
		maxparams = maximum(first, tiled_svd_results)
		for (nparam_tile,tile_approx_level) in tiled_svd_results
			if tile_approx_level < approx_level
				push!(relapprox, (approx_level, nparam/nparam_tile))
				break
			end
		end
	end

	plot(relapprox,marker=:dot, size=(600,500), label="Image as $(tilesz)x$(tilesz) Tiles" )
	xlabel!("Error %")
	hline!([1.0], label="", linewidth=2)
	annotate!(maximum(relapprox[1]), 1.0,
			("↑ Tiled Approx is Better", :right, :bottom))
	annotate!(maximum(relapprox[1]), 1.0,
			("↓ Standard Approx is Better", :right, :top))

	ylabel!("Image-as-Matrix Approx Params / Tiled Approx Params ")
end

# ╔═╡ 5b99d2df-fb5a-4fcc-b5f9-f44ca5328814
md"""
# Interpreting the Tiled SVD basis.

One of the best reasons to look at the tiled approximation is that the left singular vectors `tU` are themselves `image patches` that we can look at!

For this section, you get the best results with the `cat` picture or `earth`, but the others aren't too far off. So I suggest switching to the cat picture.
"""

# ╔═╡ e869b69b-b79f-4bd6-bd2f-6015e22028d9
md""" > The following is a little function to upscale images to make them bigger to see. It's a good exercise in multidimensional array programming in Julia as this will work for vectors, matrices, tensors, etc."""

# ╔═╡ 00c0d888-f192-4ddc-9a94-5cf9624ddd26
begin
	function upscale(X::AbstractArray, factor::Int)
		sz2 = factor.*size(X)
		Y = similar(X, sz2...)
		I1 = oneunit(first(CartesianIndices(X)))

		for xi = CartesianIndices(X)
			start = factor*(xi-I1)+I1
			last = factor*(xi)
			for yi in start:last
				Y[yi] = X[xi]
			end
		end
		return Y
	end
	upscale([1 2;3 4],3)
end

# ╔═╡ def9d9d3-2185-4d16-a63f-e415bb91eac3
md""" Below are the first 16 singular vectors of `tU` resized to be little image tiles themselves. (column-wise organization! so 1-4 are in first column...)

The colors are "green" for negative and "purple" for positive and "white" for near zero. These are scaled so that "white" always indicates 0.

> Remember that singular vectors have an ambiguous sign. So you may see different colors, but the patterns should be the same.

The first vector should be constant, which should be either pure green or purple. (Up to what your eyeball/monitor combo can show.)

But the second and third vector for `cat` show a roughly. vertical and horizontal split. You may see slight angles for other images.

Then we see other `sin`/`cos` like before in the next few. In particular, the vectors seem to alternate more frequently as the index gets higher."""

# ╔═╡ 5a491365-da8e-43bb-987d-24fccb9a2e6d
Z = mosaicview(map(i->colorview(RGB,
                      permutedims(applydivergingcolourmap(
                        reshape(tU[:,i], (tilesz,tilesz)), cmap("D2"), 0),
                      (3,1,2))),
               1:16)..., nrow=4, npad=5, fillvalue=127) |> Z -> map(clamp01nan, Z) |> Z -> upscale(Z,8)

# ╔═╡ ade6e8f2-bd1b-4569-8efb-aae3d42ee340
md""" ### Introducing the Discrete Cosine Transform

At some point in your life, you either learned about, or will learn more about, the Fourier transform and the discrete cosine transform (DCT). These are mainstays of signal processing and enable much of the modern world to work.

This isn't the place to introduce them. You have wikipedia and google... and youtube... there are excellent explainations here.

In particular, the DCT underlies the JPG file format.

### The DCT in JPG

As a rough guide, the way the JPG files are stored is that the image is broken up into little 8x8 tiles (like we did here!) and then they apply a DCT transform to each block to determine how much information to keep. This is basically the same as our low-rank approximation, except we use the SVD basis instead of the DCT.

### The DCT and SVD

The reason I mention the DCT here is that the discrete cosine transform "looks a lot like" the picture above. (Not exactly the same, but very similar.)

Here's the same picture for the DCT.
"""

# ╔═╡ 5c53bc42-5f09-482e-81a8-43c1cfa71be6
begin
	function idct_id(offset,sz)
	  A = zeros(sz,sz)
	  # find the triangular index
	  for j=1:sz
		# jth diagonal...
		if offset > j
		  offset -= j
		else
		  # in this diagonal...
		  A[j-offset+1,offset] = 1
		  break
		end
	  end
	  return A
	end
	Zdct = mosaicview(map(i->colorview(RGB,
                      permutedims(applydivergingcolourmap(
                        idct(idct_id(i,tilesz)), cmap("D2"), 0),
                      (3,1,2))),
               1:16)..., nrow=4, npad=5, fillvalue=127) |> Z -> map(clamp01nan, Z) |> Z->upscale(Z,8)
end

# ╔═╡ 73e63f45-1458-40ca-9ec4-bb1c00932300
md"""
Yeah, not the same, but pretty darn close!  If you want it to be closer, try a bigger image. (Upload your own image that's around 1500x1000 and use 8x8 blocks... for instance.) Using a high-res Windows XP wallpaper `bliss` image (at some time availble from this URL <https://i.imgur.com/uGRFZEs.jpg>) showed excellent results.
> In fact, one time when I was looking at this notebook, I thought I had deleted the SVD-based result because I couldn't quickly tell the result on `bliss` apart from the exact result as I had been expected.
"""

# ╔═╡ 211f12bd-6d77-4d78-b486-b8f381bdde8d
md"""
### Quantifying the SVD as DCT approximation.

The DCT is a matrix transformation. That means we an ask the DCT to tell us how much those SVD-like images are `DCT-like`.

The following picture quantifies this and `it's really good` is the quick summary.  The way you can see this is because each `image` is highly localized. You basically only see a single pixel or handful of pixels with `large values`. Also, when you  ee multiple pixels with large values, they are in the same diagonal band. (It turns out that the diagonal bands are basically energy modes determined by the number of sin/cos oscillations, so we are mixing up nearly indistinguishable energy modes.) This picture would be an excellent starting point to talk more about the DCT.
"""

# ╔═╡ 224dc71b-abdf-4ad4-9821-63b1ec2021ff
dct_of_Z = mosaicview(map(i->colorview(RGB,
                      permutedims(applydivergingcolourmap(
                        dct(reshape(tU[:,i], (tilesz,tilesz))), cmap("D2"), 0),
                      (3,1,2))),
               1:16)..., nrow=4, npad=5, fillvalue=127) |> Z -> map(clamp01nan, Z) |> Z->upscale(Z,8)

# ╔═╡ 4e1f7928-3f7e-4140-a314-a56b727b5de0
md""" Just for reference (and because it's easier) here is the same picture if the result were `exact` and we used the DCT basis itself. Remember again that the sign/color is ambiguous so we can change that."""

# ╔═╡ 28279f1f-d0bc-48d5-ad95-ef0600c774f4
dct_of_Zdct = mosaicview(map(i->colorview(RGB,
                      permutedims(applydivergingcolourmap(
                        dct(idct(idct_id(i,tilesz))), cmap("D2"), 0),
                      (3,1,2))),
               1:16)..., nrow=4, npad=5, fillvalue=127) |> Z -> map(clamp01nan, Z) |> Z->upscale(Z,8)

# ╔═╡ fed9cdd4-4262-4a60-aee6-12b611c373b2
md""" ## (Under construction) Interpreting the right tiled SVD basis.

The matrix `tV` gives the expansion."""

# ╔═╡ ba5cbe49-05da-438e-b9aa-36e77843e778
heatmap(tV';colors=cmap("D2"))

# ╔═╡ fa381b5f-fee2-4974-9b0d-0b317c6fe921
histogram(vec(tV))

# ╔═╡ 23b74436-7020-4a33-82d6-5b2cb50c2a39
histogram(vec(V))

# ╔═╡ f82f8ba0-c1ae-42d8-a63a-bf605b158fad
plot(abs.(tV[1:5,:]'), linecolor=:black, alpha=0.05, legend=false)

# ╔═╡ fd45332e-2980-474a-9237-31031c4a2c07
begin
	thresh = 0.02
	sum(x->abs(x) > thresh, tV), sum(x->abs(x) > thresh, V)
end

# ╔═╡ ddf50450-5582-4ab0-9e1f-8fe0b9739edb
md""" # Studying failure modes and towards a model.

When you see results like this, you should ask _why_ is it the case that we get the DCT result from the SVD on tiled matrices. Maybe that's the just the general answer? Here is what happens on a matrix of random values instead.

I wish I had an image that showed this instead of just random noise, but we'll work with that!
"""

# ╔═╡ e26f91f9-82a5-461c-a26d-2671bb17fff4
begin
	# generate random values in a matrix with the same size as the tiled matrix X
	Random.seed!(0)
	Xr = rand(size(X)...)
	rU,rS,rV = svd(Xr)
end;

# ╔═╡ a0ab28cd-ff89-4678-a8c0-a75c384273db
mosaicview(map(i->colorview(RGB,
                      permutedims(applydivergingcolourmap(
                        reshape(rU[:,i], (tilesz,tilesz)), cmap("D2"), 0),
                      (3,1,2))),
               1:16)..., nrow=4, npad=5, fillvalue=127) |> Z -> map(clamp01nan, Z) |> Z -> upscale(Z,8)

# ╔═╡ 3f072542-81b2-445b-abeb-c5d7686e53a8
md""" What goes wrong is that the random data _only_ has high frequency information. So there is no low frequency inforamtion (nearly constant patches) to drive it towards a DCT-like answer."""

# ╔═╡ 0cd0586e-2efc-4f88-89e5-99f813435e3f
md"""We can create low-frequency information by bluring a little tile or patch we generate.

- `blur_amount` $(@bind blur_amount Slider(0:3, default=1))
"""

# ╔═╡ 2e6dbdfa-71da-4eef-847c-f78879121822
begin
	Random.seed!(0)
	rand_tile = Gray.(rand(tilesz, tilesz))
	blur_tile = imfilter(rand_tile, Kernel.gaussian(blur_amount))
	mosaicview(rand_tile, blur_tile, nrow=1, npad=5, fillvalue=127) |> Z -> upscale(Z,8)
end

# ╔═╡ b7156110-cf28-4df4-8068-9b9131122076
md""" Now what we'll do is generate a matrix like X, but where each column simulates a random patch with a random amount of blurring. """

# ╔═╡ 96587e36-0ded-480d-aa51-3b1177306bd7
begin
function random_mixed_frequency(m::Int, n::Int)
	@assert tilesz*tilesz == m # need multiple of tile-size
	R = zeros(m,n)
	for j=1:n
		T = rand(tilesz,tilesz) # generate a random tile
		# randomly blur, 0 = no blur, 3 = lots!), see above
		T = imfilter(T, Kernel.gaussian(rand(0:10)))
		R[:,j] = vec(T)
	end
	return R
end

	Xr2 = random_mixed_frequency(size(X)...)
	rU2,rS2,rV2 = svd(Xr2)
end;

# ╔═╡ 5098d01f-bb09-4155-9994-c77e3e742e72
mosaicview(map(i->colorview(RGB,
                      permutedims(applydivergingcolourmap(
                        reshape(rU2[:,i], (tilesz,tilesz)), cmap("D2"), 0),
                      (3,1,2))),
               1:16)..., nrow=4, npad=5, fillvalue=127) |> Z -> map(clamp01nan, Z) |> Z -> upscale(Z,8)

# ╔═╡ c687100b-96cd-4788-8bd5-e53e952a2762
md"""
While this looks like we get the first few modes/vectors correct. We can check with the idct transform as before.
"""

# ╔═╡ afcee1f2-f5b8-4346-87de-21afc5bbce16
mosaicview(map(i->colorview(RGB,
                      permutedims(applydivergingcolourmap(
                        dct(reshape(rU2[:,i], (tilesz,tilesz))), cmap("D2"), 0),
                      (3,1,2))),
               1:16)..., nrow=4, npad=5, fillvalue=127) |> Z -> map(clamp01nan, Z) |> Z->upscale(Z,8)

# ╔═╡ 10602dbb-d960-4488-bfff-8e4c8c8a46a7
md"""
> You can play around with this a lot. The key line above is `T = imfilter(T, Kernel.gaussian(rand(0:3)))`. See what happens if you use `T = imfilter(T, Kernel.gaussian(rand(0:1)))` to look at only a tiny bit of blurring, or `T = imfilter(T, Kernel.gaussian(rand(5:10)))` to look at tons of blurring...
"""

# ╔═╡ 655d163f-7845-430e-b784-fb0bfa7550c3
md""" # Varimax and sparse bases from SVD

> This is a more advanced concept. Please skip it unless you want to learn a lot more yourself or you already known about things like Sparse PCA and 1-norm approximations and the modern sparsity concept. I mention it only to show how this idea can be used as a jumpstart for all sorts of ideas in applied math and signal processing and stats.

There is a method in statistics to `find` more interpretable or `sparse` bases from an orthgonal basis called "varimax" or "orthomax". Karl Rohe pointed out some interesting comparisons to "sparse PCA" (Rohe and Zeng, arXiv 2020 <https://arxiv.org/abs/2004.05387>). Sparse PCA is designed around convex programming, 1-norms, and the like. But varimax is just a fast heurstic, but can be applied to any SVD output. So let's see what the varimax basis here is.

"""

# ╔═╡ 32064bb3-b50a-42b0-aa2d-5cef91a25398
md""" Compute a rank-16 basis transformation via orthomax. """

# ╔═╡ 351840fb-81a9-47c1-8635-1c1887d7f04a
begin
	# Ahh, I stole this from a github issue/pull request for ...
	# https://github.com/JuliaStats/MultivariateStats.jl/pull/130/files
	## Comparison to other implementations
	# The implementation of varimax in R row-normlises the input matrix before
	# application of the algorithm and rescales the rows afterwards.
	## Reference
	# Mikkel B. Stegmann, Karl Sjöstrand, Rasmus Larsen, "Sparse modeling of
	# landmark and texture variability using the orthomax criterion,"
	# Proc. SPIE 6144, Medical Imaging 2006: Image Processing, 61441G
	# (10 March 2006); doi: 10.1117/12.651293
	function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ)
		n, p = size(F)
		if n < 2
			return (F, Matrix{eltype(F)}(I, p, p))
		end

		# Warm up step
		# Do one step. If that first step did not lead away from the identity
		# matrix enough use a random orthogonal matrix as a starting point.
		M = svd(F' * (F .^ 3 - γ / n * F * Diagonal(vec(sum(F .^ 2, dims=1)))))
		R = M.U * M.Vt
		if norm(R - Matrix{eltype(R)}(I, p, p)) < ϵ
			R = qr(randn(p, p)).Q
		end

		# Main iterations
		d = 0
		lastchange = NaN
		converged = false
		for i in 1:maxiter
			dold  = d
			B = F * R
			M = svd(F' * (B .^ 3 - γ / n * B * Diagonal(vec(sum(B .^ 2, dims=1)))))
			R = M.U * M.Vt
			d = sum(M.S)
			lastchange = abs(d - dold) / d
			if lastchange < ϵ && i >= miniter
				converged = true
				break
			end
		end

		converged || throw(ConvergenceException(maxiter, lastchange, ϵ))

		(F * R, R)
	 end
	tU2 = orthomax(tU[:,1:16], 1.0, 20, 1000, 1e-12)[1]
end

# ╔═╡ c5f90874-bf6e-4179-ac74-d57d2253b3a2
md"""
The varimax basis is below. We see sparsity in terms of only having small regions of nonzero values.

The result is typically little "blobs" of nonzeros. This is reasonable. If you had to represent 8x8 images with only 16 _sparse_ images (those with only a few nonzero entries). Then having all 16 4x4 activation blocks is probably the best you can do!

The `coffee` image looks pretty good here.
"""

# ╔═╡ b617fb5b-6f20-4195-a13c-93809bf77539
Z2 = mosaicview(map(i->colorview(RGB,
                      permutedims(applydivergingcolourmap(
                        reshape(tU2[:,i], (tilesz,tilesz)), cmap("D2"), 0),
                      (3,1,2))),
               1:16)..., nrow=4, npad=5, fillvalue=127) |> Z -> map(clamp01nan, Z) |> Z-> upscale(Z,8)

# ╔═╡ 6ae2aaf4-29a3-4da8-947b-044ba3baed26
md"""
This "blob" analysis suggests that if we ask varimax/Orthomax for size-9 basis, we should get bigger blobs.

Which we do! See below. Again, this will depend on the image you are looking at. """

# ╔═╡ 5fe6838b-e244-4135-be25-1f2575b289b5
begin
	tU3 = orthomax(tU[:,1:9], 1.0, 20, 1000, 1e-12)[1]
	Z3 = mosaicview(map(i->colorview(RGB,
                      permutedims(applydivergingcolourmap(
                        reshape(tU3[:,i], (tilesz,tilesz)), cmap("D2"), 0),
                      (3,1,2))),
               1:9)..., nrow=3, npad=5, fillvalue=127) |> Z -> map(clamp01nan, Z) |> Z-> upscale(Z,8)
end

# ╔═╡ Cell order:
# ╟─b9473b23-957b-4d23-8ad5-893c6c4063cf
# ╟─49bd0281-151c-4306-bcf8-2eb8c618e04d
# ╟─cbd2ea36-7ef6-4ea8-82ac-34b3e96ed8d5
# ╟─59c37536-fff2-480a-80bc-f97125e117af
# ╠═ccbc7c67-e16a-48f1-a927-c002257fcc54
# ╟─2c8a2738-77b1-4803-91d3-9b5172e4a94a
# ╟─2ac8cb5d-ffff-4561-bc08-2ffc368338e2
# ╟─3112233a-543d-48f5-823f-5b1e1cafcb61
# ╟─03604f42-66ce-4a39-b1b1-f5fd49bfccdb
# ╟─e6e10ee8-2bcb-4743-86be-ffac6b591b96
# ╟─2d4125ac-5c10-4150-bc00-e2cfcbfaf5ad
# ╟─52d8a212-0b34-4306-b410-e3d3dce7de3f
# ╠═fd4b3fef-9df2-4253-81b4-98a9c59a0176
# ╟─13b28401-d752-4dae-9ddf-bfd59f3ab995
# ╟─08c9f59e-4408-409a-b324-36dedc005c65
# ╟─4632bc52-9106-4d05-8f6f-a189e577388d
# ╟─6a13f71e-983b-4d48-8787-8d58abbebbf5
# ╟─17dae8a4-010a-4816-b5a4-c9448d3cf8b1
# ╟─fddb68e0-adc2-4ff5-908d-e0204bf074e3
# ╠═c372842c-a529-4ad1-8edb-47d755018dbc
# ╠═dbace9f6-6b02-4e91-88a4-92fe27ec063d
# ╟─037bc839-97e0-44d9-a000-c3cb4a5dcf99
# ╟─36dccd6c-0024-4f49-88b0-f6b5b2901f5d
# ╠═822c2c32-a3e1-40a3-a98f-6fa11a660035
# ╠═74cc32bb-dc19-4a1e-9b6a-79fb782e00b3
# ╟─11a12eb5-2928-47c1-87c5-2faed7e1ba02
# ╟─eae7193c-f287-4f62-bb50-cfc998638ba9
# ╟─c3fd4712-7802-4685-aa5b-4b9d8f81aeab
# ╟─cf8ba3f6-3347-49e5-ad5a-92518a4755b7
# ╟─e953db96-42a3-4506-ad68-3e84cabc47bd
# ╠═606d5856-7051-4633-893a-c33f691dd8ac
# ╟─46d589b2-87a1-48a3-86a3-f784a17d6724
# ╟─90233a0d-c04b-4216-8b0b-9d7c1912150f
# ╟─02811d5e-7948-4cf1-9330-dd95e7bcd55a
# ╠═3632bf49-d6bf-4de0-9450-5233517db4ef
# ╟─342cc50a-81e9-4688-829a-e7018cff30d6
# ╠═25215211-7c33-41db-a327-7c2f45500b3c
# ╟─2ae1b465-2f67-483c-bd56-af0b859940f7
# ╟─c2233159-59f4-4516-859c-5c25170d53aa
# ╠═3e7902ab-4f75-4d14-901b-d2e380996c8f
# ╟─8c059291-0564-4ec7-80a0-13432962254a
# ╟─b688c31c-5394-40e8-9e6c-2f3fb64fffaf
# ╠═3d9f7a23-b5d3-49ef-ae6f-d3644ba0bcf8
# ╟─860b01bd-47e0-41f2-bd05-f5f3957b8370
# ╟─9fed5b5e-0106-484f-898c-1d664033b339
# ╟─5b99d2df-fb5a-4fcc-b5f9-f44ca5328814
# ╟─e869b69b-b79f-4bd6-bd2f-6015e22028d9
# ╠═00c0d888-f192-4ddc-9a94-5cf9624ddd26
# ╟─def9d9d3-2185-4d16-a63f-e415bb91eac3
# ╠═5a491365-da8e-43bb-987d-24fccb9a2e6d
# ╟─ade6e8f2-bd1b-4569-8efb-aae3d42ee340
# ╠═5c53bc42-5f09-482e-81a8-43c1cfa71be6
# ╟─73e63f45-1458-40ca-9ec4-bb1c00932300
# ╟─211f12bd-6d77-4d78-b486-b8f381bdde8d
# ╠═224dc71b-abdf-4ad4-9821-63b1ec2021ff
# ╟─4e1f7928-3f7e-4140-a314-a56b727b5de0
# ╠═28279f1f-d0bc-48d5-ad95-ef0600c774f4
# ╟─fed9cdd4-4262-4a60-aee6-12b611c373b2
# ╠═ba5cbe49-05da-438e-b9aa-36e77843e778
# ╠═fa381b5f-fee2-4974-9b0d-0b317c6fe921
# ╠═23b74436-7020-4a33-82d6-5b2cb50c2a39
# ╠═f82f8ba0-c1ae-42d8-a63a-bf605b158fad
# ╠═fd45332e-2980-474a-9237-31031c4a2c07
# ╟─ddf50450-5582-4ab0-9e1f-8fe0b9739edb
# ╠═e26f91f9-82a5-461c-a26d-2671bb17fff4
# ╠═a0ab28cd-ff89-4678-a8c0-a75c384273db
# ╟─3f072542-81b2-445b-abeb-c5d7686e53a8
# ╟─0cd0586e-2efc-4f88-89e5-99f813435e3f
# ╟─2e6dbdfa-71da-4eef-847c-f78879121822
# ╟─b7156110-cf28-4df4-8068-9b9131122076
# ╠═96587e36-0ded-480d-aa51-3b1177306bd7
# ╠═5098d01f-bb09-4155-9994-c77e3e742e72
# ╟─c687100b-96cd-4788-8bd5-e53e952a2762
# ╠═afcee1f2-f5b8-4346-87de-21afc5bbce16
# ╟─10602dbb-d960-4488-bfff-8e4c8c8a46a7
# ╟─655d163f-7845-430e-b784-fb0bfa7550c3
# ╠═32064bb3-b50a-42b0-aa2d-5cef91a25398
# ╠═351840fb-81a9-47c1-8635-1c1887d7f04a
# ╟─c5f90874-bf6e-4179-ac74-d57d2253b3a2
# ╠═b617fb5b-6f20-4195-a13c-93809bf77539
# ╟─6ae2aaf4-29a3-4da8-947b-044ba3baed26
# ╠═5fe6838b-e244-4135-be25-1f2575b289b5

max_stars_repo_name: dgleich/better-svd-for-images

max_issues_repo_name: dgleich/better-svd-for-images

max_forks_repo_name: dgleich/better-svd-for-images


text:
mutable struct viewer_load_robot_t <: LCMType
    num_links::Int32
    link::Vector{viewer_link_data_t}
end

@lcmtypesetup(viewer_load_robot_t,
    link => (num_links,)
)

ext: jl


text:
module GridGraphsTests

using Test
using Gridap
using Gridap.CellValuesGallery
using Gridap.GridGraphs: GridGraphFromData
using Gridap.GridGraphs: FullGridGraphFromData

c2v = [[1,2,4,5],[2,3,5,6]]
v2c = [[1,],[1,2],[2,],[1,],[1,2],[2,]]
c2e = [[1,2,3,4],[5,6,4,7]]
e2c = [[1,],[1,],[1,],[1,2],[2,],[2,],[2,]]
e2v = [[1,2],[4,5],[1,4],[2,5],[2,3],[5,6],[3,7]]
v2e = [[1,3],[1,4,5],[5,7],[2,3],[2,4,6],[6,7]]
v2v = [[1,],[2,],[3,],[4,],[5,],[6,]]
e2e = [[1,],[2,],[3,],[4,],[5,],[6,],[7,]]
c2c = [[1,],[2,]]

c2v = CellValueFromArray(c2v)
v2c = CellValueFromArray(v2c)
c2e = CellValueFromArray(c2e)
e2c = CellValueFromArray(e2c)
e2v = CellValueFromArray(e2v)
v2e = CellValueFromArray(v2e)
v2v = CellValueFromArray(v2v)
e2e = CellValueFromArray(e2e)
c2c = CellValueFromArray(c2c)

primal = [c2v,c2e,c2c]
dual = [v2c,e2c,c2c]
graph = GridGraphFromData(primal,dual)
nd = 2
test_grid_graph(graph,nd)

data = Matrix{IndexCellVector}(undef,(3,3))
data[1,:] = [v2v,v2e,v2c]
data[2,:] = [e2v,e2e,e2c]
data[3,:] = [c2v,c2e,c2c]
graph = FullGridGraphFromData(data)
test_full_grid_graph(graph,nd)

end # module


text:
"""
	struct MissingNode{D} <: AbstractNode
		data::D
		present::Vector{Bool}
	end

	Missing node adds a generic support for missing values to nodes. It wraps the node,
	such that `data` field contain only non-missing values, which are indicated in the
	`present` mask.

"""
struct MissingNode{D} <:AbstractNode
	data::D
	present::Vector{Bool}
end

MissingNode(d) = MissingNode(d, fill(true, nobs(d)))

Base.ndims(x::MissingNode) = Colon()
LearnBase.nobs(a::MissingNode) = length(a.present)
LearnBase.nobs(a::MissingNode, ::Type{ObsDim.Last}) = nobs(a.present)

function Base.reduce(::typeof(Mill.catobs), as::Vector{T}) where {T<:MissingNode}
    data = reduce(Mill.catobs, [x.data for x in as])
    present = reduce(vcat, [x.present for x in as])
    MissingNode(data, present)
end

get_present_index(present, i::Int) = sum(view(present, 1:i))
get_present_index(present, ii::Vector{Int}) = map(i -> get_present_index(present, i), ii)

function Base.getindex(x::MissingNode, i::VecOrRange)
	p = x.present[i]
	!any(p) && return(MissingNode(x.data[1:0], p))
	ii = get_present_index(x.present, i[p])
	MissingNode(x.data[ii], p)
end

Base.hash(e::MissingNode{D}, h::UInt) where {D} = hash((D, e.data, e.present), h)
Base.:(==)(e1::MissingNode{D}, e2::MissingNode{D}) where {D} = e1.data == e2.data && e1.present == e2.present


text:
"""
    fulfill!(ts, stt, acc, libs)
1. create fulfillment flow given `stt.current_stock` and \
`libs.demand_realization` at current timestep `ts` + 1
2. update `stt.current_stock` given the fulfillment flow and \
add it to the `acc.demand_fulfillment`
3. append `stt.current_stock` to `acc.inventory_levels`
"""
function fulfill!(ts::Int, stt::States, acc::Accumulators, libs::Libraries)
    to_fulfill = filter(p -> p.periode == ts + 1, libs.demand_realization)
    fulfillment = similar(to_fulfill, 0)

    for r in eachrow(to_fulfill)
        v = min(
            r.value,
            stt.current_stock[
                stt.current_stock.id.==r.id.&&stt.current_stock.pecahan.==r.pecahan,
                :value
            ] |> first
        )
        append!(fulfillment,
            DataFrame(
                id=r.id,
                periode=r.periode,
                pecahan=r.pecahan,
                value=v
            )
        )
    end

    for r in eachrow(fulfillment)
        stt.current_stock[
            stt.current_stock.id.==r.id.&&stt.current_stock.pecahan.==r.pecahan,
            :value
        ] .-= r.value
    end

    append!(acc.demand_fulfillment, fulfillment)
    append!(acc.inventory_levels,
        DataFrame(
            id=stt.current_stock.id,
            periode=ts + 1,
            pecahan=stt.current_stock.pecahan,
            value=stt.current_stock.value
        )
    )

    return nothing
end

fulfill!(sim::Simulation) = fulfill!(sim.t, sim.stt, sim.acc, sim.libs)


text:
export Chebyshev

immutable Chebyshev <: ParametricFamily
    n::Int
end

domain(::Chebyshev) = -1..1

points(c::Chebyshev, T = Float64) = cospi.(((c.n:-1:1)-T(0.5))/c.n)

degf(c::Chebyshev) = c.n

"Recurrence relation for Chebyshev polynomials."
@inline _chebyshev_recurrence(x, T₋₁, T₋₂) = 2*x*T₋₁ - T₋₂, T₋₁

"Initialize T₋₁, T₋₂ for the Chebyshev recurrence."
@inline _chebyshev_recurrence_init{T}(x::T) = x, one(T)

macro _chebyshev_recurrence!(x, T₋₁, T₋₂)
    T₋₁ = esc(T₋₁)
    T₋₂ = esc(T₋₂)
    :(($(T₋₁), $(T₋₂)) = _chebyshev_recurrence($(esc(x)), $(T₋₁), $(T₋₂)))
end

function basis!{T}(c::Chebyshev, x::T, b::AbstractVector{T})
    @argcheck c.n == length(b)
    T₋₁, T₋₂ = _chebyshev_recurrence_init(x)
    for i in 1:c.n
        if i == 1
            b[i] = T₋₂
        else
            if i > 2
                @_chebyshev_recurrence! x T₋₁ T₋₂
            end
            b[i] = T₋₁
        end
    end
    b
end

# note: after benchmarking, it was found that this is faster than Clenshaw,
# and at least as accurate.
function evaluate{T}(c::Chebyshev, θ, x::T)
    T₋₁, T₋₂ = _chebyshev_recurrence_init(x)
    value = zero(T)
    for i in 1:c.n
        if i == 1
            value += T₋₂*θ[i]
        else
            if i > 2
                @_chebyshev_recurrence! x T₋₁ T₋₂
            end
            value += T₋₁*θ[i]
        end
    end
    value
end

function evaluate{T}(c::Chebyshev, θ, vp::ValuePartial{T})
    ## if T(x) = 2xT₋₁(x) - T₋₂(x)
    ## then T′(x) = 2T₋₁(x) + 2xT′₋₁(x) - T′₋₂(x), with T′₂ = 1 and T′₁ = 0
    x = vp.x
    value = partial = zero(T)
    T₋₁ = x
    T₋₂ = T′₋₁ = one(T)
    T′₋₂ = zero(T)
    for i in 1:c.n
        if i == 1
            value += T₋₂*θ[i]
        else
            if i > 2
                @_chebyshev_recurrence! x T′₋₁ T′₋₂
                T′₋₁ += 2*T₋₁
                @_chebyshev_recurrence! x T₋₁ T₋₂
            end
            value += T₋₁*θ[i]
            partial += T′₋₁*θ[i]
        end
    end
    value, partial
end

evaluate{T}(c::Chebyshev, θ, x::Partial{T}) = evaluate(c, θ, ValuePartial(x.x))[2]

function fit!{T}(c::Chebyshev, ys::AbstractVector{T}, θ::AbstractVector{T})
    @argcheck c.n == length(ys) == length(θ)
    x = points(c, T)
    T₋₂ = ones(T, c.n)
    T₋₁ = copy(x)
    for i in 1:c.n
        if i == 1
            θ[1] = sum(ys) / c.n
        else
            if i > 2
                for j in 1:c.n
                    @_chebyshev_recurrence! x[j] T₋₁[j] T₋₂[j]
                end
            end
            θ[i] = dot(T₋₁, ys) *2 / c.n
        end
    end
    θ
end


text:
# TODO: Test this
"""
    stage(gl::GitLink, files::Vector{String}; root::String = "", tout = 60.0, flag = true)

Stage the `files` (by copying them to the GitLink stage folder).
The GitLink Server will upload them in its next iter.
To prevent this set the acknowledge (`flag`) to false.
Staging is usuful (instead of `writewdir`) for acumulating changes and reduce the push rate of the link.
This method will sleep till (or timeout `tout`) the GitLink lock is free (which must by must of the time, but...).
If a `root` (as a kwarg) is provided, the tree structure relative to it will be respected.
Ex: if root = "A" and a file path is "A/B/C.txt" it will be staged as "gl-stage/B/C.txt".
Otherwise (root = ""), the file will be copied in the repo's stage "gl-stage/C.txt".
Returns `true` if the action was successful.
"""
function stage(gl::GitLink, files::Vector{String};
        root::String = "", tout = 60.0, flag = true
    )

    ok_flag = false
    sdir = stage_dir(gl)
    mkpath(sdir)

    lock(gl; tout) do
        if isempty(root)
            for src_file in files
                dest_file = joinpath(sdir, basename(src_file))
                cp(src_file, dest_file; force = true)
            end
        else
            root = abspath(root)
            sdir = abspath(sdir)
            for src_file in files
                src_file = abspath(src_file)
                dest_file = replace(src_file, root => sdir)
                _mkdir(dest_file)
                cp(src_file, dest_file; force = true)
            end
        end
        flag && _set_stage_token(gl)
        ok_flag = true
    end
    return ok_flag
end

"""
    stage(upfun::Function, gl::GitLink; tout = 60.0, flag = true)

Allow to modify the stage directory of the GitLink.
The function `upfun(stage_dir)` will be executed and it should copy/create/modify the files into `stage_dir`.
It is recommended that `upfun` not to be an expensive function.
The GitLink Server will upload the staged files in its next iter.
To prevent this set the acknowledge (`flag`) to false.
Staging is usuful (instead of `writewdir`) for acumulating changes and reduce the push rate of the link.
This method will sleep till (or timeout `tout`) the GitLink lock is free (which must by must of the time, but...).
Returns `true` if the action was successful.
"""
function stage(upfun::Function, gl::GitLink;
        tout = 60.0, flag = true
    )

    ok_flag = false
    sdir = stage_dir(gl)
    mkpath(sdir)

    lock(gl; tout) do
        upfun(sdir)
        flag && _set_stage_token(gl)
        ok_flag = true
    end

    return ok_flag
end
stage(gl::GitLink; kwargs) = stage((sdir) -> nothing, gl; kwargs)


text:
const max_nbit = Ref(20)
const max_nbatch = Ref(128)
const max_nshots = Ref(128)

_checkbit(x::Int) = x <= max_nbit[] ? x : error("maximum number of qubits is $(max_nbit[]), got $x.")
_checkbit(x::String) = (_checkbit(length(x)); x)
_checkbatch(x::Int) = x <= max_nbatch[] ? x : error("maximum number of nbatch is $(max_nbatch[]), got $x.")
_checkshots(x::Int) = x <= max_nshots[] ? x : error("maximum number of nshots is $(max_nshots[]), got $x.")

_checkreg(a::Symbol) = a
_checkop(a::Symbol) = a

function inspect(io::IO, b::AbstractBlock)
    Base.show(io, "text/plain", mat(b))
end

function inspect(io::IO, reg::AbstractRegister)
    Base.show(io, "text/plain", statevec(reg))
end

function print2str(obj)
    io = IOBuffer()
    print(io, obj)
    String(take!(io))
end

function inspect2str(obj)
    io = IOBuffer()
    inspect(io, obj)
    String(take!(io))
end

function setnbits!(n::Int, state)
    state.parseinfo.nbit = Int(n)
    "qubit number is set to $n"
end

function getnbits(state)
    "qubit number is to $(state.parseinfo.nbit)"
end


text:
using LabView0mqJl
include(get_script_path(raw"C:\_LabView_projects\ZMQ\LabView0mqJl.jl\src\Julia\Examples-UserFn.jl"));
setglobals(isOK=true, extant=true);
server_0mq4lv(fns);


text:
"""
    eigsolve(A::AbstractMatrix, [howmany = 1, which = :LM, T = eltype(A)]; kwargs...)
    eigsolve(f, n::Int, [howmany = 1, which = :LM, T = Float64]; kwargs...)
    eigsolve(f, x₀, [howmany = 1, which = :LM]; kwargs...)
    eigsolve(f, x₀, howmany, which, algorithm)

Compute `howmany` eigenvalues from the linear map encoded in the matrix `A` or by the
function `f`. Return eigenvalues, eigenvectors and a `ConvergenceInfo` structure.

### Arguments:
The linear map can be an `AbstractMatrix` (dense or sparse) or a general function or callable
object. If an `AbstractMatrix` is used, a starting vector `x₀` does not need to be provided,
it is then chosen as `rand(T, size(A,1))`. If the linear map is encoded more generally as a
a callable function or method, the best approach is to provide an explicit starting guess `x₀`.
Note that `x₀` does not need to be of type `AbstractVector`, any type that behaves as a vector
and supports the required methods (see KrylovKit docs) is accepted. If instead of `x₀` an integer
`n` is specified, it is assumed that `x₀` is a regular vector and it is initialized to `rand(T,n)`,
where the default value of `T` is `Float64`, unless specified differently.

The next arguments are optional, but should typically be specified. `howmany` specifies how many
eigenvalues should be computed; `which` specifies which eigenvalues should be targetted. Valid
specifications of `which` are given by
*   `LM`: eigenvalues of largest magnitude
*   `LR`: eigenvalues with largest (most positive) real part
*   `SR`: eigenvalues with smallest (most negative) real part
*   `LI`: eigenvalues with largest (most positive) imaginary part, only if `T <: Complex`
*   `SI`: eigenvalues with smallest (most negative) imaginary part, only if `T <: Complex`
*   [`ClosestTo(λ)`](@ref): eigenvalues closest to some number `λ`
!!! note "Note about selecting `which` eigenvalues"
    Krylov methods work well for extremal eigenvalues, i.e. eigenvalues on the periphery of
    the spectrum of the linear map. Even with `ClosestTo`, no shift and invert is performed.
    This is useful if, e.g., you know the spectrum to be within the unit circle in the complex
    plane, and want to target the eigenvalues closest to the value `λ = 1`.

The argument `T` acts as a hint in which `Number` type the computation should be performed, but
is not restrictive. If the linear map automatically produces complex values, complex arithmetic
will be used even though `T<:Real` was specified.

### Return values:
The return value is always of the form `vals, vecs, info = eigsolve(...)` with
*   `vals`: a `Vector` containing the eigenvalues, of length at least `howmany`, but could be
    longer if more eigenvalues were converged at the same cost. Eigenvalues will be real if
    [`Lanczos`](@ref) was used and complex if [`Arnoldi`](@ref) was used (see below).
*   `vecs`: a `Vector` of corresponding eigenvectors, of the same length as `vals`. Note that
    eigenvectors are not returned as a matrix, as the linear map could act on any custom Julia
    type with vector like behavior, i.e. the elements of the list `vecs` are objects that are
    typically similar to the starting guess `x₀`, up to a possibly different `eltype`. In particular,
    for a general matrix (i.e. with `Arnoldi`) the eigenvectors are generally complex and are
    therefore always returned in a complex number format.
    When the linear map is a simple `AbstractMatrix`, `vecs` will be `Vector{Vector{<:Number}}`.
*   `info`: an object of type [`ConvergenceInfo`], which has the following fields
    -   `info.converged::Int`: indicates how many eigenvalues and eigenvectors were actually
        converged to the specified tolerance `tol` (see below under keyword arguments)
    -   `info.residual::Vector`: a list of the same length as `vals` containing the residuals
        `info.residual[i] = f(vecs[i]) - vals[i] * vecs[i]`
    -   `info.normres::Vector{<:Real}`: list of the same length as `vals` containing the norm
        of the residual `info.normres[i] = norm(info.residual[i])`
    -   `info.numops::Int`: number of times the linear map was applied, i.e. number of times
        `f` was called, or a vector was multiplied with `A`
    -   `info.numiter::Int`: number of times the Krylov subspace was restarted (see below)
!!! warning "Check for convergence"
    No warning is printed if not all requested eigenvalues were converged, so always check
    if `info.converged >= howmany`.

### Keyword arguments:
Keyword arguments and their default values are given by:
*   `tol::Real`: the requested accuracy (corresponding to the 2-norm of the residual for
    Schur vectors, not the eigenvectors). If you work in e.g. single precision (`Float32`),
    you should definitely change the default value.
*   `krylovdim::Integer`: the maximum dimension of the Krylov subspace that will be constructed.
    Note that the dimension of the vector space is not known or checked, e.g. `x₀` should not
    necessarily support the `Base.length` function. If you know the actual problem dimension
    is smaller than the default value, it is useful to reduce the value of `krylovdim`, though
    in principle this should be detected.
*   `maxiter::Integer`: the number of times the Krylov subspace can be rebuilt; see below for
    further details on the algorithms.
*   `orth::Orthogonalizer`: the orthogonalization method to be used, see [`Orthogonalizer`](@ref)
*   `issymmetric::Bool`: if the linear map is symmetric, only meaningful if `T<:Real`
*   `ishermitian::Bool`: if the linear map is hermitian
The default values are given by `tol = KrylovDefaults.tol`, `krylovdim = KrylovDefaults.krylovdim`,
`maxiter = KrylovDefaults.maxiter`, `orth = KrylovDefaults.orth`; see [`KrylovDefaults`](@ref) for details.

The default value for the last two parameters depends on the method. If an `AbstractMatrix`
is used, `issymmetric` and `ishermitian` are checked for that matrix, ortherwise the default
values are `issymmetric = false` and `ishermitian = T <: Real && issymmetric`.

### Algorithm
The last method, without default values and keyword arguments, is the one that is finally called,
and can also be used directly. Here, one specifies the algorithm explicitly as either [`Lanczos`](@ref),
for real symmetric or complex hermitian problems, or [`Arnoldi`](@ref), for general problems.
Note that these names refer to the process for building the Krylov subspace, but the actual
algorithm is an implementation of the Krylov-Schur algorithm, which can dynamically shrink and
grow the Krylov subspace, i.e. the restarts are so-called thick restarts where a part of the
current Krylov subspace is kept.

!!! note "Note about convergence"
    In case of a general problem, where the `Arnoldi` method is used, convergence of an eigenvalue
    is not based on the norm of the residual `norm(f(vecs[i]) - vals[i]*vecs[i])` for the eigenvectors
    but rather on the norm of the residual for the corresponding Schur vectors.

    See also [`schursolve`](@ref) if you want to use the partial Schur decomposition directly,
    or if you are not interested in computing the eigenvectors, and want to work in real arithmetic
    all the way true (if the linear map and starting guess are real).
"""
function eigsolve end

"""
    ClosestTo(λ)

A simple `struct` to be used in combination with [`eigsolve`](@ref) or [`schursolve`](@ref) to
indicate which eigenvalues need to be targetted, namely those closest to `λ`.
"""
struct ClosestTo{T}
    λ::T
end

const Selector = Union{ClosestTo, Symbol}

function eigsolve(A::AbstractMatrix, howmany::Int = 1, which::Selector = :LM, T::Type = eltype(A); kwargs...)
    eigsolve(A, rand(T, size(A,1)), howmany, which; kwargs...)
end
function eigsolve(f, n::Int, howmany::Int = 1, which::Selector = :LM, T::Type = Float64; kwargs...)
    eigsolve(f, rand(T, n), howmany, which; kwargs...)
end
function eigsolve(f, x₀, howmany::Int = 1, which::Selector = :LM; kwargs...)
    alg = eigselector(f, eltype(x₀); kwargs...)
    checkwhich(which) || error("Unknown eigenvalue selector: which = $which")
    if alg isa Lanczos
        if which == :LI || which == :SI || which isa ClosestTo{<:Complex}
            error("Eigenvalue selector which = $which invalid: real eigenvalues expected with Lanczos algorithm")
        end
    elseif eltype(x₀) <: Real
        if which == :LI || which == :SI || which isa ClosestTo{<:Complex}
            error("Eigenvalue selector which = $which invalid: work in complex arithmetic by providing a complex starting vector `x₀`")
        end
    end
    eigsolve(f, x₀, howmany, which, alg)
end

function eigselector(f, T::Type;
    issymmetric::Bool = false, ishermitian::Bool = T<:Real && issymmetric,
    krylovdim::Int = KrylovDefaults.krylovdim, maxiter::Int = KrylovDefaults.maxiter,
    tol::Real = KrylovDefaults.tol, orth::Orthogonalizer = KrylovDefaults.orth)
    if (T<:Real && issymmetric) || ishermitian
        return Lanczos(krylovdim = krylovdim, maxiter = maxiter, tol=tol, orth = orth)
    else
        return Arnoldi(krylovdim = krylovdim, maxiter = maxiter, tol=tol, orth = orth)
    end
end
function eigselector(A::AbstractMatrix, T::Type;
    issymmetric::Bool = issymmetric(A), ishermitian::Bool = ishermitian(A),
    krylovdim::Int = KrylovDefaults.krylovdim, maxiter::Int = KrylovDefaults.maxiter,
    tol::Real = KrylovDefaults.tol, orth::Orthogonalizer = KrylovDefaults.orth)
    if (T<:Real && issymmetric) || ishermitian
        return Lanczos(krylovdim = krylovdim, maxiter = maxiter, tol=tol, orth = orth)
    else
        return Arnoldi(krylovdim = krylovdim, maxiter = maxiter, tol=tol, orth = orth)
    end
end

checkwhich(::ClosestTo) = true
checkwhich(s::Symbol) = s in (:LM, :LR, :SR, :LI, :SI)

function eigsort(which::ClosestTo)
    by = x->abs(x-which.λ)
    rev = false
    return by, rev
end
function eigsort(which::Symbol)
    if which == :LM
        by = abs
        rev = true
    elseif which == :LR
        by = real
        rev = true
    elseif which == :SR
        by = real
        rev = false
    elseif which == :LI
        by = imag
        rev = true
    elseif which == :SI
        by = imag
        rev = false
    else
        error("invalid specification of which eigenvalues to target: which = $which")
    end
    return by, rev
end


text:
"""
    hostbyteorder() -> ByteOrder

Get the byte order for this machine.
"""
hostbyteorder() = ifelse(Base.ENDIAN_BOM == 0x04030201, LittleEndian, BigEndian)

"""
    machstfrombyteorder(order = hostbyteorder()) -> Base.CodeUnits{UInt8,String}

Get the default machine stamp for the passed byte order.
"""
function machstfrombyteorder(order::ByteOrder=hostbyteorder())
    return ifelse(order == LittleEndian, MACHINE_STAMP_LITTLE, MACHINE_STAMP_BIG)
end

"""
    byteorderfrommachst(machst) -> ByteOrder

Get the byte order from machine stamp.
"""
function byteorderfrommachst(machst)
    return if machst[1] == MACHINE_STAMP_LITTLE[1] && (
        machst[2] == MACHINE_STAMP_LITTLE[2] || machst[2] == MACHINE_STAMP_LITTLE_ALT[2]
    )
        LittleEndian
    elseif machst[1] == MACHINE_STAMP_BIG[1] && machst[2] == MACHINE_STAMP_BIG[2]
        BigEndian
    else
        @warn "Unrecognized machine stamp $(machst). Assuming little endian." maxlog = 1
        LittleEndian
    end
end

function bswaptoh(order::ByteOrder)
    return ifelse(order == LittleEndian, maybeswap(ltoh), maybeswap(ntoh))
end
bswaptoh(machst) = bswaptoh(byteorderfrommachst(machst))

function bswapfromh(order::ByteOrder)
    return ifelse(order == LittleEndian, maybeswap(htol), maybeswap(hton))
end
bswapfromh(machst) = bswapfromh(byteorderfrommachst(machst))

maybeswap(fswap, x) = x
maybeswap(fswap, x::Number) = fswap(x)
maybeswap(fswap, x::NTuple{N,Number}) where {N} = ntuple(i -> fswap(x[i]), Val{N}())
maybeswap(fswap) = x -> maybeswap(fswap, x)

"""
    padtruncto!(x::AbstractVector, n; value = zero(eltype(x)))

Pad `x` with `value` or truncate it until its length is exactly `n`.
"""
function padtruncto!(x, n; value=zero(eltype(x)))
    l = length(x)
    if l < n
        append!(x, fill(eltype(x)(value), n - l))
    elseif l > n
        resize!(x, n)
    end
    return x
end

function checkmagic(io)
    magic = read(io, 6)
    seek(io, 0)
    return get(COMPRESSOR_MAGICS, magic[1:3], get(COMPRESSOR_MAGICS, magic, :none))
end

checkextension(path) = get(COMPRESSOR_EXTENSIONS, splitext(path)[end], :none)

function compressstream(io, type)
    return get(COMPRESSIONS, type, (compressor=NoopStream,)).compressor(io)
end

function decompressstream(io, type)
    return get(COMPRESSIONS, type, (decompressor=NoopStream,)).decompressor(io)
end

function _unsafe_string(pointer, size)
    return rstrip(unsafe_string(convert(Ptr{UInt8}, pointer), size), '\0')
end


text:
# Basic tests that check that no plot crashes
function t_plotting()

	X = rand(2,10)			# random 2-D data
	y = round.(rand(10)) 		# random labels

	Du = j4pr.datacell(X)		# unlabeled datacell
	Dl = j4pr.datacell(X,y)		# labeled datacell

	pl = j4pr.lineplot(1) 		# plots first variable when piped data
	pl2 = j4pr.lineplot(sin,1) 	# plots the sine of first variable when piped data
	psc = j4pr.scatterplot() 	# plots scatterplot
	pd1 = j4pr.densityplot1d() 	# plots densityplot 1d
	pd2 = j4pr.densityplot2d() 	# plots densityplot 2d

	for p in [pl, pl2, psc, pd1, pd2]
		Test.@test try
			Du |> p;
			Dl |> p;
			true
		catch
			false
		end
	end

	if v"0.6" <= VERSION < v"0.7-"
	Test.@test try
		j4pr.rocplot((X./sum(X,1),y)|>j4pr.findop(1,j4pr.ROC.TPr(),1.));
		true
	catch
		false
	end
	end

end


text:
using ImageFeatures

img = fill(false,5,5)
img[3,:] .= true

println(hough_transform_standard(img))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using SynchrotronKernel

"""
    Constants for synchrotron
"""

const c_light    = 2.9979e10
const m_e        = 9.10953e-28
const q_e        = 1.602176487e-20 * c_light
const C_crit     = 3q_e / ( 4π * m_e * c_light ) # Donnert+16, MNRAS 462, 2014–2032 (2016), Eg. 20
                                                        #  -> converted to dimensionless momentum
const γ          = 5.0/3.0
const mJy_factor = 1.e26       # conversion factor from [erg/cm^3/Hz/s] to mJy/cm.

"""
    ν_over_ν_crit(p::Real, B::Real, ν::Real)

Computes the fraction ``x = \frac{ν}{ν_c}`` needed for the synchrotron Kernel.
See Donnert+16, MNRAS 462, 2014–2032 (2016), Eq. 19.
"""
ν_over_ν_crit(p::Real, B::Real, ν::Real, sinθ::Real = 1.0) = ν / (C_crit * B * sinθ * p^2)

"""
    integrate_θ(x_in::Real, θ_steps::Integer=100)

Pitch angle integration in Donnert+16, Eq. 17.
"""
function integrate_θ(x_in::Real, θ_steps::Integer = 64)

    dθ = 0.5π / θ_steps
    K = 0.0
    @inbounds for θ ∈ LinRange(0.0, 0.5π, θ_steps)
        sinθ = sin(θ)
        x = x_in / sinθ
        K += sinθ^2 * synchrotron_kernel(x) * dθ
    end
    return K
end

"""
    integrate_θ(x_in::Real, θ_steps::Integer=100)

Pitch angle integration in Donnert+16, Eq. 17 using Simpson's rule.
"""
function integrate_θ_simpson(x_in::Real, θ_steps::Integer = 50)

    dθ = 0.5π / θ_steps
    #θ = LinRange(0.0, 0.5π, θ_steps)

    F = Vector{Float64}(undef, θ_steps)
    F_mid = Vector{Float64}(undef, θ_steps)

    # sin(0.0) = 0.0
    F[1] = 0.0
    F_mid[1] = sin(0.5dθ)^2 * synchrotron_kernel(x_in / sin(0.5dθ))
    K = 0.0
    @inbounds for i = 2:θ_steps
        # boundary
        sinθ = sin((i - 1) * dθ)
        x = x_in / sinθ
        F[i] = sinθ^2 * synchrotron_kernel(x)

        # mid point
        sinθ = sin((i - 0.5) * dθ)
        x = x_in / sinθ
        F_mid[i] = sinθ^2 * synchrotron_kernel(x)

        # Simpson rule: https://en.wikipedia.org/wiki/Simpson%27s_rule
        K += dθ / 6.0 * (F[i] + F[i-1] + 4F_mid[i])
    end
    return K
end

"""
    synchrotron_emission( f_p::Vector{<:Real},
                                   q::Vector{<:Real},
                                   cut::Real,
                                   B_cgs::Real;
                                   ν0::Real=1.44e9
                                   integrate_pitch_angle::Bool=false,
                                   convert_to_mJy::Bool=false,
                                   reduce_spectrum::Bool = true)

Computes the synchrotron emission (in ``[erg/cm^3/Hz/s]``) for a CR spectrum as described in Donnert+16, MNRAS 462, 2014–2032 (2016), Eq. 17.

# Arguments
- `f_p::Vector{<:Real}`: Spectral Norm of bins in cgs units.
- `q::Vector{<:Real}`:   Slopes of bins.
- `cut::Real`:           Spectral cutoff.
- `B_cgs::Real`:         Magnetic field strength (absolute value) in Gauss.
- `par::CRMomentumDistributionConfig`: Parameters of spectrum

# Keyword Arguments
- `ν0::Real=1.44e9`:                  Observation frequency in ``Hz``.
- `integrate_pitch_angle::Bool=false`: Explicitly integrates over the pitch angle. If `false` assumes ``sin(θ) = 1``.
- `convert_to_mJy::Bool=false`:       Convert the result from ``[erg/cm^3/Hz/s]`` to ``mJy/cm``.
- `reduce_spectrum::Bool=true`:       Return result as spectrum if set to false, otherwise calculate total emission.
"""
function synchrotron_emission( f_p::Vector{<:Real},
                               q::Vector{<:Real},
                               cut::Real,
                               B_cgs::Real,
                               par::CRMomentumDistributionConfig;
                               ν0::Real = 1.44e9,
                               integrate_pitch_angle::Bool = false,
                               convert_to_mJy::Bool = false,
                               reduce_spectrum::Bool = true )

    # prefactor to Eq. 17
    j_ν_prefac = q_e * √(3) / (m_e * c_light^2)

    # include magnetic field into this
    j_ν_prefac *= B_cgs

    # if run without pitch angle integration
    # sinθ = 1.0 -> integral factor π/2
    if !integrate_pitch_angle
        j_ν_prefac *= 0.5π
    end

    # storage array for synchrotron emissivity
    jν = Vector{Float64}(undef, par.Nbins)

    # loop over spectrum
    @inbounds for i = 1:par.Nbins

        # beginning of bin
        p_start = par.pmin * 10.0^((i - 1) * par.bin_width)

        # check if bin is below cutoff
        if p_start > cut
            # set remaining bins to zero
            jν[i:end] = 0
            break
        end

        # x from Eq. 19
        x = ν_over_ν_crit(p_start, B_cgs, ν0)
        # pitch-angle integral
        if integrate_pitch_angle
            K = integrate_θ_simpson(x)
        else
            K = synchrotron_kernel(x)
        end

        # energy density at momentum p * integrated synchrotron kernel
        F_start = 4π * p_start^2 * f_p[i] * K

        # middle of bin
        p_mid = par.pmin * 10.0^((i - 1 / 2) * par.bin_width)
        f_p_mid = f_p[i] * (p_mid / p_start)^(-q[i])

        # x from Eq. 19
        x = ν_over_ν_crit(p_mid, B_cgs, ν0)

        if integrate_pitch_angle
            K = integrate_θ_simpson(x)
        else
            K = synchrotron_kernel(x)
        end

        F_mid = 4π * p_mid^2 * f_p_mid * K

        # end of bin
        p_end = par.pmin * 10.0^(i * par.bin_width)
        if p_end > cut
            p_end = cut
        end

        # construct norm at end of spectrum
        f_p_end = f_p[i] * (p_end / p_start)^(-q[i])

        # x from Eq. 19
        x = ν_over_ν_crit(p_end, B_cgs, ν0)

        if integrate_pitch_angle
            K = integrate_θ_simpson(x)
        else
            K = synchrotron_kernel(x)
        end

        F_end = 4π * p_end^2 * f_p_end * K

        # bin width
        dp = p_end - p_start

        # store total synchrotron emissivity
        # Simpson rule: https://en.wikipedia.org/wiki/Simpson%27s_rule
        jν[i] = dp / 6.0 * (F_start + F_end + 4F_mid)
    end

    if convert_to_mJy
        j_ν_prefac *= mJy_factor
    end

    if reduce_spectrum
        # return total emissivity
        return j_ν_prefac * sum(jν)
    else
        # return emissivity per bin
        return j_ν_prefac .* jν
    end

end


text:
using Test
import GenericLinearAlgebra
import GenericLinearAlgebra: Ac_mul_A_RFP, TriangularRFP

@testset "Rectuangular Full Pack Format" begin

    @testset "Core generic functionality" for n in (6, 7),
                                           uplo in (:U, :L)

        A = rand(10, n)

        @testset "Hermitian" begin
            AcA_RFP = Ac_mul_A_RFP(A, uplo)

            @test size(AcA_RFP, 1) == n
            @test size(AcA_RFP, 2) == n
            @test size(AcA_RFP, 3) == 1
            @test_throws BoundsError AcA_RFP[0, 1]
            @test_throws BoundsError AcA_RFP[1, 0]
            @test_throws BoundsError AcA_RFP[n + 1, 1]
            @test_throws BoundsError AcA_RFP[1, n + 1]

            @test AcA_RFP[2, 1]             == AcA_RFP[1, 2]
            @test AcA_RFP[end - 2, end - 1] == AcA_RFP[end - 1, end - 2]
        end

        @testset "Triangular" begin
            Atr_RFP = TriangularRFP(triu(A'A), uplo)
            @test size(Atr_RFP, 1) == n
            @test size(Atr_RFP, 2) == n
            @test size(Atr_RFP, 3) == 1

            # Indexing not implemented yet for Atr_RFP
            # @test_throws BoundsError Atr_RFP[0, 1]
            # @test_throws BoundsError Atr_RFP[1, 0]
            # @test_throws BoundsError Atr_RFP[n + 1, 1]
            # @test_throws BoundsError Atr_RFP[1, n + 1]

            @test_broken Atr_RFP[2, 1]             == Atr_RFP[1, 2]
            @test_broken Atr_RFP[end - 2, end - 1] == Atr_RFP[end - 1, end - 2]
        end
    end

    @testset "Hermitian with element type: $elty. Problem size: $n" for
        elty in (Float32, Float64, Complex{Float32}, Complex{Float64}),
            n in (6, 7),
                uplo in (:L, :U)

        A       = rand(elty, 10, n)
        AcA     = A'A
        AcA_RFP = Ac_mul_A_RFP(A, uplo)
        o       = ones(elty, n)

        @test AcA   ≈ A'A
        @test AcA\o ≈ AcA_RFP\o
        @test inv(AcA) ≈ inv(AcA_RFP)
        @test inv(cholesky(AcA)) ≈ inv(factorize(AcA_RFP))
    end

    @testset "Hermitian with element type: $elty. Problem size: $n" for
        elty in (Float32, Float64, Complex{Float32}, Complex{Float64}),
            n in (6, 7),
                uplo in (:L, :U)

        A     = lu(rand(elty, n, n)).U
        A     = uplo == :U ? A : copy(A')
        A_RFP = TriangularRFP(A, uplo)
        o     = ones(elty, n)

        @test_broken A ≈ A_RFP
        @test        A ≈ Array(A_RFP)
        @test      A\o ≈ A_RFP\o
        @test   inv(A) ≈ Array(inv(A_RFP))
    end
end


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

x = [1, 2, 3, 4, 5]
y = [1.0, 2.0, NaN, 4.0, 5.0]
plot(x=x, y=y, ymin=y.-0.5, ymax=y.+0.5, Geom.line, Geom.point, Geom.ribbon)


text:
function read_map(path)
    lines = open(path, "r") do file
        readlines(file)
    end
    height_line = split(lines[2])
    width_line = split(lines[3])
    height = parse(Int, height_line[2])
    width = parse(Int, width_line[2])
    map_matrix = Matrix{Char}(undef, height, width)
    for i in 1:height
        line = lines[4 + i]
        for j in 1:width
            map_matrix[i, j] = line[j]
        end
    end
    return map_matrix
end

function read_scenario(path)
    lines = open(path, "r") do file
        readlines(file)
    end
    scenario = DataFrame(;
        bucket=Int[],
        map=String[],
        width=Int[],
        height=Int[],
        start_i=Int[],
        start_j=Int[],
        goal_i=Int[],
        goal_j=Int[],
        optimal_length=Float64[],
    )
    for line in @view lines[2:end]
        line_split = split(line, "\t")
        bucket = parse(Int, line_split[1]) + 1
        map = line_split[2]
        width = parse(Int, line_split[3])
        height = parse(Int, line_split[4])
        start_x = parse(Int, line_split[5]) + 1
        start_y = parse(Int, line_split[6]) + 1
        goal_x = parse(Int, line_split[7]) + 1
        goal_y = parse(Int, line_split[8]) + 1
        optimal_length = parse(Float64, line_split[9])

        start_i = height - start_y + 1
        start_j = start_x
        goal_i = height - goal_y + 1
        goal_j = goal_x

        line_tup = (
            bucket=bucket,
            map=map,
            width=width,
            height=height,
            start_i=start_i,
            start_j=start_j,
            goal_i=goal_i,
            goal_j=goal_j,
            optimal_length=optimal_length,
        )
        push!(scenario, line_tup)
    end
    return scenario
end


text:
# module ComparisonDataset

using Images, Random, Revise
import Base: size, rand, ==
import Random: rand!

includet("rect.jl")
includet("data_generation.jl")
includet("patch.jl")
includet("sample.jl")

# export Size, Patch, Sample
# export render
# end # module


text:
"""
    convert_to_changes(s, M; drop_identities = false) # todo: change function-name

Given the initial state of the solution `s` and a list of accepted moves `M`, compute for each move in `M` which changes it induced on `s`, in terms of the index set of `s`. Put differently, this function computes what parts of a solution are changed throughout a search process defined by a sequence of moves `M`.

# Usage
```julia
julia> A
2×3 Matrix{Int64}:
 1  4  4
 5  4  1

julia> M = [s -> swap!(s, i, j) for i in 1:6, j in 1:6];

julia> convert_to_changes(A, M, drop_identities = false)
36-element Vector{Any}:
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}()
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}(2 => (5 => 1), 1 => (1 => 5))
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}(3 => (4 => 5), 1 => (5 => 4))
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}()
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}()
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}(6 => (1 => 4), 1 => (4 => 1))
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}()
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}()
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}(2 => (1 => 5), 3 => (5 => 1))
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}(4 => (4 => 5), 2 => (5 => 4))
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}()
 ⋮
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}()
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}(5 => (1 => 4), 4 => (4 => 1))
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}()
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}(5 => (4 => 5), 6 => (5 => 4))
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}()
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}()
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}(6 => (4 => 1), 3 => (1 => 4))
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}()
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}(5 => (5 => 1), 6 => (1 => 5))
 Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}}()
```

"""
function convert_to_changes(s, M; drop_identities = true)
    s_ = copy(s)
    D = Dict{eltype(M), Dict{eltype(eachindex(s)),  Pair{Union{eltype(s),Nothing},Union{eltype(s),Nothing}}}}()
    for m in M
        s_ = m(s_) # mutating call
        d = soldiff(s, s_)
        (!isempty(d) || !drop_identities) && push!(D, m => d)
        s  = m(s)
    end
    return D
end

"""
    soldiff(s1::S, s2::S) where {S}

Compute for two solutions `s1` and `s2` the set of indices whose values changed, and how they changed. Works for non-equal index sets (e.g. Matrices of different dimensions).

# Usage

```julia
julia> A = rand(1:5, 2, 3)
2×3 Matrix{Int64}:
 1  1  1
 2  5  2

julia> B = rand(1:5, 2, 2)
2×2 Matrix{Int64}:
 2  1
 3  5

julia> d = soldiff(A,B)
Dict{Int64, Pair{Union{Nothing, Int64}, Union{Nothing, Int64}}} with 4 entries:
  5 => 1=>nothing
  6 => 2=>nothing
  2 => 2=>3
  1 => 1=>2
```
"""
function soldiff(s1::S, s2::S) where {S}
    #I1 = collect(eachindex(s1))
    #I2 = collect(eachindex(s2))
    I1 = eachindex(s1)
    I2 = eachindex(s2)
    I = I1 ∩ I2  # smaller set
    I1_ = Iterators.filter(i -> i ∉ I ,I1)
    I2_ = Iterators.filter(i -> i ∉ I, I2)
    idx_to_Δ = Dict{eltype(I), Pair{Union{eltype(s1),Nothing},Union{eltype(s2),Nothing}}}()
    for i in I
        if s1[i] ≠ s2[i]
            idx_to_Δ[i] = s1[i] => s2[i]
        end
    end
    for i in I1_
        idx_to_Δ[i] = s1[i] => nothing
    end

    for i in I2_
        idx_to_Δ[i] = nothing => s2[i]
    end
    return idx_to_Δ
end

max_stars_repo_name: Michiel-VL/PolyaViz

max_issues_repo_name: Michiel-VL/PolyaViz


text:
module Postprocessing

using LinearAlgebra
using FFTW
using OffsetArrays
using ToeplitzMatrices

"""
    reference(algorithm)

Returns a citable reference for `algorithm`, e.g. a bibtex entry
for a scientific article or a book.
"""
function reference end

include("total_variation.jl")
include("fourier_pade.jl")

export reference
export total_variation, total_variation_denoising, total_variation_denoising!,
       group_sparse_total_variation_denoising, group_sparse_total_variation_denoising!
export fourier_pade

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AzureFilesIdentityBasedAuthentication <: SwaggerModel
    directoryServiceOptions::Any # spec type: Union{ Nothing, String } # spec name: directoryServiceOptions
    activeDirectoryProperties::Any # spec type: Union{ Nothing, ActiveDirectoryProperties } # spec name: activeDirectoryProperties

    function AzureFilesIdentityBasedAuthentication(;directoryServiceOptions=nothing, activeDirectoryProperties=nothing)
        o = new()
        validate_property(AzureFilesIdentityBasedAuthentication, Symbol("directoryServiceOptions"), directoryServiceOptions)
        setfield!(o, Symbol("directoryServiceOptions"), directoryServiceOptions)
        validate_property(AzureFilesIdentityBasedAuthentication, Symbol("activeDirectoryProperties"), activeDirectoryProperties)
        setfield!(o, Symbol("activeDirectoryProperties"), activeDirectoryProperties)
        o
    end
end # type AzureFilesIdentityBasedAuthentication

const _property_map_AzureFilesIdentityBasedAuthentication = Dict{Symbol,Symbol}(Symbol("directoryServiceOptions")=>Symbol("directoryServiceOptions"), Symbol("activeDirectoryProperties")=>Symbol("activeDirectoryProperties"))
const _property_types_AzureFilesIdentityBasedAuthentication = Dict{Symbol,String}(Symbol("directoryServiceOptions")=>"String", Symbol("activeDirectoryProperties")=>"ActiveDirectoryProperties")
Base.propertynames(::Type{ AzureFilesIdentityBasedAuthentication }) = collect(keys(_property_map_AzureFilesIdentityBasedAuthentication))
Swagger.property_type(::Type{ AzureFilesIdentityBasedAuthentication }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AzureFilesIdentityBasedAuthentication[name]))}
Swagger.field_name(::Type{ AzureFilesIdentityBasedAuthentication }, property_name::Symbol) =  _property_map_AzureFilesIdentityBasedAuthentication[property_name]

const _allowed_AzureFilesIdentityBasedAuthentication_directoryServiceOptions = ["None", "AADDS", "AD"]

function check_required(o::AzureFilesIdentityBasedAuthentication)
    true
end

function validate_property(::Type{ AzureFilesIdentityBasedAuthentication }, name::Symbol, val)
    if name === Symbol("directoryServiceOptions")
        Swagger.validate_param(name, "AzureFilesIdentityBasedAuthentication", :enum, val, _allowed_AzureFilesIdentityBasedAuthentication_directoryServiceOptions)
    end
end


text:
######################################################
#
# Conservative Binomial Cascade Function
#
# Date:   2014-09-28
# Author: Hessel Tuinhof
#
# Function to construct a conservative binomial cascade.
# The conservative measure preserves mass exactly.
#
# Reference:
# - Riedi, Rudolf H.: Multifractal Processes
# (http://www.stat.rice.edu/~riedi/Publ/PDF/MP.pdf)
# - Calvet, L. E.; Fisher, A. J.: Multifractal Volatility
#
# Required packages: Distributions
#
# Input:
#   m0 ... submass, float in (0,1)
#   k ... number of cascades, positive integer
#   startmass ... total mass, default: 1
#
# Output:
#   mu ... vector contains all measures
#
######################################################

using Distributions

function BinomialConservative(m0, k, startmass = 1)
	rand(Bernoulli()) == 1 ? M = m0 : M = 1-m0
	mu = [M*startmass, (1-M)*startmass]
	if k > 1
		tmp = Float64[]
		mu = append!(tmp, [BinomialConservative(m0, k-1, mu[1]), BinomialConservative(m0, k-1, mu[2])])
	end
	return mu
end

max_stars_repo_name: hesseltuinhof/masterthesis

max_issues_repo_name: hesseltuinhof/masterthesis

max_forks_repo_name: hesseltuinhof/masterthesis


text:
function parse_ev_technologies!(dtr::DieterModel, df::DataFrame)
# function parse_ev_technologies!(dtr::DieterModel, path::AbstractString)
    # df = CSV.read(path)
    dtr.sets[:ElectricVehicles] = disallowmissing(unique(df[!,:EV]))

    params = map_idcol(df)
    for (k,v) in params update_dict!(dtr.parameters, k, v) end

    return nothing
end

function parse_ev_demand!(dtr::DieterModel, df::DataFrame)
# function parse_ev_demand!(dtr::DieterModel, path::AbstractString)
    # params = CSV.read(path) |> map_dfheader_to_col
    params = map_dfheader_to_col(df)
    update_dict!(dtr.parameters, :EvDemand, params)

    return nothing
end

function parse_ev_power!(dtr::DieterModel, df::DataFrame)
# function parse_ev_power!(dtr::DieterModel, path::AbstractString)
    # params = CSV.read(path) |> map_dfheader_to_col
    params = map_dfheader_to_col(df)
    update_dict!(dtr.parameters, :EvPower, params)

    return nothing
end

function calc_ev_quantity!(dtr::DieterModel)
    p = dtr.parameters
    EV = dtr.sets[:ElectricVehicles]
    s = p[:EvShare]
    q = dtr.settings[:ev] * 1e6

    p[:AbsoluteEvDemand] = Dict(ev => s[ev] * q .* p[:EvDemand][ev] for ev in EV)
    p[:AbsoluteEvPower] = Dict(ev => s[ev] * q .* p[:EvPower][ev] for ev in EV)
    p[:AbsoluteEvCapacity] = Dict(ev => s[ev] * q * p[:EvCapacity][ev] for ev in EV)

    return nothing
end


text:
module Cbc

# This 'using' is required to suppress a warning about Clp not having Libdl in
# its dependencies (Libdl is used by BinaryProvider), e.g.:
# bicycle1885/CodecZlib.jl#26.
using Libdl

if isfile(joinpath(@__DIR__, "..", "deps", "deps.jl"))
    include(joinpath(@__DIR__, "..", "deps", "deps.jl"))
else
    error("Cbc not properly installed. Please run Pkg.build(\"Cbc\")")
end

include("CbcCInterface.jl")
include("MPB_wrapper.jl")
include("MOI_wrapper.jl")

using Cbc.CbcMathProgSolverInterface
export CbcSolver

end


text:
# ---
# title: 1488. Avoid Flood in The City
# id: problem1488
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Hash Table
# link: <https://leetcode.com/problems/avoid-flood-in-the-city/description/>
# hidden: true
# ---
#
# Your country has an infinite number of lakes. Initially, all the lakes are
# empty, but when it rains over the `nth` lake, the `nth` lake becomes full of
# water. If it rains over a lake which is **full of water** , there will be a
# **flood**. Your goal is to avoid the flood in any lake.
#
# Given an integer array `rains` where:
#
#   * `rains[i] > 0` means there will be rains over the `rains[i]` lake.
#   * `rains[i] == 0` means there are no rains this day and you can choose **one lake** this day and **dry it**.
#
# Return _an array`ans`_ where:
#
#   * `ans.length == rains.length`
#   * `ans[i] == -1` if `rains[i] > 0`.
#   * `ans[i]` is the lake you choose to dry in the `ith` day if `rains[i] == 0`.
#
# If there are multiple valid answers return **any** of them. If it is
# impossible to avoid flood return **an empty array**.
#
# Notice that if you chose to dry a full lake, it becomes empty, but if you
# chose to dry an empty lake, nothing changes. (see example 4)
#
#
#
# **Example 1:**
#
#
#
#     Input: rains = [1,2,3,4]
#     Output: [-1,-1,-1,-1]
#     Explanation: After the first day full lakes are [1]
#     After the second day full lakes are [1,2]
#     After the third day full lakes are [1,2,3]
#     After the fourth day full lakes are [1,2,3,4]
#     There's no day to dry any lake and there is no flood in any lake.
#
#
# **Example 2:**
#
#
#
#     Input: rains = [1,2,0,0,2,1]
#     Output: [-1,-1,2,1,-1,-1]
#     Explanation: After the first day full lakes are [1]
#     After the second day full lakes are [1,2]
#     After the third day, we dry lake 2. Full lakes are [1]
#     After the fourth day, we dry lake 1. There is no full lakes.
#     After the fifth day, full lakes are [2].
#     After the sixth day, full lakes are [1,2].
#     It is easy that this scenario is flood-free. [-1,-1,1,2,-1,-1] is another acceptable scenario.
#
#
# **Example 3:**
#
#
#
#     Input: rains = [1,2,0,1,2]
#     Output: []
#     Explanation: After the second day, full lakes are  [1,2]. We have to dry one lake in the third day.
#     After that, it will rain over lakes [1,2]. It's easy to prove that no matter which lake you choose to dry in the 3rd day, the other one will flood.
#
#
# **Example 4:**
#
#
#
#     Input: rains = [69,0,0,0,69]
#     Output: [-1,69,1,1,-1]
#     Explanation: Any solution on one of the forms [-1,69,x,y,-1], [-1,x,69,y,-1] or [-1,x,y,69,-1] is acceptable where 1 <= x,y <= 10^9
#
#
# **Example 5:**
#
#
#
#     Input: rains = [10,20,20]
#     Output: []
#     Explanation: It will rain over lake 20 two consecutive days. There is no chance to dry any lake.
#
#
#
#
# **Constraints:**
#
#   * `1 <= rains.length <= 10^5`
#   * `0 <= rains[i] <= 10^9`
#
#
## @lc code=start
using LeetCode

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


text:
using BenchmarkTools

using Setfield: @lens

using Bifurcations
using Bifurcations: Codim1, Codim2, special_intervals
using Bifurcations.Continuations: as, pre_solve!, sweep!
using Bifurcations.Codim2LimitCycle: FoldLimitCycleProblem
using Bifurcations.Examples: MorrisLecar

function make_morris_lecar_hopf_solver()
    solver1 = init(MorrisLecar.make_prob())
    solve!(solver1)
    hopf_point, = special_intervals(solver1, Codim1.PointTypes.hopf)

    hopf_prob = BifurcationProblem(
        hopf_point,
        solver1,
        (@lens _.z),
        (-1.0, 1.0),
    )
    return init(
        hopf_prob;
        nominal_angle_rad = 0.01,
    )
end

function make_morris_lecar_flc_solver()
    hopf_solver = make_morris_lecar_hopf_solver()
    solve!(hopf_solver)
    bautin_point, = special_intervals(hopf_solver, Codim2.PointTypes.bautin)

    flc_prob = FoldLimitCycleProblem(
        bautin_point,
        hopf_solver;
        num_mesh = 30,
        degree = 4,
    )
    return init(
        flc_prob;
        start_from_nearest_root = true,
        bidirectional_first_sweep = false,
        max_branches = 0,
        max_samples = 3,
        nominal_angle_rad = 2π * (5 / 360),
    )
end

suite = BenchmarkGroup()
suite["sweep!(hopf_solver)"] = @benchmarkable(
    sweep!(hopf_solver),
    setup = (hopf_solver = pre_solve!(make_morris_lecar_hopf_solver())))
suite["sweep!(flc_solver)"] = @benchmarkable(
    sweep!(flc_solver),
    setup = (flc_solver = pre_solve!(make_morris_lecar_flc_solver())))

if !isdefined(:SUITE)
    SUITE = BenchmarkGroup()
end
SUITE["morris_lecar"] = suite


text:
TEST_CASES = [
    # projectors, expected_fused, expected_transitions
    (([1 1; 0 1], [1 0; 1 1]), [1 0; 0 1], [[1 1; 0 1], [1 0; 1 1]]),
    (([1 1; 0 0], [1 0; 0 1]), [1 0; 0 1], [[1 1; 0 0], [1 0; 0 1]]),
    (
        ([1 1 1; 0 0 0; 1 1 1], [1 0 0; 0 1 0; 0 0 1]),
        [1 0 0; 0 1 0; 0 0 1],
        [[1 1 1; 0 0 0; 1 1 1], [1 0 0; 0 1 0; 0 0 1]]
    ),
    (
        ([1 1 1; 0 0 0; 1 1 0], [1 0 1; 0 1 0; 0 0 1]),
        [1 0 0; 0 1 0; 0 0 1],
        [[1 1 1; 0 0 0; 1 1 0], [1 0 1; 0 1 0; 0 0 1]]
    )
]

@testset "fuse_projectors correctly fuses projectors." for (projectors, expected_fused, expected_transactions) ∈ TEST_CASES
    fused, transitions = fuse_projectors(projectors)
    @test fused == expected_fused
    @test transitions == expected_transactions
end

@testset "normalize_probability properly deals with negative values" begin
    values = [1.0, -1.0, 15.324, -2.2134, 0.123, 0.0]
    new_values = [2.2134, 2.2134, 15.324, 2.2134, 2.2134, 2.2134]
    #  sum = 26.391
    expected_prob = [
        0.08386950096623849, 0.08386950096623849, 0.5806524951688076,
        0.08386950096623849, 0.08386950096623849, 0.08386950096623849
    ]

    @test normalize_probability(values) ≈ expected_prob ≈ new_values / sum(new_values)
end


text:
"""
    example_task(cur_input_string="")

Lorem ipsum dolor sit amet.
"""
function example_task(cur_input_string="")
  cur_output_string = "404"

  if cur_input_string != ""
    cur_output_string *= " - $(cur_input_string)"
  end

  println(cur_output_string)
end


text:
module Private
export get_most_union_all
function get_most_union_all(expr, mod :: Module)

    if isa(expr, Expr) && expr.head == :curly
        get_most_union_all(expr.args[1], mod)
    else
        @eval mod $expr
    end
end

end


text:
module TMP

using WriteVTK
using WriteVTK: VTKFile, VTKDataType, AbstractVTKDataset, XMLDocument
using WriteVTK: XMLElement, DatasetFile, AbstractFieldData, VTKUnstructuredGrid
using WriteVTK: create_root, xml_name, set_attribute, IS_LITTLE_ENDIAN
using WriteVTK: new_child, file_extension, add_extension, root
using WriteVTK: attribute, find_element, child_elements, datatype_str
using WriteVTK: save_file
import WriteVTK: vtk_save

struct ParallelDatasetFile <: VTKFile
    pvtkargs::Dict
    xdoc::XMLDocument
    dataset::DatasetFile
    path::AbstractString
    function ParallelDatasetFile(
        pvtkargs::Dict,
        xdoc::XMLDocument,
        dataset::DatasetFile,
        path::AbstractString)
        new(pvtkargs,xdoc,dataset,path)
    end
end

function _default_pvtkargs(pvtkargs)
    mandatory = [:part,:nparts]
    optional = [:ismain,:ghost_level]
    allkeys = vcat(mandatory,optional)

    d = Dict{Symbol,Any}(pvtkargs)
    for k in keys(d)
        if !(k in allkeys)
            throw(ArgumentError("$k is not a valid key in pvtkargs."))
        end
    end
    for k in mandatory
        if !(haskey(d,k))
            throw(ArgumentError("$k is a mandatory key in pvtkargs."))
        end
    end
    if ! haskey(d,:ismain)
      d[:ismain] = (d[:part] == 1)
    end
    if ! haskey(d,:ghost_level)
        d[:ghost_level] = 0
    end
    d
end

function vtk_save(vtk::ParallelDatasetFile)
  if isopen(vtk)
    _pvtk_grid_body(vtk)
    if vtk.pvtkargs[:ismain]
       save_file(vtk.xdoc, vtk.path)
    end
    vtk_save(vtk.dataset)
    close(vtk)
  end
  return [vtk.path] :: Vector{String}
end

function new_pdata_array(xParent,
                         type::Type{<:VTKDataType},
                         name::AbstractString,
                         Nc=nothing)
  xDA = new_child(xParent, "PDataArray")
  set_attribute(xDA, "type", datatype_str(type))
  set_attribute(xDA, "Name", name)
  if Nc != nothing
    set_attribute(xDA, "NumberOfComponents", Nc)
  end
end

function get_extension(filename::AbstractString)
  path, ext = splitext(filename)
  ext
end

function get_path(filename::AbstractString)
  path, ext = splitext(filename)
  path
end

function parallel_file_extension(g::AbstractVTKDataset)
  ext=file_extension(g)
  replace(ext,"."=>".p")
end

function parallel_xml_name(g::AbstractVTKDataset)
  "P"*xml_name(g)
end

function parallel_xml_write_header(pvtx::XMLDocument,dtype::AbstractVTKDataset)
  xroot = create_root(pvtx, "VTKFile")
  set_attribute(xroot, "type", parallel_xml_name(dtype))
  set_attribute(xroot, "version", "1.0")
  if IS_LITTLE_ENDIAN
    set_attribute(xroot, "byte_order", "LittleEndian")
  else
    set_attribute(xroot, "byte_order", "BigEndian")
  end
  xroot
end

function add_pieces!(grid_xml_node,
                     prefix::AbstractString,
                     extension::AbstractString,
                     num_pieces::Int)
  for i=1:num_pieces
    piece=new_child(grid_xml_node,"Piece")
    fn = _serial_filename(i,num_pieces,prefix,extension)
    set_attribute(piece,"Source",fn)
  end
  grid_xml_node
end

function add_ppoints!(grid_xml_node;
                      type::Type{<:VTKDataType}=Float64)
  ppoints=new_child(grid_xml_node,"PPoints")
  new_pdata_array(ppoints,type,"Points",3)
end

function add_ppoint_data!(pdataset::ParallelDatasetFile,
                          name::AbstractString;
                          type::Type{<:VTKDataType}=Float64,
                          Nc::Integer=1)
  xroot=root(pdataset.xdoc)
  dtype = xml_name_to_VTK_Dataset(pdataset.dataset.grid_type)
  grid_xml_node=find_element(xroot,parallel_xml_name(dtype))
  @assert grid_xml_node != nothing
  grid_xml_node_ppoint=find_element(grid_xml_node,"PPointData")
  if (grid_xml_node_ppoint==nothing)
    grid_xml_node_ppoint=new_child(grid_xml_node,"PPointData")
  end
  new_pdata_array(grid_xml_node_ppoint,type,name,Nc)
end

function add_pcell_data!(pdataset::ParallelDatasetFile,
                         name::AbstractString;
                         type::Type{<:VTKDataType}=Float64,
                         Nc=nothing)
  xroot=root(pdataset.xdoc)
  dtype = xml_name_to_VTK_Dataset(pdataset.dataset.grid_type)
  grid_xml_node=find_element(xroot,parallel_xml_name(dtype))
  @assert grid_xml_node != nothing
  grid_xml_node_pcell=find_element(grid_xml_node,"PCellData")
  if (grid_xml_node_pcell==nothing)
    grid_xml_node_pcell=new_child(grid_xml_node,"PCellData")
  end
  new_pdata_array(grid_xml_node_pcell,type,name,Nc)
end

function Base.setindex!(pvtk::ParallelDatasetFile,
                        data,
                        name::AbstractString,
                        loc::AbstractFieldData)
  pvtk.dataset[name,loc]=data
end

function Base.setindex!(vtk::ParallelDatasetFile, data, name::AbstractString)
  vtk.dataset[name]=data
end

function get_dataset_extension(dataset::DatasetFile)
  path, ext = splitext(dataset.path)
  @assert ext != ""
  ext
end

function get_dataset_path(dataset::DatasetFile)
  path, ext = splitext(dataset.path)
  @assert ext != ""
  path
end

function xml_name_to_VTK_Dataset(xml_name::AbstractString)
    if (xml_name=="ImageData")
      VTKImageData()
    elseif (xml_name=="RectilinearGrid")
      VTKRectilinearGrid()
    elseif (xml_name=="PolyData")
      VTKPolyData()
    elseif (xml_name=="StructuredGrid")
      VTKStructuredGrid()
    elseif (xml_name=="UnstructuredGrid")
      VTKUnstructuredGrid()
    else
      @assert false
    end
end

function num_children(element)
  length(collect(child_elements(element)))
end

function get_child_attribute(element,child,attr)
  children=collect(child_elements(element))
  attribute(children[child],attr)
end

function get_dataset_xml_type(dataset::DatasetFile)
  r=root(dataset.xdoc)
  attribute(r,"type")
end

function get_dataset_xml_grid_element(dataset::DatasetFile,
                                      element::AbstractString)
   r=root(dataset.xdoc)
   uns=find_element(r,get_dataset_xml_type(dataset))
   piece=find_element(uns,"Piece")
   find_element(piece,element)
end

const string_to_VTKDataType = Dict("Int8"=>Int8,
                                   "UInt8"=>UInt8,
                                   "Int16"=>Int16,
                                   "UInt16"=>UInt16,
                                   "Int32"=>Int32,
                                   "UInt32"=>UInt32,
                                   "Int64"=>Int64,
                                   "UInt64"=>UInt64,
                                   "Float32"=>Float32,
                                   "Float64"=>Float64,
                                   "String"=>String)

"""
    pvtk_grid(args...;pvtkargs,kwargs...)

Return a handler representing a parallel vtk file, which can be
eventually written to file with `vtk_save`.

Positional and keyword arguments in `args` and `kwargs`
are passed to `vtk_grid` verbatim (except file names that are augmented with the
corresponding part id).

The extra keyword argument `pvtkargs` contains
extra options (as a `Dict{Symbol,Any}` or a `Vector{Pair{Symbol,Any}}`)
that only apply for parallel vtk file formats.

Mandatory keys in `pvtkargs` are:

- `:part` current (1-based) part id
- `:nparts` total number of parts

Default key/value pairs in `pvtkargs` are
- `:ismain=>part==1` True if the current part id `part` is the main (the only one that will write the .pvtk file).
- `:ghost_level=>0` Ghost level.
"""
function pvtk_grid(filename::AbstractString,args...;pvtkargs,kwargs...)
    _pvtkargs = _default_pvtkargs(pvtkargs)
    part = _pvtkargs[:part]
    nparts = _pvtkargs[:nparts]
    bname = basename(filename)
    path = mkpath(filename)
    prefix = joinpath(path,bname)
    extension = ""
    fn = _serial_filename(part,nparts,prefix,extension)
    vtk = vtk_grid(fn,args...;kwargs...)
    _pvtk_grid_header(_pvtkargs,vtk,path)
end

function _serial_filename(part,nparts,prefix,extension)
  p = lpad(part,ceil(Int,log10(nparts)),'0')
  fn = prefix*"_$p"*extension
end

function _pvtk_grid_header(
    pvtkargs::Dict,
    dataset::DatasetFile,
    filename::AbstractString)

    pvtx  = XMLDocument()
    dtype = xml_name_to_VTK_Dataset(dataset.grid_type)
    xroot = parallel_xml_write_header(pvtx,dtype)

    # Generate parallel grid node
    grid_xml_node=new_child(xroot,"P"*dataset.grid_type)
    set_attribute(grid_xml_node, "GhostLevel", string(pvtkargs[:ghost_level]))

    prefix=joinpath(filename,basename(filename))
    extension=file_extension(dtype)
    add_pieces!(grid_xml_node,prefix,extension,pvtkargs[:nparts])

    pextension = parallel_file_extension(dtype)
    pfilename = add_extension(filename,pextension)
    pdataset=ParallelDatasetFile(pvtkargs,pvtx,dataset,pfilename)

    # Generate PPoints
    points=get_dataset_xml_grid_element(dataset,"Points")
    type=get_child_attribute(points,1,"type")
    add_ppoints!(grid_xml_node,type=string_to_VTKDataType[type])

    pdataset
  end

function _pvtk_grid_body(pdataset::ParallelDatasetFile)
    dataset=pdataset.dataset
    # Generate PPointData
    pointdata=get_dataset_xml_grid_element(dataset,"PointData")
    if pointdata != nothing
      if (num_children(pointdata)>0)
        for child=1:num_children(pointdata)
          name=get_child_attribute(pointdata,child,"Name")
          type=get_child_attribute(pointdata,child,"type")
          Nc=get_child_attribute(pointdata,child,"NumberOfComponents")
          add_ppoint_data!(pdataset,
                           name;
                           type=string_to_VTKDataType[type],
                           Nc=parse(Int64,Nc))
        end
      end
    end

    # Generate PCellData
    celldata=get_dataset_xml_grid_element(dataset,"CellData")
    if celldata!=nothing
      if (num_children(celldata)>0)
        for child=1:num_children(celldata)
          name=get_child_attribute(celldata,child,"Name")
          type=get_child_attribute(celldata,child,"type")
          Nc=get_child_attribute(celldata,child,"NumberOfComponents")
          add_pcell_data!(pdataset,
                          name;
                          type=string_to_VTKDataType[type],
                          Nc=(Nc==nothing ? nothing : parse(Int64,Nc)))
        end
      end
    end
    pdataset
  end

end #module