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

using Markdown
using InteractiveUtils

# ╔═╡ 08af84c0-2456-11eb-326c-2d675ca197af
using SpectralVis

# ╔═╡ 57a60da0-2461-11eb-07df-2b3ab48a6675
using Colors

# ╔═╡ 9ff8ea50-2461-11eb-3143-d1b38243f297
using Plots; plotly()

# ╔═╡ 14f2d870-245d-11eb-1fe3-779aa49070f1
using SSpline

# ╔═╡ f07cef7e-014a-4d59-8e31-4d394391adda
begin
	#import DarkMode
    #DarkMode.enable(theme="material-palenight", cm_config=Dict("tabSize" => 4))
	#DarkMode.enable()
    #DarkMode.Toolbox(theme="default")
	html"""<style> main {max-width: 900px;}"""
end

# ╔═╡ 78adafde-247e-11eb-12e1-ddea99d50711
md"# A quick tour of the `SpectralVis` and `SSpline` packages"

# ╔═╡ a8bb86de-29be-474d-94e6-1a03f9afe70f
md"## First, let’s get the basic stuff out of the way:"

# ╔═╡ f444ebf0-2460-11eb-2ed3-c5db8e304e95
md"We need to get a few packages to deal with spectra"

# ╔═╡ 1cb25910-2461-11eb-1f85-a9ed9362e326
md"**SpectralVis** contains (almost) everything you need for spectral calculations, spectral data, etc."

# ╔═╡ 5dacb780-2461-11eb-277d-e7226385a5ec
md"The **Colors** package is also needed"

# ╔═╡ 79763630-2461-11eb-0b0d-15b00898c512
md"We also would like to do some plots, so we also use the **Plots** package with the **PyPlot** backend."

# ╔═╡ 0d1268b0-2461-11eb-15a9-87ca552b10c5
md"Then we need the spectral spline package **SSpline**, which is a specialized spline package for spectral calculations. Apart from specialized functionality, it’s also faster than the **Interpolations** package."

# ╔═╡ 1afb7e20-245d-11eb-1ca8-9da92491cd50
md"Now we’re ready to have some fun!"

# ╔═╡ d04f4ee1-2609-4355-9846-5e155557d709
md"## Spectral environments"

# ╔═╡ 7d8d8c1e-ed75-472f-b08b-a93c29ab588c
md"""
The package `SpectralVis` automatically creates a *spectral environment* which defines the necessary parameters to work with spectra. These parameters are stored in the mutable struct `SPECENV` and can be requested and changed to suit the necessary requirements.

The fields contained in `SPECENV` are:

`λmin`: minimum wavelength, ideally identical with λmin of `cmf`, default: 390.0 nm

`Δλ`: wavelength stepwidth, default: 1.0 nm

`λmax`: maximum wavelength, ideally identical with λmax of `cmf`, default: 830.0 nm

`cmf`: Color matching function (CMF) of type `CMatch`, default: `CIE12_10` (CIE 2012 10° observer)

options: `CIE31` (CIE1931 2° observer), `CIE31_J` (CIE1931 2° observer with Judd-correction), `CIE31_JV` (Judd-Vos-correction, `CIE64` (CIE1964 10° observer), `CIE12_2` (CIE2012 2° observer), `CIE12_10` (CIE2012 10° observer)

`ex`: extrapolation type, with the options `:zero`, `:none`, `:boundary`, `:linear`, `:quadratic` (not yet implemented). Default: `:linear`
"""

# ╔═╡ f5760d98-77e3-4847-bb79-97f45b50b7da
env = SPECENV

# ╔═╡ 0881c790-2462-11eb-05e1-9113e664ec28
md"First, let’s create a light source! We have several options for these."

# ╔═╡ 2f7258b0-2462-11eb-31a2-8da7d9fb1dca
md"One option is to create a black body illuminant. Let’s say our light source is a black body radiator at a temperature of 6500K. We want to create an illuminant in the visual range of wavelength, say between 390 and 830 nm, with a resolution of 1 nm:"

# ╔═╡ 5216c132-2462-11eb-11ef-259ba512e80a
bb6500K = blackbody_illuminant(env,6500)

# ╔═╡ 60284748-a78f-4dfa-9927-18ad85bdf5d0
md"What is the type of the spectrum?"

# ╔═╡ ad2414b2-2462-11eb-3ea5-955a5a918a30
typeof(bb6500K)

# ╔═╡ 27bae2a1-2dc1-4a75-84ed-d329ce62abea
md"`LSpec` stands for Luminance Spectrum. We are dealing with a light source."

# ╔═╡ 437e5e6d-1469-4767-a0a9-8fed5ecf6250
md"Show the wavelengths at which this spectrum is sampled:"

# ╔═╡ d99b2330-2462-11eb-2bee-6b79c56cf0cf
bb6500K.λ

# ╔═╡ 7db83947-d5dd-4a3a-9f63-efe1da7eb26c
md"And the sampled luminance values:"

# ╔═╡ e34ef68e-2462-11eb-394c-5712ac9de6da
bb6500K.l

# ╔═╡ e8969630-2462-11eb-3c02-996a99d88691
md"So, the `LSpec` struct has two fields: λ and l.

`λ` contains the sampled wavelengths in *nm*, `l` contains the sampled *luminance* at the sampled wavelengths. For reflectance spectra of type `RSpec`, these are `λ` and `r` (reflectance)."

# ╔═╡ f24e2d8e-2477-11eb-3a27-dfebdfc43440
md"Let’s plot the spectrum of our light source:"

# ╔═╡ 0c6bf812-2478-11eb-14ea-57c6522324f3
plot(bb6500K.λ,bb6500K.l,labels=false)

# ╔═╡ 2c7f7c80-2478-11eb-1728-a1ab0b5ef106
md"Whoa! The luminance values are crazy high! That’s because of the units `[J/m³/s]`. It’s easy to mitigate that. Let’s normalize this spectrum to a peak value of 1.0. To normalize a spectrum, use the `normalize_spec` function as follows:"

# ╔═╡ 00a64130-2463-11eb-1ed2-29a671508315
bb65=normalize_spec(bb6500K)

# ╔═╡ 1e781300-2463-11eb-2b3e-ff33d0043b07
md"Let’s have a look at the spectrum in a graph:"

# ╔═╡ 369a82b0-2463-11eb-06c0-a79be68e16a6
plot(bb65.λ,bb65.l,labels=false)

# ╔═╡ 47f2352e-2463-11eb-3976-d77c1624e204
md"That looks much better! It’s convenient and common practice in color metrics to normalize illuminants to a peak value of 1."

# ╔═╡ 60fc2740-247a-11eb-031a-8bef9409f862
md"Another useful method is using a CIE Daylight generator, which creates daylight spectra according to the CIE daylight standard. This mimics the spectral distribution of natural daylight spectra. `D_series_illuminant` produces daylight spectra, like *D65* or *D50*, according to a so-called *correlated color temperature*, also abbreviated as *CCT*. To create an illuminant with a CCT of 6500K, similar to the blackbody radiator, use the `D_series_illuminant` function as follows:"

# ╔═╡ fdfa8220-2462-11eb-0cf8-c188ac491781
d6500=D_series_illuminant(6500)

# ╔═╡ 89f88560-2463-11eb-38be-eb669b15ab60
md"Before we plot this D65 spectrum, let’s normalize it as well. The standard CIE daylight function generates spectra that are normalized to a value 100.0 at a wavelength of 560 nm, which results in spectra with a spectral distribution that mimics the proper relative energy output of different D-series illuminants. For our purposes, let’s normalize the D65 spectrum:"

# ╔═╡ 97598240-2463-11eb-0ec9-0be99eea414d
d65=normalize_spec(d6500)

# ╔═╡ d0a02310-2463-11eb-2788-5dcb718cf9cb
md"Now we can plot both spectra for comparison:"

# ╔═╡ d811596e-2463-11eb-0de0-792feca7ffa0
begin
	plot(bb65.λ, bb65.l, linewidth=2, color= :blue, label="black body at 6500K")
	plot!(d65.λ, d65.l, linewidth=2, color= :gold2, label="D65 daylight generator")
end

# ╔═╡ e8a8a220-2463-11eb-2a70-cd5fbe573691
md"Voilá! Looking good!"

# ╔═╡ fe2f5850-2463-11eb-29a7-93f2db686897
md"But we’d like both spectra to have the *same wavelength range*. And both should also have the *same spectral resolution*. The *d65* spectrum has a resolution of only 5 nm. Let’s change that!

For demonstration purposes, we will first take a short detour to the `SSpline` package. Let’s use it to adapt the resolution and range of the *d65* spectrum.
Later, we’ll see a more convenient way to adapt spectra to arbitrary conditions, resolutions, etc.

Now, let’s do it the `SSpline` way first:"

# ╔═╡ 69e30fb0-2464-11eb-0066-d15aa8c969e3
md"## A quick look at the `SSpline` package"

# ╔═╡ 9b44a9b0-2464-11eb-138b-0d9775344d9f
md"First, we create a spline object. The Daylight illuminant generator creates linearly interpolated values, so we don’t need to bother with any fancy spline interpolations in this case. We only need a linear spline:"

# ╔═╡ c1a3a930-2464-11eb-1bef-996a6b7468f5
d65spl = linearspline(d65.λ,d65.l)

# ╔═╡ c77141b0-2464-11eb-1efa-a73f2a72b3b0
md"As a result, we get an object of type `LinearSpline`. It contains all the spline coefficients at the knots of the spline, and the knots themselves:"

# ╔═╡ c5490aa0-254d-11eb-383c-adb08ba2eedf
typeof(d65spl)

# ╔═╡ 07995c00-2465-11eb-125c-b5a167e0fc90
md"Now we want to interpolate this spline at certain values. For the black body illuminant we used the range `390:1:830`, and these are going to be our new `λ` values for the D65 illuminant as well:"

# ╔═╡ 385b81b0-2465-11eb-151f-878bb717e157
λ=collect(390.0:1.0:830.0)

# ╔═╡ a57b2320-247b-11eb-17d1-9fd1ecce6ed5
md"Now we create the new, interpolated values with the `interp` function. The function returns two vectors, one for the `λ` values, one for the luminance values:"

# ╔═╡ 492e25b0-2465-11eb-36eb-491d2c186e8b
lam,s=interp(d65spl,λ)

# ╔═╡ e54233b0-2465-11eb-1f8d-2f087aa4d31e
md"Now all that’s left is creating a new *Spectrum* type object containing these values. In our case we need to create a luminance spectrum because we want the spectrum to be a light source:"

# ╔═╡ 0c92e630-2466-11eb-20c4-9fe8bede9360
d65f=luminance_spec(lam,s)

# ╔═╡ 1d179140-2466-11eb-2d71-11e9622e355b
md"What does it look like in the graph?"

# ╔═╡ cb75ef1e-2466-11eb-26b4-d58a8167559c
begin
	plot(d65f.λ,d65f.l, linestyle=:dash, labels=false)
	plot!(bb65.λ,bb65.l,linestyle=:solid, labels=false)
end

# ╔═╡ ddff1cc0-2466-11eb-0709-1baba5e3fca6
md"Excellent!"

# ╔═╡ 0e233640-247e-11eb-193a-ed5ab087c968
md"## Back to Spectra"

# ╔═╡ f9023e30-2466-11eb-346f-8380306e6ab3
md"There are other ways to generate spectra. The function `led_spec` simulates the spectrum of a light emitting diode:"

# ╔═╡ b9028080-2473-11eb-3bc2-91cb08806a45
led=led_spec(390,830,1,455,30)

# ╔═╡ dde6edf0-2473-11eb-2671-2b1f152e2ced
md"Creates an LED spectrum in the range 390 to 830 nm, with a wavelength step width of 1 nm. The *peak wavelength* is at 455 nm, and the *half spectral width* (bandwidth at 50% of peak intensity) is 30 nm. Let’s have a look at the spectrum:"

# ╔═╡ 62d4c690-2474-11eb-2665-c724b32a74f4
plot!(led.λ,led.l,linestyle=:solid, labels=false)

# ╔═╡ 7bd22cf0-2474-11eb-0266-8bb6aeb08b46
md"And last, but not least, there is a method to create a special kind of spectrum, the so-called *block spectrum*. This type of spectrum is the theoretical extreme of spectra that can be produced. Physically, these do not exist, but they are very useful, if not essential, in researching e.g. the gamuts, which means the range of all theoretically producible colors of color spaces:"

# ╔═╡ f95f7100-2474-11eb-3072-f1d2fdb0ba57
bspec=block_spec(390,830,1,501,533,1.0,0)

# ╔═╡ 23b4ff10-2475-11eb-3324-176b6bc48b81
plot!(bspec.λ,bspec.r,color=:orange, labels=false)

# ╔═╡ 3c46d260-2475-11eb-0204-c3da17762579
md"As you can see, this kind of spectrum only has the values 0.0 or 1.0. This creates the theoretically most saturated colors of any given bandwidth. Block spectra exist as two different types. **Type 1**, which you can see above, creates the purest possible colors with a single center wavelength. Type 1 block spectra are zero everywhere except in the block defined by min und max values.

### Usage:

`block_spec(λmin,λmax,Δλ,λs,λe,str=1.0,mode=0)`

`λmin, λmax, Δλ`: range and resolution of spectrum

`λs,λe`: range of block peak

`str=1.0`: strength of block, default=1.0. Setting this to other values can be useful in certain cases.

`mode=0`: mode for single peak block spectrum, default = 0
"

# ╔═╡ b3ee6870-247d-11eb-0587-ef4e7989be87
md"Type 1 block spectra create the inverse of type 0 block spectra. They are 1 (or a defined maximum value) everywhere except between the defined min and max values:"

# ╔═╡ bbd38ff0-2475-11eb-30e6-939c1fae753f
bspec0=block_spec(390,830,1,410,570,1.0,1)

# ╔═╡ 0c572270-2476-11eb-0972-89e9f7f99f3d
plot!(bspec0.λ,bspec0.r,color=:black,linestyle=:dot, labels=false)

# ╔═╡ b2a5b300-2505-11eb-1ecc-098dda235b36
md"## Spectral environments, or how to adapt spectral properties"

# ╔═╡ e880d120-255b-11eb-12d4-d945762279a1
md"Let’s define a light source and a reflective surface. We use the spectrum of the Martian sky, as measured by the Opportunity rover, as light source:"

# ╔═╡ 10fbabc0-255c-11eb-0225-fb05d87dd7f2
light = mars_sky_opportunity()

# ╔═╡ 0e81b72c-1f09-4145-bd54-058ae1b7ded9
md"And our reflective surface is the spectrum of the blue patch of the calibration target on the Mars Exploration Rover:"

# ╔═╡ 23f09380-255c-11eb-1d49-89943d828c8c
blue = mer_caltarget(:blue)

# ╔═╡ 3e1e3c80-255c-11eb-23dc-75a55d31786e
begin
	plot(blue.λ,blue.r,color=:blue,label="MER caltarget “blue”")
	plot!(light.λ,light.l,color=:orange,label="Mars Sky Opportunity measurement")
end

# ╔═╡ a9cb5ab0-8b16-11eb-3d10-dfc3953fbe2b
blue

# ╔═╡ f77e81a0-8a72-11eb-1e00-f9266532685f
light

# ╔═╡ 4c3765d0-255c-11eb-3e16-c98892838d0e
md"It’s obvious that both spectra have different ranges and also different and irregular sample points. The surface spectrum even reaches pretty far into the IR region. This can easily be mitigated by adapting both spectra to a set of rules defined by the *spectral environment*.
For that, we need to define a **spectral environment**, which contains all information about the overall valid traits of the spectra we want to work with. This is important to ensure and enable compatibility between spectra etc. which might have different wavelength ranges, resolutions, etc. The function to set up a spectral environment is:

`set_specenv(λmin,Δλ,λmax,extrapolation)`

possible values for extrapolation:

`:zero`, `:constant`, `:linear`, `:quadratic`"

# ╔═╡ 15adccb0-255d-11eb-2b83-c9564c07b62b
md"We want to work in a spectral range that satisfies the latest color matching functions defined by the CIE in 2012, which are defined in a range of 390 nm to 830 nm in steps of 1 nm. Furthermore, for the sake of simplicity, we define that all extrapolated values should be set to zero. Extrapolation is necessary if a spectrum is not defined up to the boundaries (390 / 830 nm)."

# ╔═╡ 1a26c720-255c-11eb-09fe-cb4cd9cafd46
senv=set_specenv();senv.ex=:boundary

# ╔═╡ a2061550-8a70-11eb-1d1c-1b884f18ae48
senv

# ╔═╡ a386f890-255d-11eb-0fae-2326d5d88d34
md"To adapt a spectrum to the environment, use the function

`adapt_spec(spec::Spectrum,specenv::SpecEnvironment,interporder...)`

Valid arguments for `interporder` (order of spline interpolation):

1st order interpolation:

`:linear` for linear splines

2nd order interpolation, not yet implemented:

`:quadratic` for quadratic splines

3rd order interpolation:

`:cubic, :natural` for natural cubic splines (zero curvature at boundaries)

`:cubic, :periodic` for periodic cubic splines (both endpoints must have same value)

`:cubic, :clamped, slope_start, slope_end` for clamped cubic splines with defined slopes at boundaries

`:cubic, :constrained` for cubic constrained splines"

# ╔═╡ 0c5a9a70-255e-11eb-26e0-8f1b434ed385
l=adapt_spec(light,senv,:cubic,:natural)

# ╔═╡ 42a8468e-255e-11eb-0495-5b80974193c4
b=adapt_spec(blue,senv,:cubic,:natural)

# ╔═╡ 59847aa0-255e-11eb-36a6-8f10abd28695
md"Now we can check and compare the original and the adapted versions of the splines. The boxes mark the original irregular spectral samples, the smaller circles mark the resampled spectra as defined in the environment. You can see that for the **blue** spectrum the lower range was extrapolated according to the cubic spline defined by the original knots. Here is a zoomed in version of the two spectra:"

# ╔═╡ cf17c5be-255d-11eb-0cb9-3703acb4ae96
begin
	scatter(blue.λ,blue.r,color=:blue,marker=:square,markersize=6,label="MER caltarget “blue”",legend=:topleft,xlim=(385,430),ylim=(0.1,0.55))
	scatter!(light.λ,light.l,color=:orange,marker=:square,markersize=6,label="Mars Sky Opportunity measurement")
	scatter!(b.λ,b.r,color=:lightblue,linestyle=:dash,markersize=2,label="MER caltarget resampled")
	scatter!(l.λ,l.l,color=:gold,linestyle=:dash,markersize=2,label="Mars Sky Opportunity resampled")
end

# ╔═╡ d5d44030-255f-11eb-1bf2-fbc628979e90
md"Let’s have a look at the opposite end:"

# ╔═╡ c8bd44a0-255f-11eb-06d0-c911969126bf
begin
	scatter(blue.λ,blue.r,color=:blue,marker=:cross,markersize=6,label="MER caltarget “blue”",legend=:bottomleft,xlim=(800,850),ylim=(-0.1,1.0))
	scatter!(light.λ,light.l,color=:orange,marker=:cross,markersize=6,label="Mars Sky Opportunity measurement")
	scatter!(b.λ,b.r,color=:lightblue,linestyle=:dash,markersize=2,label="MER caltarget spec_env")
	scatter!(l.λ,l.l,color=:yellow,linestyle=:dash,markersize=2,label="Mars Sky Opportunity spec_env")
end

# ╔═╡ 9016a8f8-561f-4f0c-bd5e-8e3437a36801
begin
	scatter(blue.λ,blue.r,color=:blue,marker=:square,markersize=3,label="MER caltarget “blue”",legend=:topleft,xlim=(380,840))
	scatter!(light.λ,light.l,color=:orange,marker=:square,markersize=3,label="Mars Sky Opportunity measurement")
	scatter!(b.λ,b.r,color=:lightblue,linestyle=:dash,markersize=1,label="MER caltarget resampled")
	scatter!(l.λ,l.l,color=:gold,linestyle=:dash,markersize=1,label="Mars Sky Opportunity resampled")
end

# ╔═╡ e74f5600-255b-11eb-3a98-718912d58c6c
md"# Spectral operators"

# ╔═╡ d3d88a20-2505-11eb-3fe8-a7183ff72a46
md"We’re actually interested the resulting color and spectrum with a given light source and reflectance. `SpectralVis` contains a few overloaded operators that make this kind of calculation easy. Let’s define a light source and a reflecting surface:"

# ╔═╡ 3d0f5aa0-2506-11eb-3072-3d7840704c25
lightsource = mars_sky_opportunity()

# ╔═╡ 1f18b440-85f3-4e23-bb8d-1c7889c1c7c7
surface = mer_caltarget(:green)

# ╔═╡ c5ff81a0-2506-11eb-305c-8721c7eb5987
green_spl = spline3(surface.λ,surface.r,:natural)

# ╔═╡ 1497e433-4a73-4074-8b7e-228d3f029b50
green_extrap = extrap(green_spl,390:830,:linear)

# ╔═╡ 0b7b056d-0f00-4f82-873f-53825874e76f
green_int = interp(green_extrap,390:1:830)

# ╔═╡ cd636a24-b558-4c33-adca-ece8652b1640
light_spl = spline3(lightsource.λ,lightsource.l,:natural)

# ╔═╡ ad8c32f7-0a8f-4ab0-a088-c688945f96be
light_extrap = extrap(light_spl,390:830,:linear)

# ╔═╡ 9cc68e94-7008-4af2-bd56-fdaa922de29e
light_int = interp(light_spl,390:830)

# ╔═╡ 924c6ba0-2508-11eb-25d9-8fe5faf11bad
begin
	plot(surface.λ,surface.r,color=:green,label="MER caltarget “green”")
	plot!(lightsource.λ,lightsource.l,color=:orange,label="Mars Sky Opportunity measurement")
	scatter!(green_int[1],green_int[2],color=:cyan,markersize=1)
	scatter!(light_int[1],light_int[2],color=:orange,markerkize=2)
end

# ╔═╡ 32c456b0-2ed2-11eb-2c69-f3989cde9b93
surface.λ

# ╔═╡ 3e9c69a0-2ed2-11eb-3dbf-e7a0ecb05d4c
surface.r

# ╔═╡ 67776940-8b22-11eb-0fa2-3de9f03d6679
lightsource.λ[end]

# ╔═╡ af8d15e0-8b22-11eb-31c9-2be22127eb1d
begin
	senv.λmax=840.0
	senv.ex=:boundary
	adapt_spec(lightsource,senv,:cubic,:natural)
end

# ╔═╡ 6126a56a-681f-4cce-b1f5-5cf6138d51bf
senv

# ╔═╡ 899288cb-714a-4399-a9ab-911b2e0c2c37
extrap(linearspline(lightsource.λ,lightsource.l),390.0:890.0,senv.ex)

# ╔═╡ 8a955ace-2e99-11eb-25dd-e9a7242b4f6c
begin
	lightspec=adapt_spec(lightsource,senv,:linear)
	greenspec=adapt_spec(surface,senv,:cubic,:natural)
end

# ╔═╡ 0f23dce0-89d5-11eb-0b0b-7193920bef8b
begin
	plot(lightspec.λ,lightspec.l,color=:red,style=:dashdot,thickness=3,label="Martian sky (light source)")
	plot!(greenspec.λ,greenspec.r,color=:green,label="MER caltarget green (refl)")
	#plot!(refl.λ,refl.r,color=:black,style=:dashdot,label="reflected light")
	#plot!(greenfilt.λ,greenfilt.t)
end

# ╔═╡ afcd7890-89d5-11eb-30db-4f57ca767fef
begin
	refl=*(greenspec,greenspec)
	#greenf=TSpec(greenspec.λ,greenspec.r,.05)
	#greenfilt=greenf*greenf
end

# ╔═╡ 2c64c8da-779e-4a2d-b1fd-cb739b316a17
begin
	→(s::SpectralVis.Spectrum,Δλ) = shift_spec(s,Δλ)
	←(s::SpectralVis.Spectrum,Δλ) = shift_spec(s,-Δλ)
end

# ╔═╡ 66307660-89da-11eb-0576-8d571268e6dd
plot(led.λ,led.l)

# ╔═╡ c5e0d2c5-1f93-4475-b8ae-8a40d35f3db5
begin
	d65s=→(d65,23.3)
	plot(d65s.λ,d65s.l)
	plot!(d65.λ,d65.l)
end

# ╔═╡ 983b7c23-29fb-452d-beba-2aec2b8abe22
cm=cmf(LMS2006_10)

# ╔═╡ 0213f1eb-4cbe-4aa1-b8c1-4a469836723a
begin
	long=cm.l
	mid=cm.m
	short=cm.s
end

# ╔═╡ 70626740-a681-4da3-956c-1ee4c4badcdb
begin
	plot(cm.λ,cm.l,yrange=(0.000001,1.1),color=:red)
	plot!(cm.λ,cm.m,color=:green)
	plot!(cm.λ,cm.s,color=:blue)
end

# ╔═╡ 7fd0091a-640d-4df0-aabe-646c5f8742e0
source=D_series_luminance(6500)

# ╔═╡ 37dac578-8fc2-4c04-87d6-81b4822fc7ed
Macbeth_color_chart

# ╔═╡ 5e6acbc7-f548-4b5b-8774-443a92e06db8
begin
	cc = env.cmf
	#dd = shift_colormatch(cc, :l, -20)
	plot(cc.λ, cc.l, color = :red)
	plot!(cc.λ, cc.m, color = :green)
	plot!(cc.λ, cc.s, color = :blue)
end

# ╔═╡ Cell order:
# ╠═f07cef7e-014a-4d59-8e31-4d394391adda
# ╟─78adafde-247e-11eb-12e1-ddea99d50711
# ╟─a8bb86de-29be-474d-94e6-1a03f9afe70f
# ╟─f444ebf0-2460-11eb-2ed3-c5db8e304e95
# ╟─1cb25910-2461-11eb-1f85-a9ed9362e326
# ╠═08af84c0-2456-11eb-326c-2d675ca197af
# ╟─5dacb780-2461-11eb-277d-e7226385a5ec
# ╠═57a60da0-2461-11eb-07df-2b3ab48a6675
# ╟─79763630-2461-11eb-0b0d-15b00898c512
# ╠═9ff8ea50-2461-11eb-3143-d1b38243f297
# ╟─0d1268b0-2461-11eb-15a9-87ca552b10c5
# ╠═14f2d870-245d-11eb-1fe3-779aa49070f1
# ╟─1afb7e20-245d-11eb-1ca8-9da92491cd50
# ╟─d04f4ee1-2609-4355-9846-5e155557d709
# ╟─7d8d8c1e-ed75-472f-b08b-a93c29ab588c
# ╠═f5760d98-77e3-4847-bb79-97f45b50b7da
# ╟─0881c790-2462-11eb-05e1-9113e664ec28
# ╟─2f7258b0-2462-11eb-31a2-8da7d9fb1dca
# ╠═5216c132-2462-11eb-11ef-259ba512e80a
# ╟─60284748-a78f-4dfa-9927-18ad85bdf5d0
# ╠═ad2414b2-2462-11eb-3ea5-955a5a918a30
# ╟─27bae2a1-2dc1-4a75-84ed-d329ce62abea
# ╟─437e5e6d-1469-4767-a0a9-8fed5ecf6250
# ╠═d99b2330-2462-11eb-2bee-6b79c56cf0cf
# ╟─7db83947-d5dd-4a3a-9f63-efe1da7eb26c
# ╠═e34ef68e-2462-11eb-394c-5712ac9de6da
# ╟─e8969630-2462-11eb-3c02-996a99d88691
# ╟─f24e2d8e-2477-11eb-3a27-dfebdfc43440
# ╠═0c6bf812-2478-11eb-14ea-57c6522324f3
# ╟─2c7f7c80-2478-11eb-1728-a1ab0b5ef106
# ╠═00a64130-2463-11eb-1ed2-29a671508315
# ╟─1e781300-2463-11eb-2b3e-ff33d0043b07
# ╠═369a82b0-2463-11eb-06c0-a79be68e16a6
# ╟─47f2352e-2463-11eb-3976-d77c1624e204
# ╟─60fc2740-247a-11eb-031a-8bef9409f862
# ╠═fdfa8220-2462-11eb-0cf8-c188ac491781
# ╟─89f88560-2463-11eb-38be-eb669b15ab60
# ╠═97598240-2463-11eb-0ec9-0be99eea414d
# ╟─d0a02310-2463-11eb-2788-5dcb718cf9cb
# ╠═d811596e-2463-11eb-0de0-792feca7ffa0
# ╟─e8a8a220-2463-11eb-2a70-cd5fbe573691
# ╟─fe2f5850-2463-11eb-29a7-93f2db686897
# ╟─69e30fb0-2464-11eb-0066-d15aa8c969e3
# ╟─9b44a9b0-2464-11eb-138b-0d9775344d9f
# ╠═c1a3a930-2464-11eb-1bef-996a6b7468f5
# ╟─c77141b0-2464-11eb-1efa-a73f2a72b3b0
# ╠═c5490aa0-254d-11eb-383c-adb08ba2eedf
# ╟─07995c00-2465-11eb-125c-b5a167e0fc90
# ╠═385b81b0-2465-11eb-151f-878bb717e157
# ╟─a57b2320-247b-11eb-17d1-9fd1ecce6ed5
# ╠═492e25b0-2465-11eb-36eb-491d2c186e8b
# ╟─e54233b0-2465-11eb-1f8d-2f087aa4d31e
# ╠═0c92e630-2466-11eb-20c4-9fe8bede9360
# ╟─1d179140-2466-11eb-2d71-11e9622e355b
# ╠═cb75ef1e-2466-11eb-26b4-d58a8167559c
# ╟─ddff1cc0-2466-11eb-0709-1baba5e3fca6
# ╟─0e233640-247e-11eb-193a-ed5ab087c968
# ╟─f9023e30-2466-11eb-346f-8380306e6ab3
# ╠═b9028080-2473-11eb-3bc2-91cb08806a45
# ╟─dde6edf0-2473-11eb-2671-2b1f152e2ced
# ╠═62d4c690-2474-11eb-2665-c724b32a74f4
# ╟─7bd22cf0-2474-11eb-0266-8bb6aeb08b46
# ╠═f95f7100-2474-11eb-3072-f1d2fdb0ba57
# ╠═23b4ff10-2475-11eb-3324-176b6bc48b81
# ╟─3c46d260-2475-11eb-0204-c3da17762579
# ╟─b3ee6870-247d-11eb-0587-ef4e7989be87
# ╠═bbd38ff0-2475-11eb-30e6-939c1fae753f
# ╠═0c572270-2476-11eb-0972-89e9f7f99f3d
# ╟─b2a5b300-2505-11eb-1ecc-098dda235b36
# ╠═e880d120-255b-11eb-12d4-d945762279a1
# ╠═10fbabc0-255c-11eb-0225-fb05d87dd7f2
# ╟─0e81b72c-1f09-4145-bd54-058ae1b7ded9
# ╠═23f09380-255c-11eb-1d49-89943d828c8c
# ╠═3e1e3c80-255c-11eb-23dc-75a55d31786e
# ╠═a9cb5ab0-8b16-11eb-3d10-dfc3953fbe2b
# ╠═f77e81a0-8a72-11eb-1e00-f9266532685f
# ╟─4c3765d0-255c-11eb-3e16-c98892838d0e
# ╟─15adccb0-255d-11eb-2b83-c9564c07b62b
# ╠═1a26c720-255c-11eb-09fe-cb4cd9cafd46
# ╠═a2061550-8a70-11eb-1d1c-1b884f18ae48
# ╟─a386f890-255d-11eb-0fae-2326d5d88d34
# ╠═0c5a9a70-255e-11eb-26e0-8f1b434ed385
# ╠═42a8468e-255e-11eb-0495-5b80974193c4
# ╟─59847aa0-255e-11eb-36a6-8f10abd28695
# ╠═cf17c5be-255d-11eb-0cb9-3703acb4ae96
# ╟─d5d44030-255f-11eb-1bf2-fbc628979e90
# ╠═c8bd44a0-255f-11eb-06d0-c911969126bf
# ╠═9016a8f8-561f-4f0c-bd5e-8e3437a36801
# ╟─e74f5600-255b-11eb-3a98-718912d58c6c
# ╟─d3d88a20-2505-11eb-3fe8-a7183ff72a46
# ╠═3d0f5aa0-2506-11eb-3072-3d7840704c25
# ╠═1f18b440-85f3-4e23-bb8d-1c7889c1c7c7
# ╠═c5ff81a0-2506-11eb-305c-8721c7eb5987
# ╠═1497e433-4a73-4074-8b7e-228d3f029b50
# ╠═0b7b056d-0f00-4f82-873f-53825874e76f
# ╠═cd636a24-b558-4c33-adca-ece8652b1640
# ╠═ad8c32f7-0a8f-4ab0-a088-c688945f96be
# ╠═9cc68e94-7008-4af2-bd56-fdaa922de29e
# ╠═924c6ba0-2508-11eb-25d9-8fe5faf11bad
# ╠═32c456b0-2ed2-11eb-2c69-f3989cde9b93
# ╠═3e9c69a0-2ed2-11eb-3dbf-e7a0ecb05d4c
# ╠═67776940-8b22-11eb-0fa2-3de9f03d6679
# ╠═af8d15e0-8b22-11eb-31c9-2be22127eb1d
# ╠═6126a56a-681f-4cce-b1f5-5cf6138d51bf
# ╠═899288cb-714a-4399-a9ab-911b2e0c2c37
# ╠═8a955ace-2e99-11eb-25dd-e9a7242b4f6c
# ╠═0f23dce0-89d5-11eb-0b0b-7193920bef8b
# ╠═afcd7890-89d5-11eb-30db-4f57ca767fef
# ╠═2c64c8da-779e-4a2d-b1fd-cb739b316a17
# ╠═66307660-89da-11eb-0576-8d571268e6dd
# ╠═c5e0d2c5-1f93-4475-b8ae-8a40d35f3db5
# ╠═983b7c23-29fb-452d-beba-2aec2b8abe22
# ╠═0213f1eb-4cbe-4aa1-b8c1-4a469836723a
# ╠═70626740-a681-4da3-956c-1ee4c4badcdb
# ╠═7fd0091a-640d-4df0-aabe-646c5f8742e0
# ╠═37dac578-8fc2-4c04-87d6-81b4822fc7ed
# ╠═5e6acbc7-f548-4b5b-8774-443a92e06db8


text:
module LessOLS

include("sample.jl")
include("lm.jl")

end # module

ext: jl


text:
mutable struct Player <: AbstractPlayer
  game::Game{Vector{AbstractPlayer}}
  grid::Union{AbstractGrid, Nullable{AbstractGrid}}

  bag::Union{AbstractBag, Nullable{AbstractBag}}

  piece::Union{AbstractPiece, Nullable{AbstractPiece}}
  shadow::Union{AbstractPiece, Nullable{AbstractPiece}}
  hold::Union{AbstractPiece, Nullable{AbstractPiece}}

  round::Union{AbstractRound, Nullable{AbstractRound}}
  state::Union{AbstractState, Nullable{AbstractState}}
  clock::Union{AbstractClock, Nullable{AbstractClock}}

  action::AbstractString
  chain::AbstractString

  level::Int
  score::Int

  lines::Int
  combo::Int

  stash::Int
  glues::Int
end

function Player(cur_game::AbstractGame=Game())
  cur_player = Player(
    cur_game,
    Nullable{AbstractGrid}(),
    Nullable{AbstractBag}(),
    Nullable{AbstractPiece}(),
    Nullable{AbstractPiece}(),
    Nullable{AbstractPiece}(),
    Nullable{AbstractRound}(),
    Nullable{AbstractState}(),
    Nullable{AbstractClock}(),
    "",
    "",
    1,
    0,
    0,
    0,
    0,
    0
  )

  push!(cur_game.players, cur_player)

  State(cur_player)
  Clock(cur_player)

  Grid(cur_player)
  Bag(cur_player)

  cur_player.shadow = Piece(cur_player)

  cur_player
end


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 250252f0-6e61-11ec-2ffe-ebddaa9f0fb3
begin
	using Pkg
	Pkg.activate("../../")
	using Revise
end

# ╔═╡ 0eeb5bb1-bb41-43cc-8c14-eb9c1ebf6fad
using ImgProcMic, FourierTools, IndexFunArrays, ImageShow, FFTW, TestImages, Colors

# ╔═╡ 1b2ebe1a-f3a7-494c-bcca-18df5403ee65
using PlutoTest, Noise, Tullio, PlutoUI, Images, Plots

# ╔═╡ 4aeceb65-f226-41cf-988b-273f42c75059
using SpecialFunctions

# ╔═╡ c11d25d8-8fe9-476e-9a6d-66f625f29348
md"# 1 Convolution with Small Kernels

When the kernel (e.g. the PSF) of a convolution is spatially very small (e.g. $3 \times 3$), it is more efficient to perform the convolution with a sum approach (_for loops_) in real space.

In this part, you implement the convolution via the sum approach.

You will notice that the question arises, how to treat the boundary pixels. We follow the approach, that the output image is smaller than the input image. This is due to the reason that we only convolve inside the image where the kernel fits completely in.

If you need a graphical visualization, see [here](https://kharshit.github.io/img/edge_detection.png).
"

# ╔═╡ 6f87bb34-ae2b-4c8e-9c21-621644445df3
"""
	conv_s(img, kernel)

Convolution of `img` with `kernel`.

## Procedure
* Determine the correct output size
* Move the kernel over the image like [here](https://kharshit.github.io/img/edge_detection.png).
"""
function conv_s(img, kernel)
	out = similar(img, size(img) .- size(kernel) .+ 1)

	# runtime complexity is:
	for i in 1:size(out, 1)
		for j in 1:size(out, 2)
			# value we are going to write into out
			# that value accumulates the elementwise multiplication
			# and summation
			v = zero(eltype(img))
			for ki in 1:size(kernel, 1)
				for kj in 1:size(kernel, 2)
					# -1 because of Julia's one-based indexing style
					v += kernel[ki, kj] * img[i + ki - 1, j + kj - 1]
				end
			end
			out[i, j] = v
		end
	end
	return out
end

# ╔═╡ b030ccf9-9a48-4410-8ac1-f3a2a02b5b5c
md"runtime complexity is:

$\mathcal{O}(L^2 \cdot K^2)$

where $L$ is the image size, and $K$ the kernel size.
"

# ╔═╡ a879425f-a152-4d4d-b6b8-f2a03b7a452f
img = Float32.(testimage("mandril_gray"));

# ╔═╡ 012dd588-350b-4c2f-99c9-6371e443e32c
similar(img, (2,2))

# ╔═╡ 61ccfb00-513e-4bca-80eb-133aa58a1def
begin
	# todo create a kernel of size (5,5) filled with 1s
	# but afterwards normalize the sum to 1
	kernel = ones(Float32, (5,5))
	kernel ./= sum(kernel)
end

# ╔═╡ d51fadc0-3227-4e80-96f2-269a8fd2bac2
img_b = conv_s(img, kernel);

# ╔═╡ caab13b5-12a6-4058-8931-1f8a9abe57e8
[gray_show(FourierTools.select_region(img, new_size=size(img_b))) gray_show(img_b)]

# ╔═╡ cc090365-274f-4e56-9be5-cbc1a32295ad
kernel_h_filter = repeat([1 0 -1], inner=(3, 1))

# ╔═╡ b63022e3-3197-4dfe-9b7f-3ee8babce73a
md"
Those examples show the vertical and horizontal edges of the mandril
"

# ╔═╡ 7ec68244-fe43-4418-b9ee-0cb7e67572c9
[gray_show(conv_s(img, kernel_h_filter)) gray_show(conv_s(img, kernel_h_filter'))]

# ╔═╡ 2bd635be-6f07-40df-a25e-cfdf2dd22d6f
md"## Benchmark with FFT based
Is your method faster than FFT based conv?
Rerun below cell more often!

I guess, you won't beat Tullio!
"

# ╔═╡ 3765b001-deb1-49b2-87a1-0100d36d3087
function conv_tullio(img, k)
	@tullio out_t[i, j] := (img[i + ik - 1, j + jk - 1] *
							kernel[ik, jk]) (i in 1:508) (j in 1:508)
end

# ╔═╡ 57427d81-961c-488e-9a9d-e625a2e536c3
# run me multiple times!
with_terminal() do
	print("Real space conv\t\t")
	@time conv_s(img, kernel)
	print("FFT based conv\t\t")
	kernel_big = ifftshift(FourierTools.select_region(kernel, new_size=size(img)))
	@time conv(img, kernel_big)
	print("Tullio based conv\t")
	@time conv_tullio(img, kernel)
end

# ╔═╡ 6b4028de-e605-455f-947d-62649f2e26fa
md"## 1 Test"

# ╔═╡ 53ae1016-55d7-4545-a29d-b9a4fc90b854
PlutoTest.@test size(conv_s(ones((3,3)), ones((2,2)))) == ((2,2))

# ╔═╡ 01428578-4258-413d-8b8b-27203527fd2b
img_tt = zeros((6,6));

# ╔═╡ 81605617-d091-47b2-bd85-1cd8ae0fef18
img_tt[:, 1:3] .= 10;

# ╔═╡ f3fb21e1-990a-42ac-8a5d-52154be23cfb
# error disappears if size is correct
PlutoTest.@test conv_s(img_tt, kernel_h_filter) ≈ [0.0 30.0 30.0 0.0; 0.0 30.0 30.0 0.0; 0.0 30.0 30.0 0.0; 0.0 30.0 30.0 0.0]

# ╔═╡ 94694b2a-715a-450b-a024-fc474ce19975
md"# 2 Cell Detection

Based on correlation, we can come up with a simple cell detection algorithm. For the detection we need a template cell.

However, such correlation based methods are not very robust with respect to noise and distortions.

Btw, what you implemented above is mathematically not a convolution but a correlation (we didn't reverse the kernel).

The equation we want to implement in this part, is called _normalized cross-correlation_ since it divides by the local standard deviation of intensities:

$\text{NCC}\{I, K\}[i, j] = \frac{\sum_{m, n} I[i + m - 1, j + n - 1] \cdot K[m, n]}{\sqrt{\sum_{m, n} I[i + m - 1, j + n - 1]^2} \cdot \sqrt{\sum_{m, n} K[m, n]^2}}$
"

# ╔═╡ f43bd606-d202-4f74-8e46-28bab5af6d45
# load cells
cells = 1 .- Float32.(Gray.(testimage("blobs")));

# ╔═╡ 4ed29eeb-0cbc-4810-8a6c-1b7aaa0014be
gray_show(cells)

# ╔═╡ 1664be9b-08bc-4144-b191-b82bc9e2cee2
"""
	get_template(cells[, t_size=(20, 20)])

Extracts a template cell at position `pos` with template size `t_size`.

# Procedure
* Start the extract at `pos`.
* End the extraction at `pos .+ t_size`.
* Return the extracted part
"""
function get_template(cells, pos, t_size=(20, 20))
	# we need the -1 because of the one based indexing in julia
	return cells[pos[1]:pos[1] + t_size[1] - 1, pos[2]:pos[2] + t_size[2] - 1]
end

# ╔═╡ d27183d2-0d2b-4868-b34b-a90b2b5a0636
# you can change to different cells as well
cell_pos = (74, 53)

# ╔═╡ 0be4d313-6a6e-4eb6-bdc9-006ef6a726e6
# you can also change size
cell_size = (23, 23)

# ╔═╡ fe0632ef-e1b2-4a6f-ba53-44865b7971d8
cell = get_template(cells, cell_pos, cell_size);

# ╔═╡ 5f9744fa-a818-46bf-b83e-942ad1baac12
gray_show(cell)

# ╔═╡ cd861f28-8923-4ff7-ba24-6875cfa99f39
"""
	NCC(img, kernel)

Normalized cross-correlation of `img` with `kernel`.

## Procedure
* Determine the correct output size
* Move the kernel over the image like [here](https://kharshit.github.io/img/edge_detection.png).
* Calculate the local standard deivations and divide by them
"""
function NCC(img, kernel)
	out = similar(img, size(img) .- size(kernel) .+ 1)

	# runtime complexity is:
	for i in 1:size(out, 1)
		for j in 1:size(out, 2)
			# value we are going to write into out
			# that value accumulates the elementwise multiplication
			# and summation
			numerator = zero(eltype(img))
			denominator1 = zero(eltype(img))
			denominator2 = zero(eltype(img))

			for m in 1:size(kernel, 1)
				for n in 1:size(kernel, 2)
					# -1 because of Julia's one-based indexing style
					numerator += kernel[m, n] * img[i + m - 1, j + n - 1]
					denominator1 += img[i + m - 1, j + n - 1]^2
					denominator2 += kernel[m, n]^2
				end
			end
			out[i, j] = numerator / (√(denominator1) * √(denominator2))
		end
	end
	return out
end

# ╔═╡ 3fa412f2-d107-4caa-abf6-bf47c8333027
corr = NCC(cells, cell);

# ╔═╡ ff16b0cb-b8ba-4d3a-bcc8-da5c120ea898
gray_show(corr)

# ╔═╡ 0c108a6b-599a-4581-8d8d-68b9ae5865df
"""
	find_max(arr)

Returns the value and the index of the maximum intensity.

### Procedure
* Iterate over the array
* Update the current value if you find a larger one

### Examples
```julia-repl
julia> find_max([1 2 3; 4 5 6; 100 8 9])
(100, (3, 1))
```
"""
function find_max(arr)
	# TODO: initate inds and value meaningfully
	inds = (0, 0)
	value = eltype(arr)(-Inf)

	#for cind in CartesianIndices(arr)
	#	v = arr[cind]
	#	if v > value
	#		value = v
	#		inds = Tuple(cind)
	#	end
	#end

	for i in 1:size(arr, 1)
		for j in 1:size(arr, 2)
			v = arr[i, j]
			if v > value
				value = v
				inds = (i, j)
			end
		end
	end
	# TODO: return value and inds
	return value, inds
end

# ╔═╡ 2d7d48ac-c2e1-4acd-8e21-3a9314e19390
collect(eachindex(img))[2]

# ╔═╡ d1660777-a3a2-4abb-a394-4872da699c4e
Float64(-Inf)

# ╔═╡ 49a583a5-ef6b-4e5b-8ea2-0e9a5088b3cc
md"Now check whether `find_max` finds your template in the image again"

# ╔═╡ 64af243c-77f9-4c87-8f59-a4e506e3cb46
find_max(corr)

# ╔═╡ 5610563c-ea7f-40c6-ba8c-5a816bc71a2e
PlutoTest.@test find_max(corr)[2] == cell_pos

# ╔═╡ 37f43e35-015a-4596-96b0-13ccb720c767
md"## 2 Test"

# ╔═╡ bef2f5e8-e6ef-49df-9125-1d814feb31d7
kernel_ncc = [9 10; 13 14];

# ╔═╡ cf18b86d-6c82-4335-b46c-60cd161504b5
img_ncc = [Float64((i - 1) * 4 + j) for i = 1:4, j = 1:4];

# ╔═╡ 2cdd244d-33cb-487e-9fa2-069f8a845ef9
arr_rand = randn((32,32));

# ╔═╡ 8e5f5bf2-5b93-4c52-9960-ea1fb60888ce
PlutoTest.@test size(get_template(cells, (10, 10), (3, 7))) == (3,7)

# ╔═╡ 15626656-af92-44d3-83dd-9775172caed8
PlutoTest.@test get_template(cells, (10, 10), (3, 3)) ≈ Float32[0.15686274 0.18823528 0.31372547; 0.18823528 0.25098038 0.4078431; 0.25098038 0.34509802 0.43921566]

# ╔═╡ a9db1ffb-cc86-4d78-8f24-8b847443aae8
PlutoTest.@test Tuple.(find_max(arr_rand)) == Tuple.(findmax(arr_rand))

# ╔═╡ e144fe4f-cb8e-4158-a59a-8a93faca1f0c
PlutoTest.@test NCC(img_ncc, kernel_ncc) ≈ [0.9376736528408373 0.9683640522700839 0.9836212432229419; 0.9958743946641203 0.9981668178901743 0.9993408297073572; 0.9999999999999999 0.9999029998713893 0.9996660511437935]

# ╔═╡ b34bb590-0fd4-491f-a392-a538d16a86b6

# ╔═╡ 1fc3e23d-e283-43a1-a645-cfbaa3821870
md"# 3 Watershed
Based on the correlation map, we can use that as seed for a watershed algorithm.

What we need is an image with _seeds_. By that we mean an image full of 0s.
At the rough positions of the cells we put a 1.
"

# ╔═╡ 2fc20959-426a-410f-900b-bdde889426fd
"""
	get_seeds(cells, cell)

The `cells` you want to segment with the template `cell`.

# Procedure
* calculate the `NCC` arr (called `corr`)
* use `findlocalmaxima(corr)` to get the local peaks
* return an array which is 0 everywhere, except at the local peaks, where it is 1
"""
function get_seeds(cells, cell)
	# TODO: calc corr (call method from above)
	corr = NCC(cells, cell)

	# TODO: create an arr_n filled with zeros
	arr_n = zeros(eltype(corr), size(corr))
	# TODO: use findlocalmaxima
	maxs = findlocalmaxima(corr)

	# TODO: iterate over local maxima and set the corresponding
	# values of the arr_n to 1
	for cinds in maxs
		arr_n[cinds] = 1
	end

	return arr_n
end

# ╔═╡ e0f00df4-7c63-448c-b005-eda95cb4cdb8
findlocalmaxima(corr)

# ╔═╡ 7d95732e-94f4-448c-8e19-7969271059a9
10 * gray_show(get_seeds(cells, cell)) + gray_show(corr)

# ╔═╡ 86a8bf70-3c82-4100-8b63-5d010523a196
# we call watershed algorithm
segments = watershed(.- corr, label_components(get_seeds(cells, cell)));

# ╔═╡ 3641d445-1fff-416d-849b-f9d5715fef0f
heatmap(segments.image_indexmap .* (FourierTools.select_region(cells, new_size=size(segments.image_indexmap)) .> 0.3), c=:glasbey_hv_n256)

# ╔═╡ cc688fb0-0da2-4d5c-a0e8-bcff8af05ad1
gray_show(reverse(cells, dims=(1)))

# ╔═╡ e559b4b6-e26c-446e-b67e-b196e8ac3021
md"# 3 Test"

# ╔═╡ 6b966b4c-4b5c-4d62-a842-efe204b1c0fb
PlutoTest.@test get_seeds(Float32[1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16], [3 4; 5 6]) ≈ Float32[0.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.0]

# ╔═╡ 37c5f515-74a3-4bd0-b7e6-6e498b375e23

# ╔═╡ cfe15c1b-7b29-4514-82c4-e6ffe86456b0
md"# 4 NCC starts to fail
For which values does the NCC start to fail?
"

# ╔═╡ 814d4261-c433-47a9-a35a-16d83be3d4f4
"""
	draw(img, template, pos)

Draw the `template` into the `img`.
The top left corner aligns each with `pos`.

So we basically set the `template` into the `img`.
"""
function draw(img, template, pos)
	img[pos[1]:pos[1] + size(template, 1) - 1,
		pos[2]:pos[2] + size(template, 2) - 1] .= template
	return img
end

# ╔═╡ 3e30b081-e4df-4518-98aa-96e6e0880a4a
md"σ = $(@bind σ Slider(0:0.01:0.4, show_value=true))"

# ╔═╡ d45bcd94-6adf-40eb-baee-e8edab966f07
md"N = $(@bind N Slider(1:300, show_value=true))"

# ╔═╡ eb005f06-1b2f-45fe-b6b2-abfc6b4dbb1b
md"
First, we shear the image, then we add Poisson noise and finally Gauss noise.
"

# ╔═╡ a01fab74-b77b-4b6f-8415-b0f60530593b
cells_n = add_gauss(poisson(shear(cells, 20), N), σ);

# ╔═╡ de0edf52-cf26-4265-8f10-d8720dfee53f
corr_n = NCC(cells_n, cell);

# ╔═╡ 15dc3591-8141-4e8d-9778-8290de7627ad
[ gray_show(cells_n) gray_show(draw(cells_n, cell .* 3, find_max(corr_n)[2]))]

# ╔═╡ 5ef8d15b-b407-4f18-8ed3-b0d4c69b6345
gray_show(corr_n)

# ╔═╡ 5e32a3b3-30a9-4cfe-8005-b2eebfd5a4b4
md"Those values should roughly agree. Remember we sheared the array, so they won't exactly be the same"

# ╔═╡ b94655ba-382e-402e-a0c0-89d9d4d91bb6
find_max(corr_n)

# ╔═╡ 3ee0c65f-02b2-41b9-a95b-9932203223c8
find_max(corr)

# ╔═╡ 84a06e2e-13b9-44ea-a1e3-809791017656
md"## 4 Test"

# ╔═╡ 0d4763cc-f47a-4de0-9cde-12a73556b5e3
PlutoTest.@test draw(zeros((5, 5)), ones((2, 2)), (2, 2)) ≈ [0.0 0.0 0.0 0.0 0.0; 0.0 1.0 1.0 0.0 0.0; 0.0 1.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0]

# ╔═╡ 15b6196d-0b87-45fd-8303-770969386d10
md"# 5 Simulating a Scalar High NA PSF
We often discussed how you can simulate a simple PSF:
* start with an array full of 1s in Fourier space
* low pass filter it
* take the Fourier transform
* take `abs2`

However, that procedure suffers from hard edge artifacts in Fourier space.
Fortunately, in 2D we know that the solution of a perfect PSF is a

$\text{jinc}^2(x) = \left(\frac{J_1(x)}{x}\right)^2$

where $J_1$ is the Bessel function of the first kind.

We now want to simulate PSF but with an apodization function multiplied in Fourier space (e.g. the aplantic factor).

"

# ╔═╡ 75852dd4-6cc9-4040-84e2-cdb0ce812ec4
"""
	simulate_PSF(s, r, θ, [T=Float64])

Simulate the incoherent 2D PSF with radius `r` and total output size `s`.
This functions returns the PSF centered around the center pixel.
`θ` is the angle over the field size needed for the aplanatic factor.

The procedure is as following:
* create a disc with values of 1 and radius `r` (use `rr2`)
* multiply `cos.(θ)`
* take the `ift` of it.
* normalize the sum to 1
* center pixcel should be located in the center
"""
function simulate_PSF(s, r, θ, T=Float64)
	disc = T.(rr2(T, s) .< r^2)

	# a modification of the perfect pupil
	disc .*= cos.(θ)
	apsf = ift(disc)
	psf = abs2.(apsf)
	psf ./= sum(psf)

    return psf
end

# ╔═╡ bfa26eb6-d9ea-4a5d-80a5-84f048ed3ea0
"""
	simulate_PSF_jinc(s, scale, θ, [T=Float64])

`s` is the size. `scale` is the argument passed to `rr`.
`θ` is the angle over the field size needed for the aplanatic factor.
This functions returns the PSF centered around the center pixel.

## Procedure
* Call `jinc.` of the `rr` function which is scaled by `scale`
* Take `ft`.
* multiply aplanatic factor
* take `ift`.
* take `abs2.` and normalize sum
"""
function simulate_PSF_jinc(s, scale, θ, T=Float64)
	# ideal APSF
	apsf = jinc.(rr(s, scale=scale))

	# go to Fourier space and apply cos factor
	pupil = ft(apsf)
	pupil .*= cos.(θ)

	# go to real space again
	apsf_n = ift(pupil)

	# calculate psf
	psf = abs2.(apsf_n)
	psf ./= sum(psf)

	return psf
end

# ╔═╡ 23cc65b9-76d8-43c9-a398-22212b00c76a
θ = abs.(asin.(min.(1, rr(Float32, (100, 100), scale=0.05))));

# ╔═╡ 88b9c633-01ef-4bf2-9180-0bb2b59f4766
# todo
# What is the correct scaling such that both methods return a same scaled PSF
# as a hint, look at the cell below
# factor of 5 comes from: 10 / 100 * 2
# 10 is the radius of the PSF in Fourier space

scale = 1 / 5 # fix value

# ╔═╡ 170b5e91-dfaa-4f13-ac45-01b875ba02d8
# hint, look at this output
gray_show(abs.(ft(jinc.(scale .* rr((100, 100))))))

# ╔═╡ 19e2a53c-7716-4fd3-8716-87586c60ec73
gray_show(rr((100, 100)) .< 10)

# ╔═╡ 497f6e35-4b9a-4349-8366-9f54aa1476c5
md"### Can you spot differences?"

# ╔═╡ 71510434-8685-49ee-b514-31d1c7c411d8
psf_rr = simulate_PSF((100, 100), 10, θ);

# ╔═╡ 195bb3ce-5ade-4c31-9385-12cc4c119762
psf_jinc = simulate_PSF_jinc((100, 100), scale, θ);

# ╔═╡ 8d2c5280-d9c1-414d-a64c-8d8ef5a88b10
gray_show(log.(psf_rr), set_zero=true)

# ╔═╡ abf9ed24-ea46-4067-b786-3a6b23be4c80
gray_show(log.(psf_jinc), set_zero=true)

# ╔═╡ e3d0bb6e-e9c9-4a2d-97bf-73041958771b
md"## 5 Test"

# ╔═╡ f4b0f15f-016e-41ba-a999-e9f792a534fc
PlutoTest.@test simulate_PSF((3, 3), 10, 1) ≈ [0.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.0]

# ╔═╡ db517950-7031-4495-8198-f054eb833e16
PlutoTest.@test simulate_PSF_jinc((3, 3), 1000000, 1) ≈ [0.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.0]

# ╔═╡ 9a09fa1a-8045-4e93-9140-7f57915146ac
PlutoTest.@test simulate_PSF((10, 10), 3, 0.3) ≈ Float32[0.00040000005 0.0 0.0006111457 0.0 0.0041888547 0.01 0.0041888547 0.0 0.0006111457 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0006111457 0.0 0.0009337476 0.0 0.006400001 0.015278643 0.006400001 0.0 0.0009337476 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0041888547 0.0 0.006400001 0.0 0.043866254 0.10472137 0.043866254 0.0 0.006400001 0.0; 0.01 0.0 0.015278643 0.0 0.10472137 0.25 0.10472137 0.0 0.015278643 0.0; 0.0041888547 0.0 0.006400001 0.0 0.043866254 0.10472137 0.043866254 0.0 0.006400001 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0006111457 0.0 0.0009337476 0.0 0.006400001 0.015278643 0.006400001 0.0 0.0009337476 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]

# ╔═╡ ea3b65bb-2e41-4c17-9df8-35eaf69975f6
PlutoTest.@test simulate_PSF_jinc((4, 4), 0.1, 0.3, Float64) ≈ [0.05502252115995306 0.059372295904787806 0.060884063254034754 0.059372295904787806; 0.059372295904787806 0.0640036567471649 0.0656126102678531 0.0640036567471649; 0.06088406325403477 0.06561261026785312 0.06725510065275418 0.06561261026785312; 0.059372295904787806 0.0640036567471649 0.0656126102678531 0.0640036567471649]

# ╔═╡ Cell order:
# ╠═250252f0-6e61-11ec-2ffe-ebddaa9f0fb3
# ╠═0eeb5bb1-bb41-43cc-8c14-eb9c1ebf6fad
# ╠═1b2ebe1a-f3a7-494c-bcca-18df5403ee65
# ╟─c11d25d8-8fe9-476e-9a6d-66f625f29348
# ╠═6f87bb34-ae2b-4c8e-9c21-621644445df3
# ╟─b030ccf9-9a48-4410-8ac1-f3a2a02b5b5c
# ╠═012dd588-350b-4c2f-99c9-6371e443e32c
# ╠═a879425f-a152-4d4d-b6b8-f2a03b7a452f
# ╠═61ccfb00-513e-4bca-80eb-133aa58a1def
# ╠═d51fadc0-3227-4e80-96f2-269a8fd2bac2
# ╠═caab13b5-12a6-4058-8931-1f8a9abe57e8
# ╠═cc090365-274f-4e56-9be5-cbc1a32295ad
# ╟─b63022e3-3197-4dfe-9b7f-3ee8babce73a
# ╠═7ec68244-fe43-4418-b9ee-0cb7e67572c9
# ╟─2bd635be-6f07-40df-a25e-cfdf2dd22d6f
# ╠═3765b001-deb1-49b2-87a1-0100d36d3087
# ╠═57427d81-961c-488e-9a9d-e625a2e536c3
# ╟─6b4028de-e605-455f-947d-62649f2e26fa
# ╠═53ae1016-55d7-4545-a29d-b9a4fc90b854
# ╠═01428578-4258-413d-8b8b-27203527fd2b
# ╠═81605617-d091-47b2-bd85-1cd8ae0fef18
# ╠═f3fb21e1-990a-42ac-8a5d-52154be23cfb
# ╟─94694b2a-715a-450b-a024-fc474ce19975
# ╠═f43bd606-d202-4f74-8e46-28bab5af6d45
# ╠═4ed29eeb-0cbc-4810-8a6c-1b7aaa0014be
# ╠═1664be9b-08bc-4144-b191-b82bc9e2cee2
# ╠═d27183d2-0d2b-4868-b34b-a90b2b5a0636
# ╠═0be4d313-6a6e-4eb6-bdc9-006ef6a726e6
# ╠═fe0632ef-e1b2-4a6f-ba53-44865b7971d8
# ╠═5f9744fa-a818-46bf-b83e-942ad1baac12
# ╠═cd861f28-8923-4ff7-ba24-6875cfa99f39
# ╠═3fa412f2-d107-4caa-abf6-bf47c8333027
# ╠═ff16b0cb-b8ba-4d3a-bcc8-da5c120ea898
# ╠═0c108a6b-599a-4581-8d8d-68b9ae5865df
# ╠═2d7d48ac-c2e1-4acd-8e21-3a9314e19390
# ╠═d1660777-a3a2-4abb-a394-4872da699c4e
# ╟─49a583a5-ef6b-4e5b-8ea2-0e9a5088b3cc
# ╠═64af243c-77f9-4c87-8f59-a4e506e3cb46
# ╠═5610563c-ea7f-40c6-ba8c-5a816bc71a2e
# ╟─37f43e35-015a-4596-96b0-13ccb720c767
# ╠═bef2f5e8-e6ef-49df-9125-1d814feb31d7
# ╠═cf18b86d-6c82-4335-b46c-60cd161504b5
# ╠═2cdd244d-33cb-487e-9fa2-069f8a845ef9
# ╠═8e5f5bf2-5b93-4c52-9960-ea1fb60888ce
# ╠═15626656-af92-44d3-83dd-9775172caed8
# ╠═a9db1ffb-cc86-4d78-8f24-8b847443aae8
# ╠═e144fe4f-cb8e-4158-a59a-8a93faca1f0c
# ╠═b34bb590-0fd4-491f-a392-a538d16a86b6
# ╟─1fc3e23d-e283-43a1-a645-cfbaa3821870
# ╠═2fc20959-426a-410f-900b-bdde889426fd
# ╠═e0f00df4-7c63-448c-b005-eda95cb4cdb8
# ╠═7d95732e-94f4-448c-8e19-7969271059a9
# ╠═86a8bf70-3c82-4100-8b63-5d010523a196
# ╠═3641d445-1fff-416d-849b-f9d5715fef0f
# ╠═cc688fb0-0da2-4d5c-a0e8-bcff8af05ad1
# ╟─e559b4b6-e26c-446e-b67e-b196e8ac3021
# ╠═6b966b4c-4b5c-4d62-a842-efe204b1c0fb
# ╠═37c5f515-74a3-4bd0-b7e6-6e498b375e23
# ╟─cfe15c1b-7b29-4514-82c4-e6ffe86456b0
# ╠═814d4261-c433-47a9-a35a-16d83be3d4f4
# ╠═3e30b081-e4df-4518-98aa-96e6e0880a4a
# ╠═d45bcd94-6adf-40eb-baee-e8edab966f07
# ╟─eb005f06-1b2f-45fe-b6b2-abfc6b4dbb1b
# ╠═a01fab74-b77b-4b6f-8415-b0f60530593b
# ╠═15dc3591-8141-4e8d-9778-8290de7627ad
# ╠═5ef8d15b-b407-4f18-8ed3-b0d4c69b6345
# ╠═de0edf52-cf26-4265-8f10-d8720dfee53f
# ╟─5e32a3b3-30a9-4cfe-8005-b2eebfd5a4b4
# ╠═b94655ba-382e-402e-a0c0-89d9d4d91bb6
# ╠═3ee0c65f-02b2-41b9-a95b-9932203223c8
# ╟─84a06e2e-13b9-44ea-a1e3-809791017656
# ╠═0d4763cc-f47a-4de0-9cde-12a73556b5e3
# ╟─15b6196d-0b87-45fd-8303-770969386d10
# ╠═4aeceb65-f226-41cf-988b-273f42c75059
# ╠═75852dd4-6cc9-4040-84e2-cdb0ce812ec4
# ╠═bfa26eb6-d9ea-4a5d-80a5-84f048ed3ea0
# ╠═23cc65b9-76d8-43c9-a398-22212b00c76a
# ╠═88b9c633-01ef-4bf2-9180-0bb2b59f4766
# ╠═170b5e91-dfaa-4f13-ac45-01b875ba02d8
# ╠═19e2a53c-7716-4fd3-8716-87586c60ec73
# ╟─497f6e35-4b9a-4349-8366-9f54aa1476c5
# ╠═71510434-8685-49ee-b514-31d1c7c411d8
# ╠═195bb3ce-5ade-4c31-9385-12cc4c119762
# ╠═8d2c5280-d9c1-414d-a64c-8d8ef5a88b10
# ╠═abf9ed24-ea46-4067-b786-3a6b23be4c80
# ╟─e3d0bb6e-e9c9-4a2d-97bf-73041958771b
# ╠═f4b0f15f-016e-41ba-a999-e9f792a534fc
# ╠═db517950-7031-4495-8198-f054eb833e16
# ╠═9a09fa1a-8045-4e93-9140-7f57915146ac
# ╠═ea3b65bb-2e41-4c17-9df8-35eaf69975f6


text:
const EARTH_IMG = ImageMagick.load(joinpath(dirname(@__DIR__), "assets", "raster", "50-natural-earth-1-downsampled.png"))

const COASTLINES_LINEVEC = begin
    to_nansep_vec(coordinates.(geometry.(GeoJSON.read(Base.read(joinpath(dirname(@__DIR__), "assets", "vector", "110m_coastline.geojson"), String)).features))) do feature
        Point2f0.(feature)
    end
end


text:
istraining() = false

@adjoint istraining() = true, _ -> nothing

_isactive(m) = isnothing(m.active) ? istraining() : m.active

_dropout_shape(s, ::Colon) = size(s)
_dropout_shape(s, dims) = tuple((i ∉ dims ? 1 : si for (i, si) ∈ enumerate(size(s)))...)

_dropout_kernel(y::T, p, q) where {T} = y > p ? T(1 / q) : T(0)

"""
    dropout(x, p; dims=:, active=true)

The dropout function. If `active` is `true`,
for each input, either sets that input to `0` (with probability
`p`) or scales it by `1 / (1 - p)`. `dims` specifies the unbroadcasted dimensions,
e.g. `dims=1` applies dropout along columns and `dims=2` along rows.
This is used as a regularisation, i.e. it reduces overfitting during training.

If `active` is `false`, it just returns the input `x`

Warning: when using this function, you have to manually manage the activation
state. Usually in fact, dropout is used while training
but is deactivated in the inference phase. This can be
automatically managed using the [`Dropout`](@ref) layer instead of the
`dropout` function.

The [`Dropout`](@ref) layer is what you should use in most scenarios.
"""
function dropout(x, p; dims=:, active::Bool=true)
  active || return x
  y = dropout_mask(x, p, dims=dims)
  return x .* y
end

@adjoint function dropout(x, p; dims=:, active::Bool=true)
  active || return x, Δ -> (Δ, nothing)
  y = dropout_mask(x, p, dims=dims)
  return x .* y, Δ -> (Δ .* y, nothing)
end

function dropout_mask(x, p; dims=:)
  y = rand!(similar(x, _dropout_shape(x, dims)))
  y .= _dropout_kernel.(y, p, 1 - p)
  return y
end

"""
    Dropout(p, dims=:)

Dropout layer. In the forward pass, apply the [`Flux.dropout`](@ref) function on the input.

Does nothing to the input once [`Flux.testmode!`](@ref) is `true`.
"""
mutable struct Dropout{F,D}
  p::F
  dims::D
  active::Union{Bool, Nothing}
end

function Dropout(p; dims=:)
  @assert 0 ≤ p ≤ 1
  Dropout(p, dims, nothing)
end

function (a::Dropout)(x)
  _isactive(a) || return x
  return dropout(x, a.p; dims = a.dims, active=true)
end

testmode!(m::Dropout, mode=true) =
  (m.active = (isnothing(mode) || mode == :auto) ? nothing : !mode; m)

function Base.show(io::IO, d::Dropout)
  print(io, "Dropout(", d.p)
  d.dims != (:) && print(io, ", dims = $(repr(d.dims))")
  print(io, ")")
end

"""
    AlphaDropout(p)

A dropout layer. Used in
[Self-Normalizing Neural Networks](https://arxiv.org/abs/1706.02515).
The AlphaDropout layer ensures that mean and variance of activations
remain the same as before.

Does nothing to the input once [`testmode!`](@ref) is true.
"""
mutable struct AlphaDropout{F}
  p::F
  active::Union{Bool, Nothing}
  function AlphaDropout(p, active = nothing)
    @assert 0 ≤ p ≤ 1
    new{typeof(p)}(p, active)
  end
end

function (a::AlphaDropout)(x)
  _isactive(a) || return x
  λ = eltype(x)(1.0507009873554804934193349852946)
  α = eltype(x)(1.6732632423543772848170429916717)
  α1 = eltype(x)(-λ*α)
  noise = randn(eltype(x), size(x))
  x = @. x*(noise > (1 - a.p)) + α1 * (noise < (1 - a.p))
  A = (a.p + a.p * (1 - a.p) * α1 ^ 2)^0.5
  B = -A * α1 * (1 - a.p)
  x = @. A * x + B
  return x
end

testmode!(m::AlphaDropout, mode = true) =
  (m.active = (isnothing(mode) || mode == :auto) ? nothing : !mode; m)

"""
    LayerNorm(h::Integer)

A [normalisation layer](https://arxiv.org/abs/1607.06450) designed to be
used with recurrent hidden states of size `h`. Normalises the mean and standard
deviation of each input before applying a per-neuron gain/bias.
"""
struct LayerNorm{T}
  diag::Diagonal{T}
end

LayerNorm(h::Integer) =
  LayerNorm(Diagonal(h))

@functor LayerNorm

(a::LayerNorm)(x) = a.diag(normalise(x, dims=1))

function Base.show(io::IO, l::LayerNorm)
  print(io, "LayerNorm(", length(l.diag.α), ")")
end

"""
    BatchNorm(channels::Integer, σ = identity;
              initβ = zeros, initγ = ones,
              ϵ = 1e-8, momentum = .1)

[Batch Normalization](https://arxiv.org/abs/1502.03167) layer.
`channels` should be the size of the channel dimension in your data (see below).

Given an array with `N` dimensions, call the `N-1`th the channel dimension. (For
a batch of feature vectors this is just the data dimension, for `WHCN` images
it's the usual channel dimension.)

`BatchNorm` computes the mean and variance for each each `W×H×1×N` slice and
shifts them to have a new mean and variance (corresponding to the learnable,
per-channel `bias` and `scale` parameters).

Use [`testmode!`](@ref) during inference.

# Examples
```julia
m = Chain(
  Dense(28^2, 64),
  BatchNorm(64, relu),
  Dense(64, 10),
  BatchNorm(10),
  softmax)
```
"""
mutable struct BatchNorm{F,V,W,N}
  λ::F  # activation function
  β::V  # bias
  γ::V  # scale
  μ::W  # moving mean
  σ²::W  # moving std
  ϵ::N
  momentum::N
  active::Union{Bool, Nothing}
end

BatchNorm(chs::Integer, λ = identity;
          initβ = (i) -> zeros(Float32, i), initγ = (i) -> ones(Float32, i), ϵ = 1f-5, momentum = 0.1f0) =
  BatchNorm(λ, initβ(chs), initγ(chs),
            zeros(chs), ones(chs), ϵ, momentum, nothing)

trainable(bn::BatchNorm) = (bn.β, bn.γ)

function (BN::BatchNorm)(x)
  size(x, ndims(x)-1) == length(BN.β) ||
    error("BatchNorm expected $(length(BN.β)) channels, got $(size(x, ndims(x)-1))")
  dims = length(size(x))
  channels = size(x, dims-1)
  affine_shape = ntuple(i->i == ndims(x) - 1 ? size(x, i) : 1, ndims(x))
  m = div(prod(size(x)), channels)
  γ = reshape(BN.γ, affine_shape...)
  β = reshape(BN.β, affine_shape...)
  if !_isactive(BN)
    μ = reshape(BN.μ, affine_shape...)
    σ² = reshape(BN.σ², affine_shape...)
    ϵ = BN.ϵ
  else
    T = eltype(x)
    axes = [1:dims-2; dims] # axes to reduce along (all but channels axis)
    μ = mean(x, dims = axes)
    σ² = sum((x .- μ) .^ 2, dims = axes) ./ m
    ϵ = convert(T, BN.ϵ)
    # update moving mean/std
    mtm = BN.momentum
    S = eltype(BN.μ)
    BN.μ  = (1 - mtm) .* BN.μ .+ mtm .* S.(reshape(μ, :))
    BN.σ² = (1 - mtm) .* BN.σ² .+ (mtm * m / (m - 1)) .* S.(reshape(σ², :))
  end

  let λ = BN.λ
    x̂ = (x .- μ) ./ sqrt.(σ² .+ ϵ)
    λ.(γ .* x̂ .+ β)
  end
end

@functor BatchNorm

testmode!(m::BatchNorm, mode = true) =
  (m.active = (isnothing(mode) || mode == :auto) ? nothing : !mode; m)

function Base.show(io::IO, l::BatchNorm)
  print(io, "BatchNorm($(join(size(l.β), ", "))")
  (l.λ == identity) || print(io, ", λ = $(l.λ)")
  print(io, ")")
end

expand_inst = (x, as) -> reshape(repeat(x, outer=[1, as[length(as)]]), as...)

mutable struct InstanceNorm{F,V,W,N}
  λ::F  # activation function
  β::V  # bias
  γ::V  # scale
  μ::W  # moving mean
  σ²::W  # moving std
  ϵ::N
  momentum::N
  active::Union{Bool, Nothing}
end

"""
    InstanceNorm(channels::Integer, σ = identity;
                 initβ = zeros, initγ = ones,
                 ϵ = 1e-8, momentum = .1)

[Instance Normalization](https://arxiv.org/abs/1607.08022) layer.
`channels` should be the size of the channel dimension in your data (see below).

Given an array with `N` dimensions, call the `N-1`th the channel dimension. (For
a batch of feature vectors this is just the data dimension, for `WHCN` images
it's the usual channel dimension.)

`InstanceNorm` computes the mean and variance for each each `W×H×1×1` slice and
shifts them to have a new mean and variance (corresponding to the learnable,
per-channel `bias` and `scale` parameters).

Use [`testmode!`](@ref) during inference.

# Examples
```julia
m = Chain(
  Dense(28^2, 64),
  InstanceNorm(64, relu),
  Dense(64, 10),
  InstanceNorm(10),
  softmax)
```
"""
InstanceNorm(chs::Integer, λ = identity;
          initβ = (i) -> zeros(Float32, i), initγ = (i) -> ones(Float32, i), ϵ = 1f-5, momentum = 0.1f0) =
  InstanceNorm(λ, initβ(chs), initγ(chs),
            zeros(chs), ones(chs), ϵ, momentum, nothing)

trainable(in::InstanceNorm) = (in.β, in.γ)

function (in::InstanceNorm)(x)
  size(x, ndims(x)-1) == length(in.β) ||
    error("InstanceNorm expected $(length(in.β)) channels, got $(size(x, ndims(x)-1))")
  ndims(x) > 2 ||
    error("InstanceNorm requires at least 3 dimensions. With 2 dimensions an array of zeros would be returned")
  # these are repeated later on depending on the batch size
  dims = length(size(x))
  c = size(x, dims-1)
  bs = size(x, dims)
  affine_shape = ntuple(i->i == ndims(x) - 1 || i == ndims(x) ? size(x, i) : 1, ndims(x))
  m = div(prod(size(x)), c*bs)
  γ, β = expand_inst(in.γ, affine_shape), expand_inst(in.β, affine_shape)

  if !_isactive(in)
    μ = expand_inst(in.μ, affine_shape)
    σ² = expand_inst(in.σ², affine_shape)
    ϵ = in.ϵ
  else
    T = eltype(x)

    ϵ = convert(T, in.ϵ)
    axes = 1:dims-2 # axes to reduce along (all but channels and batch size axes)
    μ = mean(x, dims = axes)
    σ² = mean((x .- μ) .^ 2, dims = axes)
    S = eltype(in.μ)
    # update moving mean/std
    mtm = in.momentum
    in.μ = dropdims(mean(repeat((1 - mtm) .* in.μ, outer=[1, bs]) .+ mtm .* S.(reshape(μ, (c, bs))), dims = 2), dims=2)
    in.σ² = dropdims(mean((repeat((1 - mtm) .* in.σ², outer=[1, bs]) .+ (mtm * m / (m - 1)) .* S.(reshape(σ², (c, bs)))), dims = 2), dims=2)
  end

  let λ = in.λ
    x̂ = (x .- μ) ./ sqrt.(σ² .+ ϵ)
    λ.(γ .* x̂ .+ β)
  end
end

@functor InstanceNorm

testmode!(m::InstanceNorm, mode = true) =
  (m.active = (isnothing(mode) || mode == :auto) ? nothing : !mode; m)

function Base.show(io::IO, l::InstanceNorm)
  print(io, "InstanceNorm($(join(size(l.β), ", "))")
  (l.λ == identity) || print(io, ", λ = $(l.λ)")
  print(io, ")")
end

"""
    GroupNorm(chs::Integer, G::Integer, λ = identity;
              initβ = (i) -> zeros(Float32, i), initγ = (i) -> ones(Float32, i),
              ϵ = 1f-5, momentum = 0.1f0)

[Group Normalization](https://arxiv.org/abs/1803.08494) layer.
This layer can outperform Batch Normalization and Instance Normalization.

`chs` is the number of channels, the channel dimension of your input.
For an array of N dimensions, the `N-1`th index is the channel dimension.

`G` is the number of groups along which the statistics are computed.
The number of channels must be an integer multiple of the number of groups.

Use [`testmode!`](@ref) during inference.

# Examples
```julia
m = Chain(Conv((3,3), 1=>32, leakyrelu;pad = 1),
          GroupNorm(32,16))
          # 32 channels, 16 groups (G = 16), thus 2 channels per group used
```
"""
mutable struct GroupNorm{F,V,W,N,T}
  G::T # number of groups
  λ::F  # activation function
  β::V  # bias
  γ::V  # scale
  μ::W  # moving mean
  σ²::W  # moving std
  ϵ::N
  momentum::N
  active::Union{Bool, Nothing}
end

GroupNorm(chs::Integer, G::Integer, λ = identity;
          initβ = (i) -> zeros(Float32, i), initγ = (i) -> ones(Float32, i), ϵ = 1f-5, momentum = 0.1f0) =
  GroupNorm(G, λ, initβ(chs), initγ(chs),
            zeros(G,1), ones(G,1), ϵ, momentum, nothing)

trainable(gn::GroupNorm) = (gn.β, gn.γ)

function(gn::GroupNorm)(x)
  size(x,ndims(x)-1) == length(gn.β) || error("Group Norm expected $(length(gn.β)) channels, but got $(size(x,ndims(x)-1)) channels")
  ndims(x) > 2 || error("Need to pass at least 3 channels for Group Norm to work")
  (size(x,ndims(x) -1))%gn.G == 0 || error("The number of groups ($(gn.G)) must divide the number of channels ($(size(x,ndims(x) -1)))")

  dims = length(size(x))
  groups = gn.G
  channels = size(x, dims-1)
  batches = size(x,dims)
  channels_per_group = div(channels,groups)
  affine_shape = ntuple(i->i == ndims(x) - 1 ? size(x, i) : 1, ndims(x))

  # Output reshaped to (W,H...,C/G,G,N)
  μ_affine_shape = ntuple(i->i == ndims(x) ? groups : 1, ndims(x) + 1)

  m = prod(size(x)[1:end-2]) * channels_per_group
  γ = reshape(gn.γ, affine_shape...)
  β = reshape(gn.β, affine_shape...)

  y = reshape(x,((size(x))[1:end-2]...,channels_per_group,groups,batches))
  if !_isactive(gn)
    og_shape = size(x)
    μ = reshape(gn.μ, μ_affine_shape...) # Shape : (1,1,...C/G,G,1)
    σ² = reshape(gn.σ², μ_affine_shape...) # Shape : (1,1,...C/G,G,1)
    ϵ = gn.ϵ
  else
    T = eltype(x)
    og_shape = size(x)
    axes = [(1:ndims(y)-2)...] # axes to reduce along (all but channels axis)
    μ = mean(y, dims = axes)
    σ² = mean((y .- μ) .^ 2, dims = axes)

    ϵ = convert(T, gn.ϵ)
    # update moving mean/std
    mtm = gn.momentum
    S = eltype(gn.μ)
    gn.μ = mean((1 - mtm) .* gn.μ .+ mtm .* S.(reshape(μ, (groups,batches))),dims=2)
    gn.σ² = mean((1 - mtm) .* gn.σ² .+ (mtm * m / (m - 1)) .* S.(reshape(σ², (groups,batches))),dims=2)
  end

  let λ = gn.λ
    x̂ = (y .- μ) ./ sqrt.(σ² .+ ϵ)

    # Reshape x̂
    x̂ = reshape(x̂,og_shape)
    λ.(γ .* x̂ .+ β)
  end
end

@functor GroupNorm

testmode!(m::GroupNorm, mode = true) =
  (m.active = (isnothing(mode) || mode == :auto) ? nothing : !mode; m)

function Base.show(io::IO, l::GroupNorm)
  print(io, "GroupNorm($(join(size(l.β), ", "))")
  (l.λ == identity) || print(io, ", λ = $(l.λ)")
  print(io, ")")
end


text: # map-reduction


text:
using Packing, GeometryBasics
using Test

# write your own tests here
include("test_rectangle.jl")
include("guillotine.jl")


text:
struct OLSUCBGreedy <: OLSUCBOptimisationAlgorithm end

function optimise_linear_sqrtquadratic(instance::CombinatorialInstance{T}, ::OLSUCBGreedy,
                                       linear::Dict{T, Float64}, sqrtquadratic::Dict{Tuple{T, T}, Float64},
                                       ε::Float64, d::Int, verbose::Int; with_trace::Bool=false) where T
  directions = Set(keys(linear))
  sol = T[]

  # As long as possible, add new arms to the solution to improve the expected reward.
  t0 = time_ns()
  n_iter = 0
  while true
    n_iter += 1

    # Pick the best direction to improve the reward.
    best_reward = 0.0
    best_direction = nothing
    for d in directions
      new_sol = copy(sol)
      push!(new_sol, d)
      if ! is_partially_acceptable(instance, new_sol)
        continue
      end

      reward = sum(linear[arm] for arm in keys(linear) if arm in new_sol) +
               sqrt(sum(sqrtquadratic[arms] for arms in keys(sqrtquadratic) if arms[1] in new_sol && arms[2] in new_sol))

      if reward > best_reward
        best_reward = reward
        best_direction = d
      end
    end

    # If no direction could be found, done!
    if best_direction == nothing
      break
    end

    # Update the solution and ensure this direction will never be tried again.
    push!(sol, best_direction)
    delete!(directions, best_direction)

    # If there is no more a direction to explore, done!
    if isempty(directions)
      break
    end
  end
  t1 = time_ns()

  if with_trace
    run_details = ESCB2Details()
    run_details.nIterations = n_iter
    run_details.bestLinearObjective = sum(linear[arm] + sqrtlinear[arm] for arm in keys(linear) if in(arm, sol))
    run_details.bestNonlinearObjective = sum(linear[arm] for arm in keys(linear) if in(arm, new_sol)) +
                                        sqrt(sum(sqrtquadratic[arms] for arms in keys(sqrtquadratic) if in(arms[1], new_sol) && in(arms[2], new_sol)))
    run_details.solverTimes = Float64[(t1 - t0) / 1_000_000_000]
    return sol, run_details
  else
    return sol
  end
end


text:
export LayerState

export HDF5DataLayer, MemoryDataLayer, AsyncHDF5DataLayer
export InnerProductLayer, ConvolutionLayer, PoolingLayer, SoftmaxLayer
export PowerLayer, SplitLayer, ElementWiseLayer, ChannelPoolingLayer
export LRNLayer, DropoutLayer, ReshapeLayer, ArgmaxLayer, HDF5OutputLayer
export CropLayer, ConcatLayer, RandomMaskLayer, TiedInnerProductLayer
export IdentityLayer, Index2OnehotLayer, MemoryOutputLayer
export SquareLossLayer, SoftmaxLossLayer, MultinomialLogisticLossLayer
export SoftlabelSoftmaxLossLayer, WassersteinLossLayer, HingeLossLayer
export AccuracyLayer, BinaryAccuracyLayer, BinaryCrossEntropyLossLayer

export RandomNormalLayer
export GaussianKLLossLayer
export setup, forward, backward, shutdown

export reset_statistics
export freeze!, unfreeze!, is_frozen

############################################################
# Implementing a Layer
#
# A Layer object is a configuration of how a layer is going
# to behave. The following fields are needed by the neural
# network engine:
#
# - tops: An array of strings, as the name of output blobs.
#       Note tops is optional for LossLayer and StatLayer.
# - bottoms: An array of strings for the name of input blobs.
#       Data Layers can omit this field.
#
# To add functional operation for the layer, a layer state
# type should be defined as a sub-type of LayerState. A
# LayerState object stores all the states and intermediate
# computation results of each layer. The following fields
# are needed by the neural network engine:
#
# - blobs: An array of blobs. Corresponding to the output
#       of this layer in the forward pass.
#
# If the layer needs back propagated gradient from upper
# layer, then the following fields are required:
#
# - blobs_diff: An array of blobs. Corresponding to the
#       gradient of the objective function with respect
#       to the output of this layer in the backward pass.
#       Note this value is computed by the upper layer.
#
# Optional each layer can have an activation function
# by defining a field called neuron. The activation
# function computation (both forward and backward)
# are handled automatically by the engine.
#
# If the layer has its own parameters that need to be updated
# during optimization, the following field should be defined
#
# - parameters: vector of Parameters
#
# Note the layer needs to compute the gradient of the parameter
# during the backward pass, but everything else will be taken
# care of automatically by the engine, including the followings:
#
# - initialize the parameters before optimization starts
# - compute the regularization and its gradient during forward
#   and backward pass, respectively
# - update the parameter after one forward-backward pass
#
# Then the following functions need to be defined
#
# - setup(backend::Backend,layer::MyLayer,inputs::Vector{Blob})
#   This function construct the layer state object and do
#   necessary initialization. The inputs are initialized
#   with proper shape, but not necessarily with valid data
#   values. The constructed layer state should be returned.
#
# - forward(backend::Backend,state::MyLayerState,inputs::Vector{Blob})
#   This function do the forward computation: inputs are
#   forward computed output from bottom layers.
#
# - backward(backend::Backend,state::MyLayerState,inputs::Vector{Blob},diffs::Vector{Blob})
#   This function do the backward computation: inputs are
#   the same as in forward, diffs contains blobs to hold
#   gradient with respect to the bottom layer input. Some
#   blob in this vector might be "undefined", meaning that
#   blob do not want to get back propagated gradient. This
#   procedure also compute gradient with respect to layer
#   parameters if any.
############################################################

############################################################
# Default layer characterizations
############################################################
@characterize_layer(Layer,
  is_source  => false, # data layer, takes no bottom blobs
  is_sink    => false, # top layer, produces no top blobs (loss, accuracy, etc.)
  has_param  => false, # contains trainable parameters
  has_neuron => false, # has a neuron
  can_do_bp  => false, # can do back-propagation
  is_inplace => false, # do inplace computation, does not has own top blobs
  has_loss   => false, # produce a loss
  has_stats  => false, # produce statistics
)

function setup(backend::Backend, layer::Layer, inputs::Vector{Blob}, diffs::Vector{Blob})
  error("Not implemented, should setup layer state")
end

# Default overload when there is no shared_parameters
function setup(backend::Backend, layer::Layer, shared_parameters, inputs::Vector{Blob}, diffs::Vector{Blob})
  setup(backend, layer, inputs, diffs)
end

function shutdown(backend::Backend, state::LayerState)
  error("Not implemented, please define an empty function explicitly if not needed")
end

function forward(backend::Backend, state::LayerState, inputs::Vector{Blob})
  error("Not implemented, should do forward computing")
end

function backward(backend::Backend, state::LayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  error("Not implemented, please define an empty function explicitly if not needed")
end

function param_key(layer::Layer)
  @assert has_param(layer)
  key = layer.param_key
  return isempty(key) ? layer.name : key
end

# If a layer state is frozen, its parameters will not be trained, also there is no
# need to compute the gradients for the parameters
function freeze!(state::LayerState)
  @assert !has_param(state.layer) "Layers with parameters should implement their own freeze function"
  # freeze has no effects for layers without parameters
end
function unfreeze!(state::LayerState)
  @assert !has_param(state.layer) "Layers with parameters should implement their own unfreeze function"
end
function is_frozen(state::LayerState)
  @assert !has_param(state.layer) "Layers with parameters should implement their own is_frozen function"
  false # layers without parameters are never frozen
end

#############################################################
# Display layers
#############################################################
import Base.show
export show, show_layer

function show(io::IO, layer::Layer)
  print(io, "$(typeof(layer))($(layer.name))")
end
function show_layer_details(io::IO, state::LayerState)
  # overload this function to add layer details
end
function show_layer(io::IO, state::LayerState, inputs::Vector{Blob})
  println(io, "............................................................")
  layer_title = " *** $(state.layer)"
  if isa(io, Base.TTY)
    print_with_color(:blue, io, layer_title)
  else
    print(io, layer_title)
  end
  println(io)

  show_layer_details(io, state)
  if !is_source(state.layer)
    println(io, "    Inputs ----------------------------")
    for i = 1:length(inputs)
      println(io, "     $(@sprintf("%9s", state.layer.bottoms[i])): $(inputs[i])")
    end
  end
  if !is_sink(state.layer) && !is_inplace(state.layer)
    println(io, "    Outputs ---------------------------")
    for i = 1:length(state.blobs)
      println(io, "     $(@sprintf("%9s", state.layer.tops[i])): $(state.blobs[i])")
    end
  end
end

#############################################################
# Data Layers
#############################################################
include("layers/hdf5-data.jl")
include("layers/memory-data.jl")
include("layers/async-hdf5-data.jl")

#############################################################
# General Computation Layers
#############################################################
include("layers/inner-product.jl")
include("layers/convolution.jl")
include("layers/pooling.jl")
include("layers/softmax.jl")
include("layers/power.jl")
include("layers/split.jl")
include("layers/element-wise.jl")
include("layers/channel-pooling.jl")
include("layers/lrn.jl")
include("layers/dropout.jl")
include("layers/reshape.jl")
include("layers/crop.jl")
include("layers/concat.jl")
include("layers/random-mask.jl")
include("layers/tied-inner-product.jl")
include("layers/identity.jl")
include("layers/random-normal.jl")

#############################################################
# Utility layers
#############################################################
include("layers/argmax.jl")
include("layers/hdf5-output.jl")
include("layers/index2onehot.jl")
include("layers/memory-output.jl")

#############################################################
# Loss Layers
#############################################################
include("layers/square-loss.jl")
include("layers/multinomial-logistic-loss.jl")
include("layers/softmax-loss.jl")
include("layers/softlabel-softmax-loss.jl")
include("layers/wasserstein-loss.jl")
include("layers/binary-cross-entropy-loss.jl")
include("layers/gaussian-kl-loss.jl")
include("layers/hinge-loss.jl")

#############################################################
# Statistics Layers
#############################################################
include("layers/accuracy.jl")
include("layers/binary-accuracy.jl")


text:
using StaticArrays

const triangle1an = [
    SVector{3}([0.333333333333333 0.333333333333333 0.333333333333333])
]

const triangle1wn = [
  1.0,
]

const triangle3an = [
    SVector{3}([0.666666666666667 0.166666666666667 0.166666666666667]),
    SVector{3}([0.166666666666667 0.666666666666667 0.166666666666667]),
    SVector{3}([0.166666666666667 0.166666666666667 0.666666666666667])
]

const triangle3wn = [
    0.333333333333333,
    0.333333333333333,
    0.333333333333333
]

const triangle6an = [
    SVector{3}([0.816847572980459 0.091576213509771 0.091576213509771]),
    SVector{3}([0.091576213509771 0.816847572980459 0.091576213509771]),
    SVector{3}([0.091576213509771 0.091576213509771 0.816847572980459]),
    SVector{3}([0.445948490915965 0.445948490915965 0.108103018168071]),
    SVector{3}([0.445948490915965 0.108103018168071 0.445948490915965]),
    SVector{3}([0.108103018168071 0.445948490915965 0.445948490915965])
]

const triangle6wn = [
    0.109951743655322,
    0.109951743655322,
    0.109951743655322,
    0.223381589678012,
    0.223381589678012,
    0.223381589678012
]

const triangle10an = [
    SVector{3}([0.3333333333333333 0.3333333333333333 0.333333333333333]),
    SVector{3}([0.0762564021238483 0.0762564021238483 0.847487195752303]),
    SVector{3}([0.0762564021238483 0.847487195752303 0.0762564021238483]),
    SVector{3}([0.847487195752303 0.0762564021238483 0.0762564021238483]),
    SVector{3}([0.322653005725289 0.0709487609446569 0.606398233330054] ),
    SVector{3}([0.322653005725289 0.606398233330054 0.0709487609446569] ),
    SVector{3}([0.0709487609446569 0.322653005725289 0.606398233330054] ),
    SVector{3}([0.0709487609446569 0.606398233330054 0.322653005725289] ),
    SVector{3}([0.606398233330054 0.322653005725289 0.0709487609446569] ),
    SVector{3}([0.606398233330054 0.0709487609446569 0.322653005725289] )
]

const triangle10wn = [
    0.19937905503631695,
    0.06819071348353711,
    0.06819071348353711,
    0.06819071348353711,
    0.09934146741884528,
    0.09934146741884528,
    0.09934146741884528,
    0.09934146741884528,
    0.09934146741884528,
    0.09934146741884528
]

#x = [0.06302787662265326, 0.24063968500834432, 0.4664324392487761, 0.29444203981813605, 0.0443521797062553]
#w = [0.049548433357090374, 0.1282255662128493, 0.03160255077585674, 0.06197839149376846]
const triangle15an = [
    SVector{3}([0.0630278766226533 0.0630278766226533 0.873944246754694]),
    SVector{3}([0.0630278766226533 0.873944246754694 0.0630278766226533]),
    SVector{3}([0.873944246754694 0.0630278766226533 0.0630278766226533]),
    SVector{3}([0.240639685008344 0.240639685008344 0.518720629983311]  ),
    SVector{3}([0.240639685008344 0.518720629983311 0.240639685008344]  ),
    SVector{3}([0.518720629983311 0.240639685008344 0.240639685008344]  ),
    SVector{3}([0.466432439248776 0.466432439248776 0.0671351215024478] ),
    SVector{3}([0.466432439248776 0.0671351215024478 0.466432439248776] ),
    SVector{3}([0.0671351215024478 0.466432439248776 0.466432439248776] ),
    SVector{3}([0.294442039818136 0.0443521797062553 0.661205780475609] ),
    SVector{3}([0.294442039818136 0.661205780475609 0.0443521797062553] ),
    SVector{3}([0.0443521797062553 0.294442039818136 0.661205780475609] ),
    SVector{3}([0.0443521797062553 0.661205780475609 0.294442039818136] ),
    SVector{3}([0.661205780475609 0.294442039818136 0.0443521797062553] ),
    SVector{3}([0.661205780475609 0.0443521797062553 0.294442039818136] )
]

const triangle15wn = [
    0.049548433357090374,
    0.049548433357090374,
    0.049548433357090374,
    0.1282255662128493,
    0.1282255662128493,
    0.1282255662128493,
    0.03160255077585674,
    0.03160255077585674,
    0.03160255077585674,
    0.06197839149376846,
    0.06197839149376846,
    0.06197839149376846,
    0.06197839149376846,
    0.06197839149376846,
    0.06197839149376846
]

#x = [0.045189009784135412, 0.40360397981780544, 0.48151983478264854, 0.22206316553736064, 0.030424361731352323, 0.1369912012735261, 0.21829007096767486]
#w = [0.025993571032099375, 0.094080073456857274, 0.051617202570398299, 0.035351705092115611, 0.045469538044873581]
const triangle21an = [
    SVector{3}([0.0451890097841354 0.0451890097841354 0.909621980431729]),
    SVector{3}([0.0451890097841354 0.909621980431729 0.0451890097841354]),
    SVector{3}([0.909621980431729 0.0451890097841354 0.0451890097841354]),
    SVector{3}([0.403603979817805 0.403603979817805 0.192792040364389]  ),
    SVector{3}([0.403603979817805 0.192792040364389 0.403603979817805]  ),
    SVector{3}([0.192792040364389 0.403603979817805 0.403603979817805]  ),
    SVector{3}([0.481519834782649 0.481519834782649 0.0369603304347029] ),
    SVector{3}([0.481519834782649 0.0369603304347029 0.481519834782649] ),
    SVector{3}([0.0369603304347029 0.481519834782649 0.481519834782649] ),
    SVector{3}([0.222063165537361 0.0304243617313523 0.747512472731287] ),
    SVector{3}([0.222063165537361 0.747512472731287 0.0304243617313523] ),
    SVector{3}([0.0304243617313523 0.222063165537361 0.747512472731287] ),
    SVector{3}([0.0304243617313523 0.747512472731287 0.222063165537361] ),
    SVector{3}([0.747512472731287 0.222063165537361 0.0304243617313523] ),
    SVector{3}([0.747512472731287 0.0304243617313523 0.222063165537361] ),
    SVector{3}([0.136991201273526 0.218290070967675 0.644718727758799]  ),
    SVector{3}([0.136991201273526 0.644718727758799 0.218290070967675]  ),
    SVector{3}([0.218290070967675 0.136991201273526 0.644718727758799]  ),
    SVector{3}([0.218290070967675 0.644718727758799 0.136991201273526]  ),
    SVector{3}([0.644718727758799 0.136991201273526 0.218290070967675]  ),
    SVector{3}([0.644718727758799 0.218290070967675 0.136991201273526]  )
]

const triangle21wn = [
    0.025993571032099375,
    0.025993571032099375,
    0.025993571032099375,
    0.09408007345685727,
    0.09408007345685727,
    0.09408007345685727,
    0.0516172025703983,
    0.0516172025703983,
    0.0516172025703983,
    0.03535170509211561,
    0.03535170509211561,
    0.03535170509211561,
    0.03535170509211561,
    0.03535170509211561,
    0.03535170509211561,
    0.04546953804487358,
    0.04546953804487358,
    0.04546953804487358,
    0.04546953804487358,
    0.04546953804487358,
    0.04546953804487358
]

#x = [0.032653640853443702, 0.14196281952419448, 0.28318116029651036, 0.02759039318927288, 0.16543435604751824, 0.60187823586421474, 0.026774393775283171, 0.13185966517985423, 0.32750749987987121]
#w = [ 0.00089863287010520762, 0.013780183679627252, 0.047407203691024753, 0.047279917201778536, 0.024833733334694335, 0.030636740264071611, 0.056812768636667919]
const triangle28an = [
    SVector{3}([0.333333333333333 0.333333333333333 0.333333333333333]    ),
    SVector{3}([0.0326536408534437 0.0326536408534437 0.934692718293113]),
    SVector{3}([0.0326536408534437 0.934692718293113 0.0326536408534437]),
    SVector{3}([0.934692718293113 0.0326536408534437 0.0326536408534437]),
    SVector{3}([0.141962819524194 0.141962819524194 0.716074360951611]  ),
    SVector{3}([0.141962819524194 0.716074360951611 0.141962819524194]  ),
    SVector{3}([0.716074360951611 0.141962819524194 0.141962819524194]  ),
    SVector{3}([0.283181160296510 0.283181160296510 0.433637679406979]  ),
    SVector{3}([0.283181160296510 0.433637679406979 0.283181160296510]  ),
    SVector{3}([0.433637679406979 0.283181160296510 0.283181160296510]  ),
    SVector{3}([0.0275903931892729 0.165434356047518 0.806975250763209] ),
    SVector{3}([0.0275903931892729 0.806975250763209 0.165434356047518] ),
    SVector{3}([0.165434356047518 0.0275903931892729 0.806975250763209] ),
    SVector{3}([0.165434356047518 0.806975250763209 0.0275903931892729] ),
    SVector{3}([0.806975250763209 0.0275903931892729 0.165434356047518] ),
    SVector{3}([0.806975250763209 0.165434356047518 0.0275903931892729] ),
    SVector{3}([0.601878235864215 0.0267743937752832 0.371347370360502] ),
    SVector{3}([0.601878235864215 0.371347370360502 0.0267743937752832] ),
    SVector{3}([0.0267743937752832 0.601878235864215 0.371347370360502] ),
    SVector{3}([0.0267743937752832 0.371347370360502 0.601878235864215] ),
    SVector{3}([0.371347370360502 0.601878235864215 0.0267743937752832] ),
    SVector{3}([0.371347370360502 0.0267743937752832 0.601878235864215] ),
    SVector{3}([0.131859665179854 0.327507499879871 0.540632834940274]  ),
    SVector{3}([0.131859665179854 0.540632834940274 0.327507499879871]  ),
    SVector{3}([0.327507499879871 0.131859665179854 0.540632834940274]  ),
    SVector{3}([0.327507499879871 0.540632834940274 0.131859665179854]  ),
    SVector{3}([0.540632834940274 0.131859665179854 0.327507499879871]  ),
    SVector{3}([0.540632834940274 0.327507499879871 0.131859665179854]  ),
]

const triangle28wn = [
    0.0008986328701052076,
    0.013780183679627252,
    0.013780183679627252,
    0.013780183679627252,
    0.04740720369102475,
    0.04740720369102475,
    0.04740720369102475,
    0.047279917201778536,
    0.047279917201778536,
    0.047279917201778536,
    0.024833733334694335,
    0.024833733334694335,
    0.024833733334694335,
    0.024833733334694335,
    0.024833733334694335,
    0.024833733334694335,
    0.03063674026407161,
    0.03063674026407161,
    0.03063674026407161,
    0.03063674026407161,
    0.03063674026407161,
    0.03063674026407161,
    0.05681276863666792,
    0.05681276863666792,
    0.05681276863666792,
    0.05681276863666792,
    0.05681276863666792,
    0.05681276863666792
]

const triangleAn = [
    triangle1an,
    triangle3an,
    triangle6an,
    triangle10an,
    triangle15an,
    triangle21an,
    triangle28an

]

const triangleWn = [
    triangle1wn,
    triangle3wn,
    triangle6wn,
    triangle10wn,
    triangle15wn,
    triangle21wn,
    triangle28wn

]


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: sagemaker_a2i_runtime
using AWS.Compat
using AWS.UUIDs

"""
    delete_human_loop(human_loop_name)
    delete_human_loop(human_loop_name, params::Dict{String,<:Any})

Deletes the specified human loop for a flow definition. If the human loop was deleted, this
operation will return a ResourceNotFoundException.

# Arguments
- `human_loop_name`: The name of the human loop that you want to delete.

"""
function delete_human_loop(HumanLoopName; aws_config::AbstractAWSConfig=global_aws_config())
    return sagemaker_a2i_runtime(
        "DELETE",
        "/human-loops/$(HumanLoopName)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_human_loop(
    HumanLoopName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sagemaker_a2i_runtime(
        "DELETE",
        "/human-loops/$(HumanLoopName)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_human_loop(human_loop_name)
    describe_human_loop(human_loop_name, params::Dict{String,<:Any})

Returns information about the specified human loop. If the human loop was deleted, this
operation will return a ResourceNotFoundException error.

# Arguments
- `human_loop_name`: The name of the human loop that you want information about.

"""
function describe_human_loop(
    HumanLoopName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return sagemaker_a2i_runtime(
        "GET",
        "/human-loops/$(HumanLoopName)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_human_loop(
    HumanLoopName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sagemaker_a2i_runtime(
        "GET",
        "/human-loops/$(HumanLoopName)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_human_loops(flow_definition_arn)
    list_human_loops(flow_definition_arn, params::Dict{String,<:Any})

Returns information about human loops, given the specified parameters. If a human loop was
deleted, it will not be included.

# Arguments
- `flow_definition_arn`: The Amazon Resource Name (ARN) of a flow definition.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"CreationTimeAfter"`: (Optional) The timestamp of the date when you want the human loops
  to begin in ISO 8601 format. For example, 2020-02-24.
- `"CreationTimeBefore"`: (Optional) The timestamp of the date before which you want the
  human loops to begin in ISO 8601 format. For example, 2020-02-24.
- `"MaxResults"`: The total number of items to return. If the total number of available
  items is more than the value specified in MaxResults, then a NextToken is returned in the
  output. You can use this token to display the next page of results.
- `"NextToken"`: A token to display the next page of results.
- `"SortOrder"`: Optional. The order for displaying results. Valid values: Ascending and
  Descending.
"""
function list_human_loops(
    FlowDefinitionArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return sagemaker_a2i_runtime(
        "GET",
        "/human-loops",
        Dict{String,Any}("FlowDefinitionArn" => FlowDefinitionArn);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_human_loops(
    FlowDefinitionArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sagemaker_a2i_runtime(
        "GET",
        "/human-loops",
        Dict{String,Any}(
            mergewith(
                _merge, Dict{String,Any}("FlowDefinitionArn" => FlowDefinitionArn), params
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    start_human_loop(flow_definition_arn, human_loop_input, human_loop_name)
    start_human_loop(flow_definition_arn, human_loop_input, human_loop_name, params::Dict{String,<:Any})

Starts a human loop, provided that at least one activation condition is met.

# Arguments
- `flow_definition_arn`: The Amazon Resource Name (ARN) of the flow definition associated
  with this human loop.
- `human_loop_input`: An object that contains information about the human loop.
- `human_loop_name`: The name of the human loop.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"DataAttributes"`: Attributes of the specified data. Use DataAttributes to specify if
  your data is free of personally identifiable information and/or free of adult content.
"""
function start_human_loop(
    FlowDefinitionArn,
    HumanLoopInput,
    HumanLoopName;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sagemaker_a2i_runtime(
        "POST",
        "/human-loops",
        Dict{String,Any}(
            "FlowDefinitionArn" => FlowDefinitionArn,
            "HumanLoopInput" => HumanLoopInput,
            "HumanLoopName" => HumanLoopName,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function start_human_loop(
    FlowDefinitionArn,
    HumanLoopInput,
    HumanLoopName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sagemaker_a2i_runtime(
        "POST",
        "/human-loops",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "FlowDefinitionArn" => FlowDefinitionArn,
                    "HumanLoopInput" => HumanLoopInput,
                    "HumanLoopName" => HumanLoopName,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    stop_human_loop(human_loop_name)
    stop_human_loop(human_loop_name, params::Dict{String,<:Any})

Stops the specified human loop.

# Arguments
- `human_loop_name`: The name of the human loop that you want to stop.

"""
function stop_human_loop(HumanLoopName; aws_config::AbstractAWSConfig=global_aws_config())
    return sagemaker_a2i_runtime(
        "POST",
        "/human-loops/stop",
        Dict{String,Any}("HumanLoopName" => HumanLoopName);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function stop_human_loop(
    HumanLoopName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sagemaker_a2i_runtime(
        "POST",
        "/human-loops/stop",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("HumanLoopName" => HumanLoopName), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
using Documenter
using PowerSimulations

makedocs(
    sitename = "PowerSimulations.jl",
    format = Documenter.HTML(mathengine = Documenter.MathJax()),
    modules = [PowerSimulations],
    strict = true,
    authors = "Jose Daniel Lara, Clayton Barrows and Dheepak Krishnamurthy",
    pages = Any[
        "Introduction" => "index.md",
        #"Quick Start Guide" => "qs_guide.md",
        "Logging" => "man/logging.md",
        "Operation Model" => "man/op_problem.md",
        "Simulation Recorder" => "man/simulation_recorder.md",
        "API" => Any["PowerSimulations" => "api/PowerSimulations.md"],
    ],
)

deploydocs(
    repo = "github.com/NREL-SIIP/PowerSimulations.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "master",
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#"],
)


text:
using Test

@testset "CHSH facet regression using relabels" begin

using BellScenario: Degeneracy, LocalPolytope, ConvexPolytope

function chsh_inequalities()
    # positivity
    I_ineq = [0 0 0 0 0 -1 0 0 0]
    # non-trivial CH inequality
    II_ineq = [0 -1 0 -1 0 1 1 1 -1]

    (α_out_relabels, β_out_relabels) = Degeneracy.bipartite_output_relabels((2,2),(2,2),"non-signaling")
    in_relabels = Degeneracy.bipartite_input_relabels((2,2),(2,2),"non-signaling")

    ineqs = []
    # TODO: resolve the 3x redundancy in relabeled inequalities
    for α_out_relabel in α_out_relabels
        for β_out_relabel in β_out_relabels
            for in_relabel in in_relabels
                push!(ineqs, I_ineq*in_relabel*α_out_relabel*β_out_relabel)
                push!(ineqs, II_ineq*in_relabel*α_out_relabel*β_out_relabel)
            end
        end
    end

    unique(ineqs)
end

@testset "verifying polytope" begin

    dir = "./test/regression/files/"

    vertices = LocalPolytope.vertices((2,2),(2,2))
    @test size(vertices) == (16,)
    @test size(vertices[1]) == (8,1)

    constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)
    @test constraints[2] == []
    @test constraints[3] == []

    ineq_matches = chsh_inequalities()

    inequalities = constraints[1]
    issetequal(inequalities, ineq_matches)
end

end


text:
Base.eval(:(have_color = true))

using Documenter
using Glob

using ClickHouse

makedocs(
    sitename="ClickHouse.jl",
    pages = [
        "index.md",
        "Manual" => ["client.md"],
        "_index.md"
    ]
)

deploydocs(
    repo = "github.com/maxmouchet/ClickHouse.jl.git",
)


text:
module SymmetricSparseMatricesTests

using Test

@time @testset "SymmetricSparseMatrix" begin include("SymmetricSparseMatrix/runtests.jl") end

end # module


text:
"""
    FeatureDescriptor

A feature, in this context, refers to a distinct attribute manifested
by an atomic system.
A FeatureDescriptor describes a feature/class of features - i.e., its name,
possible values, etc.

A FeatureDescriptor does NOT store any actual instances of the value(s) of the
feature it describes.
Simply put, it can be understood to be "features of a feature".

All FeatureDescriptors MUST also describe an encoding and decoding scheme.
This can (and should) be easily done using a Codec.
"""
module FeatureDescriptor

using Base: Int16
using ..ChemistryFeaturization.AbstractType:
    AbstractAtoms, AbstractFeatureDescriptor, AbstractCodec

import ..ChemistryFeaturization.encodable_elements
encodable_elements(fd::AbstractFeatureDescriptor) =
    throw(MethodError(encodable_elements, fd))
export encodable_elements

import ..ChemistryFeaturization.encode
"""
    encode(fd::AbstractAtomFeatureDescriptor, atoms::AbstractAtoms)
Encode `atoms` using the feature descriptor `fd`.
"""
encode(fd::AbstractFeatureDescriptor, atoms::AbstractAtoms) = fd(atoms)
export encode

import ..ChemistryFeaturization.decode
"""
    decode(fd::AbstractAtomFeatureDescriptor, encoded_feature)
Decode `encoded_feature` using the feature descriptor `fd`.
"""
decode(fd::AbstractFeatureDescriptor, encoded_feature) =
    fd.encoder_decoder(fd, encoded_feature)
export decode

(codec::AbstractCodec)(fd::AbstractFeatureDescriptor, a::AbstractAtoms) = error(
    "Logic specifying how $(typeof(codec))'s encoding mechanism actually encodes $(typeof(fd)) is undefined.",
)
(codec::AbstractCodec)(fd::AbstractFeatureDescriptor, encoded_feature) = error(
    "Logic specifying how $(typeof(codec))'s decoding mechanism actually decodes $(typeof(fd)) is undefined.",
)

output_shape(efd::AbstractFeatureDescriptor) = output_shape(efd, efd.encoder_decoder)
export output_shape

include("abstractfeatures.jl")

include("bondfeatures.jl")

include("elementfeature.jl")
export ElementFeatureDescriptor, encode, decode, output_shape

include("orbitalfeature.jl")
export OrbitalFeatureDescriptor

include("pairfeature.jl")
export PairFeatureDescriptor

include("speciesfeature.jl")
export SpeciesFeatureDescriptor

end


text:
# # @__NAME__

# PREAMBLE

# PKG_SETUP

# ## Floating-point error

# In computers, real numbers are commonly approximated using floating-point numbers, such as Julia's `Float64`. Unfortunately, not all real numbers can be exactly represented as a finite-size floating-point number, and the results of operations on floating-point numbers can only approximate the results of applying the operation to a true real number. This results in peculiarities like:

2.6 - 0.7 - 1.9

# IntervalArithmetic.jl can be used to quantify floating point error, by computing _rigorous_ worst-case bounds on floating point error, within which the true result is _guaranteed_ to lie.

using Pkg # hide
Pkg.activate(@__DIR__) # hide
Pkg.instantiate() # hide
using IntervalArithmetic

# IntervalArithmetic.jl provides the `Interval` type, which stores an upper and a lower bound:

i = Interval(1.0, 2.0)

#-

dump(i)

# IntervalArithmetic.jl provides overloads for most common Julia functions that take these bounds into account. For example:

i + i

#-

sin(i)

# Note that the bounds computed by IntervalArithmetic.jl take floating point error into account. Also note that a given real number, once converted to (approximated by) a floating-point number may not be equal to the original real number. To rigorously construct an `Interval` that contains a given real number as an input, IntervalArithmetic.jl provides the `@interval` macro:

i = @interval(2.9)
i.lo === i.hi

#-

dump(i)

# Compare this to

i = Interval(2.9)
i.lo === i.hi

#-

dump(i)

# As an example, consider again the peculiar result from before, now using interval arithmetic:

i = @interval(2.6) - @interval(0.7) - @interval(1.9)

# showing that the true result, `0`, is indeed in the guaranteed interval, and indeed:

using Test
@test (2.6 - 0.7 - 1.9) ∈ i

# ## Accuracy of RigidBodyDynamics.jl's `mass_matrix`

# Let's use IntervalArithmetic.jl to establish rigorous bounds on the accuracy of the accuracy of the `mass_matrix` algorithm for the Acrobot (double pendulum) in a certain configuration. Let's get started.

using RigidBodyDynamics

# We'll create a `Mechanism` by parsing the Acrobot URDF, passing in `Interval{Float64}` as the type used to store the parameters (inertias, link lengths, etc.) of the mechanism. Note that the parameters parsed from the URDF are treated as floating point numbers (i.e., like `Interval(2.9)` instead of `@interval(2.9)` above).

const T = Interval{Float64}
srcdir = dirname(pathof(RigidBodyDynamics))
urdf = joinpath(srcdir, "..", "test", "urdf", "Acrobot.urdf")
const mechanism = parse_urdf(urdf; scalar_type=T)
state = MechanismState(mechanism)

# Let's set the initial joint angle of the shoulder joint to the smallest `Interval{Float64}` containing the real number $1$, and similarly for the elbow joint:

shoulder, elbow = joints(mechanism)
set_configuration!(state, shoulder, @interval(1))
set_configuration!(state, elbow, @interval(2));

# And now we can compute the mass matrix as normal:

M = mass_matrix(state)

# Woah, those bounds look pretty big. RigidBodyDynamics.jl must not be very accurate! Actually, things aren't so bad; the issue is just that IntervalArithmetic.jl isn't kidding when it comes to guaranteed bounds, and that includes printing the numbers in shortened form. Here are the lengths of the intervals:

err = map(x -> x.hi - x.lo, M)

#-

@test maximum(abs, err) ≈ 0 atol = 1e-14

# ## Rigorous (worst-case) uncertainty propagation

# IntervalArithmetic.jl can also be applied to propagate uncertainty in a rigorous way when the inputs themselves are uncertain. Consider for example the case that we only know the joint angles up to $\pm 0.05$ radians:

set_configuration!(state, shoulder, @interval(0.95, 1.05))
set_configuration!(state, elbow, @interval(1.95, 2.05));

# and let's compute bounds on the center of mass position:

center_of_mass(state)

# Note that the bounds on the $y$-coordinate are very tight, since our mechanism only lives in the $x$-$z$ plane.

max_stars_repo_path: examples/7. Rigorous error bounds using IntervalArithmetic/7. Rigorous error bounds using IntervalArithmetic.jl