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

using Markdown
using InteractiveUtils

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

# ╔═╡ c3e52bf2-ca9a-11ea-13aa-03a4335f2906
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
			Pkg.PackageSpec(name="Plots", version="1.6-1"),
			Pkg.PackageSpec(name="PlutoUI", version="0.6.8-0.6"),
			])
	using Plots
	using PlutoUI
	using LinearAlgebra
end

# ╔═╡ 1df32310-19c4-11eb-0824-6766cd21aaf4
md"_homework 7, version 1_"

# ╔═╡ 1e01c912-19c4-11eb-269a-9796cccdf274
# WARNING FOR OLD PLUTO VERSIONS, DONT DELETE ME

html"""

"""

# ╔═╡ 1e109620-19c4-11eb-013e-1bc95c14c2ba
md"""

# **Homework 7**: _Raytracing in 2D_
`18.S191`, fall 2020

This notebook contains _built-in, live answer checks_! In some exercises you will see a coloured box, which runs a test case on your code, and provides feedback based on the result. Simply edit the code, run it, and the check runs again.

_For MIT students:_ there will also be some additional (secret) test cases that will be run as part of the grading process, and we will look at your notebook and write comments.

Feel free to ask questions!
"""

# ╔═╡ 1e202680-19c4-11eb-29a7-99061b886b3c
# edit the code below to set your name and kerberos ID (i.e. email without @mit.edu)

student = (name = "Alison Davey", kerberos_id = "none")

# you might need to wait until all other cells in this notebook have completed running.
# scroll around the page to see what's up

# ╔═╡ 1df82c20-19c4-11eb-0959-8543a0d5630d
md"""

Submission by: **_$(student.name)_** ($(student.kerberos_id)@mit.edu)
"""

# ╔═╡ 1e2cd0b0-19c4-11eb-3583-0b82092139aa
md"_Let's create a package environment:_"

# ╔═╡ 92290e54-1940-11eb-1a24-5d1eaee9f6ca
md"""
## **Exercise 1:** _Walls_

As discussed the lecture, event-driven simulations are the traditional method used for raytracing. Here, we look for any objects in our path and _analytically_ determine how far away they are. From there, we take one big timestep all the way to the surface boundary, calculate refraction or reflection to see what direction we are moving in, and then seek out any other object we could potentially run into.

So let's start simple with determining when a ray of light could intersect with a wall.

#### Exercise 1.1 - _what is a wall?_

To start, let's create the concept of a wall. For our purposes, walls will be infinitely long, so we only need to create an object that has a position and a normal vector at that position:
"""

# ╔═╡ d851a202-1ca0-11eb-3da0-51fcb656783c
abstract type Object end

# ╔═╡ 99c61b74-1941-11eb-2323-2bdb7c120a28
struct Wall <: Object
	"Position"
	position::Vector{Float64}

	"Normal vector"
	normal::Vector{Float64}
end

# ╔═╡ 0906b340-19d3-11eb-112c-e568f69deb5d
test_wall = Wall(
	[8,-1],
	normalize([-3,1]),
)

# ╔═╡ 6de1bafc-1a01-11eb-3d67-c9d9b6c3cea8
function plot_object!(p, wall::Wall)
	# old_xlims = xlims(p)
	# old_ylims = ylims(p)

	adjacent = [wall.normal[2], -wall.normal[1]]

	a = wall.position + adjacent * 20
	b = wall.position - adjacent * 20

	line = [a, b]

	plot!(p, first.(line), last.(line), label="Wall")
	# xlims!(p, old_xlims)
	# xlims!(p, old_xlims)
end

# ╔═╡ 5f551588-1ac4-11eb-1f86-197442f1ef1d
md"""
In our simulations, we will enclose our scene in a box of **four walls**, to make sure that no rays can escape the scene. We have written this box (i.e. vector of walls) below, but we are still missing the roof.
"""

# ╔═╡ ac9bafaa-1ac4-11eb-16c4-0df8133f9c98
box_scene = [
	Wall(
		[10,0],
		[-1,0]
		),
	Wall(
		[-10,0],
		[1,0]
		),
	Wall(
		[0,-10],
		[0,1]
		),
	Wall(
		[0,10],
		[0,-1]
		)
	]

# ╔═╡ 293776f8-1ac4-11eb-21db-9d023c09e89f
md"""
👉 Modify the definition of `box_scene` to be a vector of 4 walls, instead of 3. The fourth wall should be positioned at `[0,10]`, and point downwards.
"""

# ╔═╡ e5ed6098-1c70-11eb-0b58-31d1830b9a10
md"""
In the next exercise, we will find the intersection of a ray of light and a wall. To represent light, we create a `struct` called **`Photon`**, holding the position and travel direction of a single particle of light. We also include the _index of refraction_ of the medium it is currently traveling in, we will use this later.
"""

# ╔═╡ 24b0d4ba-192c-11eb-0f66-e77b544b0510
struct Photon
	"Position vector"
	p::Vector{Float64}

	"Direction vector"
	l::Vector{Float64}

    "Current Index of Refraction"
	ior::Real
end

# ╔═╡ 925e98d4-1c78-11eb-230d-994518f0060e
test_photon = Photon([-1, 2], normalize([1,-.8]), 1.0)

# ╔═╡ eabca8ce-1c73-11eb-26ad-271f6eba889b
function plot_photon_arrow!(p, photon::Photon, length=2; kwargs...)
	line = [photon.p, photon.p .+ length*photon.l]

	plot!(p, first.(line), last.(line); lw=2, arrow=true, color=:darkred, kwargs..., label=nothing)
	scatter!(p, photon.p[1:1], photon.p[2:2]; color=:darkred, markersize=3, label=nothing, size=(150,150), tickfont = (4, :orange), kwargs...)
end

# ╔═╡ aa43ef1c-1941-11eb-04de-552719a08da0
md"""
$(html"

")
#### Exercise 1.2 - _how far is the wall?_
We will write a function that finds the location where a photon hits the wall. Instead of moving the photon forward in small timesteps until we reach the wall, we will compute the intersection directly, making use of the fact that the wall is a geometrically simple object.

Our function will return one of two possible types: a `Miss` or a `Intersection`. We define these types below, and both definitions need some elaboration.
"""

# ╔═╡ 8acef4b0-1a09-11eb-068d-79a259244ed1
struct Miss end

# ╔═╡ 8018fbf0-1a05-11eb-3032-95aae07ca78f
struct Intersection{T<:Object}
	object::T
	distance::Float64
	point::Vector{Float64}
end

# ╔═╡ e9c5d68c-1ac2-11eb-04ec-3b72eb133239
md"""
##### `Miss`
is a struct with _no fields_. It does not contain any information, except the fact that it is a `Miss`. You create a new `Miss` object like so:
"""

# ╔═╡ 5a9d00f6-1ac3-11eb-01fb-53c35796e766
a_miss = Miss()

# ╔═╡ 5aa7c4e8-1ac3-11eb-23f3-03bd58e75c4b
md"""
##### `Intersection`
is a **parametric type**. The first field (`object`) is of type `T`, and `T` is a subtype of `Object`. Have a look at the definition above, and take note of how we write such statements in Julia syntax.

We also could have used `Object` directly as the type for the field `object`, but what's special about parametric types is that `T` becomes "part of the type". Let's have a look at an example:
"""

# ╔═╡ 9df1d0f2-1ac3-11eb-0eac-d90eccca669c
test_intersection_1 = Intersection(test_wall, 3.0, [1.0,2.0])

# ╔═╡ bc10541e-1ac3-11eb-0b5f-916922f1a8e8
typeof(test_intersection_1)

# ╔═╡ d39f149e-1ac3-11eb-39a2-41c2030d7d49
md"""
You see that `Wall` is **included in the type**. This will be very useful later, when we want to do something different _depending on the intersected object_ (wall, sphere, etc.) using multiple dispatch. We can write one method for `::Intersection{Sphere}`, and one for `::Intersection{Wall}`.
"""

# ╔═╡ e135d490-1ac2-11eb-053e-914051f16e31
md"""
 $(html"

")

##### Wall geometry
So, how do we find the location where it hits the wall? Well, because our walls are infinitely long, we are essentially trying to find the point at which 2 lines intersect.

To do this, we can combine a few dot products: one to find how far away we are, and another to scale that distance. Mathematically, it would look like:

$D = -\frac{(p_{\text{ray}} - p_{\text{wall}})\cdot \hat n}{\hat \ell \cdot \hat n},$

where $p$ is the position, $\hat \ell$ is the direction of the light, and $\hat n$ is the normal vector for the wall. subscripts $i$, $r$, and $w$ represent the intersection point, ray, and wall respectively. The result is $D$, the amount that the photon needs to travel until it hits the wall.

👉 Write a function `intersection_distance` that implements this formula, and returns $D$. You can use `dot(a,b)` to compute the vector dot product ``a \cdot b``.
"""

# ╔═╡ abe3de54-1ca0-11eb-01cd-11fe798bfb97
function intersection_distance(photon::Photon, wall::Wall)
	return -dot((photon.p - wall.position),wall.normal)/dot(photon.l,wall.normal)
end

# ╔═╡ 42d65f56-1aca-11eb-1079-e32f85554349
md"""
 $(html"

")
#### Exercise 1.3 - _hitting the wall_

👉 Write a function `intersection` that takes a `photon` and a `wall`, and returns either a `Miss` or an `Intersection`, based on the result of `intersection_distance(photon, wall)` ``= D``.

If $D$ is _positive_, then the photon will hit the wall, and we should return an `Intersection`. We already have the intersected object, and we have $D$, our intersection distance. To find the intersection _point_, we use the photon's position and velocity.

$p_{\text{intersection}} = p_{\text{ray}} + D\hat \ell$

If $D$ is _negative_ (or zero), then the wall is _behind_ the photon - we should return a `Miss`.

##### Floating points
We are using _floating points_ (`Float64`) to store positions, distances, etc., which means that we need to account for small errors. Like in the lecture, we will not check for `D > 0`, but `D > ϵ` with `ϵ = 1e-3`.
"""

# ╔═╡ a5847264-1ca0-11eb-0b45-eb5388f6e688
function intersection(photon::Photon, wall::Wall; ϵ=1e-3)
	D = intersection_distance(photon, wall)
	if D > ϵ
		return Intersection(wall, D, photon.p+D*photon.l)
	else
		return Miss()
	end
end

# ╔═╡ 7f286ccc-1c75-11eb-1270-95a87840b300
@bind dizzy_angle Slider(0:0.0001:2π, default=2.2)

# ╔═╡ 6544be90-19d3-11eb-153c-218025f738c6
dizzy = Photon([0, 1], normalize([cos(dizzy_angle + π),sin(dizzy_angle + π)]), 1.0)

# ╔═╡ d70380a4-1ad0-11eb-1184-f7e9b84a83ad
md"""
 $(html"

")
#### Exercise 1.4 - _which wall?_

We are now able to find the `Intersection` of a single photon with a single wall (or detect a `Miss`). Great! To make our simulation more interesting, we will combine **multiple walls** into a single scene.
"""

# ╔═╡ 55187168-1c78-11eb-1182-ab4336b577a4
philip = Photon([3, 0], normalize([.5,-1]), 1.0)

# ╔═╡ 2158a356-1a05-11eb-3f5b-4dfa810fc602
ex_1_scene = [box_scene..., test_wall]

# ╔═╡ 87a8e280-1c7c-11eb-2bb0-034011f6c10f
md"""
When we shoot a photon at the scene, we compute the intersections between the photon and every object in the scene. Click on the vector below to see all elements:
"""

# ╔═╡ 4d69c36a-1c73-11eb-3ae3-23900db09c27
md"""
There are two misses and three intersections. Just what we hoped!
"""

# ╔═╡ 5342430e-1c79-11eb-261c-15abd0f8cfc1
md"""
So which of these **five** results should we use to determine what the photon does next? It should be the _closest intersection_.

Because we used two different types for hits and misses, we can express this in a charming way. We define what it means for one to be better than the other:
"""

# ╔═╡ 6c37c5f4-1a09-11eb-08ae-9dce752f29cb
begin
	Base.isless(a::Miss, b::Miss) = false
	Base.isless(a::Miss, b::Intersection) = false
	Base.isless(a::Intersection, b::Miss) = true
	Base.isless(a::Intersection, b::Intersection) = a.distance  Note that we did not define the `sort` and `minimum` methods ourselves! We only added methods for `Base.isless`.

By taking the minimum, we have found our closest hit! Let's turn this into a function.

👉 Write a function `closest_hit` that takes a `photon` and a vector of objects. Calculate the vector of `Intersection`s/`Miss`es, and return the `minimum`.
"""

# ╔═╡ e9c6a0b8-1ad0-11eb-1606-0319caf0948a
md"""
 $(html"

")
## **Exercise 2:** _Mirrors_
"""

# ╔═╡ 522e6b22-194d-11eb-167c-052e65f6b703
md"""
Now we're going to make a bold claim: All walls in this simulation are mirrors. This is just for simplicity so we don't need to worry about rays stopping at the boundaries.

We are already able to find the intersection of a light ray with a mirror, but we still need to tell our friendly computer what a _reflection_ is.
"""

# ╔═╡ dad5acfa-194c-11eb-27f9-01f40342a681
md"
#### Exercise 2.1 - _reflect_

For this one, we need to implement a reflection function. This one is way easier than refraction. All we need to do is find how much of the light is moving in the direction of the surface's normal and subtract that twice.

$\ell_2 = \ell_1 - 2(\ell_1\cdot \hat n)\hat n$

Where $\ell_1$ and $\ell_2$ are the photon directions before and after the reflection off a surface with normal ``\hat{n}``. Let's write that in code:
"

# ╔═╡ 43306bd4-194d-11eb-2e30-07eabb8b29ef
reflect(ℓ₁::Vector, n̂::Vector)::Vector = ℓ₁ - 2 * dot(ℓ₁, n̂) * n̂

# ╔═╡ 70b8401e-1c7e-11eb-16b2-d54d8f66d71a
md"""
👉 Verify that the function `reflect` works by writing a simple test case:
"""

# ╔═╡ 48668080-22b4-11eb-3c4f-eb6aa1169382
ex_1_scene[1]

# ╔═╡ 79532662-1c7e-11eb-2edf-57e7cfbc1eda
reflect(philip.l, ex_1_scene[1].normal )

# ╔═╡ b6614d80-194b-11eb-1edb-dba3c29672f8
md"""
#### Exercise 2.2 - _step_

Our event-driven simulation is a stepping method, but instead of taking small steps in time, we take large steps from one collision event to the next.

👉 Write a function `interact` that takes a photon and a `hit::Intersection{Wall}` and returns a new `Photon` at the next step. The new photon is located at the hit point, its direction is reflected off the wall's normal and the `photon.ior` is reused.
"""

# ╔═╡ 2c6defd0-1ca1-11eb-17db-d5cb498f3265
function interact(photon::Photon, hit::Intersection{Wall})
	return Photon(hit.point, reflect(photon.l,hit.object.normal), photon.ior)
end

# ╔═╡ 3f727a2c-1c80-11eb-3608-e55ccb9786d9
md"""
For convenience, we define a function `step_ray` that combines these two actions: it finds the closest hit and computes the interaction.
"""

# ╔═╡ a45e1012-194d-11eb-3252-bb89daed3c8d
md"
Great! Next, we will repeat this action to trace the path of a photon.
"

# ╔═╡ 7ba5dda0-1ad1-11eb-1c4e-2391c11f54b3
md"""
#### Exercise 2.3 - _accumulate_

👉 Write a function `trace` that takes an initial `Photon`, a vector of `Object`s and `N`, the number of steps to make. Return a vector of `Photon`s. Try to use `accumulate`.
"""

# ╔═╡ 3cd36ac0-1a09-11eb-1818-75b36e67594a
@bind mirror_test_ray_N Slider(1:30; default=4)

# ╔═╡ 7478330a-1c81-11eb-2f9f-099f1111032c
md"""
#### Recap
In Exercise 3 and 4, we will add a `Sphere` type, and our scene will consist of `Wall`s (mirrors) and `Sphere`s (lenses). But before we move on, let's review what we have done so far.

Our main character is a `Photon`, which bounces around a scene made up of `Wall`s.

1. Using `intersection(photon, wall::Wall)` we can find either an `Intersection` (containing the `wall`, the `distance` and the `point`) or a `Miss`.
2. Our scene is just a `Vector` or objects, and we compute the intersection between the photon and every object.
3. By adding `Base.isless` methods we have told Julia how to compare hit/miss results, and we get the closest one using `minimum(all_intersections)`.
4. We wrote a function `interact(photon, hit::Intersection{Wall})` that returns a new photon after interacting with a wall collision.

We repeat these four steps to trace a ray through the scene.

---

In the next two exercises we will reuse some of the functionality that we have already written, using multiple dispatch! For example, we add a method `intersection(photon, sphere::Sphere)`, and steps 2 and 3 magically also work with spheres!

"""
# We have a type `Photon` and a type `Wall`, and u

# ╔═╡ ba0a869a-1ad1-11eb-091f-916e9151f052
md"""
 $(html"

")

## **Exercise 3:** _Spheres_
Now that we know how to bounce light around mirrors, we want to simulate a _spherical lens_ to make things more interesting. Let's define a `Sphere`.
"""

# ╔═╡ 3aa539ce-193f-11eb-2a0f-bbc6b83528b7
struct Sphere <: Object
	center::Vector{Float64}
	radius::Real
	ior::Real
end

# ╔═╡ caa98732-19cd-11eb-04ce-2f018275cf01
function plot_object!(p::Plots.Plot, sphere::Sphere)
	points = [
		sphere.center .+ sphere.radius .* [cos(ϕ), sin(ϕ)]
		for ϕ in LinRange(0,2π,50)
	]

	plot!(p, points .|> first, points .|> last, seriestype=:shape, label="Sphere", fillopacity=0.2)
	p
end

# ╔═╡ eff9329e-1a05-11eb-261f-734127d36750
function plot_scene(objects::Vector{<:Object}; kwargs...)
	p = plot(aspect_ratio=:equal; kwargs...)

	for o in objects
		plot_object!(p, o)
	end
	p
end

# ╔═╡ e45e1d36-1a12-11eb-2720-294c4be6e9fd
plot_scene([test_wall], size=(150,100), tickfont = (4, :orange))

# ╔═╡ 0393dd3a-1a06-11eb-18a9-494ae7a26bc0
plot_scene(box_scene, legend=false, size=(150,150), tickfont = (4, :orange))

# ╔═╡ 76d4351c-1c78-11eb-243f-5f6f5e485d5d
let
	p = plot_scene(box_scene, legend=false, size=(400,200))
	plot_photon_arrow!(p, test_photon, 7)
	p
end

# ╔═╡ 5501a700-19ec-11eb-0ded-53e41f7f821a
let
	p = plot_scene(ex_1_scene, legend=false, size=(400,200))
	plot_photon_arrow!(p, philip, 5, size=(150,150), legendfontsize = 2, tickfont = (4, :orange))
end

# ╔═╡ e5c0e960-19cc-11eb-107d-39b397a783ab
example_sphere = Sphere(
	[7, -6],
	2,
	1.5,
)

# ╔═╡ 2a2b7284-1ade-11eb-3b71-d17fe2ca638a
plot_scene([example_sphere], legend=false, xlim=(-15,15), ylim=(-10,10), size=(150,150), tickfont = (4, :orange))

# ╔═╡ e2a8d1d6-1add-11eb-0da1-cda1492a950c
md"
#### Exercise 3.1
Just like with the `Wall`, our first step is to be able to find the intersection point of a ray of light and a sphere.

This one is a bit more challenging than the intersection with the wall, in particular because there are 3 potential outcomes of a line interacting with a sphere:
- No intersection
- 1 intersection
- 2 intersections

As shown below:
"

# ╔═╡ 337918f4-194f-11eb-0b45-b13fef3b23bf
PlutoUI.Resource("https://upload.wikimedia.org/wikipedia/commons/6/67/Line-Sphere_Intersection_Cropped.png")

# ╔═╡ 492b257a-194f-11eb-17fb-f770b4d3da2e
md"""
So we need a way of finding all of these.

To start, let's look at the intersection of a **point** and a sphere. So long as the relative distance between the photon and the sphere's center satisfies the sphere equation, we can be considered inside of the sphere. More specifically, we are inside the sphere if:

$(x_s-x_p)^2+(y_s-y_p) < r^2.$

where the $s$ and $p$ subscripts represent the sphere and photon, respectively. We know we are *on* the sphere if

$(x_s-x_p)^2+(y_s-y_p) = r^2.$

Let's rewrite this in vector notation as:

$(\mathbf{R} - \mathbf{S})\cdot(\mathbf{R} - \mathbf{S}) = r^2,$

where $\mathbf{R}$ and $\mathbf{S}$ are the $x$, $y$, and $z$ location of the photon and sphere, respectively.

Returning to the timestepping example from above, we know that our ray is moving forward with time such that $\mathbf{R} = \mathbf{R}_0 + v dt = \mathbf{R}_0 + \ell t$. We now need to ask ourselves if there is any time when our ray interacts with the sphere. Plugging this in to the dot product from above, we get

$(\mathbf{R}_0 + \ell t - \mathbf{S})\cdot(\mathbf{R}_0 + \ell t - \mathbf{S}) = r^2$

To solve this for $t$, we first need to reorder everything into the form of a polynomial, such that:

$t^2(\ell\cdot\ell)+2t\ell\cdot(\mathbf{R_0}-\mathbf{S})+(\mathbf{R}_0-\mathbf{S})\cdot(\mathbf{R}_0-\mathbf{S}) - r^2=0.$

This can be solved with the good ol' fashioned quadratic equation:

$\frac{-b\pm\sqrt{b^2-4ac}}{2a},$

where $a = \ell\cdot\ell$, $b = 2\ell\cdot(\mathbf{R}_0-\mathbf{S})$, and $c=(\mathbf{R}_0-\mathbf{S})\cdot(\mathbf{R}_0-\mathbf{S}) - r^2$

If the quadratic equation returns no roots, there is no intersection. If it returns 1 root, the ray just barely hits the edge of the sphere. If it returns 2 roots, it goes right through!

The easiest way to check this is by looking at the discriminant $d = b^2-4ac$.

```math
\text{Number of roots} = \left\{
    \begin{align}
       &0, \qquad \text{if } d < 0 \\
       &1, \qquad \text{if } d = 0 \\
       &2, \qquad \text{if } d > 0 \\
    \end{align}
\right.
```

In the case that there are 2 roots, the second root corresponds to when the ray would interact with the far edge of the sphere *if there were no refraction or reflection!*; therefore, we only care about returning the closest point.

With all this said, we are ready to write some code.

👉 Write a new method `intersection` that takes a `Photon` and a `Sphere`, and returns either a `Miss` or an `Intersection`, using the method described above. Go back to Exercise 1.3 where we defined the first method, and see how we adapt it to a sphere.
"""

# ╔═╡ 392fe192-1ca1-11eb-36c4-f9bd2b01a5e5
function intersection(photon::Photon, sphere::Sphere; ϵ=1e-3)

	a = dot(photon.l,photon.l)
	b = 2*dot(photon.l,(photon.p-sphere.center))
	c = dot((photon.p-sphere.center),(photon.p-sphere.center))-sphere.radius^2
	d = b^2-4a*c

	if d < ϵ
		return Miss()

	else
		t1 = (-b-sqrt(d))/2a
		t2 = (-b+sqrt(d))/2a

		t = if t1 > ϵ
			t1
		elseif t2 > ϵ
			t2
		else
			return Miss()
		end

		return Intersection(sphere, t, photon.p+t*photon.l)
	end
end

# ╔═╡ a306e880-19eb-11eb-0ff1-d7ef49777f63
test_intersection = intersection(dizzy, test_wall)

# ╔═╡ 3663bf80-1a06-11eb-3596-8fbbed28cc38
let
	p = plot_scene([test_wall])
	plot_photon_arrow!(p, dizzy, 4; label="Philip")

	try
		scatter!(p, test_intersection.point[1:1], test_intersection.point[2:2], label="Intersection point", size=(150,150), legendfontsize = 2, tickfont = (4, :orange))
	catch
	end

	p
end

# ╔═╡ 1b0c0e4c-1c73-11eb-225d-23c731455755
all_intersections = [intersection(philip, o) for o in ex_1_scene]

# ╔═╡ e055262c-1c73-11eb-14de-2f537a19b012
let
	p = plot_scene(ex_1_scene)

	plot_photon_arrow!(p, philip, 4; label="Philip")
	for (i,hit) in enumerate(all_intersections)
		if hit isa Intersection
			scatter!(p, hit.point[1:1], hit.point[2:2], label="intersection $i",size=(150,150), legendfontsize = 2, tickfont = (4, :orange))
		end
	end
	p |> as_svg
end

# ╔═╡ c3090e4a-1a09-11eb-0f32-d3bbfd9992e0
sort(all_intersections)

# ╔═╡ 63ef21c6-1c7a-11eb-2f3c-c5ac16bc289f
minimum(all_intersections)

# ╔═╡ 19cf420e-1c7c-11eb-1cb8-dd939fee1276
function closest_hit(photon::Photon, objects::Vector{<:Object})
	return minimum([intersection(photon, o) for o in objects])
end

# ╔═╡ b8cd4112-1c7c-11eb-3b2d-29170ad9beb5
test_closest = closest_hit(philip, ex_1_scene)

# ╔═╡ a99c40bc-1c7c-11eb-036b-7fe6e9b937e5
let
	p = plot_scene(ex_1_scene)

	plot_photon_arrow!(p, philip, 4; label="Philip")

	scatter!(p, test_closest.point[1:1], test_closest.point[2:2], label="Closest hit", size=(150,150), legendfontsize = 2, tickfont = (4, :orange))

	p |> as_svg
end

# ╔═╡ af5c6bea-1c9c-11eb-35ae-250337e4fc86
test_sphere = Sphere(
	[3, 1],
	2,
	1.5,
)

# ╔═╡ 251f0262-1a0c-11eb-39a3-09be67091dc8
sphere_intersection = intersection(philip, test_sphere)

# ╔═╡ b3ab93d2-1a0b-11eb-0f5a-cdca19af3d89
ex_3_scene = [test_sphere, box_scene...]

# ╔═╡ 83aa9cea-1a0c-11eb-281d-699665da2b4f
let
	p = plot_scene(ex_3_scene)

	plot_photon_arrow!(p, philip, 4; label="Philip")
	if sphere_intersection isa Intersection
		scatter!(p, sphere_intersection.point[1:1], sphere_intersection.point[2:2], label="Intersection point", size=(150,150), legendfontsize = 2, tickfont = (4, :orange))
	end

	p |> as_svg
end

# ╔═╡ 71dc652e-1c9c-11eb-396c-cfd9ee2261fe
md"""
👉 Change the definition of `test_sphere` to test different situations:
- Hit the circle
- Miss the circle
- Start inside the circle (you should hit the exit boundary)
"""

# ╔═╡ 584ce620-1935-11eb-177a-f75d9ad8a399
md"""
 $(html"

")
## **Exercise 4:** _Lenses_

For this, we will start with refraction from the surface of water and then move on to a spherical lens.

So, how does refraction work? Well, every time light enters a new medium that is more dense than air, it will bend towards the normal to the surface, like so:

$(RemoteResource("https://upload.wikimedia.org/wikipedia/commons/thumb/3/3f/Snells_law2.svg/800px-Snells_law2.svg.png", :width=>200, :style=>"display: block; margin: auto;"))
"""

# ╔═╡ 78915326-1937-11eb-014f-fff29b3660a0
md"""
This can be described by Snell's law:

$\frac{n_1}{n_2} = \frac{v_2}{v_1} = \frac{\sin(\theta_2)}{\sin(\theta_1)}$

Here, $n$ is the index of refraction, $v$ is the speed (not velocity (sorry for the notation!)), and $\theta$ is the angle with respect to the surface normal. Any variables with an subscript of 1 are in the outer medium (air), and any variables with a subscript 2 are in the inner medium (water).

This means that we can find the angle of the new ray of light as

$\sin(\theta_2) = \frac{n_1}{n_2}\sin(\theta_1)$

The problem is that $\sin$ is slow, so we typically want to rewrite this in terms of vector operations. This means that we want to rewrite everything to be in terms of dot products, but because $A\cdot B = |A||B|cos(\theta)$, we really want to rewrite everything in terms of cosines first. So, using the fact that $\sin(\theta)^2 + \cos(\theta)^2 = 1$, we can rewrite the above equation to be:

$\sin(\theta_2) = \frac{n_1}{n_2}\sqrt{1-\cos(\theta_1)^2}$

We also know that

$\cos(\theta_2) = \sqrt{1-\sin(\theta_2)^2} = \sqrt{1-\left(\frac{n_1}{n_2}\right)^2\left(1-\cos(\theta_1)^2\right)}.$

Finally, we know that the new light direction should be the same as the old one, but shifted towards (or away) from the normal according to the new refractive index.
In particular:

$n_2 \ell _2 = {n_1} \ell _1 + (n_1\cos(\theta_1)-n_2\cos(\theta_2))\hat n,$

where $\hat n$ is the normal from the water's surface. Rewriting this, we find:

$\ell _2 = \left(\frac{n_1}{n_2}\right) \ell _1 + \left(\left(\frac{n_1}{n_2}\right)\cos(\theta_1)-\cos(\theta_2)\right)\hat n.$

Now, we already know $\cos(\theta_2)$ in terms of $\cos(\theta_1)$, so we can just plug that in... But first, let's do some simplifications, such that

$r = \frac{n_1}{n_2}$

and

$c = -\hat n \cdot \ell_1.$

Now, we can rewrite everything such that

$\ell_2 = r\ell_1 + \left(rc-\sqrt{1-r^2(1-c^2)}\right)\hat n.$

The last step is to write this in code with a function that takes the light direction, the normal, and old and new indices of refraction:
"""

# ╔═╡ 14dc73d2-1a0d-11eb-1a3c-0f793e74da9b
function refract(
		ℓ₁::Vector, n̂::Vector,
		old_ior, new_ior
	)

	r = old_ior / new_ior

	n̂_oriented = if -dot(ℓ₁, n̂) < 0
		-n̂
	else
		n̂
	end

	c = -dot(ℓ₁, n̂_oriented)

	normalize(r * ℓ₁ + (r*c - sqrt(1 - r^2 * (1 - c^2))) * n̂_oriented)
end

# ╔═╡ 71b70da6-193e-11eb-0bc4-f309d24fd4ef
md"

Now to move on to lenses. Like in lecture, we will focus exclusively on spherical lenses. Ultimately, there isn't a big difference between a lens and a spherical drop of water. It just has a slightly different refractive index and it's normal is defined slightly differently.
"

# ╔═╡ 54b81de0-193f-11eb-004d-f90ec43588f8
md"""
We need a helper functions to find the normal of the sphere's surface at any position. Remember that the normal will always be pointing perpendicularly from the surface of the sphere. This means that no matter what point you are at, the normal will just be a normalized vector of your current location minus the sphere's position:
"""

# ╔═╡ 6fdf613c-193f-11eb-0029-957541d2ed4d
function sphere_normal_at(p::Vector{Float64}, s::Sphere)
	normalize(p - s.center)
end

# ╔═╡ 392c25b8-1add-11eb-225d-49cfca27bef4
md"""
👉 Write a new method for `interact` that takes a `photon` and a `hit` of type `Intersection{Sphere}`, that implements refraction. It returns a new `Photon` positioned at the hit point, with the refracted velocity and the new index of refraction.
"""

# ╔═╡ 427747d6-1ca1-11eb-28ae-ff50728c84fe
function interact(photon::Photon, hit::Intersection{Sphere})

	normal 	= sphere_normal_at(hit.point, hit.object)

	if photon.ior == 1.0
		new_ior = hit.object.ior
	else
		new_ior = 1.0
	end

	return Photon(hit.point, refract(photon.l, normal, photon.ior, new_ior), new_ior)
end

# ╔═╡ 0b03316c-1c80-11eb-347c-1b5c9a0ae379
test_new_photon = interact(philip, test_closest)

# ╔═╡ fb70cc0c-1c7f-11eb-31b5-87b168a66e19
let
	p = plot_scene(ex_1_scene)

	plot_photon_arrow!(p, philip, 4; label="Philip")
	plot_photon_arrow!(p, test_new_photon, 4; label="Philip after interaction", size=(150,150), legendfontsize = 2, tickfont = (4, :orange))

	p |> as_svg
end

# ╔═╡ 76ef6e46-1a06-11eb-03e3-9f40a86dc9aa
function step_ray(photon::Photon, objects::Vector{<:Object})
	hit = closest_hit(photon, objects)
	interact(photon, hit)
end

# ╔═╡ 8549346e-22c0-11eb-20a3-b9964fff65e3
function trace(photon::Photon, scene::Vector{<:Object}, N)
	accumulate(1:N; init=photon) do old_photon, _
		step_ray(old_photon, scene)
	end
end

# ╔═╡ 1ee0787e-1a08-11eb-233b-43a654f70117
let
	p = plot_scene(ex_1_scene, legend=false, xlim=(-11,11), ylim=(-11,11))

	path = trace(philip, ex_1_scene, mirror_test_ray_N)

	line = [philip.p, [r.p for r in path]...]
	plot!(p, first.(line), last.(line), lw=5, color=:pink)

	plot_photon_arrow!(p, philip)
	plot_photon_arrow!.([p], path, size=(150,150), legendfontsize = 2, tickfont = (4, :orange))

	p
end |> as_svg

# ╔═╡ dced1fd0-1c9e-11eb-3226-17dc1e09e018
md"""
To test your code, modify the definition of `test_lens_photon` and `test_lens` below.
"""

# ╔═╡ 65aec4fc-1c9e-11eb-1c5a-6dd7c533d3b8
test_lens_photon = Photon([0,0], [1,0], 1.00)

# ╔═╡ 5895d9ae-1c9e-11eb-2f4e-671f2a7a0150
test_lens = Sphere([5, -2.5], 3, 1.5)

# ╔═╡ 83acf10e-1c9e-11eb-3426-bb28e7bc6c79
let
	scene = [test_lens, box_scene...]
	N = 3

	p = plot_scene(scene, legend=false, xlim=(-11,11), ylim=(-11,11))

	path = accumulate(1:N; init=test_lens_photon) do old_photon, _
		step_ray(old_photon, scene)
	end

	line = [test_lens_photon.p, [r.p for r in path]...]
	plot!(p, first.(line), last.(line), lw=5, color=:red)

	p
end |> as_svg

# ╔═╡ 13fef49c-1c9e-11eb-2aa3-d3aa2bfd0d57
md"""
By defining a method for `interact` that takes a sphere intersection, we are now able to use the machinery developed in Exercise 2 to simulate a scene with both lenses and mirrors. Let's try it out!
"""

# ╔═╡ c492a1f8-1a0c-11eb-2c38-5921c39cf5f8
@bind sphere_test_ray_N Slider(1:30; default=4)

# ╔═╡ 155aa5e0-22c3-11eb-2596-5b77fbd02fe0
let
	scene = [test_lens, box_scene...]
	p = plot_scene(scene, legend=false, xlim=(-11,11), ylim=(-11,11))

	path = accumulate(1:sphere_test_ray_N; init=test_lens_photon) do old_photon, i
		step_ray(old_photon, scene)
	end

	line = [test_lens_photon.p, [r.p for r in path]...]
	plot!(p, first.(line), last.(line), lw=5, color=:red)

	p
end |> as_svg

# ╔═╡ c00eb0a6-cab2-11ea-3887-070ebd8d56e2
md"
#### Spherical aberration
Now we can put it all together into an image of spherical aberration!
"

# ╔═╡ 3dd0a48c-1ca3-11eb-1127-e7c43b5d1666
md"""
👉 Recreate the spherical aberration figure from [the lecture](https://www.youtube.com/watch?v=MkkZb5V6HqM) (around the end of the video), and make the index of refraction interactive using a `Slider`. _Or make something else!_
"""

# ╔═╡ b314f676-233d-11eb-1858-7b5e101bfffc
@bind spher_ab_ior Slider(1:.1:2; default=1.5)

# ╔═╡ 270762e4-1ca4-11eb-2fb4-392e5c3b3e04
let
	spher_ab_scene = [
	Wall(
		[20,0],
		[-1,0]
		),
	Wall(
		[-20,0],
		[1,0]
		),
	Wall(
		[0,-10],
		[0,1]
		),
	Wall(
		[0,10],
		[0,-1]
		)
	]

	spher_ab_lens = Sphere([0.,0.], 9, spher_ab_ior)
	scene = [spher_ab_lens, spher_ab_scene...]
	N = 3
	p = plot_scene(scene, legend=false, xlim=(-21,21), ylim=(-11,11))
	spher_ab_lens_photon = Photon([-20.0,-10.],[1.,0.],1.0)

	for _ in 1:20
		spher_ab_lens_photon = Photon(spher_ab_lens_photon.p + [0.,1.],
			spher_ab_lens_photon.l,
			spher_ab_lens_photon.ior)

		path = accumulate(1:N; init=spher_ab_lens_photon) do old_photon, _
			step_ray(old_photon, scene)
		end

		line = [spher_ab_lens_photon.p, [r.p for r in path]...]
		plot!(p, first.(line), last.(line), lw=.5, color=:red,
			size=(300,150), tickfont = (2, :black))
	end

	p
end |> as_svg

# ╔═╡ bbf730c8-1ca6-11eb-3bb0-1188046339ac
md"""
## **Exercise XX:** _Lecture transcript_
(MIT students only)
Please see the link for hw 7 transcript document on [Canvas](https://canvas.mit.edu/courses/5637).
We want each of you to correct about 500 lines, but don’t spend more than 20 minutes on it.
See the the beginning of the document for more instructions.

👉 Please mention the name of the video(s) and the line ranges you edited:
"""

# ╔═╡ cbd8f164-1ca6-11eb-1440-bdaabf73a6c7
lines_i_edited = md"""
Abstraction, lines 1-219; Array Basics, lines 1-137; Course Intro, lines 1-144 (_for example_)
"""

# ╔═╡ ebd05bf0-19c3-11eb-2559-7d0745a84025
if student.name == "Jazzy Doe" || student.kerberos_id == "jazz"
	md"""
	!!! danger "Before you submit"
	    Remember to fill in your **name** and **Kerberos ID** at the top of this notebook.
	"""
end

# ╔═╡ ec275590-19c3-11eb-23d0-cb3d9f62ba92
md"## Function library

Just some helper functions used in the notebook."

# ╔═╡ ec31dce0-19c3-11eb-1487-23cc20cd5277
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]))

# ╔═╡ ad5a7420-1c7f-11eb-042f-115a9ef4c676
hint(md"`Intersection` contains the intersected object, so you can retrieve the wall using `hit.object`, and the normal using `hit.object.normal`.")

# ╔═╡ c25caf08-1a13-11eb-3c4d-0567faf4e662
md"""
You can use `ray.ior == 1.0` to check whether this is a ray _entering_ or _leaving_ the sphere.
""" |> hint

# ╔═╡ ec3ed530-19c3-11eb-10bb-a55e77550d1f
almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]))

# ╔═╡ ec4abc12-19c3-11eb-1ca4-b5e9d3cd100b
still_missing(text=md"Replace `missing` with your answer.") = Markdown.MD(Markdown.Admonition("warning", "Here we go!", [text]))

# ╔═╡ ec57b460-19c3-11eb-2142-07cf28dcf02b
keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]))

# ╔═╡ ec5d59b0-19c3-11eb-0206-cbd1a5415c28
yays = [md"Fantastic!", md"Splendid!", md"Great!", md"Yay ❤", md"Great! 🎉", md"Well done!", md"Keep it up!", md"Good job!", md"Awesome!", md"You got the right answer!", md"Let's move on to the next section."]

# ╔═╡ ec698eb0-19c3-11eb-340a-e319abb8ebb5
correct(text=rand(yays)) = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]))

# ╔═╡ 0e9a240c-1ac5-11eb-1a7e-b3c43c459484
let
	if length(box_scene) != 4
		keep_working()
	elseif !(box_scene isa Vector{Wall})
		keep_working(md"`box_scene` should be a Vector of `Wall` objects.")
	else
		w = last(box_scene)

		if w.position != [0,10]
			keep_working(md"The wall's position is not correct.")
		elseif w.normal != [0,-1]
			keep_working(md"The wall's direction is not correct.")
		else
			correct()
		end
	end
end

# ╔═╡ 0787f130-1aca-11eb-24b4-2ff2ddd0bc48
let
	p = Photon([5,0], [1,0], 1.0)
	w = Wall([10,10], normalize([-1,-1]))

	result = intersection_distance(p, w)

	if result isa Missing
		still_missing()
	elseif !(result isa Real)
		keep_working(md"You need to return a number.")
	else
		if abs(result - (20 - 5)) > 0.1
			if abs(-result - (20 - 5)) > 0.1
				keep_working(md"The returned distance is not correct.")
			else
				keep_working(md"Did you forget the minus sign?")
			end
		else
			correct()
		end
	end
end

# ╔═╡ 038d5e88-1ac7-11eb-2020-a9d7e19feebc
let
	p = Photon([5,0], [1,0], 1.0)
	w = Wall([10,10], normalize([-1,-1]))

	result = intersection(p, w)

	if result isa Missing
		still_missing()
	elseif !(result isa Miss || result isa Intersection)
		keep_working(md"You need to return a `Miss` or a `Intersection`.")
	else
		if result isa Miss
			keep_working(md"You returned a `Miss` for a photon that hit the wall.")
		else
			if abs(result.distance - (20 - 5)) > 0.1
				keep_working(md"The returned distance is not correct.")
			else
				p = Photon([5,0], [-1,0], 1.0)
				result = intersection(p, w)

				if !(result isa Miss)
					almost(md"What should happen when ``D < 0``?")
				else
					p = Photon([10,10], [1,0], 1.0)
					w = Wall([10,10], normalize([-1,-1]))

					result = intersection(p, w)
					if result isa Miss
						correct()
					else
						almost(md"Remember to use ``\epsilon``.")
					end
				end
			end
		end
	end
end

# ╔═╡ ec7638e0-19c3-11eb-1ca1-0b3aa3b40240
not_defined(variable_name) = Markdown.MD(Markdown.Admonition("danger", "Oopsie!", [md"Make sure that you define a variable called **$(Markdown.Code(string(variable_name)))**"]))

# ╔═╡ ec85c940-19c3-11eb-3375-a90735beaec1
TODO = html"TODO"

# ╔═╡ 8cfa4902-1ad3-11eb-03a1-736898ff9cef
TODO_note(text) = Markdown.MD(Markdown.Admonition("warning", "TODO note", [text]))

# ╔═╡ Cell order:
# ╟─1df32310-19c4-11eb-0824-6766cd21aaf4
# ╟─1df82c20-19c4-11eb-0959-8543a0d5630d
# ╟─1e01c912-19c4-11eb-269a-9796cccdf274
# ╟─1e109620-19c4-11eb-013e-1bc95c14c2ba
# ╟─1e202680-19c4-11eb-29a7-99061b886b3c
# ╟─1e2cd0b0-19c4-11eb-3583-0b82092139aa
# ╟─c3e52bf2-ca9a-11ea-13aa-03a4335f2906
# ╟─92290e54-1940-11eb-1a24-5d1eaee9f6ca
# ╠═d851a202-1ca0-11eb-3da0-51fcb656783c
# ╠═99c61b74-1941-11eb-2323-2bdb7c120a28
# ╠═0906b340-19d3-11eb-112c-e568f69deb5d
# ╠═e45e1d36-1a12-11eb-2720-294c4be6e9fd
# ╟─6de1bafc-1a01-11eb-3d67-c9d9b6c3cea8
# ╟─eff9329e-1a05-11eb-261f-734127d36750
# ╟─5f551588-1ac4-11eb-1f86-197442f1ef1d
# ╠═ac9bafaa-1ac4-11eb-16c4-0df8133f9c98
# ╠═0393dd3a-1a06-11eb-18a9-494ae7a26bc0
# ╟─293776f8-1ac4-11eb-21db-9d023c09e89f
# ╟─0e9a240c-1ac5-11eb-1a7e-b3c43c459484
# ╟─e5ed6098-1c70-11eb-0b58-31d1830b9a10
# ╠═24b0d4ba-192c-11eb-0f66-e77b544b0510
# ╠═925e98d4-1c78-11eb-230d-994518f0060e
# ╟─76d4351c-1c78-11eb-243f-5f6f5e485d5d
# ╠═eabca8ce-1c73-11eb-26ad-271f6eba889b
# ╟─aa43ef1c-1941-11eb-04de-552719a08da0
# ╠═8acef4b0-1a09-11eb-068d-79a259244ed1
# ╠═8018fbf0-1a05-11eb-3032-95aae07ca78f
# ╟─e9c5d68c-1ac2-11eb-04ec-3b72eb133239
# ╠═5a9d00f6-1ac3-11eb-01fb-53c35796e766
# ╟─5aa7c4e8-1ac3-11eb-23f3-03bd58e75c4b
# ╠═9df1d0f2-1ac3-11eb-0eac-d90eccca669c
# ╠═bc10541e-1ac3-11eb-0b5f-916922f1a8e8
# ╟─d39f149e-1ac3-11eb-39a2-41c2030d7d49
# ╟─e135d490-1ac2-11eb-053e-914051f16e31
# ╠═abe3de54-1ca0-11eb-01cd-11fe798bfb97
# ╟─0787f130-1aca-11eb-24b4-2ff2ddd0bc48
# ╟─42d65f56-1aca-11eb-1079-e32f85554349
# ╠═a5847264-1ca0-11eb-0b45-eb5388f6e688
# ╟─038d5e88-1ac7-11eb-2020-a9d7e19feebc
# ╟─6544be90-19d3-11eb-153c-218025f738c6
# ╟─a306e880-19eb-11eb-0ff1-d7ef49777f63
# ╠═3663bf80-1a06-11eb-3596-8fbbed28cc38
# ╟─7f286ccc-1c75-11eb-1270-95a87840b300
# ╟─d70380a4-1ad0-11eb-1184-f7e9b84a83ad
# ╠═55187168-1c78-11eb-1182-ab4336b577a4
# ╟─2158a356-1a05-11eb-3f5b-4dfa810fc602
# ╠═5501a700-19ec-11eb-0ded-53e41f7f821a
# ╟─87a8e280-1c7c-11eb-2bb0-034011f6c10f
# ╠═1b0c0e4c-1c73-11eb-225d-23c731455755
# ╟─4d69c36a-1c73-11eb-3ae3-23900db09c27
# ╠═e055262c-1c73-11eb-14de-2f537a19b012
# ╟─5342430e-1c79-11eb-261c-15abd0f8cfc1
# ╠═6c37c5f4-1a09-11eb-08ae-9dce752f29cb
# ╟─052dc502-1c7a-11eb-2316-d3a1eef2af94
# ╠═c3090e4a-1a09-11eb-0f32-d3bbfd9992e0
# ╟─55f475a8-1c7a-11eb-377e-91d07fa0bdb6
# ╠═63ef21c6-1c7a-11eb-2f3c-c5ac16bc289f
# ╟─6cf7df1a-1c7a-11eb-230b-df1333f191c7
# ╠═19cf420e-1c7c-11eb-1cb8-dd939fee1276
# ╠═b8cd4112-1c7c-11eb-3b2d-29170ad9beb5
# ╠═a99c40bc-1c7c-11eb-036b-7fe6e9b937e5
# ╟─e9c6a0b8-1ad0-11eb-1606-0319caf0948a
# ╟─522e6b22-194d-11eb-167c-052e65f6b703
# ╟─dad5acfa-194c-11eb-27f9-01f40342a681
# ╠═43306bd4-194d-11eb-2e30-07eabb8b29ef
# ╟─70b8401e-1c7e-11eb-16b2-d54d8f66d71a
# ╠═48668080-22b4-11eb-3c4f-eb6aa1169382
# ╠═79532662-1c7e-11eb-2edf-57e7cfbc1eda
# ╟─b6614d80-194b-11eb-1edb-dba3c29672f8
# ╠═2c6defd0-1ca1-11eb-17db-d5cb498f3265
# ╟─ad5a7420-1c7f-11eb-042f-115a9ef4c676
# ╠═0b03316c-1c80-11eb-347c-1b5c9a0ae379
# ╠═fb70cc0c-1c7f-11eb-31b5-87b168a66e19
# ╟─3f727a2c-1c80-11eb-3608-e55ccb9786d9
# ╠═76ef6e46-1a06-11eb-03e3-9f40a86dc9aa
# ╟─a45e1012-194d-11eb-3252-bb89daed3c8d
# ╟─7ba5dda0-1ad1-11eb-1c4e-2391c11f54b3
# ╠═8549346e-22c0-11eb-20a3-b9964fff65e3
# ╠═3cd36ac0-1a09-11eb-1818-75b36e67594a
# ╠═1ee0787e-1a08-11eb-233b-43a654f70117
# ╟─7478330a-1c81-11eb-2f9f-099f1111032c
# ╟─ba0a869a-1ad1-11eb-091f-916e9151f052
# ╠═3aa539ce-193f-11eb-2a0f-bbc6b83528b7
# ╟─caa98732-19cd-11eb-04ce-2f018275cf01
# ╠═e5c0e960-19cc-11eb-107d-39b397a783ab
# ╠═2a2b7284-1ade-11eb-3b71-d17fe2ca638a
# ╟─e2a8d1d6-1add-11eb-0da1-cda1492a950c
# ╟─337918f4-194f-11eb-0b45-b13fef3b23bf
# ╟─492b257a-194f-11eb-17fb-f770b4d3da2e
# ╠═392fe192-1ca1-11eb-36c4-f9bd2b01a5e5
# ╠═251f0262-1a0c-11eb-39a3-09be67091dc8
# ╠═83aa9cea-1a0c-11eb-281d-699665da2b4f
# ╠═af5c6bea-1c9c-11eb-35ae-250337e4fc86
# ╠═b3ab93d2-1a0b-11eb-0f5a-cdca19af3d89
# ╟─71dc652e-1c9c-11eb-396c-cfd9ee2261fe
# ╟─584ce620-1935-11eb-177a-f75d9ad8a399
# ╟─78915326-1937-11eb-014f-fff29b3660a0
# ╠═14dc73d2-1a0d-11eb-1a3c-0f793e74da9b
# ╟─71b70da6-193e-11eb-0bc4-f309d24fd4ef
# ╟─54b81de0-193f-11eb-004d-f90ec43588f8
# ╠═6fdf613c-193f-11eb-0029-957541d2ed4d
# ╟─392c25b8-1add-11eb-225d-49cfca27bef4
# ╟─c25caf08-1a13-11eb-3c4d-0567faf4e662
# ╠═427747d6-1ca1-11eb-28ae-ff50728c84fe
# ╟─dced1fd0-1c9e-11eb-3226-17dc1e09e018
# ╠═65aec4fc-1c9e-11eb-1c5a-6dd7c533d3b8
# ╠═5895d9ae-1c9e-11eb-2f4e-671f2a7a0150
# ╠═83acf10e-1c9e-11eb-3426-bb28e7bc6c79
# ╟─13fef49c-1c9e-11eb-2aa3-d3aa2bfd0d57
# ╠═c492a1f8-1a0c-11eb-2c38-5921c39cf5f8
# ╟─155aa5e0-22c3-11eb-2596-5b77fbd02fe0
# ╟─c00eb0a6-cab2-11ea-3887-070ebd8d56e2
# ╟─3dd0a48c-1ca3-11eb-1127-e7c43b5d1666
# ╠═b314f676-233d-11eb-1858-7b5e101bfffc
# ╠═270762e4-1ca4-11eb-2fb4-392e5c3b3e04
# ╟─bbf730c8-1ca6-11eb-3bb0-1188046339ac
# ╟─cbd8f164-1ca6-11eb-1440-bdaabf73a6c7
# ╟─ebd05bf0-19c3-11eb-2559-7d0745a84025
# ╟─ec275590-19c3-11eb-23d0-cb3d9f62ba92
# ╟─ec31dce0-19c3-11eb-1487-23cc20cd5277
# ╟─ec3ed530-19c3-11eb-10bb-a55e77550d1f
# ╟─ec4abc12-19c3-11eb-1ca4-b5e9d3cd100b
# ╟─ec57b460-19c3-11eb-2142-07cf28dcf02b
# ╟─ec5d59b0-19c3-11eb-0206-cbd1a5415c28
# ╟─ec698eb0-19c3-11eb-340a-e319abb8ebb5
# ╟─ec7638e0-19c3-11eb-1ca1-0b3aa3b40240
# ╟─ec85c940-19c3-11eb-3375-a90735beaec1
# ╟─8cfa4902-1ad3-11eb-03a1-736898ff9cef


text:
include("broken_spines.jl")
include("density.jl")
include("clustering.jl")
include("visualize.jl")

#SBMs
include("sbms.jl")
include("poissonsbm.jl")
include("sbmopt.jl")

ext: jl


text:
@testset "Paraboloid" begin

P = Paraboloid()

@test P.id == 1
@test P.h0 == 0.0
@test P.h1 == 1.0
@test P.radius == 1.0
@test GRAYTR.can_intersect(P)

# Create with transformation
P = Paraboloid(GRAYTR.scaling(2.0, 1.0, 1.0))
@test GRAYTR.obj_bounds(P) == GRAYTR.BoundingBox(Point3(-1, -1, 0), Point3(1, 1, 1))
@test GRAYTR.world_bounds(P) == GRAYTR.BoundingBox(Point3(-2, -1, 0), Point3(2, 1, 1))

# Transform after creation
T = GRAYTR.scaling(1.0, 2.0, 1.0)
P = T(Paraboloid())
@test GRAYTR.obj_bounds(P) == GRAYTR.BoundingBox(Point3(-1, -1, 0), Point3(1, 1, 1))
@test GRAYTR.world_bounds(P) == GRAYTR.BoundingBox(Point3(-1, -2, 0), Point3(1, 2, 1))

# intersectP
D = Paraboloid(translation(0, 0, -1))
R = GRAYTR.Ray(Point3(0, 0, 1), Vector3(0, 0, -1))
@test GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0, 0, 1), Vector3(0, 0, -1), 10, Inf, 1)
@test !GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0, 2, 1), Vector3(0, 0, -1))
@test !GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0.9*sqrt(2), 0.9*sqrt(2), 1), Vector3(0, 0, -1))
@test !GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0, 0, -2), Vector3(0, 0, 1))
@test GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0, 0, 0), Vector3(0, 1, 0))
@test GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(3, 0, 1), Vector3(-1, 0, 0))
@test !GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(3, 0, -2), Vector3(-1, 0, 0))
@test !GRAYTR.intersectP(R, D)

# intersect

R = GRAYTR.Ray(Point3(0, 0, 1), Vector3(0, 0, -1))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg.p ≈ Point3(0, 0, -1)
@test dg.n ≈ Normal3(0, 0, 1)
@test t ≈ 2.0

R = GRAYTR.Ray(Point3(0, 0, -3), Vector3(0, 0, 1))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg.p ≈ Point3(0, 0, -1)
@test dg.n ≈ Normal3(0, 0, -1)
@test t ≈ 2.0

R = GRAYTR.Ray(Point3(0.9*sqrt(2), 0.9*sqrt(2), 1), Vector3(0, 0, -1))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg == nothing
@test isnan(t)
@test isnan(e)

R = GRAYTR.Ray(Point3(1, 0, -0.5), Vector3(-1, 0, 0))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg.p ≈ Point3(1/sqrt(2), 0.0, -0.5)
@test dg.n ≈ Normal3(sqrt(2), 0.0, -1.0) / sqrt(3)
@test t ≈ 1 - 1/sqrt(2)

R = GRAYTR.Ray(Point3(0, 0, 1), Vector3(0, 0, -1), 10, Inf, 1)
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg == nothing
@test isnan(t)
@test isnan(e)

R = GRAYTR.Ray(Point3(3, 0, 1), Vector3(-1, 0, 0))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg == nothing
@test isnan(t)
@test isnan(e)

end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
module TestUtil

using FactCheck

importall Orchestra.Util

include("fixture_learners.jl")
using .FixtureLearners
nfcp = NumericFeatureClassification()

facts("Orchestra util functions") do
  context("holdout returns proportional partitions") do
    n = 10
    right_prop = 0.3
    (left, right) = holdout(n, right_prop)

    @fact size(left, 1) => n - (n * right_prop)
    @fact size(right, 1) => n * right_prop
    @fact intersect(left, right) => isempty
    @fact size(union(left, right), 1) => n
  end

  context("kfold returns k partitions") do
    num_instances = 10
    num_partitions = 3
    partitions = kfold(num_instances, num_partitions)

    @fact size(partitions, 1) => num_partitions
    [@fact length(partition) >= 6 => true for partition in partitions]
  end

  context("score calculates accuracy") do
    learner = PerfectScoreLearner({:problem => nfcp})
    predictions = fit_and_transform!(learner, nfcp)

    @fact score(
      :accuracy, nfcp.test_labels, predictions
    ) => 100.0
  end
  context("score throws exception on unknown metric") do
    learner = PerfectScoreLearner({:problem => nfcp})
    predictions = fit_and_transform!(learner, nfcp)

    @fact_throws score(
      :fake, nfcp.test_labels, predictions
    )
  end

  context("infer_eltype returns inferred elements type") do
    vector = [1,2,3,"a"]
    @fact infer_eltype(vector[1:3]) => Int
  end

  context("nested_dict_to_tuples produces list of tuples") do
    nested_dict = {
      :a => [1,2],
      :b => {
        :c => [3,4,5]
      }
    }
    expected_set = Set({
      ([:a], [1,2]),
      ([:b,:c], [3,4,5])
    })
    set = nested_dict_to_tuples(nested_dict)

    @fact set => expected_set
  end

  context("nested_dict_set! assigns values") do
    nested_dict = {
      :a => 1,
      :b => {
        :c => 2
      }
    }
    expected_dict = {
      :a => 1,
      :b => {
        :c => 3
      }
    }
    nested_dict_set!(nested_dict, [:b,:c], 3)

    @fact nested_dict => expected_dict
  end

  context("nested_dict_merge merges two nested dictionaries") do
    first = {
      :a => 1,
      :b => {
        :c => 2,
        :d => 3
      }
    }
    second = {
      :a => 4,
      :b => {
        :d => 5
      }
    }
    expected = {
      :a => 4,
      :b => {
        :c => 2,
        :d => 5
      }
    }
    actual = nested_dict_merge(first, second)

    @fact actual => expected
  end

  context("create_transformer produces new transformer") do
    learner = AlwaysSameLabelLearner({:label => :a})
    new_options = {:label => :b}
    new_learner = create_transformer(learner, new_options)

    @fact learner.options[:label] => :a
    @fact new_learner.options[:label] => :b
    @fact true => !isequal(learner, new_learner)
  end
end

end # module


text:
function AUKFTUIntermediate(T::Type, num_x::Number)
    xi_temp = Vector{T}(undef, num_x)
    UKFTUIntermediate(
        Augmented(Matrix{T}(undef, num_x, num_x), Matrix{T}(undef, num_x, num_x)),
        Augmented(xi_temp, xi_temp),
        TransformedSigmaPoints(Vector{T}(undef, num_x), Matrix{T}(undef, num_x, 4 * num_x), MeanSetWeightingParameters(0.0)), # Weighting parameters will be reset
        TransformedSigmaPoints(Vector{T}(undef, num_x), Matrix{T}(undef, num_x, 4 * num_x), MeanSetWeightingParameters(0.0)),
        Vector{T}(undef, num_x),
        Matrix{T}(undef, num_x, num_x)
    )
end

AUKFTUIntermediate(num_x::Number) = AUKFTUIntermediate(Float64, num_x)

function AUKFMUIntermediate(T::Type, num_x::Number, num_y::Number)
    UKFMUIntermediate(
        Augmented(Matrix{T}(undef, num_x, num_x), Matrix{T}(undef, num_y, num_y)),
        Augmented(Vector{T}(undef, num_x), Vector{T}(undef, num_y)),
        Vector{T}(undef, num_y),
        TransformedSigmaPoints(Vector{T}(undef, num_y), Matrix{T}(undef, num_y, 2 * num_x + 2 * num_y), MeanSetWeightingParameters(0.0)), # Weighting parameters will be reset
        TransformedSigmaPoints(Vector{T}(undef, num_y), Matrix{T}(undef, num_y, 2 * num_x + 2 * num_y), MeanSetWeightingParameters(0.0)),
        Vector{T}(undef, num_y),
        Matrix{T}(undef, num_y, num_y),
        Matrix{T}(undef, num_x, num_y),
        Matrix{T}(undef, num_y, num_y),
        Matrix{T}(undef, num_x, num_y),
        Vector{T}(undef, num_x),
        Matrix{T}(undef, num_x, num_x)
    )
end

AUKFMUIntermediate(num_x::Number, num_y::Number) = AUKFMUIntermediate(Float64, num_x, num_y)

function calc_kalman_gain_and_posterior_covariance(P::Augmented, Pᵪᵧ, S, consider)
    K = calc_kalman_gain(Pᵪᵧ, S, consider)
    P_posterior = calc_posterior_covariance(P.P, Pᵪᵧ, K, consider)
    K, P_posterior
end

function calc_kalman_gain_and_posterior_covariance!(s_chol, kalman_gain, p_post, P::Augmented, Pᵪᵧ, S)
    K = calc_kalman_gain!(s_chol, kalman_gain, Pᵪᵧ, S)
    P_posterior = calc_posterior_covariance!(p_post, P.P, Pᵪᵧ, K)
    K, P_posterior
end

function time_update(x, P::Union{<:AbstractMatrix, <:Cholesky}, f, Q::Augment; weight_params::AbstractWeightingParameters = WanMerweWeightingParameters())
    time_update(x, Augmented(P, Q), f, Q, weight_params)
end

function time_update!(tu::UKFTUIntermediate, x, P::Union{<:AbstractMatrix, <:Cholesky}, f!, Q::Augment; weight_params::AbstractWeightingParameters = WanMerweWeightingParameters())
    time_update!(tu, x, Augmented(P, Q), f!, Q, weight_params)
end

function measurement_update(x, P::Union{<:AbstractMatrix, <:Cholesky}, y, h, R::Augment; weight_params::AbstractWeightingParameters = WanMerweWeightingParameters(), consider = nothing)
    measurement_update(x, Augmented(P, R), y, h, R, weight_params = weight_params, consider = consider)
end

function measurement_update!(mu::UKFMUIntermediate, x, P::Union{<:AbstractMatrix, <:Cholesky}, y, h!, R::Augment; weight_params::AbstractWeightingParameters = WanMerweWeightingParameters())
    measurement_update!(mu, x, Augmented(P, R), y, h!, R, weight_params = weight_params)
end


text:
using DaylightPlots, Plots
using Base.Test

include("sunrise.jl")
include("plots.jl")
include("tlims.jl")


text:
using Test
using Clustering
using Distances
include("test_helpers.jl")

@testset "dbscan() (DBSCAN clustering)" begin

@testset "Argument checks" begin
    Random.seed!(34568)
    @test_throws ArgumentError dbscan(randn(2, 3), 1.0, 1)
    @test_throws ArgumentError dbscan(randn(1, 1), 1.0, 1)
    @test_throws ArgumentError dbscan(randn(2, 2), 0.0, 1)
    @test_throws ArgumentError dbscan(randn(2, 2), 1.0, 0)
    @test dbscan(randn(2, 2), 0.5, 1) isa DbscanResult
end

Random.seed!(34568)

X1 = randn(2, 200) .+ [0., 5.]
X2 = randn(2, 200) .+ [-5., 0.]
X3 = randn(2, 200) .+ [5., 0.]
X = hcat(X1, X2, X3)
n = size(X,2)

D = pairwise(Euclidean(), X, dims=2)

R = dbscan(D, 1.0, 10)
@test isa(R, DbscanResult)
k = length(R.seeds)
# println("k = $k")
@test k == 3
@test all(R.assignments .<= k)
@test length(R.assignments) == n
@test length(R.counts) == k
for c = 1:k
    @test count(==(c), R.assignments) == R.counts[c]
end
@test all(R.counts .>= 180)

@testset "Support for arrays other than Matrix{T}" begin
    @testset "$(typeof(M))" for M in equivalent_matrices(D)
        R2 = dbscan(M, 1.0, 10)  # run on complete subarray
        @test R2.assignments == R.assignments
    end
end

@testset "normal points" begin
    p0 = randn(StableRNG(0), 3, 1000)
    p1 = randn(StableRNG(1), 3, 1000) .+ [3.0, 3.0, 0.0]
    p2 = randn(StableRNG(2), 3, 1000) .+ [-3.0, -3.0, 0.0]

    points = [p0 p1 p2]

    # FIXME Current tests depend too much on a specific random sequence
    #       We need better tests, that check point coordinates rather their indices
    inds_1 = [1, 3, 4, 5, 6, 9, 10, 12, 18, 22, 26, 29, 33, 35, 36, 39, 40, 42, 43, 46, 48, 50, 51, 52, 56, 57, 58, 60, 62, 63, 65, 70, 71, 72, 73, 74, 76, 80, 81, 84, 85, 86, 90, 91, 94, 95, 97, 100, 101, 102, 104, 107, 108, 112, 113, 114, 116, 117, 118, 119, 123, 124, 125, 126, 128, 129, 130, 131, 133, 134, 135, 136, 137, 138, 142, 145, 155, 157, 159, 160, 161, 162, 167, 168, 169, 170, 172, 174, 175, 177, 180, 181, 182, 184, 185, 187, 189, 190, 191, 197, 199, 204, 205, 208, 209, 212, 215, 217, 218, 219, 221, 223, 225, 227, 228, 229, 230, 231, 237, 239, 240, 241, 247, 248, 249, 251, 254, 256, 259, 261, 264, 265, 266, 268, 274, 277, 282, 283, 284, 285, 287, 288, 289, 290, 293, 294, 295, 298, 304, 305, 307, 308, 311, 312, 316, 317, 319, 320, 321, 323, 325, 330, 335, 339, 340, 343, 344, 345, 347, 363, 364, 365, 366, 367, 370, 371, 373, 378, 383, 385, 388, 391, 393, 396, 400, 402, 403, 404, 405, 406, 409, 411, 412, 415, 416, 417, 418, 421, 422, 423, 425, 426, 427, 430, 433, 435, 440, 441, 442, 444, 448, 450, 451, 453, 462, 464, 467, 472, 473, 474, 476, 482, 484, 485, 488, 489, 490, 492, 494, 496, 497, 498, 499, 500, 501, 503, 504, 505, 506, 508, 515, 519, 520, 526, 529, 530, 531, 532, 533, 536, 537, 542, 548, 556, 559, 562, 563, 565, 566, 567, 570, 574, 575, 576, 582, 584, 587, 588, 590, 591, 598, 600, 601, 602, 603, 604, 605, 608, 609, 612, 613, 614, 617, 621, 622, 623, 625, 627, 628, 629, 635, 636, 639, 641, 647, 650, 653, 655, 657, 659, 660, 661, 662, 665, 666, 667, 670, 671, 673, 674, 675, 676, 677, 679, 681, 683, 686, 688, 691, 694, 695, 696, 699, 701, 704, 706, 708, 711, 712, 713, 715, 717, 719, 720, 723, 724, 727, 729, 730, 731, 735, 739, 740, 741, 742, 743, 744, 746, 747, 750, 751, 755, 756, 761, 770, 772, 773, 774, 775, 780, 784, 787, 788, 790, 792, 794, 797, 800, 801, 802, 805, 806, 808, 809, 813, 814, 815, 816, 817, 821, 822, 824, 826, 827, 828, 830, 832, 833, 834, 835, 837, 843, 846, 847, 848, 850, 851, 854, 855, 857, 859, 862, 863, 864, 867, 869, 870, 872, 873, 875, 876, 878, 879, 880, 881, 884, 886, 887, 888, 889, 890, 892, 894, 901, 902, 908, 909, 913, 914, 917, 918, 919, 920, 922, 924, 928, 933, 934, 935, 936, 938, 940, 941, 943, 944, 948, 949, 950, 952, 953, 954, 960, 961, 965, 966, 970, 971, 979, 980, 983, 985, 986, 987, 990, 991, 993, 996, 1000, 1339, 2143]
    inds_2 = [132, 1001, 1003, 1006, 1008, 1011, 1014, 1015, 1017, 1018, 1019, 1020, 1023, 1024, 1027, 1028, 1034, 1036, 1039, 1042, 1044, 1045, 1047, 1049, 1051, 1052, 1056, 1057, 1058, 1059, 1064, 1065, 1068, 1070, 1071, 1076, 1081, 1084, 1087, 1089, 1090, 1093, 1094, 1095, 1096, 1097, 1099, 1100, 1102, 1103, 1108, 1110, 1111, 1112, 1113, 1119, 1120, 1123, 1124, 1125, 1130, 1131, 1136, 1140, 1142, 1143, 1146, 1147, 1156, 1158, 1161, 1162, 1167, 1168, 1172, 1174, 1176, 1177, 1178, 1179, 1183, 1186, 1187, 1190, 1191, 1192, 1193, 1200, 1201, 1202, 1203, 1206, 1209, 1210, 1212, 1213, 1215, 1217, 1219, 1222, 1226, 1229, 1230, 1231, 1232, 1233, 1239, 1241, 1242, 1244, 1246, 1247, 1249, 1250, 1251, 1256, 1257, 1258, 1260, 1261, 1263, 1264, 1265, 1266, 1268, 1275, 1276, 1282, 1285, 1286, 1287, 1291, 1293, 1294, 1295, 1300, 1303, 1307, 1308, 1313, 1315, 1318, 1320, 1325, 1331, 1333, 1336, 1337, 1341, 1345, 1346, 1347, 1348, 1350, 1351, 1355, 1358, 1360, 1361, 1362, 1364, 1365, 1368, 1370, 1372, 1373, 1374, 1378, 1379, 1381, 1382, 1383, 1386, 1392, 1393, 1394, 1396, 1397, 1398, 1400, 1401, 1405, 1406, 1408, 1410, 1413, 1415, 1416, 1418, 1419, 1420, 1421, 1426, 1431, 1433, 1434, 1437, 1441, 1445, 1446, 1447, 1448, 1452, 1453, 1454, 1455, 1459, 1462, 1463, 1464, 1466, 1467, 1468, 1473, 1474, 1476, 1477, 1478, 1480, 1484, 1485, 1487, 1489, 1490, 1492, 1493, 1499, 1501, 1502, 1503, 1504, 1505, 1507, 1514, 1515, 1517, 1519, 1521, 1522, 1524, 1526, 1528, 1529, 1534, 1541, 1542, 1544, 1545, 1546, 1551, 1552, 1553, 1555, 1556, 1561, 1564, 1566, 1567, 1568, 1569, 1571, 1574, 1575, 1576, 1583, 1586, 1588, 1589, 1590, 1592, 1594, 1596, 1597, 1598, 1599, 1601, 1602, 1603, 1604, 1606, 1607, 1608, 1609, 1610, 1612, 1615, 1618, 1619, 1620, 1621, 1622, 1623, 1624, 1627, 1629, 1633, 1635, 1641, 1643, 1646, 1647, 1648, 1649, 1650, 1651, 1652, 1654, 1656, 1658, 1659, 1661, 1662, 1663, 1664, 1665, 1666, 1667, 1668, 1669, 1670, 1673, 1678, 1683, 1684, 1685, 1690, 1692, 1696, 1700, 1701, 1703, 1705, 1706, 1708, 1709, 1710, 1712, 1713, 1716, 1718, 1719, 1720, 1722, 1723, 1725, 1726, 1727, 1729, 1730, 1736, 1737, 1738, 1739, 1740, 1742, 1743, 1744, 1747, 1748, 1749, 1752, 1755, 1758, 1761, 1769, 1771, 1775, 1776, 1777, 1785, 1787, 1791, 1793, 1794, 1795, 1796, 1797, 1798, 1799, 1803, 1805, 1806, 1808, 1811, 1816, 1818, 1821, 1822, 1827, 1828, 1829, 1830, 1831, 1834, 1838, 1839, 1845, 1849, 1850, 1851, 1852, 1853, 1857, 1859, 1864, 1867, 1869, 1870, 1871, 1872, 1878, 1886, 1888, 1889, 1898, 1900, 1901, 1904, 1908, 1912, 1913, 1914, 1915, 1916, 1917, 1919, 1921, 1924, 1929, 1932, 1933, 1935, 1936, 1938, 1940, 1941, 1942, 1948, 1949, 1951, 1952, 1954, 1955, 1957, 1962, 1964, 1965, 1966, 1973, 1976, 1977, 1978, 1979, 1984, 1985, 1988, 1993, 1994, 1996, 1998, 1999]
    inds_3 = [589, 703, 2002, 2004, 2005, 2006, 2008, 2010, 2014, 2015, 2016, 2017, 2019, 2022, 2023, 2024, 2025, 2031, 2032, 2035, 2036, 2038, 2041, 2042, 2044, 2046, 2048, 2052, 2053, 2056, 2057, 2058, 2059, 2060, 2063, 2066, 2070, 2071, 2072, 2073, 2075, 2076, 2078, 2080, 2081, 2083, 2085, 2088, 2089, 2093, 2096, 2097, 2098, 2099, 2101, 2103, 2105, 2106, 2107, 2108, 2109, 2110, 2111, 2112, 2113, 2114, 2115, 2116, 2117, 2120, 2124, 2125, 2126, 2127, 2128, 2129, 2135, 2136, 2138, 2142, 2144, 2146, 2147, 2151, 2152, 2155, 2163, 2164, 2165, 2166, 2172, 2173, 2176, 2177, 2178, 2185, 2186, 2187, 2189, 2190, 2191, 2193, 2195, 2196, 2197, 2198, 2200, 2201, 2203, 2204, 2205, 2211, 2213, 2214, 2215, 2218, 2219, 2221, 2228, 2231, 2233, 2236, 2237, 2238, 2239, 2240, 2241, 2242, 2245, 2249, 2250, 2251, 2252, 2253, 2258, 2259, 2260, 2265, 2270, 2273, 2274, 2275, 2277, 2279, 2280, 2281, 2282, 2283, 2284, 2285, 2286, 2287, 2290, 2291, 2292, 2294, 2296, 2297, 2299, 2301, 2303, 2304, 2305, 2306, 2310, 2311, 2312, 2316, 2317, 2322, 2324, 2325, 2329, 2331, 2332, 2333, 2334, 2335, 2336, 2339, 2341, 2342, 2343, 2344, 2349, 2350, 2351, 2355, 2356, 2359, 2360, 2362, 2364, 2365, 2368, 2370, 2372, 2377, 2380, 2384, 2385, 2386, 2388, 2393, 2403, 2404, 2405, 2406, 2407, 2410, 2411, 2412, 2415, 2418, 2420, 2421, 2422, 2423, 2427, 2428, 2430, 2433, 2434, 2437, 2439, 2440, 2441, 2442, 2443, 2445, 2449, 2452, 2453, 2455, 2458, 2459, 2462, 2464, 2473, 2474, 2479, 2481, 2482, 2484, 2485, 2487, 2488, 2489, 2493, 2494, 2495, 2499, 2504, 2508, 2511, 2513, 2515, 2517, 2521, 2524, 2528, 2530, 2534, 2535, 2536, 2537, 2539, 2540, 2541, 2543, 2545, 2547, 2548, 2549, 2550, 2551, 2555, 2560, 2561, 2562, 2564, 2570, 2572, 2574, 2576, 2578, 2583, 2584, 2585, 2587, 2594, 2595, 2597, 2598, 2599, 2602, 2603, 2605, 2608, 2610, 2611, 2612, 2614, 2618, 2620, 2621, 2623, 2625, 2626, 2627, 2629, 2631, 2632, 2634, 2636, 2637, 2641, 2643, 2647, 2648, 2649, 2652, 2655, 2656, 2657, 2663, 2670, 2672, 2674, 2675, 2676, 2677, 2679, 2680, 2685, 2687, 2691, 2693, 2695, 2696, 2697, 2698, 2700, 2702, 2703, 2706, 2707, 2711, 2713, 2715, 2716, 2717, 2718, 2719, 2721, 2722, 2723, 2724, 2726, 2728, 2730, 2736, 2737, 2739, 2740, 2741, 2745, 2747, 2750, 2752, 2754, 2755, 2758, 2759, 2763, 2764, 2765, 2767, 2770, 2771, 2772, 2774, 2777, 2783, 2784, 2786, 2787, 2790, 2794, 2797, 2800, 2801, 2802, 2803, 2804, 2806, 2807, 2808, 2811, 2817, 2818, 2819, 2822, 2823, 2827, 2830, 2833, 2838, 2839, 2842, 2843, 2844, 2845, 2846, 2850, 2851, 2852, 2857, 2861, 2862, 2863, 2866, 2876, 2877, 2878, 2880, 2881, 2882, 2884, 2885, 2888, 2890, 2891, 2893, 2894, 2895, 2897, 2902, 2904, 2905, 2906, 2909, 2910, 2911, 2915, 2918, 2919, 2922, 2924, 2925, 2926, 2928, 2929, 2930, 2931, 2932, 2933, 2934, 2935, 2936, 2940, 2941, 2942, 2947, 2948, 2950, 2951, 2952, 2960, 2966, 2972, 2973, 2974, 2976, 2977, 2978, 2983, 2985, 2987, 2990, 2991, 2992, 2994, 2996]

    clusters = dbscan(points, 0.3, min_neighbors=2, min_cluster_size=100, leafsize=20)
    @test length(clusters) == 3
    @test clusters[1].core_indices == inds_1
    @test clusters[2].core_indices == inds_2
    @test clusters[3].core_indices == inds_3

    @testset "Issue #84" begin
        clusters = dbscan(convert(Matrix{Float32}, points), 0.3f0, min_neighbors=2, min_cluster_size=100, leafsize=20)
        @test length(clusters) == 3
        clusters = dbscan(convert(Matrix{Float32}, points), 0.3, min_neighbors=2, min_cluster_size=100, leafsize=20)
        @test length(clusters) == 3
    end
end

end


text:
export Fun,evaluate,values,points,extrapolate
export coefficients,ncoefficients
export integrate,differentiate,domain,space,linesum,linenorm

include("Domain.jl")
include("Space.jl")

##  Constructors

type Fun{S,T}
    coefficients::Vector{T}
    space::S
    function Fun(coeff::Vector{T},sp::S)
        @assert length(coeff)≤dimension(sp)
        new(coeff,sp)
    end
end

"""
    Fun(coefficients,space)

returns a Fun with coefficients in the space
"""
Fun(coeff::Vector,sp::Space) = Fun{typeof(sp),eltype(coeff)}(coeff,sp)
Fun{T<:Integer}(coeff::Vector{T},sp::Space) = Fun(1.0coeff,sp)

function Fun(v::Vector{Any},sp::Space)
    if isempty(v)  || all(x->isa(x,Number) && x==0,v)
        Fun(Float64[],sp)
    else
        error("Cannot convert $v to a Fun of type $sp")
    end
end

##Coefficient routines
#TODO: domainscompatible?

"""
    coefficients(fun,space)

returns the coefficients of a fun in a possibly different space
"""
function coefficients(f::Fun,msp::Space)
    #zero can always be converted
    if ncoefficients(f)==1 && f.coefficients[1]==0
        f.coefficients
    else
        coefficients(f.coefficients,space(f),msp)
    end
end
coefficients{T<:Space}(f::Fun,::Type{T})=coefficients(f,T(domain(f)))
coefficients(f::Fun)=f.coefficients
coefficients(c::Number,sp::Space)=Fun(c,sp).coefficients

##Convert routines

Base.convert{T,S}(::Type{Fun{S,T}},f::Fun{S})=Fun(convert(Vector{T},f.coefficients),f.space)
Base.convert{T,S}(::Type{Fun{S,T}},f::Fun)=Fun(Fun(convert(Vector{T},f.coefficients),f.space),S(domain(f)))  #TODO: this line is incompatible with space conversion

Base.convert{T,S}(::Type{Fun{S,T}},x::Number) =
    x==0?zeros(T,S(AnyDomain())):x*ones(T,S(AnyDomain()))
Base.convert{S}(::Type{Fun{S}},x::Number) =
    x==0?zeros(S(AnyDomain())):x*ones(S(AnyDomain()))
Base.convert{IF<:Fun}(::Type{IF},x::Number)=Fun(x)
Base.promote_rule{T,V,S}(::Type{Fun{S,T}},::Type{Fun{S,V}})=Fun{S,promote_type(T,V)}

# promotion of * to fix 0.5 bug
if VERSION ≥ v"0.5.0-rc1+1"
    Base.promote_op{N,S,T}(::typeof(*),::Type{N},::Type{Fun{S,T}}) = Fun{S,promote_type(N,T)}
    Base.promote_op{N,S,T}(::typeof(*),::Type{Matrix{N}},::Type{Matrix{Fun{S,T}}}) =
        Matrix{Fun{S,promote_type(N,T)}}
end

Base.zero(::Type{Fun})=Fun(0.)
Base.zero{T,S<:Space}(::Type{Fun{S,T}})=zeros(T,S(AnyDomain()))
Base.one{T,S<:Space}(::Type{Fun{S,T}})=ones(T,S(AnyDomain()))
for op in (:(Base.zeros),:(Base.ones))
    @eval ($op){S,T}(f::Fun{S,T})=$op(T,f.space)
end

Base.zero(f::Fun)=zeros(f)
Base.one(f::Fun)=ones(f)

Base.eltype{S,T}(::Fun{S,T})=T

setspace(v::AbstractVector,s::Space) = Fun(v,s)
setspace(f::Fun,s::Space) = Fun(f.coefficients,s)

## domain

## General routines

domain(f::Fun) = domain(f.space)
domain{T<:Fun}(v::AbstractMatrix{T}) = map(domain,v)

setdomain(f::Fun,d::Domain) = Fun(f.coefficients,setdomain(space(f),d))

for op in (:tocanonical,:tocanonicalD,:fromcanonical,:fromcanonicalD,:invfromcanonicalD)
    @eval $op(f::Fun,x...) = $op(domain(f),x...)
end

for op in (:tocanonical,:tocanonicalD)
    @eval $op(d::Domain) = $op(d,Fun(identity,d))
end
for op in (:fromcanonical,:fromcanonicalD,:invfromcanonicalD)
    @eval $op(d::Domain) = $op(d,Fun(identity,canonicaldomain(d)))
end

space(f::Fun) = f.space
spacescompatible(f::Fun,g::Fun) = spacescompatible(space(f),space(g))
canonicalspace(f::Fun) = canonicalspace(space(f))
canonicaldomain(f::Fun) = canonicaldomain(domain(f))

##Evaluation

function evaluate(f::AbstractVector,S::Space,x...)
    csp=canonicalspace(S)
    if spacescompatible(csp,S)
        error("Override evaluate for " * string(typeof(csp)))
    else
        evaluate(coefficients(f,S,csp),csp,x...)
    end
end

evaluate(f::Fun,x) = evaluate(f.coefficients,f.space,x)
evaluate(f::Fun,x,y,z...) = evaluate(f.coefficients,f.space,Vec(x,y,z...))

@compat (f::Fun)(x...)=evaluate(f,x...)

for op in (:(Base.first),:(Base.last))
    @eval $op{S,T}(f::Fun{S,T}) = f($op(domain(f)))
end

## Extrapolation

# Default extrapolation is evaluation. Override this function for extrapolation enabled spaces.
extrapolate(f::AbstractVector,S::Space,x...)=evaluate(f,S,x...)

# Do not override these
extrapolate(f::Fun,x) = extrapolate(f.coefficients,f.space,x)
extrapolate(f::Fun,x,y,z...) = extrapolate(f.coefficients,f.space,Vec(x,y,z...))

##Data routines

values(f::Fun,dat...)=itransform(f.space,f.coefficients,dat...)
points(f::Fun)=points(f.space,ncoefficients(f))
ncoefficients(f::Fun)=length(f.coefficients)

function Base.stride(f::Fun)
    # Check only for stride 2 at the moment
    # as higher stride is very rare anyways
    M=maxabs(f.coefficients)
    for k=2:2:ncoefficients(f)
        if abs(f.coefficients[k])>40*M*eps()
            return 1
        end
    end

    2
end

## Manipulate length

pad!(f::Fun,n::Integer) = (pad!(f.coefficients,n);f)
pad(f::Fun,n::Integer) = Fun(pad(f.coefficients,n),f.space)

function chop!{S,T}(f::Fun{S,T},tol::Real)
    chop!(f.coefficients,tol)
    if length(f.coefficients) == 0
        f.coefficients = [zero(T)]
    end

    f
end
chop(f::Fun,tol)=chop!(Fun(copy(f.coefficients),f.space),tol)
chop!(f::Fun)=chop!(f,eps(eltype(f.coefficients)))

## Addition and multiplication

for op in (:+,:-,:(.+),:(.-))
    @eval begin
        function $op(f::Fun,g::Fun)
            if spacescompatible(f,g)
                n = max(ncoefficients(f),ncoefficients(g))
                f2 = pad(f,n); g2 = pad(g,n)

                Fun(($op)(f2.coefficients,g2.coefficients),isambiguous(domain(f))?g.space:f.space)
            else
                m=union(f.space,g.space)
                if isa(m,NoSpace)
                    error("Cannot "*string($op)*" because no space is the union of "*string(typeof(f.space))*" and "*string(typeof(g.space)))
                end
                $op(Fun(f,m),Fun(g,m)) # convert to same space
            end
        end
        $op{S,T<:Number}(f::Fun{S,T},c::T)=c==0?f:$op(f,Fun(c))
        $op(f::Fun,c::Number)=$op(f,Fun(c))
        $op(f::Fun,c::UniformScaling)=$op(f,c.λ)
        $op(c::UniformScaling,f::Fun)=$op(c.λ,f)
    end
end

# equivalent to Y+=a*X
axpy!(a,X::Fun,Y::Fun)=axpy!(a,coefficients(X,space(Y)),Y)
function axpy!(a,xcfs::Vector,Y::Fun)
    if a!=0
        n=ncoefficients(Y); m=length(xcfs)

        if n≤m
            resize!(Y.coefficients,m)
            for k=1:n
                @inbounds Y.coefficients[k]+=a*xcfs[k]
            end
            for k=n+1:m
                @inbounds Y.coefficients[k]=a*xcfs[k]
            end
        else #X is smaller
            for k=1:m
                @inbounds Y.coefficients[k]+=a*xcfs[k]
            end
        end
    end

    Y
end

-(f::Fun)=Fun(-f.coefficients,f.space)
for op in (:-,:(.-))
    @eval $op(c::Number,f::Fun)=-$op(f,c)
end

for op = (:*,:.*,:./,:/)
    @eval $op(f::Fun,c::Number) = Fun($op(f.coefficients,c),f.space)
end

for op = (:*,:.*,:+,:(.+))
    @eval $op(c::Number,f::Fun)=$op(f,c)
end

function .^{S,T}(f::Fun{S,T},k::Integer)
    if k == 0
        ones(space(f))
    elseif k==1
        f
    elseif k > 1
        f.*f.^(k-1)
    else
        1./f.^(-k)
    end
end

Base.inv{S,T}(f::Fun{S,T})=1./f

# Integrals over two Funs, which are fast with the orthogonal weight.

export bilinearform, linebilinearform, innerproduct, lineinnerproduct

# Having fallbacks allow for the fast implementations.

defaultbilinearform(f::Fun,g::Fun)=sum(f.*g)
defaultlinebilinearform(f::Fun,g::Fun)=linesum(f.*g)

bilinearform(f::Fun,g::Fun)=defaultbilinearform(f,g)
bilinearform(c::Number,g::Fun)=sum(c*g)
bilinearform(g::Fun,c::Number)=sum(g*c)

linebilinearform(f::Fun,g::Fun)=defaultbilinearform(f,g)
linebilinearform(c::Number,g::Fun)=linesum(c*g)
linebilinearform(g::Fun,c::Number)=linesum(g*c)

# Conjugations

innerproduct(f::Fun,g::Fun)=bilinearform(conj(f),g)
innerproduct(c::Number,g::Fun)=bilinearform(conj(c),g)
innerproduct(g::Fun,c::Number)=bilinearform(conj(g),c)

lineinnerproduct(f::Fun,g::Fun)=linebilinearform(conj(f),g)
lineinnerproduct(c::Number,g::Fun)=linebilinearform(conj(c),g)
lineinnerproduct(g::Fun,c::Number)=linebilinearform(conj(g),c)

## Norm

for (OP,SUM) in ((:(Base.norm),:(Base.sum)),(:linenorm,:linesum))
    @eval begin
        $OP(f::Fun) = $OP(f,2)

        function $OP(f::Fun,p::Number)
            if p < 1
                return error("p should be 1 ≤ p ≤ ∞")
            elseif 1 ≤ p < Inf
                return abs($SUM(abs2(f)^(p/2)))^(1/p)
            else
                return maxabs(f)
            end
        end

        function $OP(f::Fun,p::Int)
            if 1 ≤ p < Inf
                return iseven(p) ? abs($SUM(abs2(f)^div(p,2)))^(1/p) : abs($SUM(abs2(f)^(p/2)))^(1/p)
            else
                return error("p should be 1 ≤ p ≤ ∞")
            end
        end
    end
end

## Mapped functions

Base.transpose(f::Fun) = f  # default no-op

for op = (:(Base.real),:(Base.imag),:(Base.conj))
    @eval ($op){S<:Space{RealBasis}}(f::Fun{S}) = Fun(($op)(f.coefficients),f.space)
end

Base.conj(f::Fun)=error("Override conj for $(typeof(f))")

Base.abs2{S<:Space{RealBasis},T<:Real}(f::Fun{S,T})=f^2
Base.abs2{S<:Space{RealBasis},T<:Complex}(f::Fun{S,T})=real(f)^2+imag(f)^2
Base.abs2(f::Fun)=f*conj(f)

##  integration

function Base.cumsum(f::Fun)
    cf = integrate(f)
    cf - first(cf)
end

Base.cumsum(f::Fun,d::Domain)=cumsum(Fun(f,d))
Base.cumsum(f::Fun,d)=cumsum(f,Domain(d))

function differentiate(f::Fun,k::Integer)
    @assert k >= 0
    (k==0)?f:differentiate(differentiate(f),k-1)
end

==(f::Fun,g::Fun) =  (f.coefficients == g.coefficients && f.space == g.space)
function Base.isapprox(f::Fun,g::Fun)
    if spacescompatible(f,g)
        m=min(ncoefficients(f),ncoefficients(g))
        tol=100eps()  # TODO: normalize by norm of f/g

        for k=1:m
            if !isapprox(f.coefficients[k],g.coefficients[k])
                return false
            end
        end
        for k=m+1:ncoefficients(f)
            if abs(f.coefficients[k])>tol
                return false
            end
        end
        for k=m+1:ncoefficients(g)
            if abs(g.coefficients[k])>tol
                return false
            end
        end

        true
    else
        sp=union(f.space,g.space)
        if isa(sp,NoSpace)
            false
        else
            isapprox(Fun(f,sp),Fun(g,sp))
        end
    end
end

Base.isapprox(f::Fun,g::Number)=isapprox(f,g*ones(space(f)))
Base.isapprox(g::Number,f::Fun)=isapprox(g*ones(space(f)),f)

Base.isreal{S,T<:Real}(f::Fun{S,T})=basistype(S)<:RealBasis
Base.isreal(f::Fun)=false

# sum, integrate, and idfferentiate are in CalculusOperator

function reverseorientation(f::Fun)
    csp=canonicalspace(f)
    if spacescompatible(csp,space(f))
        error("Implement reverseorientation for $(typeof(f))")
    else
        reverseorientation(Fun(f,csp))
    end
end

## non-vector notation

*(f::Fun,g::Fun)=f.*g
^(f::Fun,k::Integer)=f.^k
^(f::Fun,k::Union{Number,Fun})=f.^k
/(c::Union{Number,Fun},g::Fun)=c./g

include("constructors.jl")


text:
### SOURCE BEGIN

function Data.schema(jl_result::Result)
    Data.schema(jl_result, Data.Field)
end

function Data.schema(jl_result::Result, ::Type{Data.Field})
    types = map(jl_result.not_null, column_types(jl_result)) do not_null, col_type
        not_null ? col_type : Union{col_type, Missing}
    end

    return Data.Schema(
        types,
        column_names(jl_result),
        num_rows(jl_result),
    )
end

function Data.isdone(jl_result::Result, row, col)
    Data.isdone(jl_result, row, col, num_rows(jl_result), num_columns(jl_result))
end

function Data.isdone(jl_result::Result, row, col, rows, cols)
    row > rows || col > cols
end

Data.streamtype(::Type{Result}, ::Type{Data.Field}) = true
Data.accesspattern(jl_result::Result) = RandomAccess()

# function Data.streamfrom(
#     jl_result::Result,
#     ::Type{Data.Field},
#     ::Type{Union{T, Missing}},
#     row::Int,
#     col::Int,
# )::Union{T, Missing} where T
#     if libpq_c.PQgetisnull(jl_result.result, row - 1, col - 1) == 1
#         return missing
#     else
#         oid = jl_result.column_oids[col]
#         return jl_result.column_funcs[col](PQValue{oid}(jl_result, row, col))::T
#     end
# end

@inline _non_null_type(::Type{Union{T, Missing}}) where {T} = T
@inline _non_null_type(::Type{T}) where {T} = T

# allow types that aren't just unions to handle nulls
function Data.streamfrom(
    jl_result::Result,
    ::Type{Data.Field},
    ::Type{T},
    row::Int,
    col::Int,
)::T where T>:Missing
    NNT = _non_null_type(T)

    if libpq_c.PQgetisnull(jl_result.result, row - 1, col - 1) == 1
        return missing
    else
        oid = jl_result.column_oids[col]
        return jl_result.column_funcs[col](PQValue{oid}(jl_result, row, col))::NNT
    end
end

# if a user says they don't want Missing, error on NULL
function Data.streamfrom(
    jl_result::Result,
    ::Type{Data.Field},
    ::Type{T},
    row::Int,
    col::Int,
)::T where T
    if libpq_c.PQgetisnull(jl_result.result, row - 1, col - 1) == 1
        error("Unexpected NULL at column $col row $row")
    end

    oid = jl_result.column_oids[col]
    return jl_result.column_funcs[col](PQValue{oid}(jl_result, row, col))
end

### SOURCE END

### SINK BEGIN

"""
    Statement(sch::Data.Schema, ::Type{Data.Row}, append, connection::Connection, query::AbstractString) -> Statement

Construct a `Statement` for use in streaming with DataStreams.
This function is called by `Data.stream!(source, Statement, connection, query)`.
"""
function Statement(
    sch::Data.Schema,
    ::Type{Data.Row},
    append::Bool,  # ignored
    connection::Connection,
    query::AbstractString,
)
    return prepare(connection, query)
end

Data.weakrefstrings(::Type{<:Statement}) = false
Data.streamtypes(::Type{<:Statement}) = [Data.Row]

function Data.streamto!(sink::Statement, ::Type{Data.Row}, row, row_num, col_num)
    parameters = Vector{Parameter}(length(row))

    # this should change to be whatever custom pgtype conversion function we invent
    map!(parameters, values(row)) do val
        if ismissing(val)
            missing
        elseif val isa AbstractString
            convert(String, val)
        else
            string(val)
        end
    end

    clear!(execute(sink, parameters; throw_error=true))
end

### SINK END

### FETCH BEGIN

# fetch! is not part of the DataSreams API

"""
    fetch!(sink::Union{T, Type{T}}, result::Result, args...; kwargs...) where {T} -> T

Stream data to `sink` or a new structure of type T using [`Data.stream!`](https://juliadata.github.io/DataStreams.jl/stable/#Data.stream!-1).
Any trailing `args` or `kwargs` are passed to `Data.stream!`.
`result` is cleared upon completion.
"""
function fetch!(sink, result::Result, args...; kwargs...)
    if result.cleared
        error(LOGGER, "Cannot fetch a cleared Result")
    end

    data = Data.stream!(result, sink, args...; kwargs...)
    clear!(result)
    return Data.close!(data)
end

### FETCH END


text:
module LuxurySparse

using LinearAlgebra, SparseArrays, Random
using StaticArrays: SVector, SMatrix, SDiagonal, SArray
using Base: @propagate_inbounds

import Base: copyto!, *, kron, -
import LinearAlgebra: ishermitian
import SparseArrays: SparseMatrixCSC, nnz, nonzeros, dropzeros!, findnz, issparse
import Base: getindex, size, similar, copy, show

export I, fast_invperm, isdense, allocated_coo

include("Core.jl")
include("IMatrix.jl")
include("PermMatrix.jl")
include("SparseMatrixCOO.jl")
include("SSparseMatrixCSC.jl")

include("conversions.jl")
include("promotions.jl")
include("staticize.jl")
include("arraymath.jl")
include("linalg.jl")
include("kronecker.jl")

end


text:
using BenchmarkTools
using KernelFunctions

X = rand(10, 20)
Xc = ColVecs(X)
Xr = RowVecs(X)
Xv = collect.(eachcol(X))
Y = rand(10, 20)
Yc = ColVecs(Y)
Yr = RowVecs(Y)
Yv = collect.(eachcol(Y))
SUITE = BenchmarkGroup()

kernels = Dict(
    "SqExponential" => SqExponentialKernel(),
    "Exponential" => ExponentialKernel(),
)

inputtypes = Dict(
    "ColVecs" => (Xc, Yc),
    "RowVecs" => (Xr, Yr),
    "Vecs" => (Xv, Yv),
)

functions = Dict(
    "kernelmatrixX" => (kernel, X, Y) -> kernelmatrix(kernel, X),
    "kernelmatrixXY" => (kernel, X, Y) -> kernelmatrix(kernel, X, Y),
    "kernelmatrix_diagX" => (kernel, X, Y) -> kernelmatrix_diag(kernel, X),
    "kernelmatrix_diagXY" => (kernel, X, Y) -> kernelmatrix_diag(kernel, X, Y)
)

for (kname, kernel) in kernels
    SUITE[kname] = sk = BenchmarkGroup()
    for (inputname, (X, Y)) in inputtypes
        sk[inputname] = si = BenchmarkGroup()
        for (fname, f) in functions
            si[fname] = @benchmarkable $f($kernel, $X, $Y)
        end
    end
end

# Uncomment the following to run benchmark locally

# tune!(SUITE)

# results = run(SUITE, verbose=true)


text:
export preprocess_binarize

using mlpack._Internal.io

import mlpack_jll
const preprocess_binarizeLibrary = mlpack_jll.libmlpack_julia_preprocess_binarize

# Call the C binding of the mlpack preprocess_binarize binding.
function preprocess_binarize_mlpackMain()
  success = ccall((:preprocess_binarize, preprocess_binarizeLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module preprocess_binarize_internal
  import ..preprocess_binarizeLibrary

end # module

"""
    preprocess_binarize(input; [dimension, threshold, verbose])

This utility takes a dataset and binarizes the variables into either 0 or 1
given threshold. User can apply binarization on a dimension or the whole
dataset.  The dimension to apply binarization to can be specified using the
`dimension` parameter; if left unspecified, every dimension will be binarized.
The threshold for binarization can also be specified with the `threshold`
parameter; the default threshold is 0.0.

The binarized matrix may be saved with the `output` output parameter.

For example, if we want to set all variables greater than 5 in the dataset `X`
to 1 and variables less than or equal to 5.0 to 0, and save the result to `Y`,
we could run

```julia
julia> using CSV
julia> X = CSV.read("X.csv")
julia> Y = preprocess_binarize(X; threshold=5)
```

But if we want to apply this to only the first (0th) dimension of `X`,  we could
instead run

```julia
julia> using CSV
julia> X = CSV.read("X.csv")
julia> Y = preprocess_binarize(X; dimension=0, threshold=5)
```

# Arguments

 - `input::Array{Float64, 2}`: Input data matrix.
 - `dimension::Int`: Dimension to apply the binarization. If not set, the
      program will binarize every dimension by default.  Default value `0`.

 - `threshold::Float64`: Threshold to be applied for binarization. If not
      set, the threshold defaults to 0.0.  Default value `0`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output::Array{Float64, 2}`: Matrix in which to save the output.

"""
function preprocess_binarize(input;
                             dimension::Union{Int, Missing} = missing,
                             threshold::Union{Float64, Missing} = missing,
                             verbose::Union{Bool, Missing} = missing,
                             points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, preprocess_binarizeLibrary), Nothing, ());

  IORestoreSettings("Binarize Data")

  # Process each input argument before calling mlpackMain().
  IOSetParamMat("input", input, points_are_rows)
  if !ismissing(dimension)
    IOSetParam("dimension", convert(Int, dimension))
  end
  if !ismissing(threshold)
    IOSetParam("threshold", convert(Float64, threshold))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output")
  # Call the program.
  preprocess_binarize_mlpackMain()

  return IOGetParamMat("output", points_are_rows)
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api.concept;
#
# import typedb.client.api.TypeDBTransaction;
# import typedb.client.api.concept.thing.Attribute;
# import typedb.client.api.concept.thing.Entity;
# import typedb.client.api.concept.thing.Relation;
# import typedb.client.api.concept.thing.Thing;
# import typedb.client.api.concept.type.AttributeType;
# import typedb.client.api.concept.type.EntityType;
# import typedb.client.api.concept.type.RelationType;
# import typedb.client.api.concept.type.RoleType;
# import typedb.client.api.concept.type.ThingType;
# import typedb.client.api.concept.type.Type;
#
# import javax.annotation.CheckReturnValue;
#
# public interface Concept {
#
#     @CheckReturnValue
#     default boolean isType() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isThingType() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isEntityType() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isAttributeType() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isRelationType() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isRoleType() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isThing() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isEntity() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isAttribute() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isRelation() {
#         return false;
#     }
#
#     @CheckReturnValue
#     Type asType();
#
#     @CheckReturnValue
#     ThingType asThingType();
#
#     @CheckReturnValue
#     EntityType asEntityType();
#
#     @CheckReturnValue
#     AttributeType asAttributeType();
#
#     @CheckReturnValue
#     RelationType asRelationType();
#
#     @CheckReturnValue
#     RoleType asRoleType();
#
#     @CheckReturnValue
#     Thing asThing();
#
#     @CheckReturnValue
#     Entity asEntity();
#
#     @CheckReturnValue
#     Attribute<?> asAttribute();
#
#     @CheckReturnValue
#     Relation asRelation();
#
#     @CheckReturnValue
#     Remote asRemote(TypeDBTransaction transaction);
#
#     @CheckReturnValue
#     boolean isRemote();
#
#     interface Remote extends Concept {
#
#         void delete();
#
#         @CheckReturnValue
#         boolean isDeleted();
#
#         @Override
#         @CheckReturnValue
#         Type.Remote asType();
#
#         @Override
#         @CheckReturnValue
#         ThingType.Remote asThingType();
#
#         @Override
#         @CheckReturnValue
#         EntityType.Remote asEntityType();
#
#         @Override
#         @CheckReturnValue
#         RelationType.Remote asRelationType();
#
#         @Override
#         @CheckReturnValue
#         AttributeType.Remote asAttributeType();
#
#         @Override
#         @CheckReturnValue
#         RoleType.Remote asRoleType();
#
#         @Override
#         @CheckReturnValue
#         Thing.Remote asThing();
#
#         @Override
#         @CheckReturnValue
#         Entity.Remote asEntity();
#
#         @Override
#         @CheckReturnValue
#         Relation.Remote asRelation();
#
#         @Override
#         @CheckReturnValue
#         Attribute.Remote<?> asAttribute();
#     }
# }


text:
immutable StingerGraph
    max_nv::Int64
    max_neblock::Int64
    max_netypes::Int64
    max_nvtypes::Int64
    num_insertions::Int64
    num_deletions::Int64
    num_insertions_last_batch::Int64
    num_deletions_last_batch::Int64

    batch_time::Float64
    update_time::Float64

    queue_size::Int64
    dropped_batches::Int64
    vertices_start::Int64
    physmap_start::Int64
    etype_names_start::Int64
    vtype_names_start::Int64

    ETA_start::Int64
    ebpool_start::Int64
    size_t::UInt

    #cache_pad::Array{Int64}(5)
    #storage - Zero sized array (use `storageptr` to get a pointer to this)
end

export max_nv, max_neblock, max_netypes, max_nvtypes, num_insertions, num_deletions,
num_insertions_last_batch, num_deletions_last_batch, batch_time, update_time, queue_size,
dropped_batches, vertices_start, physmap_start, etype_names_start, vtype_names_start,
ETA_start, ebpool_start, size_t

const stingergraphfields = fieldnames(StingerGraph)

@enum StingerFields max_nv=1 max_neblock=2 max_netypes=3 max_nvtypes=4 num_insertions=5 num_deletions=6 num_insertions_last_batch=7 num_deletions_last_batch=8 batch_time=9 update_time=10 queue_size=11 dropped_batches=12 vertices_start=13 physmap_start=14 etype_names_start=15 vtype_names_start=16 ETA_start=17 ebpool_start=18 size_t=19


text:
import .TikzPictures: TikzPicture

# External library import.
import .TikZ

""" Convert our TikZ picture type to `TikzPicture`'s picture type.
"""
function TikzPicture(pic::TikZ.Picture; preamble::String="")::TikzPicture
  data = join(sprint(TikZ.pprint, stmt) for stmt in pic.stmts)
  options = join((sprint(TikZ.pprint, prop) for prop in pic.props), ",")
  preamble = join([
    preamble,
    "\\usepackage{amssymb}",
    # FIXME: These TikZ library dependencies should be stored in TikZ.Picture.
    "\\usetikzlibrary{arrows.meta}",
    "\\usetikzlibrary{calc}",
    "\\usetikzlibrary{decorations.markings}",
    "\\usetikzlibrary{positioning}",
    "\\usetikzlibrary{shapes.geometric}",
  ], "\n")
  TikzPicture(data; options=options, preamble=preamble)
end

function Base.show(io::IO, ::MIME"image/svg+xml", pic::TikZ.Picture)
  show(io, MIME"image/svg+xml"(), TikzPicture(pic))
end


text:
module UseGA_AEM
using Libdl
using Cxx

mutable struct EMforward
    TLM  :: Cxx.CxxCore.CppPtr
    THM  :: Cxx.CxxCore.CppPtr
    PX   :: Float64
    PY   :: Float64
    PZ   :: Float64
    SXLM :: Array{Float64, 1}
    SYLM :: Array{Float64, 1}
    SZLM :: Array{Float64, 1}
    SXHM :: Array{Float64, 1}
    SYHM :: Array{Float64, 1}
    SZHM :: Array{Float64, 1}
end

function init_GA_AEM()
    #CXXJL_ROOTDIR should be set at least when first time installing Cxx
    #export CXXJL_ROOTDIR=/apps/gcc/5.2.0

    # Importing shared library and header file
    path_to_lib = "/g/data1a/zk34/ga-aem/julia"
    addHeaderDir(path_to_lib, kind=C_System)

    ga_aem_path = "/home/547/ar0754/zk34/ga-aem"

    # Julia compiler also needs to find headers in these dirs
    addHeaderDir("/g/data1a/zk34/ga-aem/src", kind=C_System)
    addHeaderDir("/g/data1a/zk34/ga-aem/submodules/cpp-utils/src/", kind=C_System)
    addHeaderDir("/apps/fftw3/3.3.7-gcc/include", kind=C_System)

    cxxinclude("gatdaem1d_julia.h")
    Libdl.dlopen(path_to_lib * "/gatdaem1d_julia.so", Libdl.RTLD_GLOBAL)

    # Create the cTDEmSystem class object
    stmfile = "/g/data1a/zk34/ga-aem/examples/bhmar-skytem/stmfiles/Skytem-LM.stm"
    TLM = @cxxnew cTDEmSystem(pointer(stmfile));
    stmfile = "/g/data1a/zk34/ga-aem/examples/bhmar-skytem/stmfiles/Skytem-HM.stm"
    THM = @cxxnew cTDEmSystem(pointer(stmfile));

    # Allocate reusable storage for outputs
    PX = 0.0
    PY = 0.0
    PZ = 0.0
    nwindowsLM = @cxx TLM->NumberOfWindows;
    nwindowsHM = @cxx THM->NumberOfWindows;
    SXLM = Array{Float64, 1}(undef,nwindowsLM)
    SYLM = Array{Float64, 1}(undef,nwindowsLM)
    SZLM = Array{Float64, 1}(undef,nwindowsLM)

    SXHM = Array{Float64, 1}(undef,nwindowsHM)
    SYHM = Array{Float64, 1}(undef,nwindowsHM)
    SZHM = Array{Float64, 1}(undef,nwindowsHM)

    # call it once bug
    # Set the input arrays
    geometryLM   = [35.0,0.0,0.0,0.0,-17.0,0.0,+2.0,0.0,0.0,0.0]
    geometryHM   = [35.0,0.0,0.0,0.0,-17.0,0.0,+0.2,0.0,0.0,0.0]
    conductivity = [0.01, 1.0, 0.001]
    thickness    = [20.0, 40.0]
    nlayers      = length(conductivity)

    @cxx TLM->forwardmodel(nlayers,pointer(conductivity),pointer(thickness),
                           pointer(geometryLM),Ref(PX),Ref(PY),Ref(PZ),
                           pointer(SXLM),pointer(SYLM),pointer(SZLM))
    @cxx THM->forwardmodel(nlayers,pointer(conductivity),pointer(thickness),
                           pointer(geometryHM),Ref(PX),Ref(PY),Ref(PZ),
                           pointer(SXHM),pointer(SYHM),pointer(SZHM))

    EMforward(TLM, THM, PX, PY, PZ, SXLM, SYLM, SZLM, SXHM, SYHM, SZHM)
end

function(em::EMforward)()

    # Set the input arrays
    geometryLM   = [35.0,0.0,0.0,0.0,-17.0,0.0,+2.0,0.0,0.0,0.0]
    geometryHM   = [35.0,0.0,0.0,0.0,-17.0,0.0,+0.2,0.0,0.0,0.0]
    conductivity = [0.01, 1.0, 0.001]
    thickness    = [20.0, 40.0]
    nlayers      = length(conductivity)

    @cxx em.TLM->forwardmodel(nlayers,pointer(conductivity),pointer(thickness),
                           pointer(geometryLM),Ref(em.PX),Ref(em.PY),Ref(em.PZ),
                           pointer(em.SXLM),pointer(em.SYLM),pointer(em.SZLM))
    @cxx em.THM->forwardmodel(nlayers,pointer(conductivity),pointer(thickness),
                           pointer(geometryHM),Ref(em.PX),Ref(em.PY),Ref(em.PZ),
                           pointer(em.SXHM),pointer(em.SYHM),pointer(em.SZHM))
end

mutable struct Geometry
    geometryLM  :: Array{Float64, 1}
    geometryHM  :: Array{Float64, 1}
end

function Geometry(;ztx  =  35.0,
                   rrx  = -17.0,
                   dzrxLM =   2.0,
                   dzrxHM =   0.2
                 )
    @assert ztx  > 0.0
    @assert rrx  < 0.0
    @assert dzrxLM > 0.0
    @assert dzrxHM > 0.0

    geometryLM   = [ztx,0.0,0.0,0.0,rrx,0.0,dzrxLM,0.0,0.0,0.0]
    geometryHM   = [ztx,0.0,0.0,0.0,rrx,0.0,dzrxHM,0.0,0.0,0.0]

    Geometry(geometryLM, geometryHM)
end

function (em::EMforward)(g::Geometry, conductivity::Array{Float64, 1}, thickness::Array{Float64, 1})
    nlayers      = length(conductivity)
    @cxx em.TLM->forwardmodel(nlayers,pointer(conductivity),pointer(thickness),
                           pointer(g.geometryLM),Ref(em.PX),Ref(em.PY),Ref(em.PZ),
                           pointer(em.SXLM),pointer(em.SYLM),pointer(em.SZLM))
    @cxx em.THM->forwardmodel(nlayers,pointer(conductivity),pointer(thickness),
                           pointer(g.geometryHM),Ref(em.PX),Ref(em.PY),Ref(em.PZ),
                           pointer(em.SXHM),pointer(em.SYHM),pointer(em.SZHM))
end

mutable struct EMoperator
    em :: EMforward
    g  :: Geometry
end

function (op::EMoperator)(ztx::Float64, conductivity::Array{Float64, 1}, thickness::Array{Float64,1 })
    op.g.geometryLM[1] = ztx
    op.g.geometryHM[1] = ztx
    op.em(op.g, conductivity, thickness)
end

mutable struct Sounding
    dataLM    :: Array{Float64, 1}
    dataHM    :: Array{Float64, 1}
    sdLM      :: Array{Float64, 1}
    sdHM      :: Array{Float64, 1}
    x         :: Array{Float64, 1}
    ztx       :: Float64
    thickness :: Array{Float64, 1}
end

end


text:
using Symbolics
using SymbolicUtils, Test
using Symbolics: symtype, shape, wrap, unwrap, Unknown, Arr, arrterm
using Base: Slice
using SymbolicUtils: Sym, term, operation

@testset "arrays" begin
    @variables X[1:5, 1:5] Y[1:5, 1:5]
    @test typeof(X) <: Arr
    @test shape(X) == Slice.((1:5, 1:5))
    @test shape(Y) == Slice.((1:5, 1:5))

    A = Y[2, :]
    @test typeof(A) <: Arr{Num, 1}
    @test axes(A) == (1:5,)

    B = A[3:5]
    @test axes(B) == (Slice(1:3),)

    i = Sym{Int}(:i)
    j = Sym{Int}(:j)
    @test symtype(X[i,j]) == Real
    @test symtype(X[1,j]) == Real
end

@testset "getindex" begin
    @variables X[1:5, 1:5] Y[1:5, 1:5]

    @test isequal(X[1,1], wrap(term(getindex, unwrap(X), 1,1)))

    XX = unwrap(X)
    @test isequal(unwrap(X[1, :]), Symbolics.@arrayop(XX[1,:], (j,), XX[1, j]))
    @test isequal(unwrap(X[:, 2]), Symbolics.@arrayop(XX[:,2], (i,), XX[i, 2]))
    @test isequal(unwrap(X[:, 2:3]), Symbolics.@arrayop(XX[:,2:3], (i,j), XX[i, j], (+), (j in 2:3)))

    @variables t x[1:4](t)
    @syms i::Int
    @test isequal(x[i], operation(unwrap(x[i]))(t))
end

getdef(v) = getmetadata(v, Symbolics.VariableDefaultValue)
@testset "broadcast & scalarize" begin
    @variables A[1:5,1:3]=42 b[1:3]=[2, 3, 5] t x[1:4](t)
    AA = Symbolics.scalarize(A)
    bb = Symbolics.scalarize(b)
    @test all(isequal(42), getdef.(AA))
    @test getdef.(bb) == [2, 3, 5]
    @test isequal(Symbolics.scalarize([b.*1; b.*1]), [bb; bb])
    @test isequal(Symbolics.scalarize(b.^1), bb)
    c = A*b

    @test isequal(collect(sin.(x)),
                  sin.([x[i] for i in 1:4]))

    @test isequal(Symbolics.scalarize(sin.(A.*c)[1,1]),
                  sin(A[1, 1]*(b[1]*A[1, 1] +
                               b[2]*A[1, 2] +
                               b[3]*A[1, 3])))

    D = Differential(t)
    @test isequal(collect(D.(x) ~ x), map(i->D(x[i]) ~ x[i], eachindex(x)))
    @test_throws ArgumentError A ~ t
end


text:
const fields = ["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid", "cid"]
const required_fields = fields[begin:end-1]
const record_pattern = Regex("(?:(?:" * join(fields, "|") * raw"):\S+(?:\s|$))+", "ms")
const field_pattern = Regex("(" * join(fields, "|") * raw"):(\S+)\s?")

function validate_range(str::AbstractString, len::Int, low::Int, high::Int)::Bool
    if length(str) != len
        return false
    end
    val = parse(Int, str)
    if low <= val <= high
        return true
    end
    return false
end

function validate_byr(field::AbstractString)::Bool
    return validate_range(field, 4, 1920, 2002)
end

function validate_iyr(field::AbstractString)::Bool
    return validate_range(field, 4, 2010, 2020)
end

function validate_eyr(field::AbstractString)::Bool
    return validate_range(field, 4, 2020, 2030)
end

function validate_hgt(field::AbstractString)::Bool
    if length(field) < 3
        return false
    end
    unit = field[end-1:end]
    val = chop(field, tail = 2)
    if unit == "cm"
        return validate_range(val, 3, 150, 193)
    elseif unit == "in"
        return validate_range(val, 2, 59, 76)
    end
    return false
end

function validate_hcl(field::AbstractString)::Bool
    return occursin(r"^#[0-9a-f]{6}$", field)
end

function validate_ecl(field::AbstractString)::Bool
    return field in ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]
end

function validate_pid(field::AbstractString)::Bool
    return occursin(r"^\d{9}$", field)
end

function validate_required_fields(record)::Bool
    return all(f -> haskey(record, f), required_fields)
end

const validation_functions = Dict(f => getfield(Main, Symbol("validate_$(f)")) for f in required_fields)

function validate_passport(passport)::Bool
    if !validate_required_fields(passport)
        return false
    end

    for field in required_fields
        if !validation_functions[field](passport[field])
            return false
        end
    end

    return true
end

passports = open("passport_data.txt") do f
    map(passport -> Dict(f.captures[1] => f.captures[2] for f in eachmatch(field_pattern, passport)), map(r -> r.match, eachmatch(record_pattern, read(f, String))))
end

println(sum(p -> validate_required_fields(p), passports))
println(sum(p -> validate_passport(p), passports))


text:
"""
    MatrixReshaped(D, n, p)
```julia
D::MultivariateDistribution  base distribution
n::Integer   number of rows
p::Integer   number of columns
```
Reshapes a multivariate distribution into a matrix distribution with n rows and
p columns.

"""
struct MatrixReshaped{S<:ValueSupport,D<:MultivariateDistribution{S}} <:
       MatrixDistribution{S}
    d::D
    num_rows::Int
    num_cols::Int
    function MatrixReshaped(
        d::D,
        n::N,
        p::N,
    ) where {
        D<:MultivariateDistribution{
            S,
        },
    } where {S<:ValueSupport} where {N<:Integer}
        (n > 0 && p > 0) || throw(ArgumentError("n and p should be positive"))
        n * p == length(d) ||
        throw(ArgumentError("Dimensions provided ($n x $p) do not match source distribution of length $(length(d))"))
        return new{S,D}(d, n, p)
    end
end

MatrixReshaped(D::MultivariateDistribution, n::Integer) =
    MatrixReshaped(D, n, n)

show(io::IO, d::MatrixReshaped) =
    show_multline(io, d, [(:num_rows, d.num_rows), (:num_cols, d.num_cols)])

#  -----------------------------------------------------------------------------
#  Properties
#  -----------------------------------------------------------------------------

size(d::MatrixReshaped) = (d.num_rows, d.num_cols)

length(d::MatrixReshaped) = length(d.d)

rank(d::MatrixReshaped) = minimum(size(d))

function insupport(d::MatrixReshaped, X::AbstractMatrix)
    return isreal(X) && size(d) == size(X) && insupport(d.d, vec(X))
end

mean(d::MatrixReshaped) = reshape(mean(d.d), size(d))
mode(d::MatrixReshaped) = reshape(mode(d.d), size(d))
cov(d::MatrixReshaped, ::Val{true} = Val(true)) =
    reshape(cov(d.d), prod(size(d)), prod(size(d)))
cov(d::MatrixReshaped, ::Val{false}) =
    ((n, p) = size(d); reshape(cov(d), n, p, n, p))
var(d::MatrixReshaped) = reshape(var(d.d), size(d))

params(d::MatrixReshaped) = (d.d, d.num_rows, d.num_cols)

@inline partype(
    d::MatrixReshaped{S,<:MultivariateDistribution{S}},
) where {S<:Real} = S

_logpdf(d::MatrixReshaped, X::AbstractMatrix) = logpdf(d.d, vec(X))

function _rand!(rng::AbstractRNG, d::MatrixReshaped, Y::AbstractMatrix)
    rand!(rng, d.d, view(Y, :))
    return Y
end

vec(d::MatrixReshaped) = d.d


text:
# module LFP

# function chunk_bymarker(data::Array{D, N}, marker::Vector{T}, roi::Tuple{T, T}, freq::T, mbias::T=0) where {D, T, N}
#     gap = (roi[2]-roi[1]) * freq |> ceil |> Int64
#     result = zeros(D, size(data, 1), gap, size(data, 2))

#     for (midx, eachm) in enumerate(marker)
#         _start = (eachm + roi[1] + mbias) * freq |> floor |> Int64
#         result = [:, :, midx] = #TODO
#     end
#     #TODO
# end

function chunk_bymarker(data::Array{D, 2}, marker::Vector{T}, roi::Tuple{T,T}, freq::T; mbias::T=0.0) where {D, T}
    ngap = (roi[2] - roi[1]) * freq |> ceil |> Int
    _result = zeros(D, ngap, size(data, 2), length(marker))

    for (_idx, _marker) in enumerate(marker)
        _start = (_marker + roi[1] + mbias) * freq |> floor |> Int
        _result[:, :, _idx] = data[_start:_start+ngap-1, :]
    end
    _result
end

function chunk_bymarker(data::Vector{D}, marker::Vector{T}, roi::Tuple{T,T}, freq::T; mbias::T=0.0) where {D, T}
    ngap = (roi[2] - roi[1]) * freq |> ceil |> Int
    _result = zeros(D, ngap, length(marker))

    for (_idx, _marker) in enumerate(marker)
        _start = (_marker + roi[1] + mbias) * freq |> floor |> Int
        _result[:, _idx] = data[_start:_start+ngap-1]
    end
    _result
end

include("LFP/iSplitContainer.jl")
include("LFP/Decomposition.jl")
# for both one-dimensional and two-dimensional size.

# Container.jl
# create compact channel data (i.e. iSplit format)
# load compact channel data
# create_epoch, chunk_iSplit,

# Decomposition.jl
# morlet wavelet transformation
# cmplx_power
# cmplx_phase

# Filter.jl
# gaussian window filter

# end