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

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"),
			Pkg.PackageSpec(name="ImageMagick"),
			Pkg.PackageSpec(name="Images", version="0.23"),
			])
	using Plots
	using PlutoUI
	using LinearAlgebra
	using Images
end

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

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

# **Bonus homework 8**: _Raytracing in 3D_
`18.S191`, fall 2020

This week we will recreate the 3D raytracer from the lectures, building on what we wrote for the 2D case. Compared to the past weeks, you will notice that this homework is a little different! We have included fewer intermediate checks, it is up to you to decide how you:
- split up the problem is smaller functions
- check and visualize your progress

To guide you through the homework, you can follow along with this week's live coding session. Feel free to stay close to James's code, or to come up with your own solution.

With the raytracer done, we will be able to recreate an artwork by M.C. Escher, scroll down to have a sneak peek!

---

_For MIT students:_ this homework is **optional**.

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 = "Jazzy Doe", kerberos_id = "jazz")

# 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)
"""

# ╔═╡ e4d2f6f6-2089-11eb-0b4a-0d22f72f72ba
html"""

Lecture

Live coding session

"""

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

# ╔═╡ 4e917968-1f87-11eb-371f-e3899b76dc24
md"""
## From the last homework

Below we have included some important functions from the last homework (_Raytracing in 2D_), which we will be able to re-use for the 3D case.

There are some small changes:
1. The concept of a `Photon` now carries **color information**.
2. The `Sphere` is no longer a pure lens, it contains a `Surface` property which describes the mixture between transmission, reflection and a pure color. More on this later!
3. The `refract` function is updated to handle two edge cases, but its behaviour is generally unchanged.

Outside of these changes, all functions from the previous homework can be taken "as-is" when converting to 3D, cool!
"""

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

	"Direction vector"
	l::Vector{Float64}

	"Color associated with the photon"
	c::RGB

	ior::Real
end

# ╔═╡ c6e8d30e-205c-11eb-271c-6165a164073d
md"""
#### Intersections:
"""

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

# ╔═╡ 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

# ╔═╡ fcde90ca-2048-11eb-3e96-f9f47b6154e8
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 < b.distance
end

# ╔═╡ dc36ceaa-205c-11eb-169c-bb4c36aaec9f
md"""
#### Reflect and refract:
"""

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

# ╔═╡ 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)

	if abs(c) > 0.999
		ℓ₁
	else
		f = 1 - r^2 * (1 - c^2)
		if f < 0
			ℓ₁
		else
			normalize(r * ℓ₁ + (r*c - sqrt(f)) * n̂_oriented)
		end
	end
end

# ╔═╡ 7f0bf286-2071-11eb-0cac-6d10c93bab6c
md"""
#### Surface (new)
"""

# ╔═╡ 8a4e888c-1ef7-11eb-2a52-17db130458a5
struct Surface
	# Reflectivity
	r::Float64

	# Transmission
	t::Float64

	# Color
	c::RGBA

	# index of refraction
	ior::Float64

end

# ╔═╡ 9c3bdb62-1ef7-11eb-2204-417510bf0d72
html"""

Sphere

Aasdf

"""

# ╔═╡ cb7ed97e-1ef7-11eb-192c-abfd66238378
struct Sphere <: Object
	# Lens position
	p::Vector{Float64}

	# Lens radius
	r::Float64

	s::Surface
end

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

# ╔═╡ 452d6668-1ec7-11eb-3b0a-0b8f45b43fd5
md"""
## Camera and Skyboxes

Now we can begin looking into the 3D nature of raytracing to create visualizations similar to those in lecture.
The first step is setting up the camera and another stuct known as a *sky box* to collect all the rays of light.

Luckily, the transition from 2D to 3D for raytracing is relatively straightforward and we can use all of the functions and concepts we have built in 2D moving forward.

Firstly, the camera:

"""

# ╔═╡ 791f0bd2-1ed1-11eb-0925-13c394b901ce
md"""
### Camera

For the purposes of this homework, we will constrain ourselves to a camera pointing exclusively downward.
This is simply because camera positioning can be a bit tricky and there is no reason to make the homework more complicated than it needs to be!

So, what is the purpose of the camera?

Well, in reality, a camera is a device that collects the color information from all the rays of light that are refracting and reflecting off of various objects in some sort of scene.
Because there are a nearly infinite number of rays bouncing around the scene at any time, we will actually constrain ourselves only to rays that are entering our camera.
In poarticular, we will create a 2D screen just in front of the camera and send a ray from the camera to each pixel in the screen, as shown in the following image:

$(RemoteResource("https://upload.wikimedia.org/wikipedia/commons/thumb/8/83/Ray_trace_diagram.svg/1920px-Ray_trace_diagram.svg.png", :width=>400, :style=>"display: block; margin: auto;"))
"""

# ╔═╡ 1a446de6-1ec9-11eb-1e2f-6f4376005d24
md"""
Because we are not considering camera motion for this exercise, we will assume that the image plane is constrained to the horizontal plane, but that the camera, itself, can be some distance behind it.
This distance from the image plane to the camera is called the *focal length* and is used to determine the field of view.

From here, it's clear we need to construct:
1. A camera struct
2. A function to initialize all the rays being generated by the camera

Let's start with the struct
"""

# ╔═╡ 88576c6e-1ecb-11eb-3e34-830aeb433df1
struct Camera <: Object
	"Set of all pixels, counts as scene resolution"
	resolution::Tuple{Int64,Int64}

	"Physical size of aperture"
	aperture_width::Float64

	"Camera's distance from screen"
	focal_length::Float64

	"Camera's position"
	p::Vector{Float64}
end

# ╔═╡ e774d6a8-2058-11eb-015a-83b4b6104e6e
test_cam = Camera((400,300), 9, -10, [0,00,100])

# ╔═╡ 8f73824e-1ecb-11eb-0b28-4d1bc0eefbc3
md"""
Now we need to construct some method to create each individual ray extending from the camera to a pixel in the image plane.
"""

# ╔═╡ 4006566e-1ecd-11eb-2ce1-9d1107186784
function init_rays(cam::Camera)

	# Physical size of the aperture/image/grid
	aspect_ratio = cam.resolution[1] / cam.resolution[2]
	dim = (
		cam.aperture_width,
		cam.aperture_width / aspect_ratio
	)

	# The x, y coordinates of every pixel in our image grid
	# relative to the image center
	xs = LinRange(-0.5* dim[1], 0.5 * dim[1], cam.resolution[1])
	ys = LinRange(0.5* dim[2], -0.5 * dim[2], cam.resolution[2])

	pixel_positions = [[x, y, cam.focal_length] for y in ys, x in xs]
	directions = normalize.(pixel_positions)

	Photon.([cam.p], directions, [zero(RGB)], [1.0])
end

# ╔═╡ 156c0d7a-2071-11eb-1551-4f2d393df6c8
tiny_resolution_camera = Camera((4,3), 16, -5, [0, 20, 100])

# ╔═╡ 2838c1e4-2071-11eb-13d8-1da955fbf544
test_rays = init_rays(tiny_resolution_camera)

# ╔═╡ 595acf48-1ef6-11eb-0b46-934d17186e7b
extract_colors(rays) = map(ray -> ray.c, rays)

# ╔═╡ b7fa4512-2089-11eb-255d-4d6de9cdfb8e
extract_colors(test_rays)

# ╔═╡ 494687f6-1ecd-11eb-3ada-6f11f45aa74f
md"""
Nothing yet... time to add some objects!
### Skybox

Now that we have the concept of a camera, we can technically do a fully 3D raytracing example; however, we want to ensure that each pixel will actually *hit* something -- preferrably something with some color gradient so we can make sure our simulation is working!

For this, we will introduce the concept of a sky box, which is standard for most gaming applications.
Here, the idea is that our entire scene is held within some additional object, just like the mirrors we used in the 2D example.
The only difference here is that we will be using some texture instead of a reflective surface.
In addition, even though we are calling it a box, we'll actually be treating it as a sphere.

Because we have already worked out how to make sure we have hit the interior of a spherical lens, we will be using a similar function here.
For this part of the exercise, we will need to construct 2 things:

1. A skybox struct
2. A function that returns some color gradient to be called whenever a ray of light interacts with a sky box

So let's start with the sky box struct
"""

# ╔═╡ 9e71183c-1ef4-11eb-1802-3fc60b51ceba
struct SkyBox <: Object
	# Skybox position
	p::Vector{Float64}

	# Skybox radius
	r::Float64

	# Color function
	c::Function
end

# ╔═╡ 093b9e4a-1f8a-11eb-1d32-ad1d85ddaf42
function intersection(photon::Photon, sphere::S; ϵ=1e-3) where {S <: Union{SkyBox, Sphere}}
	a = dot(photon.l, photon.l)
	b = 2 * dot(photon.l, photon.p - sphere.p)
	c = dot(photon.p - sphere.p, photon.p - sphere.p) - sphere.r^2

	d = b^2 - 4*a*c

	if d <= 0
		Miss()
	else
		t1 = (-b-sqrt(d))/2a
		t2 = (-b+sqrt(d))/2a

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

		point = photon.p + t * photon.l

		Intersection(sphere, t, point)
	end
end

# ╔═╡ 89e98868-1fb2-11eb-078d-c9298d8a9970
function closest_hit(photon::Photon, objects::Vector{<:Object})
	hits = intersection.([photon], objects)

	minimum(hits)
end

# ╔═╡ aa9e61aa-1ef4-11eb-0b56-cd7ded52b640
md"""
Now we have the ability to create a skybox, the only thing left is to create some sort of texture function so that when the ray of light hits the sky box, we can return some form of color information.
So for this, we will basically create a function that returns back a smooth gradient in different directions depending on the position of the ray when it hits the skybox.

For the color information, we will be assigning a color to each cardinal axis.
That is to say that there will be a red gradient along $x$, a blue gradient along $y$, and a green gradient along $z$.
For this, we will need to define some extent over which the gradient will be active in 'real' units.
From there, we can say that the gradient is

$$\frac{r+D}{2D},$$

where $r$ is the ray's position when it hits the skybox, and $D$ is the extent over which the gradient is active.

So let's get to it and write the function!
"""

# ╔═╡ c947f546-1ef5-11eb-0f02-054f4e7ae871
function gradient_skybox_color(position, skybox)
	extents = skybox.r
	c = zero(RGB)

	if position[1] < extents && position[1] > -extents
		c += RGB((position[1]+extents)/(2.0*extents), 0, 0)
	end

	if position[2] < extents && position[2] > -extents
		c += RGB(0,0,(position[2]+extents)/(2.0*extents))
	end

	if position[3] < extents && position[3] > -extents
		c += RGB(0,(position[3]+extents)/(2.0*extents), 0)
	end

	return c
end

# ╔═╡ a919c880-206e-11eb-2796-55ccd9dbe619
sky = SkyBox([0.0, 0.0, 0.0], 1000, gradient_skybox_color)

# ╔═╡ 49651bc6-2071-11eb-1aa0-ff829f7b4350
md"""
Let's set up a basic scene and trace an image! Since our skybox is _spherical_ we can use **the same `intersect`** method as we use for `Sphere`s. Have a look at [the `intersect` method](#sphere-defs), we already added `SkyBox` as a possible type.
"""

# ╔═╡ daf80644-2070-11eb-3363-c577ae5846b3
basic_camera = Camera((300,200), 16, -5, [0,20,100])

# ╔═╡ df3f2178-1ef5-11eb-3098-b1c8c67cf136
md"""
To create this image, we used the ray tracing function bewlow, which takes in a camera and a set of objects / scene, and...
1. Initilializes all the rays
2. Propagates the rays forward
3. Converts everything into an image
"""

# ╔═╡ 04a86366-208b-11eb-1977-ff7e4ae6b714
md"""
## Writing a ray tracer

It's your turn! Below is the code needed to trace just the sky box, but we still need to add the ability to trace spheres.

**We recommend** that you start by just implementing _reflection_ - make every sphere reflect, regardless of its surface. Make sure that this is working well - can you see the reflection of one sphere in another sphere? Does our program get stuck in a loop?

Once you have reflections working, you can add _refraction_ and _colored spheres_. In the 2D example, we dealt specifically with spheres that could either 100% reflect or refract. In reality, it is possible for objects to either reflect or refract, something in-between.
That is to say, a ray of light can *split* when hitting an object surface, creating at least 2 more rays of light that will both return separate color values.
The color that we _perceive_ for that ray is the combination both of these colors - they are mixed.

A third possibility explored in the lecture is that the objects can also have a color associated with them and just return the color value instead of reflecting or refracting.

**You can choose!** After implementing reflection, you can implement three different spheres (you can modify the existing code, create new types, add functions, and so on), a purely reflective, purely refractive or opaquely colored sphere. You can also go straight for the more photorealistic option, which is that every sphere is a combination of these three - this is what we did in the lecture.
"""

# ╔═╡ a9754410-204d-11eb-123e-e5c5f87ae1c5
function interact(ray::Photon, hit::Intersection{SkyBox}, ::Any, ::Any)

	ray_color = hit.object.c(hit.point, hit.object)
	Photon(hit.point, ray.l, ray_color, ray.ior)
end

# ╔═╡ 086e1956-204e-11eb-2524-f719504fb95b
interact(photon::Photon, ::Miss, ::Any, ::Any) = photon

# ╔═╡ 95ca879a-204d-11eb-3473-959811aa8320
function step_ray(ray::Photon, objects::Vector{O},
			   num_intersections) where {O <: Object}

	if num_intersections == 0
		ray
	else
		hit = closest_hit(ray, objects)
		interact(ray, hit, num_intersections, objects)
	end
end

# ╔═╡ 6b91a58a-1ef6-11eb-1c36-2f44713905e1
function ray_trace(objects::Vector{O}, cam::Camera;
				   num_intersections = 10) where {O <: Object}
	rays = init_rays(cam)

	new_rays = step_ray.(rays, [objects], [num_intersections])

	extract_colors(new_rays)
end

# ╔═╡ a0b84f62-2047-11eb-348c-db83f4e6c39c
let
	scene = [sky]
	ray_trace(scene, basic_camera; num_intersections=4)
end

# ╔═╡ d1970a34-1ef7-11eb-3e1f-0fd3b8e9657f
md"""
Below, we create a scene with a number of balls inside of it.
While working on your code, work in small increments, and do frequent checks to see if your code is working. Feel free to modify this test scene, or to create a simpler one.
"""

# ╔═╡ 16f4c8e6-2051-11eb-2f23-f7300abea642
main_scene = [
	sky,
	Sphere([0,0,-25], 20,
		Surface(1.0, 0.0, RGBA(1,1,1,0.0), 1.5)),

	Sphere([0,50,-100], 20,
		Surface(0.0, 1.0, RGBA(0,0,0,0.0), 1.0)),

	Sphere([-50,0,-25], 20,
		Surface(0.0, 0.0, RGBA(0, .3, .8, 1), 1.0)),

	Sphere([30, 25, -60], 20,
		Surface(0.0, 0.75, RGBA(1,0,0,0.25), 1.5)),

	Sphere([50, 0, -25], 20,
		Surface(0.5, 0.0, RGBA(.1,.9,.1,0.5), 1.5)),

	Sphere([-30, 25, -60], 20,
		Surface(0.5, 0.5, RGBA(1,1,1,0), 1.5)),
]

# ╔═╡ 1f66ba6e-1ef8-11eb-10ba-4594f7c5ff19
let
	cam = Camera((600,360), 16, -15, [0,10,100])

	ray_trace(main_scene, cam; num_intersections=3)
end

# ╔═╡ 67c0bd70-206a-11eb-3935-83d32c67f2eb
md"""
## **Bonus:** Escher

If you managed to get through the exercises, we have a fun bonus exercise! The goal is to recreate this self-portrait by M.C. Escher:

"""

# ╔═╡ 748cbaa2-206c-11eb-2cc9-7fa74308711b
Resource("https://www.researchgate.net/profile/Madhu_Gupta22/publication/3427377/figure/fig1/AS:663019482775553@1535087571714/A-self-portrait-of-MC-Escher-1898-1972-in-spherical-mirror-dating-from-1935-titled.png", :width=>300)

# ╔═╡ 981e6bd2-206c-11eb-116d-6fad4e04ce34
md"""
It looks like M.C. Escher is a skillful raytracer, but so are we! To recreate this image, we can simplify it by having just two objects in our scene:
- A purely reflective sphere
- A skybox, containing an image of us!

Let's start with our old skybox, and set up our scene:
"""

# ╔═╡ 7a12a99a-206d-11eb-2393-bf28b881087a
escher_sphere = Sphere([0,0,0], 20,
			Surface(1.0, 0.0, RGBA(1,1,1,0.0), 1.5))

# ╔═╡ 373b6a26-206d-11eb-1e67-9debb032f69e
escher_cam = Camera((300,300), 30, -10, [0,00,30])

# ╔═╡ 5dfec31c-206d-11eb-23a2-259f2c205cb5
md"""
👆 You can modify `escher_cam` to increase or descrease the resolution!
"""

# ╔═╡ 6f1dbf48-206d-11eb-24d3-5154703e1753
let
	scene = [sky, escher_sphere]
	ray_trace(scene, escher_cam; num_intersections=3)
end

# ╔═╡ dc786ccc-206e-11eb-29e2-99882e6613af
md"""
Awesome! Next, we want to set an _image_ as our skybox, instead of a gradient. To do so, we have written a function that converts the x,y,z coordinates of the intersection point with a skybox into a latitude/longitude pair, which we can use (after scaling, rounding & clamping) as pixel coordinates to index an image!
"""

# ╔═╡ 8ebe4cd6-2061-11eb-396b-45745bd7ec55
earth = load(download("https://upload.wikimedia.org/wikipedia/commons/thumb/8/8f/Whole_world_-_land_and_oceans_12000.jpg/1280px-Whole_world_-_land_and_oceans_12000.jpg"))

# ╔═╡ 12d3b806-2062-11eb-20a8-7d1a33e4b073
function get_index_rational(A, x, y)
	a, b = size(A)

	u = clamp(floor(Int, x * (a-1)) + 1, 1, a)
	v = clamp(floor(Int, y * (b-1)) + 1, 1, b)
	A[u,v]
end

# ╔═╡ cc492966-2061-11eb-1000-d90c279c4668
function image_skybox(img)
	f = function(position, skybox)
		lon = atan(-position[1], position[3])
		lat = -atan(position[2], norm(position[[1,3]]))

		get_index_rational(img, (lat/(pi)) + .5, (lon/2pi) + .5)
	end

	SkyBox([0.0, 0.0, 0.0], 1000, f)
end

# ╔═╡ 137834d4-206d-11eb-0082-7b87bf222808
earth_skybox = image_skybox(earth)

# ╔═╡ bff27890-206e-11eb-2e40-696424a0b8be
let
	scene = [earth_skybox, escher_sphere]
	ray_trace(scene, escher_cam; num_intersections=3)
end

# ╔═╡ b0bc76f8-206d-11eb-0cad-4bde96565fed
md"""
Great! It's like the Earth, but distorted. Notice that the continents are mirrored, and that you can see both poles at the same time.

Okay, self portrait time! Let's take a picture using your webcam, and we will use it as the skybox texture:
"""

# ╔═╡ 48166866-2070-11eb-2722-556a6719c2a2
md"""
👀 wow! It's _Planet $(student.name)_, surrounded by even more $(student.name).

When you look at the drawing by Escher, you see that he only occupies a small section of the 'skybox'. Behind Escher, you can see his cozy house, and in _front_ of him (i.e. _behind_ the glass sphere, from his persective), you see a gray background.

What we need is a 360° panoramic image of our room. One option is that you make one! There are great tutorials online, and maybe you can use an app to do this with your phone.

Another option is that we approximate the panaroma by _padding_ the image of your face. You can pad the image with a solid color, you can use a gradient, you can use the earth satellite image, you can add noise, and so on. [Here](https://user-images.githubusercontent.com/6933510/98426463-31a64f00-2099-11eb-96dd-9f0f065e0da4.png) is an example of what I made.
"""

# ╔═╡ 6480b85c-2067-11eb-0262-f752d306d8ae
function padded(img)

	return missing
end

# ╔═╡ 3de614da-2091-11eb-361a-83bcf357c394
md"""
Let's put it all together!
"""

# ╔═╡ 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."

# ╔═╡ 2c36d3e8-2065-11eb-3a12-8382f8539ea6

function process_raw_camera_data(raw_camera_data)
	# the raw image data is a long byte array, we need to transform it into something
	# more "Julian" - something with more _structure_.

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

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

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

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

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

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

	RGB.(reds, greens, blues)
end

# ╔═╡ f7b5ff68-2064-11eb-3be3-554519ca4847
function camera_input(;max_size=200, default_url="https://i.imgur.com/SUmi94P.png")
"""

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

		<video playsinline autoplay></video>

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

		Enable webcam

""" |> HTML
end

# ╔═╡ 27d64432-2065-11eb-3795-e99b1d6718d2
@bind wow camera_input()

# ╔═╡ 64ce8106-2065-11eb-226c-0bcaf7e3f871
face = process_raw_camera_data(wow)

# ╔═╡ 06ac2efc-206f-11eb-1a73-9306bf5f7a9c
let
	face_skybox = image_skybox(face)
	scene = [face_skybox, escher_sphere]
	ray_trace(scene, escher_cam; num_intersections=3)
end

# ╔═╡ 7d03b258-2067-11eb-3070-1168e282b2ea
padded(face)

# ╔═╡ aa597a16-2066-11eb-35ae-3170468a90ed
@bind escher_face_data camera_input()

# ╔═╡ c68dbe1c-2066-11eb-048d-038df2c68a8b
let
	img = process_raw_camera_data(escher_face_data)
	img_padded = padded(img)

	scene = [
		image_skybox(padded(img)),

		escher_sphere,
	]

	ray_trace(scene, escher_cam; num_intersections=20)
end

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

# ╔═╡ 7c804c30-208d-11eb-307c-076f2086ae73
hint(md"""If you are getting a _"Circular Defintions"_ error - this could be because of a Pluto limitation. If two functions call each other, they need to be contained in a single cell, using a `begin end` block.""")

# ╔═╡ 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]))

# ╔═╡ 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
# ╟─1e109620-19c4-11eb-013e-1bc95c14c2ba
# ╠═1e202680-19c4-11eb-29a7-99061b886b3c
# ╟─e4d2f6f6-2089-11eb-0b4a-0d22f72f72ba
# ╟─1e2cd0b0-19c4-11eb-3583-0b82092139aa
# ╠═c3e52bf2-ca9a-11ea-13aa-03a4335f2906
# ╟─4e917968-1f87-11eb-371f-e3899b76dc24
# ╠═24b0d4ba-192c-11eb-0f66-e77b544b0510
# ╟─c6e8d30e-205c-11eb-271c-6165a164073d
# ╠═d851a202-1ca0-11eb-3da0-51fcb656783c
# ╠═8acef4b0-1a09-11eb-068d-79a259244ed1
# ╠═8018fbf0-1a05-11eb-3032-95aae07ca78f
# ╠═fcde90ca-2048-11eb-3e96-f9f47b6154e8
# ╠═89e98868-1fb2-11eb-078d-c9298d8a9970
# ╟─dc36ceaa-205c-11eb-169c-bb4c36aaec9f
# ╠═43306bd4-194d-11eb-2e30-07eabb8b29ef
# ╠═14dc73d2-1a0d-11eb-1a3c-0f793e74da9b
# ╟─7f0bf286-2071-11eb-0cac-6d10c93bab6c
# ╠═8a4e888c-1ef7-11eb-2a52-17db130458a5
# ╠═9c3bdb62-1ef7-11eb-2204-417510bf0d72
# ╠═cb7ed97e-1ef7-11eb-192c-abfd66238378
# ╠═093b9e4a-1f8a-11eb-1d32-ad1d85ddaf42
# ╠═6fdf613c-193f-11eb-0029-957541d2ed4d
# ╟─452d6668-1ec7-11eb-3b0a-0b8f45b43fd5
# ╟─791f0bd2-1ed1-11eb-0925-13c394b901ce
# ╟─1a446de6-1ec9-11eb-1e2f-6f4376005d24
# ╠═88576c6e-1ecb-11eb-3e34-830aeb433df1
# ╠═e774d6a8-2058-11eb-015a-83b4b6104e6e
# ╟─8f73824e-1ecb-11eb-0b28-4d1bc0eefbc3
# ╠═4006566e-1ecd-11eb-2ce1-9d1107186784
# ╠═156c0d7a-2071-11eb-1551-4f2d393df6c8
# ╠═2838c1e4-2071-11eb-13d8-1da955fbf544
# ╠═595acf48-1ef6-11eb-0b46-934d17186e7b
# ╠═b7fa4512-2089-11eb-255d-4d6de9cdfb8e
# ╟─494687f6-1ecd-11eb-3ada-6f11f45aa74f
# ╠═9e71183c-1ef4-11eb-1802-3fc60b51ceba
# ╟─aa9e61aa-1ef4-11eb-0b56-cd7ded52b640
# ╠═c947f546-1ef5-11eb-0f02-054f4e7ae871
# ╠═a919c880-206e-11eb-2796-55ccd9dbe619
# ╟─49651bc6-2071-11eb-1aa0-ff829f7b4350
# ╠═daf80644-2070-11eb-3363-c577ae5846b3
# ╠═a0b84f62-2047-11eb-348c-db83f4e6c39c
# ╟─df3f2178-1ef5-11eb-3098-b1c8c67cf136
# ╠═6b91a58a-1ef6-11eb-1c36-2f44713905e1
# ╟─04a86366-208b-11eb-1977-ff7e4ae6b714
# ╠═a9754410-204d-11eb-123e-e5c5f87ae1c5
# ╠═086e1956-204e-11eb-2524-f719504fb95b
# ╠═95ca879a-204d-11eb-3473-959811aa8320
# ╠═1f66ba6e-1ef8-11eb-10ba-4594f7c5ff19
# ╟─d1970a34-1ef7-11eb-3e1f-0fd3b8e9657f
# ╠═16f4c8e6-2051-11eb-2f23-f7300abea642
# ╟─7c804c30-208d-11eb-307c-076f2086ae73
# ╟─67c0bd70-206a-11eb-3935-83d32c67f2eb
# ╟─748cbaa2-206c-11eb-2cc9-7fa74308711b
# ╟─981e6bd2-206c-11eb-116d-6fad4e04ce34
# ╠═7a12a99a-206d-11eb-2393-bf28b881087a
# ╠═373b6a26-206d-11eb-1e67-9debb032f69e
# ╟─5dfec31c-206d-11eb-23a2-259f2c205cb5
# ╠═6f1dbf48-206d-11eb-24d3-5154703e1753
# ╟─dc786ccc-206e-11eb-29e2-99882e6613af
# ╟─8ebe4cd6-2061-11eb-396b-45745bd7ec55
# ╠═cc492966-2061-11eb-1000-d90c279c4668
# ╠═12d3b806-2062-11eb-20a8-7d1a33e4b073
# ╠═137834d4-206d-11eb-0082-7b87bf222808
# ╠═bff27890-206e-11eb-2e40-696424a0b8be
# ╟─b0bc76f8-206d-11eb-0cad-4bde96565fed
# ╠═27d64432-2065-11eb-3795-e99b1d6718d2
# ╟─64ce8106-2065-11eb-226c-0bcaf7e3f871
# ╠═06ac2efc-206f-11eb-1a73-9306bf5f7a9c
# ╟─48166866-2070-11eb-2722-556a6719c2a2
# ╠═6480b85c-2067-11eb-0262-f752d306d8ae
# ╠═7d03b258-2067-11eb-3070-1168e282b2ea
# ╟─3de614da-2091-11eb-361a-83bcf357c394
# ╟─aa597a16-2066-11eb-35ae-3170468a90ed
# ╠═c68dbe1c-2066-11eb-048d-038df2c68a8b
# ╟─ebd05bf0-19c3-11eb-2559-7d0745a84025
# ╟─ec275590-19c3-11eb-23d0-cb3d9f62ba92
# ╟─2c36d3e8-2065-11eb-3a12-8382f8539ea6
# ╟─f7b5ff68-2064-11eb-3be3-554519ca4847
# ╟─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

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
# SmallStrainContinuum

prop = Dict("name"=> "PlaneStrain", "rho"=> 1.0, "E"=> 2.0, "nu"=> 0.35)
ngp = 2
elnodes = [1;2;3;4]
nodes = [
    0.0 0.0
    1.0 0.0
    1.0 1.0
    0.0 1.0
]
element =  SmallStrainContinuum(nodes, elnodes, prop, ngp)

ν = 0.35
E = 2.0

H = zeros(3,3)
H[1,1] = E*(1. -ν)/((1+ν)*(1. -2. *ν));
H[1,2] = H[1,1]*ν/(1-ν);
H[2,1] = H[1,2];
H[2,2] = H[1,1];
H[3,3] = H[1,1]*0.5*(1. -2. *ν)/(1. -ν);

@testset "small strain gauss" begin
    ggp = getGaussPoints(e)
    c = e.coords
    scatter(c[:,1], c[:,2])
    scatter(ggp[:,1], ggp[:,2])
end

@testset "body force" begin
    ggp = getGaussPoints(element)
    f = @. ggp[:,1]^2 + ggp[:,2]^2
    fbody = getBodyForce(element, [f f])
    idx = [1;2;4;3]
    idx = [idx; idx .+ 4]
    @test fbody[idx] ≈ compute_fem_source_term(f, f, 1, 1, 1.)
end

m, n =  2, 2
h = 1/m

# Create a very simple mesh
elements = SmallStrainContinuum[]
prop = Dict("name"=> "PlaneStrain", "rho"=> 1.0, "E"=> 2.0, "nu"=> 0.35)
coords = zeros((m+1)*(n+1), 2)
for j = 1:n
    for i = 1:m
        idx = (m+1)*(j-1)+i
        elnodes = [idx; idx+1; idx+1+m+1; idx+m+1]
        ngp = 3
        nodes = [
            (i-1)*h (j-1)*h
            i*h (j-1)*h
            i*h j*h
            (i-1)*h j*h
        ]
        coords[elnodes, :] = nodes
        push!(elements, SmallStrainContinuum(nodes, elnodes, prop, ngp))
    end
end

# fixed on the bottom, push on the right
EBC = zeros(Int64, (m+1)*(n+1), 2)
FBC = zeros(Int64, (m+1)*(n+1), 2)
g = zeros((m+1)*(n+1), 2)
f = zeros((m+1)*(n+1), 2)

for i = 1:m+1
    for j = 1:n+1
        if i==1 || i==m+1 || j == 1|| j==n+1
            idx = (j-1)*(m+1) + i
            FBC[idx,:] .= -1
        end
    end
end

ndims = 2
domain = Domain(coords, elements, ndims, EBC, g, FBC, f)

function bt(x, y, t)
    f1 = @. x^2+y^2
    f2 = @. x^2+y^2
    [f1 f2]
end

Dstate = zeros(domain.neqs)
state = rand(domain.neqs)
velo =  zeros(domain.neqs)
acce =  zeros(domain.neqs)
globdat = GlobalData(state, Dstate, velo, acce, domain.neqs, gt, ft, bt)

assembleMassMatrix!(globdat, domain)
updateTimeDependentEssentialBoundaryCondition!(domain,globdat)
updateStates!(domain, globdat)

@testset "small strain body force full" begin
    fbody = getBodyForce(domain, globdat)
    f = eval_f_on_gauss_pts((x,y)->x^2+y^2, m, n, h)
    @test norm(fbody - compute_fem_source_term(f, f, m, n, h)) < 1e-10
end

@testset "small strain assembleMassMatrix!"
    @test maximum(abs.(compute_fem_mass_matrix(m, n, h) - globdat.M))<1e-10
end

@testset "small strain assembleStiffAndForce" begin
    fint, stiff = assembleStiffAndForce( globdat, domain)
    @test maximum(abs.(compute_fem_stiffness_matrix(H, m, n, h)-stiff))<1e-10
    @test norm(fint - stiff * state)<1e-10
end

ext: jl


text:
using CSV
using DataFrames
using DelimitedFiles
using YAML
using Query

using SLiDE

UNITS = "billions of us dollars (USD)"

# ******************************************************************************************
#   READ BLUENOTE DATA -- For benchmarking!
# ******************************************************************************************
BLUE2 = joinpath("data", "windc_output", "2z_windcdatabase")
b2_lst = [x for x in readdir(joinpath(SLIDE_DIR, BLUE2)) if occursin(".csv", x)]
# b2 = Dict(Symbol(k[1:end-4]) => sort(edit_with(
#     read_file(joinpath(BLUE2, k))), Rename(:Val, :value))) for k in b2_lst)

BLUE_DIR = joinpath("data", "windc_output", "2a_io_national_cgeparm_raw")
bluenote_lst = [x for x in readdir(joinpath(SLIDE_DIR, BLUE_DIR)) if occursin(".csv", x)]
bluenote = Dict(Symbol(k[1:end-4]) => sort(edit_with(
    read_file(joinpath(BLUE_DIR, k)), Rename(:Val, :value))) for k in bluenote_lst)

# Add supply/use info for checking.
BLUE_DIR_IN = joinpath("data", "windc_output", "1b_stream_windc_base")
[bluenote[k] = sort(edit_with(read_file(joinpath(BLUE_DIR_IN, string(k, "_units.csv"))), [
    Rename.([:Dim1,:Dim2,:Dim3,:Dim4,:Val], [:yr,:i,:j,:units,:value]);
    Replace.([:i,:j], "upper", "lower")])) for k in [:supply, :use]]

[bluenote[k][!,:value] .= round.(bluenote[k][:,:value]*1E-3, digits=3) # convert millions -> billions USD
    for k in [:supply,:use]]

# ******************************************************************************************
#   READ SETS AND SLiDE SUPPLY/USE DATA.
# ******************************************************************************************
y = read_file(joinpath("data","readfiles","list_sets.yml"));
set = Dict(Symbol(k) => sort(read_file(joinpath(y["Path"]..., ensurearray(v)...)))[:,1]
    for (k,v) in y["Input"])

# Read supply/use data.
DATA_DIR = joinpath("data", "input")
io_lst = convert_type.(Symbol, ["supply", "use"])
jo = Dict(k => read_file(joinpath(DATA_DIR, string(k, ".csv"))) for k in io_lst)

# Perform some minor edits that will apply to all parameters in order to maintain
# consistency with the data in build_national_cgeparm_raw.gdx for easier benchmarking.

for k in keys(jo)
    global jo[k] = edit_with(jo[k], [Rename(:value, k); Drop.([:value,:units], [0.,"all"], "==")])
    global jo[k] = jo[k] |> @filter(_.yr in set[:yr]) |> DataFrame

    # global jo[k][!,:value] .= round.(jo[k][:,:value]*1E-3, digits=3)
    # global jo[k][!,:units] .= UNITS
end

jo[:supply], jo[:use] = fill_zero(jo[:supply], jo[:use]; permute_keys = true);

df = jo[:supply][:,find_oftype(jo[:supply], Not(AbstractFloat))]
df = hcat(df, [edit_with(jo[k],Rename(:value,k))[:,[k]] for k in sort(collect(keys(jo)))]...)
df = fill_zero((yr = set[:yr], i = set[:i], j = set[:j]), df)

# df = copy(df)
# x = (i = set[:i], j = set[:j])
# col = [:a, :b]

function add_filter(df::DataFrame, col::Array{Symbol,1}, set)
    df = copy(df)
    cols = intersect(find_oftype(df, Not(AbstractFloat)), find_oftype(df, Not(Bool)))
    inset = [col in collect(keys(set)) for col in cols]

    vals_sets = [[set[k] for k in cols[inset]]; unique.(eachcol(df[:,cols[.!inset]]))]
    cols_sets = [cols[inset]; cols[.!inset]]

    df_sets = edit_with(DataFrame(permute(NamedTuple{Tuple(cols_sets,)}(vals_sets,))),
        Add.(col, true))
    df = leftjoin(df, df_sets, on = cols_sets)

    df = edit_with(df, Replace.(col, missing, false))
    return df
end

s = Dict()

mutable struct Filter
    set::NamedTuple
    col::Array{Symbol,1}
end

filters = []

push!(filters, Filter((yr = set[:yr], i = set[:i], j = set[:j]),     [:ys0, :id0]))
push!(filters, Filter((yr = set[:yr], i = set[:i], j = set[:fd]),    [:fd0]))
push!(filters, Filter((yr = set[:yr], i = set[:va], j = set[:j]),    [:va0]))
push!(filters, Filter((yr = set[:yr], i = set[:ts], j = set[:j]),    [:ts0]))

push!(filters, Filter((yr = set[:yr], i = set[:i], j = "exports"),   [:x0]))
push!(filters, Filter((yr = set[:yr], i = set[:i], j = "imports"),   [:m0]))
push!(filters, Filter((yr = set[:yr], i = set[:i], j = "margins"),   [:mrg0]))
push!(filters, Filter((yr = set[:yr], i = set[:i], j = "trncost"),   [:trn0]))
push!(filters, Filter((yr = set[:yr], i = set[:i], j = "ciffob"),    [:cif0]))
push!(filters, Filter((yr = set[:yr], i = set[:i], j = "duties"),    [:duty0]))
push!(filters, Filter((yr = set[:yr], i = set[:i], j = "tax"),       [:tax0]))
push!(filters, Filter((yr = set[:yr], i = set[:i], j = "subsidies"), [:sbd0]))

# Create DataFrame from filters.
df_filtered = copy(df[:,[:yr,:i,:j]])
[global df_filtered = add_filter(df_filtered, f.col, f.set) for f in filters]

# Create final DataFrames
ij = fill_zero((yr = set[:yr], i = set[:i], j = set[:j]))[:,1:end-1]
i  = fill_zero((yr = set[:yr], i = set[:i]))[:,1:end-1]
fd = fill_zero((yr = set[:yr], i = set[:i], fd = set[:fd]))[:,1:end-1]
ts = fill_zero((yr = set[:yr], ts = set[:ts], j = set[:j]))[:,1:end-1]
va = fill_zero((yr = set[:yr], va = set[:va], j = set[:j]))[:,1:end-1]

kij = [:ys0,:id0]
vij = [:supply,:use]

ki = [:x0,:m0,:mrg0,:trn0,:cif0,:duty0,:tax0,:sbd0]
vi = [:use; fill(:supply,7)]

[global i[!,k] .= df[df_filtered[:,k],v] for (k,v) in zip(ki,vi)]
[global ij[!,k] .= df[df_filtered[:,k],v] for (k,v) in zip(kij,vij)]
fd[!,:fd0] = df[df_filtered[:,:fd0],:use]
ts[!,:ts0] = df[df_filtered[:,:ts0],:use]
va[!,:va0] = df[df_filtered[:,:va0],:use]

# Treat negative inputs as outputs.
ij[!,:ys0] = ij[:,:ys0] - min.(0, ij[:,:id0])
ij[!,:id0] = max.(0, ij[:,:id0])

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using Test
using Trapz
using Unitful

@info "Started Package Testing"
vx=range(0,1,length=5)
vy=range(0,2,length=10)
vz=range(0,3,length=15)
M=[x^2+y^2+z^2 for x=vx,y=vy,z=vz]
res=28.157801083396322

@testset "Methods full integral,partial integral, 2axis, 1axis" begin
    @test trapz((vx,vy,vz), M) ≈ res

    I_xy=trapz((vx,vy,:), M)
    I_xz=trapz((vx,:,vz), M)
    I_yz=trapz((:,vy,vz), M)

    I_xy_z=trapz((vz,), I_xy)
    I_xz_y=trapz((vy,), I_xz)
    I_yz_x=trapz((vz,), I_xy)
    @test I_xy_z ≈ res
    @test I_yz_x ≈ res
    @test I_xz_y ≈ res
end

@testset "Methods only 1 axis" begin
    I_y=trapz(vy, M, Val(2))
    I_x_y=trapz(vx, I_y, Val(1))
    I_x_y_z=trapz(vz, I_x_y, Val(1))
    @test I_x_y_z ≈ res

    I_y=trapz((:,vy,:), M)
    I_x_y=trapz((vx,:), I_y)
    I_x_y_z=trapz(vz,I_x_y)
    @test I_x_y_z ≈ res
end

@testset "Method only last axis" begin
    I_z=trapz(vz, M)
    I_y=trapz(vy, I_z)
    I_x=trapz(vx, I_y)
    @test I_x ≈ res
end

@testset "Corner case of empty vector" begin
    @test trapz(zeros(Float16,()),zeros(Float16,())) == 0.0
    @test trapz(Float64[],Float64[]) == 0.0
    @test trapz(Float64[1],Float64[1]) == 0.0
end

@testset "Integrate by reshaped vector" begin
    @test trapz((1.0:10.0)',1.0:10.0)==trapz(1.0:10.0,1.0:10.0)
end

@testset "Some Inference Tests" begin
    args=((:,vy,:), M)
    @test typeof(trapz(args...))==Base.return_types(trapz,typeof.(args))[1]
    args=(vy, M, Val(2))
    @test typeof(trapz(args...))==Base.return_types(trapz,typeof.(args))[1]
    args=((vx,vy,:), M)
    @test typeof(trapz(args...))==Base.return_types(trapz,typeof.(args))[1]
end

@testset "@trapz macro" begin
    r=@trapz vx x begin
        @trapz vy y begin
            @trapz vz z begin
                x*x+y*y+z*z
            end
        end
    end
    @test r≈res
end

# test whether units work, which requires proper type information to be
# applied to the output arrays
@testset "units" begin
    x=(1:10)u"m"
    y=range(-0.1, 10, length=length(x))u"J/m"
    res = trapz(x, y)
    @test unit(res) == u"J"
    @test res ≈ 44.55u"J"
end


text:
tests = ["testini.jl",
         "testhttp.jl",
         "testsimple.jl"
]

for test in tests
    include(test)
end


text:
ToVec3(x::Vec{3}) = x
ToVec3(x::Vec{2}) = Vec(x[1], x[2], 0)

@testset "Quaternion" begin
    for T in (Float32, Float64)
        # basic constructors
        @test (@inferred Quaternion{T}((1,2,3,4)))::Quaternion{T} |> Tuple === map(T, (1,2,3,4))
        @test (@inferred Quaternion{T}(1,2,3,4))::Quaternion{T} |> Tuple === map(T, (1,2,3,4))
        @test (@inferred Quaternion{T}(Vec(1,2,3,4)))::Quaternion{T} |> Tuple === map(T, (1,2,3,4))
        @test (@inferred Quaternion{T}(Vec(1,2,3)))::Quaternion{T} |> Tuple === map(T, (0,1,2,3))
        @test (@inferred Quaternion{T}(Vec(1,2)))::Quaternion{T} |> Tuple === map(T, (0,1,2,0))
        @test (@inferred Quaternion{T}(Vec(1)))::Quaternion{T} |> Tuple === map(T, (0,1,0,0))
        @test (@inferred Quaternion{T}(4, Vec(1,2,3)))::Quaternion{T} |> Tuple === map(T, (4,1,2,3))
        @test (@inferred Quaternion{T}(4, Vec(1,2)))::Quaternion{T} |> Tuple === map(T, (4,1,2,0))
        @test (@inferred Quaternion{T}(4, Vec(1)))::Quaternion{T} |> Tuple === map(T, (4,1,0,0))
        @test (@inferred Quaternion{T}(4))::Quaternion{T} |> Tuple === map(T, (4,0,0,0))
        @test (@inferred Quaternion((T(1),2,3,4)))::Quaternion{T} |> Tuple === map(T, (1,2,3,4))
        @test (@inferred Quaternion(T(1),2,3,4))::Quaternion{T} |> Tuple === map(T, (1,2,3,4))
        @test (@inferred Quaternion(Vec(T(1),2,3,4)))::Quaternion{T} |> Tuple === map(T, (1,2,3,4))
        @test (@inferred Quaternion(Vec(T(1),2,3)))::Quaternion{T} |> Tuple === map(T, (0,1,2,3))
        @test (@inferred Quaternion(Vec(T(1),2)))::Quaternion{T} |> Tuple === map(T, (0,1,2,0))
        @test (@inferred Quaternion(Vec(T(1))))::Quaternion{T} |> Tuple === map(T, (0,1,0,0))
        @test (@inferred Quaternion(T(4), Vec(1,2,3)))::Quaternion{T} |> Tuple === map(T, (4,1,2,3))
        @test (@inferred Quaternion(T(4), Vec(1,2)))::Quaternion{T} |> Tuple === map(T, (4,1,2,0))
        @test (@inferred Quaternion(T(4), Vec(1)))::Quaternion{T} |> Tuple === map(T, (4,1,0,0))
        @test (@inferred Quaternion(T(4)))::Quaternion{T} |> Tuple === map(T, (4,0,0,0))

        # properties
        q = Quaternion{T}(1,2,3,4)
        @test propertynames(q) == (:scalar, :vector, :data)
        get_scalar = q -> q.scalar
        get_vector = q -> q.vector
        get_data = q -> q.data
        @test (@inferred get_scalar(q))::T == T(1)
        @test (@inferred get_vector(q))::Vec{3, T} == Vec{3, T}(2,3,4)
        @test (@inferred get_data(q))::NTuple{4, T} == map(T, (1,2,3,4))

        # quaternion
        @test (@inferred quaternion(T(π/4), Vec{2, T}(1,2)))::Quaternion{T} == (@inferred quaternion(T(π/4), Vec{3, T}(1,2,0)))::Quaternion{T}
        @test (@inferred quaternion(T, π/4, Vec(1,2)))::Quaternion{T} == (@inferred quaternion(T, π/4, Vec(1,2,0)))::Quaternion{T}
        q = (@inferred quaternion(T(π/4), Vec{3, T}(1,2,3)))::Quaternion{T}
        q = (@inferred quaternion(T, π/4, Vec(1,2,3)))::Quaternion{T}
        @test length(q) == 4
        @test size(q) == (4,)
        @test norm(q) ≈ 1
        @test q/q ≈ 1
        q = (@inferred quaternion(T(π/4), Vec{3, T}(1,2,3), normalize = false))::Quaternion{T}
        q = (@inferred quaternion(T, π/4, Vec(1,2,3), normalize = false))::Quaternion{T}
        @test !(norm(q) ≈ 1)
        @test q/q ≈ 1
        @test q ≈ (@inferred quaternion(T(45), Vec{3, T}(1,2,3), normalize = false, degree = true))::Quaternion{T}
        @test q ≈ (@inferred quaternion(T, 45, Vec(1,2,3), normalize = false, degree = true))::Quaternion{T}

        q = quaternion(rand(T), rand(Vec{3, T}))
        p = quaternion(rand(T), rand(Vec{3, T}), normalize = false)

        # conversion
        @test (@inferred convert(Quaternion{T}, q))::Quaternion{T} == q
        @test (@inferred convert(Quaternion{T}, 3))::Quaternion{T} == Quaternion{T}(3,0,0,0)

        # promotion
        @test (@inferred promote_rule(Quaternion{T}, T)) == Quaternion{T}
        @test (@inferred promote_rule(Quaternion{T}, Int)) == Quaternion{T}
        @test (@inferred promote_rule(Quaternion{T}, Quaternion{T})) == Quaternion{T}
        @test (@inferred promote_rule(Quaternion{T}, Quaternion{Int})) == Quaternion{T}
        @test ((@inferred promote(q, T(3)))::NTuple{2, Quaternion{T}})[2] == Quaternion(3.0,0,0,0)
        @test ((@inferred promote(q, 3))::NTuple{2, Quaternion{T}})[2] == Quaternion(3.0,0,0,0)

        # math operations
        @test (@inferred +q)::Quaternion{T} === q
        @test (@inferred +p)::Quaternion{T} === p
        @test (@inferred -q)::Quaternion{T} == -1 * q
        @test (@inferred -p)::Quaternion{T} == -1 * p
        @test (@inferred q + p)::Quaternion{T} == Quaternion(Vec(q) + Vec(p))
        @test (@inferred q - p)::Quaternion{T} == Quaternion(Vec(q) - Vec(p))
        @test (@inferred 2 * q)::Quaternion{T} == Quaternion(2 * Vec(q))
        @test (@inferred 2 * p)::Quaternion{T} == Quaternion(2 * Vec(p))
        @test (@inferred q * 2)::Quaternion{T} == Quaternion(Vec(q) * 2)
        @test (@inferred p * 2)::Quaternion{T} == Quaternion(Vec(p) * 2)
        @test (@inferred q / 2)::Quaternion{T} == Quaternion(Vec(q) / 2)
        @test (@inferred p / 2)::Quaternion{T} == Quaternion(Vec(p) / 2)
        @test (@inferred norm(q))::T ≈ (@inferred abs(q))::T
        @test (@inferred norm(p))::T ≈ (@inferred abs(p))::T

        # angleaxis
        @test (quaternion((@inferred angleaxis(q))::Tuple{T, Vec{3, T}}...)) ≈ q

        a = rand(T)
        @test (@inferred exp(log(q)))::Quaternion{T} ≈ q
        @test (@inferred exp(log(p)))::Quaternion{T} ≈ p
        @test (@inferred exp(a + q))::Quaternion{T} ≈ exp(a) * exp(q)
        @test (@inferred exp(a + p))::Quaternion{T} ≈ exp(a) * exp(p)
        @test (@inferred log(a * q))::Quaternion{T} ≈ log(a) + log(q)
        @test (@inferred log(a * p))::Quaternion{T} ≈ log(a) + log(p)
        @test (@inferred exp(Quaternion{T}(1,0,0,0)))::Quaternion{T} ≈ exp(1)

        Rq = rotmat(q)
        Rp = rotmat(p)
        r = p * q

        for dim in (2, 3)
            # check multiplications
            x = rand(Vec{dim, T})
            x3 = ToVec3(x)
            @test (q * x / q).vector ≈ Rq ⋅ x3
            @test (p * x / p).vector ≈ Rp ⋅ x3
            @test (r * x / r).vector ≈ Rp ⋅ Rq ⋅ x3
            @test (q * x * inv(q)).vector ≈ Rq ⋅ x3
            @test (p * x * inv(p)).vector ≈ Rp ⋅ x3
            @test (r * x * inv(r)).vector ≈ Rp ⋅ Rq ⋅ x3
            # inverse of rotation
            @test (inv(q) * x * q).vector ≈ inv(Rq) ⋅ x3
            # check order of multiplications
            @test ((q * x) / q).vector ≈ Rq ⋅ x3
            @test (q * (x / q)).vector ≈ Rq ⋅ x3
            # rotate
            @test rotate(x, q) ≈ rotate(x, Rq)
            @test rotate(x, p) ≈ rotate(x, Rp)
            @test rotate(x, r) ≈ rotate(x, Rp ⋅ Rq)
            @test rotate(x, inv(q)) ≈ rotate(x, inv(Rq))
            @test rotate(x, inv(p)) ≈ rotate(x, inv(Rp))
            @test rotate(x, inv(r)) ≈ rotate(x, inv(Rp ⋅ Rq))
        end
        # test with rotmat(θ, n)
        θ = rand(T)
        n = rand(Vec{3, T})
        @test rotmat(θ, n) ≈ rotmat(quaternion(θ, n))
    end
end


text:
using BibParser, BibParser.BibTeX
using Test

for file in ["test.bib"]
    println("\nstart $file")
    parsed = parse_file("../examples/$file")

    for (k,e) in pairs(parsed)
        println("|--->$(e.id)")
        for fn in fieldnames(typeof(e))
            println("|     |----$(string(fn)) -> parsed!") # TODO: better printing
        end
        for (n, v) in pairs(e.fields)
            println("|     |----$n -> $v")
        end
        println("|")
    end
end


text:
# state/control discretization module
module dis

export Discretization

struct Discretization
  pt::AbstractArray{Int, 1}      # points per
  mi::AbstractArray{Float64, 1}  # minimum
  mx::AbstractArray{Float64, 1}  # maximum

  dim::Int                       # dimension
  thr::AbstractArray{Int, 1}     # threshold
  rng::AbstractArray{Float64, 1} # range
end

function Discretization(pt, mi, mx)
  @assert length(pt) == length(mi) == length(mx)
  dim = length(pt)
  for i in 1:dim
    @assert pt[i] > 1
    @assert mx[i] > mi[i]
  end

  dim = length(pt)
  thr = circshift(cumprod(pt), 1)
  thr[1] = 1
  rng = (mx - mi)

  Discretization(pt, mi, mx, dim, thr, rng)
end

# clamp state to within the discretization region
function clamp_x!(d::Discretization, x::AbstractArray{Float64, 1})
  @assert d.dim == length(x)
  for i in 1:d.dim
    if x[i] > d.mx[i]
      x[i] = d.mx[i]
    elseif x[i] < d.mi[i]
      x[i] = d.mi[i]
    end
  end
  return x
end

function clamp_x(d::Discretization, x::AbstractArray{Float64, 1})
  nx = copy(x)
  clamp_x!(nx)
  return nx
end

# linear state and discrete state vector
function ls2xd!(d::Discretization, xd::AbstractArray{Int, 1}, ls::Int)
  @assert d.dim == length(xd)
  for i in d.dim:-1:1
    xd[i] = div(ls, d.thr[i])
    ls = rem(ls, d.thr[i])
  end
  return xd
end

function ls2xd(d::Discretization, ls::Int)
  xd = fill(0, d.dim)
  ls2xd!(d, xd, ls)
  return xd
end

function xd2ls(d::Discretization, xd::AbstractArray{Int, 1})
  ls = 0
  for i in 1:d.dim
    @assert xd[i] < d.pt[i] && xd[i] >= 0
    ls += xd[i] * d.thr[i]
  end
  return ls
end

# discrete state vector and continuous state vector
function x2xd!(d::Discretization, xd::AbstractArray{Int, 1},
               x::AbstractArray{Float64, 1})
  @assert d.dim == length(xd) == length(x)
  for i in 1:d.dim
    @assert x[i] >= d.mi[i]
    @assert x[i] <= d.mx[i]
    xd[i] = round(Int, (d.pt[i] - 1) * (x[i] - d.mi[i]) / (d.rng[i]))
  end
  return xd
end

function x2xd(d::Discretization, x::AbstractArray{Float64, 1})
  @assert d.dim == length(x)
  xd = fill(0, d.dim)
  x2xd!(d,  xd, x)
  return xd
end

function xd2x!(d::Discretization, x::AbstractArray{Float64, 1},
               xd::AbstractArray{Int, 1})
  @assert d.dim == length(x) == length(xd)
  for i in 1:d.dim
    @assert xd[i] < d.pt[i]
    @assert xd[i] >= 0
    x[i] = (xd[i] / (d.pt[i] - 1)) * d.rng[i] + d.mi[i]
  end
  return x
end

function xd2x(d::Discretization, xd::AbstractArray{Int, 1})
  @assert d.dim == length(xd)
  x = fill(0.0, d.dim)
  xd2x!(d, x, xd)
  return x
end

# continuous state vector and linear state
function ls2x!(d::Discretization, x::AbstractArray{Float64, 1}, ls::Int)
  xd = ls2xd(d, ls)
  xd2x!(d, x, xd)
  return x
end

function ls2x(d::Discretization, ls::Int)
  x = fill(0.0, d.dim)
  ls2x!(d, x, ls)
  return x
end

function x2ls(d::Discretization, x::AbstractArray{Float64, 1})
  xd = x2xd(d, x)
  ls = xd2ls(d, xd)
  return ls
end

# interpolation
function bounding_XD(d::Discretization, x::AbstractArray{Float64, 1})
  @assert length(x) == d.dim

  xd = x2xd(d, x)
  nx = xd2x(d, xd)

  XD = Array{AbstractArray{Int, 1}, 1}(undef, 2^(d.dim))
  for i in 1:length(XD)
    XD[i] = copy(xd)
  end
  #XD = [copy(xd) for i in 1:(2^(d.dim))]

  for i in 1:d.dim
    add = fill(0, d.dim)
    add[i] = 1
    use_add = 0

    if nx[i] >= x[i]
      add .*= -1
      use_add = 1
    end

    j = 0
    while j < length(XD)
      for k in 1:2^(i - 1)
        XD[j + k] += use_add * add
      end
      j += 2^(i - 1)
      use_add = use_add == 0 ? 1 : 0
    end
  end

  # account for a condition when a value is exactly at minimum
  for i in 1:d.dim
    malformed = false
    for j in 1:length(XD)
      if XD[j][i] == -1
        malformed = true
        break
      end
    end
    if malformed == true
      for j in 1:length(XD)
        XD[j][i] += 1
      end
    end
  end

  return XD
end

function bounding_X(d::Discretization, x::AbstractArray{Float64, 1})
  XD = bounding_XD(d, x)
  X = Array{AbstractArray{Float64, 1}, 1}(undef, length(XD))
  for i in 1:length(X)
    X[i] = xd2x(d, XD[i])
  end
  #X = [xd2x(xd, road) for xd in XD]
  return X
end

function interp(X::AbstractArray{AbstractArray{Float64, 1}, 1},
                x::AbstractArray{Float64, 1}, v::AbstractArray{Float64, 1})
  dim = Int(log2(length(X)))

  xi = fill(0.0, dim)
  for i in 1:dim
    xi[i] = (x[i] - X[1][i]) / (X[end][i] - X[1][i])
  end

  c = v # no need to copy
  for i in 1:dim
    len = 2^(dim - i)
    nc = fill(0.0, len)
    for j in 1:len
      nc[j] = c[2 * (j - 1) + 1] * (1.0 - xi[i]) + c[2 * (j - 1) + 2] * xi[i]
    end
    c = nc
  end

  return c[]
end

function interp(d::Discretization, x::AbstractArray{Float64, 1},
                v::AbstractArray{Float64, 1})
  X = bounding_X(d, x)
  return interp(X, x, v)
end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "SLEEF"
version = v"3.4.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/shibatch/sleef.git", "8df2bce4596ec6b86f357beebf2e428564f0a4ec"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/sleef
mkdir build-native
cd build-native
cmake \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_HOST_TOOLCHAIN} \
    -G Ninja \
    -DBUILD_TESTS=OFF \
    -DBUILD_DFT=TRUE \
    -DBUILD_QUAD=TRUE \
    ..
ninja all

cd $WORKSPACE/srcdir/sleef
if [[ $target == *w64-mingw32* ]]; then
    atomic_patch -p1 ../patches/lowercase-windows-include.patch
elif [[ $target == arm-* ]]; then
    atomic_patch -p1 ../patches/arm-neon32vfpv4.patch
fi
mkdir build-cross
cd build-cross
cmake \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -G Ninja \
    -DNATIVE_BUILD_DIR=$WORKSPACE/srcdir/sleef/build-native \
    -DSLEEF_SHOW_CONFIG=1 \
    -DBUILD_TESTS=OFF \
    -DBUILD_DFT=TRUE \
    -DBUILD_QUAD=TRUE \
    ..
ninja all
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
filter!(p -> !(arch(p) == "i686" && Sys.iswindows(p)), platforms) # i686-windows build fails

# The products that we will ensure are always built
products = [
    LibraryProduct("libsleef", :libsleef),
    LibraryProduct("libsleefquad", :libsleefquad),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    Dependency("MPFR_jll"),
    Dependency("OpenSSL_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"7.4",
    julia_compat="1.6")


text:
module TestRanges

using Test
using RLEVectors

@testset begin

# disjoin_length
@test disjoin_length([],[]) == 0
@test disjoin_length([4,5],[]) == 2
@test disjoin_length([],[9,10,12]) == 3
@test disjoin_length([1,2,3],[4,5,6]) == 6
@test disjoin_length([1,2,3],[3,5,6]) == 5
@test disjoin_length([1,3],[2,5,6]) == 5
@test disjoin_length(RLEVector([1,2], [1,3]), RLEVector([3,2,1], [2,5,6])) == 5

# disjoin
@test disjoin( [1,3,8,10],[3,4,9,10] ) == [1,3,4,8,9,10]
@test disjoin( [3,4,9,10], [1,3,8,10] ) == [1,3,4,8,9,10]

@test disjoin( RLEVector([1,2,3,4], [1,3,8,10]), RLEVector([4,5,6,7],[3,4,9,10])) == ([1,3,4,8,9,10], [1,2,3,3,4,4], [4,4,5,6,6,7])
@test disjoin( RLEVector([4,5,6,7],[3,4,9,10]), RLEVector([1,2,3,4], [1,3,8,10])) == ([1,3,4,8,9,10], [4,4,5,6,6,7], [1,2,3,3,4,4])

# rangeMeans

x = RLEVector( collect(2:2:20) )
y = RLEVector( collect(1:1:10) )

ranges = [ 1:1, 1:2, 2:3, 4:8, 1:10, 9:10, 10:10 ]
@test [ mean(x[ r ]) for r in ranges ] == rangeMeans( ranges, x )

end # testset

end # module


text:
using Test, Gtk

@testset "text" begin

@testset "GtkTextIter" begin
import Gtk: GtkTextIter, mutable

w = GtkWindow() #1.call a window
b = GtkTextBuffer() #2.call a textBuffer
b.text[String] = "test" #3.create a label on buffer's argument 'text'
v = GtkTextView(b) #4.Create a (text) view
@test v[:buffer, GtkTextBuffer] == b #

push!(w, v) #push changes: view v onto a window w
showall(w)#diplay function

its = GtkTextIter(b, 1) #q. Iter funcyiosnlity (is Sync functionality?)
ite = GtkTextIter(b, 2)

@test buffer(its) == b
@test get_gtk_property(its:ite, :text, String) == "t"
@test (its:ite).text[String] == "t"

splice!(b, its:ite)
@test b.text[String] == "est"

insert!(b, GtkTextIter(b, 1), "t")
@test b.text[String] == "test"

it = GtkTextIter(b)
@test get_gtk_property(it, :line) == 0 #lines are 0-based
@test get_gtk_property(it, :starts_line) == true

b.text[String] = "line1\nline2"
it = mutable(GtkTextIter(b))
set_gtk_property!(it, :line, 1)
@test get_gtk_property(it, :line) == 1

it1 = GtkTextIter(b, 1)
it2 = GtkTextIter(b, 1)
@test it1 == it2
it2 = GtkTextIter(b, 2)
@test (it1 == it2) == false
@test it1 < it2
it2 -= 1
@test mutable(it1) == it2

# tags
Gtk.create_tag(b, "big"; size_points = 24)
Gtk.create_tag(b, "red"; foreground = "red")
f(buffer)=Gtk.apply_tag(buffer, "big", GtkTextIter(b, 1), GtkTextIter(b, 6))
user_action(f, b)
Gtk.apply_tag(b, "red", GtkTextIter(b, 1), GtkTextIter(b, 6))
Gtk.remove_tag(b, "red", GtkTextIter(b, 1), GtkTextIter(b, 3))
Gtk.remove_all_tags(b, GtkTextIter(b, 4), GtkTextIter(b, 6))

# getproperty
@test it1.offset[Int] == 0 #Gtk indices are zero based
@test it2.offset[Int] == 0

it1 = mutable(it1)
it1.offset[Int] = 1
@test it1.offset[Int] == 1

mark = create_mark(b, it)
scroll_to(v, mark, 0, true, 0.0, 0.15)
scroll_to(v, it, 0, true, 0.0, 0.15)

# skip
skip(it2, 1, :line)
@test get_gtk_property(it2, :line) == 1
skip(it2, :backward_line)
@test get_gtk_property(it2, :line) == 0
skip(it2, :forward_line)
@test get_gtk_property(it2, :line) == 1
skip(it2, :forward_to_line_end)
it1 = GtkTextIter(b, get_gtk_property(it2, :offset)-1)
(it1:it2).text[String] == "2"

itc = convert(GtkTextIter, it2)
@test get_gtk_property(it2, :offset) == get_gtk_property(itc, :offset)

whats = [:forward_word_end, :backward_word_start, :backward_sentence_start, :forward_sentence_end]
for what in whats
    skip(mutable(it1), what)
end
whats = [:char,:line,:word,:word_cursor_position,:sentence,:visible_word,:visible_cursor_position,:visible_line,:line_end]
for what in whats
    skip(mutable(it1), 0, what)
end

# place_cursor
place_cursor(b, it2)
iter, strong, weak = Gtk.cursor_locations(v)
@test get_gtk_property(it2, :is_cursor_position) == true
@test b.cursor_position[Int] == get_gtk_property(it2, :offset)

# search
(found, its, ite) = Gtk.search(b, "line1", :backward)
@test found == true
@test (its:ite).text[String] == "line1"

place_cursor(b, ite)
(found, its, ite) = Gtk.search(b, "line2", :forward)
@test found == true
@test (its:ite).text[String] == "line2"

# GtkTextRange
range=its:ite
@test range[1] == 'l'
@test range[5] == '2'
@test_throws BoundsError range[10]

# selection
select_range(b, its, ite)
(selected, start, stop) = selection_bounds(b)
@test selected == true
@test (start:stop).text[String] == "line2"

insert!(v, start, "inserted text")

# coords
wx, wy = Gtk.buffer_to_window_coords(v, 3, 2, 2) #translate(buffer -> window )
bx, by = Gtk.window_to_buffer_coords(v, wx, wy)  #translate window -> buffer
@test bx == 3 && by == 2

iter = Gtk.text_iter_at_position(v, 3, 2) #Iter at position

destroy(w) #don't forget to destroy before the end
end

end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
function MaxEval(arg0::jint)
    return MaxEval((jint,), arg0)
end

function get_max_eval(obj::MaxEval)
    return jcall(obj, "getMaxEval", jint, ())
end

function unlimited(::Type{MaxEval})
    return jcall(MaxEval, "unlimited", MaxEval, ())
end


text:
# # Estimate Binomial draw probabilility

using DynamicHMCModels

Random.seed!(1356779)

# Define a structure to hold the data.

Base.@kwdef struct BernoulliProblem
    "Total number of draws in the data."
    n::Int
    "Number of draws ' == 1' "
    obs::Vector{Int}
end;

# Write a function to return properly dimensioned transformation.

make_transformation(model::BernoulliProblem) =
    as((p = as𝕀, ))

# Add data

model = BernoulliProblem(; n = 9, obs = rand(Binomial(9, 2/3), 3))

# Make the type callable with the parameters *as a single argument*.

function (model::BernoulliProblem)(θ)
    @unpack n, obs = model        # extract the data
    @unpack p = θ
    loglikelihood(Binomial(n, p), obs)
end

# Use a flat priors (the default, omitted) for α

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P);

# Sample chain

results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000;
  reporter = NoProgressReport()
)

posterior = P.transformation.(results.chain)

# Create Particles NamedTuple object

println()
p = as_particles(posterior)
p |> display

println()
DynamicHMC.Diagnostics.EBFMI(results.tree_statistics) |> display

println()
DynamicHMC.Diagnostics.summarize_tree_statistics(results.tree_statistics) |> display


text:
## ################################################
# Load Packages

using TimeSeriesEcon
using Plots

## ################################################
#### Part 1: Moment In Time (MIT) and TSeries

### Initialize MITs and TSeries

# `MIT` is a primitive type based on 64-bit signed integers that
# represents discrete dates. There are two ways to initialize
# `MIT`s: (1) directly `2020M8`, or (2) using the functional form `mm(2020, 8)`.

# `Unit` Frequency
mit_integer   = 2000U  # ii(2000)

# `Monthly` Frequency
mit_monthly   = 2020M8 # mm(2020, 8)

# `Quarterly` Frequency
mit_quarterly = 2020Q3 # qq(2020, 3)

# `Yearly` Frequency
mit_yearly    = 2020Y  # yy(2020)

# `TSeries` is subtype of AbstractVector and represents 1-dimensional time-series. A key
# feature of `TSeries` is the ability to use `MIT`s as indices to get and set values.

# Note: TSeries converts all values to `Float64`. The automatic conversion feature might be
# changed in the future.

series_monthly = TSeries(2020M1, rand(1:10, 6))

# There are various ways to initialize `TSeries` - here, we provided an `MIT` that represents
# the first date in the series and a vector of random values.

### Frequency

# Every instance of `MIT` and `TSeries` is equipped with `Frequency` information, stored as a parameter.

subtypes(Frequency)

frequencyof(2020M8)
frequencyof(series_monthly)

# As such, we can avoid performing operations on `TSeries` or `MIT`s of different frequencies.
# Also, haing `Frequency` parameter simplifies the retrieval of frequency specific information.

year(2020M8)
period(2020M8)

ppy(2020M8) # number of periods per year

## ################################################
#### Part 2: Indexing using MITs

### Access
# Just as integers are used to index into Julia vectors, `MIT`s are
# used to index into `TSeries`.

# Indexing using a single `MIT` returns a float value associated with that date.
# Indexing using a range of `MIT`s will return another `TSeries` instance

series_monthly

series_monthly[2020M1] # 5
series_monthly[2020M1:2020M3] # TSeries(2020M8, [5.0, 4.0, 2.0])

### Assign

series_monthly[2020M1] = -1;
series_monthly

series_monthly[2020M2:end] = -1;
series_monthly

# You can also assign values outside of the bounds that were initially declared.
series_monthly[end + 2] = -1;
series_monthly

## ################################################
#### Part 3: Plots and helpful functions

### Plotting support

# Using the Plots package, we can plot multiple `TSeries`
# with varying frequency.

plot(TSeries(2000Q1, rand(1:3, 10)),
     TSeries(2000M1, rand(4:6, 30)),
     legend=true,
     title="TSeries Plot",
     labels=["Quarterly", "Monthly"],
     size=(600, 400)
    )

png_fname = joinpath(pwd(), "TimeSeriesEcon", "tseries.png")
savefig(png_fname)

### Conversions

# !!! Note: Support for frequency conversion from low to high
# is under development.

tsmonthly = TSeries(2020M1, collect(1:12));

# Monthly -> Quarterly
tsquarterly = convert(TSeries{Quarterly}, tsmonthly)

# Monthly -> Yearly
tsyearly = convert(TSeries{Yearly}, tsmonthly)


text:
cd(@__DIR__)

@testset "reading negative delta" begin
    midi = load("negative_delta.mid")
    @test midi.tracks[1].events[3].dT == -48
end

@testset "writing negative delta" begin
    stream = IOBuffer()
    @test_throws ErrorException MIDI.writeevent(stream, MIDI.NoteOnEvent(-1, 0x90, [0x01, 0x01]))
    @test_throws ErrorException MIDI.writeevent(stream, MIDI.SetTempoEvent(-1, 0x51, [0x00, 0x00, 0x01]))
    @test_throws ErrorException MIDI.writeevent(stream, SysexEvent(-1, [0x01]))
    @test_throws ErrorException MIDI.writevariablelength(stream, -1)
    @test_throws ErrorException MIDI.writevariablelength(stream, Int(0x1FFFFFFF))
    close(stream)
end


text:
# This file was generated, do not modify it. # hide
krb_best = fitted_params(mtm).best_model
@show krb_best.knn_model.K
@show krb_best.ridge_model.lambda
@show krb_best.knn_weight

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
"""
DistributedFactorGraphs.jl provides a flexible factor graph API for use in the Caesar.jl ecosystem.

The package supplies:

- A standardized API for interacting with factor graphs
- Implementations of the API for in-memory and database-driven operation
- Visualization extensions to validate the underlying graph

"""
module DistributedFactorGraphs

##==============================================================================
## imports
##==============================================================================

using Base
using DocStringExtensions
using Requires
using Dates
using TimeZones
using Distributions
using Reexport
using JSON
using Unmarshal
using JSON2 # JSON2 requires all properties to be in correct sequence, can't guarantee that from DB.
using LinearAlgebra
using SparseArrays
using UUIDs
using Pkg
using TensorCast

# used for @defVariable
import ManifoldsBase
import ManifoldsBase: AbstractManifold, manifold_dimension
export AbstractManifold, manifold_dimension

import Base: getindex

##==============================================================================
# Exports
##==============================================================================
const DFG = DistributedFactorGraphs
export DFG
##------------------------------------------------------------------------------
## DFG
##------------------------------------------------------------------------------
export AbstractDFG
export AbstractParams, NoSolverParams
export AbstractBlobStore

# accessors & crud
export getUserId, getRobotId, getSessionId
export getDFGInfo
export getDescription, setDescription!,
       getSolverParams, setSolverParams!,
       getUserData, setUserData!,
       getRobotData, setRobotData!,
       getSessionData, setSessionData!,
       getAddHistory
export getBlobStore,
       addBlobStore!,
       updateBlobStore!,
       deleteBlobStore!,
       emptyBlobStore!,
       listBlobStores

# TODO Not sure these are needed or should work everywhere, implement in cloud?
export updateUserData!, updateRobotData!, updateSessionData!, deleteUserData!, deleteRobotData!, deleteSessionData!
export emptyUserData!, emptyRobotData!, emptySessionData!

# Graph Types: exported from modules or @reexport
export InMemoryDFGTypes, DefaultDFG

# AbstractDFG Interface
export exists,
       addVariable!, addFactor!,
       getVariable, getFactor,
       updateVariable!, updateFactor!,
       deleteVariable!, deleteFactor!,
       listVariables, listFactors,
       listSolveKeys, listSupersolves,
       getVariables, getFactors,
       isVariable, isFactor

export getindex

export isConnected

export getBiadjacencyMatrix

#summary structure
export DFGSummary

export getSummary, getSummaryGraph

# Abstract Nodes
export DFGNode, AbstractDFGVariable, AbstractDFGFactor

# Variables
export DFGVariable, DFGVariableSummary, SkeletonDFGVariable

# Factors
export DFGFactor, DFGFactorSummary, SkeletonDFGFactor

# Common
export getSolvable, setSolvable!, isSolvable
export getInternalId
export getVariableLabelNumber

# accessors
export getLabel, getTimestamp, setTimestamp, setTimestamp!, getTags, setTags!

# Node Data
export isSolveInProgress, getSolveInProgress

# CRUD & SET
export listTags, mergeTags!, removeTags!, emptyTags!

#this isn't acttually implemented. TODO remove or implement
export addTags!

##------------------------------------------------------------------------------
# Variable
##------------------------------------------------------------------------------
# Abstract Variable Data
export InferenceVariable

# accessors
export getSolverDataDict, setSolverData!
export getVariableType, getVariableTypeName

export getSolverData

export getVariableType

# VariableType functions
export getDimension, getManifold, getPointType
export getPointIdentity, getPoint, getCoordinates

export getManifolds # TODO Deprecate?

# Small Data CRUD
export SmallDataTypes, getSmallData, addSmallData!, updateSmallData!, deleteSmallData!, listSmallData, emptySmallData!
export getSmallData, setSmallData!

# CRUD & SET
export getVariableSolverData,
       addVariableSolverData!,
       updateVariableSolverData!,
       deleteVariableSolverData!,
       listVariableSolverData,
       mergeVariableSolverData!,
       deepcopySolvekeys!,
       deepcopySupersolve!

# PPE
##------------------------------------------------------------------------------
# PPE TYPES
export AbstractPointParametricEst
export MeanMaxPPE
export getPPEMax, getPPEMean, getPPESuggested, getLastUpdatedTimestamp

# accessors

export getPPEDict, getVariablePPEDict, getVariablePPE

# CRUD & SET
export getPPE,
       getVariablePPE,
       addPPE!,
       updatePPE!,
       deletePPE!,
       listPPEs,
       mergePPEs! #TODO look at rename to just mergePPE to match with cloud?

# Variable Node Data
##------------------------------------------------------------------------------
export VariableNodeData, PackedVariableNodeData

export packVariableNodeData, unpackVariableNodeData

export getSolvedCount, isSolved, setSolvedCount!, isInitialized, isMarginalized, setMarginalized!

export getNeighborhood, getNeighbors, _getDuplicatedEmptyDFG
export findFactorsBetweenNaive
export copyGraph!, deepcopyGraph, deepcopyGraph!, buildSubgraph, mergeGraph!
# Big Data
##------------------------------------------------------------------------------
export addDataEntry!, getDataEntry, updateDataEntry!, deleteDataEntry!, getDataEntries, listDataEntries, hasDataEntry, hasDataEntry
# convenience wrappers
export addDataEntry!
# aliases
export addData!

##------------------------------------------------------------------------------
# Factors
##------------------------------------------------------------------------------
# Factor Data
export GenericFunctionNodeData, PackedFunctionNodeData, FunctionNodeData
export AbstractFactor, AbstractPackedFactor
export FunctorInferenceType, PackedInferenceType # NOTE might be deprecated in the future (const aliases)
export AbstractPrior, AbstractRelative, AbstractRelativeRoots, AbstractRelativeMinimize
export FactorOperationalMemory

# accessors
export getVariableOrder
export getFactorType, getFactorFunction

# Node Data
export mergeVariableData!, mergeGraphVariableData!

# Serialization type conversion
export convertPackedType, convertStructType

##------------------------------------------------------------------------------
## Other utility functions
##------------------------------------------------------------------------------

# Sort
export natural_lt, sortDFG

# Validation
export isValidLabel

## List
export ls, lsf, ls2
export lsTypes, lsfTypes, lsTypesDict, lsfTypesDict
export lsWho, lsfWho
export isPrior, lsfPriors
export hasTags, hasTagsNeighbors

## Finding
export findClosestTimestamp, findVariableNearTimestamp

# Serialization
export packVariable, unpackVariable, packFactor, unpackFactor
export rebuildFactorMetadata!
export @defVariable

# File import and export
export saveDFG, loadDFG!, loadDFG
export toDot, toDotFile

# shortest path
export findShortestPathDijkstra
export isPathFactorsHomogeneous

# Comparisons
export
    compare,
    compareField,
    compareFields,
    compareAll,
    compareAllSpecial,
    compareVariable,
    compareFactor,
    compareAllVariables,
    compareSimilarVariables,
    compareSubsetFactorGraph,
    compareSimilarFactors,
    compareFactorGraphs

## Deprecated exports should be listed in Deprecated.jl if possible, otherwise here

## CustomPrinting.jl
export printFactor, printVariable, printNode

##==============================================================================
## Files Includes
##==============================================================================

# Entities

include("entities/AbstractDFG.jl")

include("entities/DFGFactor.jl")

include("entities/DFGVariable.jl")

include("entities/AbstractDFGSummary.jl")

include("services/AbstractDFG.jl")

# In Memory Types
# include("../attic/GraphsDFG/GraphsDFG.jl")
# @reexport using .GraphsDFGs
include("LightDFG/LightDFG.jl")
@reexport using .LightDFGs
#supported in Memory fg types
# const InMemoryDFGTypes = Union{GraphsDFG, LightDFG}
const InMemoryDFGTypes = LightDFG
const DefaultDFG = LightDFG

# Common includes
include("services/CommonAccessors.jl")
include("services/Serialization.jl")
include("services/DFGVariable.jl")
include("services/DFGFactor.jl")
include("Deprecated.jl")
include("services/CompareUtils.jl")

# Include the FilesDFG API.
include("FileDFG/FileDFG.jl")

# Custom show and printing for variable factor etc.
include("services/CustomPrinting.jl")

# To be moved as necessary.
include("Common.jl")

# Data Blob extensions
include("DataBlobs/DataBlobs.jl")

if get(ENV, "DFG_USE_CGDFG", "") == "true"
    @info "Detected ENV[\"DFG_USE_CGDFG\"]: Including optional CloudGraphsDFG (LGPL) Driver"
    include("CloudGraphsDFG/CloudGraphsDFG.jl")
    @reexport using .CloudGraphsDFGs
end

function __init__()
    @require GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231" begin
        @info "DistributedFactorGraphs.jl is adding tools using GraphPlot.jl"
        include("DFGPlots/DFGPlots.jl")
        @reexport using .DFGPlots
    end
end

end


text:
"""
    BrooksAndCorey(; n, Ks=1, l=1, α=1, θr=0, θs=1) <: UnsaturatedFlowModel
    BrooksAndCorey(; n, k, l=1, α=1, θr=0, θs=1, ρ=1e3, μ=1e-3, g=9.81) <: UnsaturatedFlowModel

Create a Brooks and Corey porous model.

# Keyword arguments
- `n`: _n_ parameter.
- `Ks=1`: saturated hydraulic conductivity.
- `k`: intrinsic permeability.
- `l=1`: _l_ parameter.
- `α=1` (`\\alpha<tab>`): _α_ parameter.
- `θr=0` (`\\theta<tab>r`): residual moisture content.
- `θs=1` (`\\theta<tab>s`): moisture content when saturated.
- `ρ=1e3` (`\\rho<tab>`): density of the fluid.
- `μ=1e-3` (`\\mu<tab>`): dynamic viscosity of the fluid.
- `g=9.81`: magnitude of the gravitational acceleration.

# References
BROOKS, R.; COREY, T. Hydraulic properties of porous media.
Hydrology Papers, Colorado State University, 1964, vol. 24, p. 37.
"""
struct BrooksAndCorey{_Tnl,_Tα,_TKs,_Tθ} <: UnsaturatedFlowModel
    n::_Tnl
    l::_Tnl
    α::_Tα
    Ks::_TKs
    θr::_Tθ
    θs::_Tθ

    function BrooksAndCorey(; n, l=1, α=1, Ks=nothing, k=nothing, θr=0, θs=1, ρ=1e3, μ=1e-3, g=9.81)

        @argcheck α>zero(α)
        @argcheck θr<θs

        Ks = _asKs(Ks=Ks, k=k, ρ=ρ, μ=μ, g=g)

        new{promote_type(typeof(n),typeof(l)),typeof(α),typeof(Ks),promote_type(typeof(θr),typeof(θs))}(n,l,α,Ks,θr,θs)
    end
end

function θh(pm::BrooksAndCorey, h)
    if h ≥ -1/pm.α
        return pm.θs
    end

    Se = 1/(pm.α*(-h))^pm.n

    return pm.θr + Se*(pm.θs - pm.θr)
end

function hθ(pm::BrooksAndCorey, θ)
    Se = (θ - pm.θr)/(pm.θs - pm.θr)

    return -1/(pm.α*Se^(1/pm.n))
end

function Ch(pm::BrooksAndCorey, h)
    if h ≥ -1/pm.α
        return zero(1/h)
    end

    return -pm.n/h*(pm.θs - pm.θr)/(pm.α*(-h))^pm.n
end

function Kθ(pm::BrooksAndCorey, θ)
    Se = (θ - pm.θr)/(pm.θs - pm.θr)

    return pm.Ks*Se^(2/pm.n + pm.l + 2)
end

function Kh(pm::BrooksAndCorey, h)
    if h ≥ -1/pm.α
        return pm.Ks
    end

    return Kθ(pm, θh(pm, h))
end

function Dθ(pm::BrooksAndCorey, θ)
    Se = (θ - pm.θr)/(pm.θs - pm.θr)
    return pm.Ks*Se^(1/pm.n + pm.l + 1)/(pm.α*pm.n*(pm.θs - pm.θr))
end


text:
"""
    print_results(
        ata_model::NoObjModel;
        results_folder = "results",
        sim_pool::DataFrames.DataFrame = DataFrame(),
    )

# Description

Print the features of the assembled tests.

# Arguments

- **`ata_model::NoObjModel`** : Required. The model built with `ATA` fuctions, `ata_model.design` matrix must be `IxT` or `nfsxT` if the items are grouped by friend sets.
- **`results_folder`** : Optional. Default: "results". The folder in which the output is stored.
- **`sim_pool::DataFrames.DataFrame`** : Optional. Default: DataFrame(). The pool with true item paramaters. For simulation studies.
"""
function print_results(
    ata_model::NoObjModel;
    results_folder = "results",
    sim_pool::DataFrames.DataFrame = DataFrame(),
)
    if !isdir(results_folder)
        mkdir(results_folder)
    else
        success!(
            ata_model,
            string(
                "There is already a folder with this name, files in ",
                results_folder,
                " will be overwritten.",
            ),
        )

    end
    if size(ata_model.output.design, 1) > 0
        T = ata_model.settings.T
        n_items = ata_model.settings.n_items
        length_max = [ata_model.constraints[t].length_max for t = 1:T]
        categories = ata_model.output.categories
        n_fs = size(ata_model.settings.fs.items, 1)
        if n_fs < ata_model.settings.n_items
            n = n_fs * T
            design = reshape(ata_model.output.design, n_fs, T)
            new_design = zeros(Float64, n_items, T)
            for t = 1:T
                new_design[
                    vcat(
                        ata_model.settings.fs.items[findall(
                            design[:, t] .== one(Float64),
                        )]...,
                    ),
                    t,
                ] .= one(Float64)
            end
            design = new_design
            DelimitedFiles.writedlm(string(results_folder, "/designItemLevel.csv"), design)
        else
            n = n_items * T
            design = reshape(ata_model.output.design, n_items, T)
        end
        #overlap
        olMatrixOUT = zeros(Int64, T, T)
        for t = 1:T
            for v = 1:T
                olMatrixOUT[t, v] = design[:, t]' * design[:, v]
            end
        end
        #DelimitedFiles.writedlm(string(results_folder,"/olMatrixOUT.csv"),olMatrixOUT)
        #TIF e ICF
        if size(ata_model.settings.irt.parameters, 1) > 0
            points = [-1.5, -1.0, -0.5, 0, 0.5, 1.0, 1.5]
            ICF = Vector{Matrix{Float64}}(undef, T)
            for t = 1:T
                ICF[t] = zeros(size(points, 1), n_items)
                for k = 1:size(points, 1)
                    ICF[t][k, :] = (item_char(
                        ata_model.settings.irt.parameters,
                        points[k],
                        model = ata_model.settings.irt.model,
                        parametrization = ata_model.settings.irt.parametrization,
                        D = ata_model.settings.irt.D,
                    )[1][
                        :,
                        :,
                        1,
                    ])
                end
            end
            IIF = Vector{Matrix{Float64}}(undef, T)

            for t = 1:T
                IIF[t] = zeros(7, n_items)
                for k = 1:7
                    IIF[t][k, :] = item_info(
                        ata_model.settings.irt.parameters,
                        points[k],
                        model = ata_model.settings.irt.model,
                        parametrization = ata_model.settings.irt.parametrization,
                        D = ata_model.settings.irt.D,
                    )
                end
            end
            #save values
            estimated = Float64[]
            if size(sim_pool, 1) > 0
                True = Float64[]
            end
            for t = 1:T
                K = size(points, 1)
                for k = 1:K
                    estimated = vcat(estimated, IIF[t][k, :]' * design[:, t])
                    if size(sim_pool, 1) > 0
                        True = vcat(
                            True,
                            item_info(
                                sim_pool,
                                points[k],
                                model = ata_model.settings.irt.model,
                                parametrization = ata_model.settings.irt.parametrization,
                                D = ata_model.settings.irt.D,
                            )' * design[:, t],
                        )
                    end
                end
            end
            tif_at_theta_pts = DataFrames.DataFrame(estimated = estimated)
            if size(sim_pool, 1) > 0
                tif_at_theta_pts[!, :True] = True
                DataFrames.rename(tif_at_theta_pts, Symbol.(["estimated", "true"]))
            end
            CSV.write(string(results_folder, "/tif_at_theta_pts.csv"), tif_at_theta_pts)
            #expected score
            es_print_irt = Vector{Vector{Float64}}(undef, T)
            for t = 1:T
                if size(ICF[t], 1) > 0
                    es_print_irt[t] = zeros(size(ICF[t], 1))
                    for k = 1:size(ICF[t], 1)
                        es_print_irt[t][k] =
                            ((ICF[t][k, :]'*design[:, t]))[1] / sum(design[:, t])
                    end
                end
            end
        end

        #number of items
        n = sum(design, dims = 1)
        design_items = zeros(Int64, Int(maximum(n)), T)
        #design_items=Array{String,2}(nothing,n_max,T)
        for t = 1:T
            design_items[collect(1:Int(n[t])), t] .= sort!(findall(design[:, t] .== 1))#sort!(CODEVar[findall(design[:,t].==1)])
        end
        open(string(results_folder, "/results.txt"), "w") do io
            write(io, "Length")
            write(io, "\r\n")
            DelimitedFiles.writedlm(io, [1:T, Int.(n)'], "\t")
            write(io, "\r\n")
            write(io, "design")
            write(io, "\r\n")
            DelimitedFiles.writedlm(io, design_items, "\t")
            write(io, "\r\n")
            if size(ata_model.settings.irt.parameters, 1) > 0
                write(io, "Expected score")
                write(io, "\r\n")
                for t = 1:T
                    DelimitedFiles.writedlm(
                        io,
                        vcat(t, round.(es_print_irt[t], digits = 3))',
                        "\t",
                    )
                    write(io, "\r\n")
                end
            end
            if size(categories, 1) > 0
                write(io, "Categorical variables")
                write(io, "\r\n")
                for cats in categories
                    write(io, cats)
                    write(io, "\r\n")
                    vals1 =
                        setdiff(unique(ata_model.settings.bank[!, Symbol(cats)]), [missing])
                    vals = copy(vals1)
                    for t = 1:T
                        valscat = Int64[]
                        for val in vals1
                            valscat_t =
                                ata_model.settings.bank[!, Symbol(cats)][design[
                                    :,
                                    t,
                                ].==1] .== val
                            valscat_t = valscat_t[.!ismissing.(valscat_t)]
                            if size(valscat_t, 1) > 0
                                valscat = vcat(valscat, sum(valscat_t))
                            else
                                valscat = vcat(valscat, 0)
                            end
                        end
                        vals = hcat(vals, valscat)
                    end
                    DelimitedFiles.writedlm(io, vals, "\t")
                    #DelimitedFiles.writedlm(io,hcat(vals,[sum(skipmissing(ata_model.settings.bank[Symbol(cats)][design[:,t].==1].== i))   for i in vals,t=1:T]),", ")
                    write(io, "\r\n")
                end

                write(io, "Overlap matrix")
                write(io, "\r\n")
                olMatrixOUT =
                    hcat(vcat("", string.(collect(1:T))), vcat(collect(1:T)', olMatrixOUT))
                DelimitedFiles.writedlm(io, olMatrixOUT, "\t")
                write(io, "\r\n")
            end
            write(io, "Item use")
            write(io, "\r\n")
            iu = [sum(design[i, :]) for i = 1:ata_model.settings.n_items]
            ids = sortperm(iu; rev = true)
            sort!(iu; rev = true)
            for i = 1:ata_model.settings.n_items
                write(io, string(ids[i]))
                write(io, "\t")
                write(io, string(iu[i]))
                write(io, "\r\n")
            end
        end
    else
        println("No solution found")
    end
end


text:
"""
"""
struct WCycle <: AbstractCycle end

"""
"""
struct WkCycle <: AbstractCycle
    k::Int
end

""
function _core_cycle!(problem, level, ::WCycle)
    process_fine!(problem, level)
    _descend!(problem, level, WCycle())
    process_coarse!(problem, level)

    process_uncoarse!(problem, level)
    _descend!(problem, level, WCycle())
    process_coarse!(problem, level)
    return nothing
end

""
@inline function _repeatcycle!(problem, level, cycle_t::WkCycle)
    for i in 1:(cycle_t.k)
        process_uncoarse!(problem, level)
        _descend!(problem, level, cycle_t)
        process_coarse!(problem, level)
    end
    return nothing
end

""
function _core_cycle!(problem, level, cycle_t::WkCycle)
    process_fine!(problem, level)
    _descend!(problem, level, cycle_t)
    process_coarse!(problem, level)

    _repeatcycle!(problem, level, cycle_t)
    return nothing
end


text:
# # Multiple Sequence Alignments
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__notebooks/09_MSA.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__notebooks/09_MSA.ipynb)
#
# For this example, we use the Multiple Sequence Alignment (MSA) of the
# catalytic domain from the *Cytosol aminopeptidase family*. We use the
# `downloadpfam` function of the MIToS' Pfam module to get the annotated
# [Stockholm](https://en.wikipedia.org/wiki/Stockholm_format) file from the
# [Pfam](pfam.xfam.org/) database:

using MIToS.Pfam
const pfam_file = downloadpfam("PF00883")

# The `read` function from MIToS' MSA allows parsing multiple sequence
# alignments. The file format is indicated by the second argument. There is
# no need for decompressing the file before reading it:

using MIToS.MSA
msa = read(pfam_file, Stockholm)

# Pfam alignments are generated using [HMMER3](http://hmmer.org/) Hidden Markov
# Model (HMM) searches, have insertion columns (indicated by dots and lowercase
# letters):

first_seq = sequencenames(msa)[1]

# **TIP:** You can use [TranscodingStreams](https://bicycle1885.github.io/TranscodingStreams.jl/stable/)
# to read compressed files. In particular, the *CodecZlib* package has the
# needed codecs for reading gzipped files.
# For example, the next cell is similar to using `zgrep` in Linux:
# ```bash
# zgrep ^A0A094YQ53_BACAO/182-489 PF00883.stockholm.gz
# ```

using CodecZlib

open(GzipDecompressorStream, pfam_file) do file
	for line in eachline(file)
		if startswith(line, first_seq)
			println(line)
		end
	end
end

# **TIP:** If you are on Linux and you have zgrep installed, you can run it as an
# [external command](https://docs.julialang.org/en/v1/manual/running-external-programs/index.html)
# and use interpolation:
# ```julia
# run(`zgrep ^$first_seq $pfam_file`)
# ```
# or using the [shell mode](https://docs.julialang.org/en/v1/stdlib/REPL/#The-different-prompt-modes-1)
# of the Julia REPL:
# ```julia
# ;zgrep ^$first_seq $pfam_file
# ```
#
# Insert columns are deleted by MIToS, therefore the first MSA column has the
# label `"69"` (column  number in the input file) instead of `"1"`:

columnnames(msa)[1]

# If we want to keep track of residue numbers in each sequence, we need to read
# the file using the keyword arguments `generatemapping` and `useidcoordinates`:

const mapped_msa = read(pfam_file,
						Stockholm,
						generatemapping=true,
						useidcoordinates=true)

#-

getsequencemapping(mapped_msa, "A0A094YQ53_BACAO/182-489")

# The mapping is stored in the MSA annotations and it is automatically updated
# when columns are deleted.
#
# The `getannot...` functions return dictionaries with the MSA annotations

getannotresidue(mapped_msa)

# Pfam has "AS" (active site) and "pAS" (predicted active site) annotations for
# some families. We use `filter` to get a dictionary with the experimentally
# annotated active site residues (indicated by `*`):

active_sites = filter(pair -> pair[1][2] == "AS", getannotresidue(mapped_msa))

# #### Exercise I
#
# Which are the UniProt residue numbers of residues at the active site of
# AMPL_BOVIN/197-508?
#
# **Hint:** You can use the functions `getsequencemapping` and `zip`.
#

## ...your solution...

#src for (pos, ann) in zip(getsequencemapping(mapped_msa, "AMPL_BOVIN/197-508"),
#src 			 		   active_sites[("AMPL_BOVIN/197-508", "AS")])
#src 	if ann == '*'
#src 		println(pos)
#src 	end
#src end

# The answer is [294 and 368](https://www.uniprot.org/uniprot/P00727).

# ## MSA statistics and plotting
#
# [Plots.jl](http://docs.juliaplots.org/latest/) is a nice library because
# it gives a single API for multiple plotting backends (packages):
# - [GR](https://github.com/jheinen/GR.jl) is fast
# - [PlotlyJS](https://github.com/sglyon/PlotlyJS.jl) is interactive
# - [UnicodePlots](https://github.com/Evizero/UnicodePlots.jl) doesn't need a graphics environment
# - [PGFPlots](https://github.com/JuliaTeX/PGFPlots.jl) has LaTeX quality
# - [PyPlot](https://github.com/JuliaPy/PyPlot.jl) for pythonists

using Plots
gr() # choose the plotting backend

# Plots.jl allows to define [plot recipes](docs.juliaplots.org/latest/recipes/)
# and MIToS uses them to define plots for common biological objects:

plot(mapped_msa)

# Plots are also really useful to visualize MSA statistics:

plotlyjs() # change to an interactive backend

#-

gap_percentage = 100.0 .* gapfraction(mapped_msa, 1) # i.e. reduce dimension 1

#-

x = names(gap_percentage, 2)

#-

y = vec(gap_percentage)

#-

plot(x, y, legend=false)
xlabel!("MSA column (input file)")
ylabel!("Gap percentage")

# Something you can easily calculate with MIToS is the percentage identity
# between all the sequences in a Multiple Sequence Alignment:

percent_id = percentidentity(mapped_msa)

# #### Exercise 1
#
# What is the sequence that has the highest identity percentage with
# AMPL_BOVIN/197-508?
#
# **Hint:** Use [NamedArrays](https://github.com/davidavdav/NamedArrays.jl)
# `Not` indexing and `names` and Julia's `findmax`:

using NamedArrays

## ...your solution...

#src identities = percent_id["AMPL_BOVIN/197-508", Not("AMPL_BOVIN/197-508")]
#src value, index = findmax(identities)
#src names(identities, 1)[index]

max_stars_repo_name: diegozea/JuliaForBioinformatics

max_issues_repo_name: diegozea/JuliaForBioinformatics

max_forks_repo_name: diegozea/JuliaForBioinformatics


text:
function featureextraction(
    data::ModelSelection.ModelSelectionData;
    fe_sqr::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_log::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_inv::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_lag::Union{Nothing, Array}=nothing,
    interaction::Union{Nothing, Array, Dict}=nothing,
    removemissings::Bool=REMOVEMISSINGS_DEFAULT
    )

    return featureextraction!(
        ModelSelection.copy_data(data),
        fe_sqr=fe_sqr,
        fe_log=fe_log,
        fe_inv=fe_inv,
        fe_lag=fe_lag,
        interaction=interaction,
        removemissings=removemissings
    )
end

function featureextraction!(
    data::ModelSelection.ModelSelectionData;
    fe_sqr::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_log::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_inv::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_lag::Union{Nothing, Array}=nothing,
    interaction::Union{Nothing, Array, Dict}=nothing,
    removemissings::Bool=REMOVEMISSINGS_DEFAULT
    )

    data = execute!(
        data,
        fe_sqr=fe_sqr,
        fe_log=fe_log,
        fe_inv=fe_inv,
        fe_lag=fe_lag,
        interaction=interaction,
        removemissings=removemissings
    )

    data = addextras(data, fe_sqr, fe_log, fe_inv, fe_lag, interaction, removemissings)

    return data
end

function execute!(
    data::ModelSelection.ModelSelectionData;
    fe_sqr::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_log::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_inv::Union{Nothing, String, Symbol, Array{String}, Array{Symbol}}=nothing,
    fe_lag::Union{Nothing, Array}=nothing,
    interaction::Union{Nothing, Array, Dict}=nothing,
    removemissings::Bool=REMOVEMISSINGS_DEFAULT
    )

    if data.intercept
        ModelSelection.remove_intercept!(data)
    end

    if fe_sqr != nothing
        data = data_add_fe_sqr(data, parse_fe_variables(fe_sqr, data.expvars))
    end

    if fe_log != nothing
        data = data_add_fe_log(data, parse_fe_variables(fe_log, data.expvars))
    end

    if fe_inv != nothing
        data = data_add_fe_inv(data, parse_fe_variables(fe_inv, data.expvars))
    end

    if fe_lag != nothing
        data = data_add_fe_lag(data, parse_fe_variables(fe_lag, data.expvars, depvar=data.depvar, is_pair=true))
    end

    if interaction != nothing
        data = data_add_interaction(data, parse_fe_variables(interaction, data.expvars))
    end

    if data.intercept
        ModelSelection.add_intercept!(data)
    end

    if removemissings
        data = ModelSelection.filter_data_by_empty_values(data)
    end

    data = ModelSelection.convert_data(data)

    return data
end


text:
# kink linear: f(x)=max(x,gamma*x) where gamma=0.25
gamma=0.25
kinklin(x)=max(x,gamma*x)

Fkinklin(x::Array{Float64,2})=(max(x,gamma*x),nothing)
Fkinklin_inplace(value,auxvalue,x::Array{Float64,2})=copy!(value,max(x,gamma*x))

Dkinklin(derivativeIDX,f_c,faux_c,grad_c,grad_n,x)=axpy!(1.0,grad_c.*(gamma+(1-gamma)*(x.>0)),grad_n)

if PROC=="GPU"
    Fkinklin_inplace(value,aux,x::CudaArray)=kinklin!(x,value)
    function Dkinklin(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray)
        tmp=CudaArray(Float64,size(grad_c)); fill!(tmp,0.0)
        A_emult_Bg0!(grad_c,x,tmp)
        scale!((1-gamma),tmp)
        axpy!(gamma,grad_c,tmp)
        axpy!(1.0,tmp,grad_n)
        free(tmp)
    end
end

Derivative[Fkinklin]=Dkinklin
Inplace[Fkinklin]=Fkinklin_inplace

kinklin(A::ADnode)=ADnode(Fkinklin,A)
export kinklin

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff