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

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

# ╔═╡ 5e688928-e939-11ea-0e16-fbc80af390ab
using LinearAlgebra

# ╔═╡ a50b5f48-e8d5-11ea-1f05-a3741b5d15ba
html"<button onclick=present()>Present</button>"

# ╔═╡ 8a6fed4c-e94b-11ea-1113-d56f56fb293b
br = HTML("
")

# ╔═╡ dc53f316-e8c8-11ea-150f-1374dbce114a
md"""# Welcome to 18.S191 -- Fall 2020!

### Introduction to Computational Thinking for Real-World Problems"""

# ╔═╡ c3f43d66-e94b-11ea-02bd-23cfeb878ff1
br

# ╔═╡ c6c77738-e94b-11ea-22f5-1dce3dbcc3ca
md"### <https://github.com/mitmath/18S191>"

# ╔═╡ cf80793a-e94b-11ea-0120-f7913ae06f22
br

# ╔═╡ d1638d96-e94b-11ea-2ff4-910e399f864d
md"##### Alan Edelman, David P. Sanders, Grant Sanderson, James Schloss"

# ╔═╡ 0117246a-e94c-11ea-1a76-c981ce8e725d
md"##### & Philip the Corgi"

# ╔═╡ 27060098-e8c9-11ea-2fe0-03b39b1ddc32
md"""# Class outline

### Data and computation

- Module 1: Analyzing images

- Module 2: Particles and ray tracing

- Module 3: Epidemic spread

- Module 4: Climate change
"""

# ╔═╡ 4fc58814-e94b-11ea-339b-cb714a63f9b6
md"## Tools

- Julia programming language: <http://www.julialang.org/learning>

- Pluto notebook environment
"

# ╔═╡ f067d3b8-e8c8-11ea-20cb-474709ffa99a
md"""# Module 1: Images"""

# ╔═╡ 37c1d012-ebc9-11ea-2dfe-8b86bb78f283
4 + 4

# ╔═╡ a0a97214-e8d2-11ea-0f46-0bfaf016ab6d
md"""## Data takes many forms

- Time series:
  - Number of infections per day
  - Stock price each minute
  - A piece for violin broadcast over the radio
$(HTML("
"))

- Video:
  - The view from a window of a self-driving car
  - A hurricane monitoring station
$(HTML("
"))

- Images:
  - Diseased versus healthy tissue in a scan
  - Deep space via the Hubble telescope
  - Can your social media account recognise your friends?
"""

# ╔═╡ 1697a756-e93d-11ea-0b6e-c9c78d527993
md"## Capture your own image!"

# ╔═╡ af28faca-ebb7-11ea-130d-0f94bf9bd836

# ╔═╡ ee1d1596-e94a-11ea-0fb4-cd05f62471d3
md"##"

# ╔═╡ 8ab9a978-e8c9-11ea-2476-f1ef4ba1b619
md"""## What is an image?"""

# ╔═╡ 38c54bfc-e8cb-11ea-3d52-0f02452f8ba1
md"Albrecht Dürer:"

# ╔═╡ 983f8270-e8c9-11ea-29d2-adeccb5a7ffc
# md"# begin
# 	using Images

# 	download("https://i.stack.imgur.com/QQL8X.jpg", "durer.jpg")

# 	load("durer.jpg")
# end

md"![](https://i.stack.imgur.com/QQL8X.jpg)"

# ╔═╡ 2fcaef88-e8ca-11ea-23f7-29c48580f43c
md"""##

An image is:

- A 2D representation of a 3D world

- An approximation

"""

# ╔═╡ 7636c4b0-e8d1-11ea-2051-757a850a9d30
begin
	image_text =
	md"""
	## What *is* an image, though?

	- A grid of coloured squares called **pixels**

	- A colour for each pair $(i, j)$ of indices

	- A **discretization**

	"""

	image_text
end

# ╔═╡ bca22176-e8ca-11ea-2004-ebeb103116b5
md"""
## How can we store an image in the computer?

- Is it a 1D array (`Vector`)?

- A 2D array (`Matrix`)?

- A 3D array (`tensor`)?
"""

# ╔═╡ 0ad91f1e-e8d2-11ea-2c18-93f66c906a8b
md"""## If in doubt: Ask Julia!

- Let's use the `Images.jl` package to load an image and see what we get
"""

# ╔═╡ de373816-ec79-11ea-2772-ebdca52246ac
begin
	import Pkg
	Pkg.activate(mktempdir())
end

# ╔═╡ 552129ae-ebca-11ea-1fa1-3f9fa00a2601
begin
	Pkg.add(["Images", "ImageIO", "ImageMagick"])
	using Images
end

# ╔═╡ fbe11200-e938-11ea-12e9-6125c1b56b25
begin
	Pkg.add("PlutoUI")
	using PlutoUI
end

# ╔═╡ 54c1ba3c-e8d2-11ea-3564-bdaca8563738
# defines a variable called `url`
# whose value is a string (written inside `"`):

url = "https://i.imgur.com/VGPeJ6s.jpg"

# ╔═╡ 6e0fefb6-e8d4-11ea-1f9b-e7a3db40df39
philip_file = download(url, "philip.jpg")  # download to a local file

# ╔═╡ 9c359212-ec79-11ea-2d7e-0124dad5f127
philip = load(philip_file)

# ╔═╡ 7703b032-ebca-11ea-3074-0b80a077078e
philip

# ╔═╡ 7eff3522-ebca-11ea-1a65-59e66a4e72ab
typeof(philip)

# ╔═╡ c9cd6c04-ebca-11ea-0990-5fa19ff7ed97
RGBX(0.9, 0.1, 0.1)

# ╔═╡ 0d873d9c-e93b-11ea-2425-1bd79677fb97
md"##"

# ╔═╡ 6b09354a-ebb9-11ea-2d5a-3b75c5ae7aa9
RGB(0.8,1,1)

# ╔═╡ 2d6c434e-e93b-11ea-2678-3b9db4975089
md"##"

# ╔═╡ 2b14e93e-e93b-11ea-25f1-5f565f80e778
typeof(philip)

# ╔═╡ 0bdc6058-e8d5-11ea-1889-3f706cea7a1f
md"""##

- According to Julia / Pluto, the variable `philip` *is* an image

- Julia always returns output

- The output can be displayed in a "rich" way

$(HTML("
"))

- Arthur C. Clarke:

> Any sufficiently advanced technology is indistinguishable from magic.
"""

# ╔═╡ e61db924-ebca-11ea-2f79-f9f1c121b7f5
size(philip)

# ╔═╡ ef60fcc4-ebca-11ea-3f69-155afffe8ea8
philip

# ╔═╡ fac550ec-ebca-11ea-337a-dbc16848c617
philip[1:1000, 1:400]

# ╔═╡ 42aa8cfe-e8d5-11ea-3cb9-c365b98e7a8c
md"
## How big is Philip?

- He's pretty big:
"

# ╔═╡ 4eea5710-e8d5-11ea-3978-af66ee2a137e
size(philip)

# ╔═╡ 57b3a0c2-e8d5-11ea-15aa-8da4549f849b
md"- Which number is which?"

# ╔═╡ 03a7c0fc-ebba-11ea-1c71-79d750c97b16
philip

# ╔═╡ e6fd68fa-e8d8-11ea-3dc4-274caceda222
md"# So, what *is* an image?"

# ╔═╡ 63a1d282-e8d5-11ea-0bba-b9cdd32a218b
typeof(philip)

# ╔═╡ fc5e1af0-e8d8-11ea-1077-07216ff96d29
md"""
- It's an `Array`

- The `2` means that it has **2 dimensions** (a **matrix**)

$(HTML("
"))

- `RGBX{Normed{UInt8,8}}` is the type of object stored in the array

- A Julia object representing a colour

- RGB = Red, Green, Blue
"""

# ╔═╡ c79dd836-e8e8-11ea-029d-57be9899979a
md"## Getting pieces of an image"

# ╔═╡ ae260168-e932-11ea-38fd-4f2c6f43e21c
begin
	(h, w) = size(philip)
	head = philip[(h ÷ 2):h, (w ÷ 10): (9w ÷ 10)]
	# `÷` is typed as \div <TAB>  -- integer division
end

# ╔═╡ 47d1bc04-ebcb-11ea-3643-d1ba8dea57c8
size(head)

# ╔═╡ 72400458-ebcb-11ea-26b6-678ae1de8e23
size(philip)

# ╔═╡ f57ea7c2-e932-11ea-0d52-4112187bcb38
md"## Manipulating matrices

- An image is just a matrix, so we can manipulate *matrices* to manipulate the *image*
"

# ╔═╡ 740ed2e2-e933-11ea-236c-f3c3f09d0f8b
[head head head]

# ╔═╡ 6128a5ba-e93b-11ea-03f5-f170c7b90b25
md"##"

# ╔═╡ 78eafe4e-e933-11ea-3539-c13feb894ef6
[
 head                   reverse(head, dims=2)
 reverse(head, dims=1)  reverse(reverse(head, dims=1), dims=2)
]

# ╔═╡ bf3f9050-e933-11ea-0df7-e5dcff6bb3ee
md"## Manipulating an image

- How can we get inside the image and change it?

- There are two possibilities:

  - **Modify** (**mutate**) numbers inside the array -- useful to change a small piece

  - Create a new **copy** of the array -- useful to alter everything together
"

# ╔═╡ 212e1f12-e934-11ea-2f35-51c7a6c8dff1
md"## Painting a piece of an image

- Let's paint a corner red

- We'll copy the image first so we don't destroy the original
"

# ╔═╡ 117a98c0-e936-11ea-3aac-8f66337cea68
new_phil = copy(head)

# ╔═╡ 8004d076-e93b-11ea-29cc-a1bfcc75e87f
md"##"

# ╔═╡ 3ac63296-e936-11ea-2144-f94bdbd60eaf
red = RGB(1, 0, 0)

# ╔═╡ 3e3f841a-e936-11ea-0a81-1b95fe0faa83
for i in 1:100
	for j in 1:300
		new_phil[i, j] = red
	end
end

# ╔═╡ 5978db50-e936-11ea-3145-059a51be2281
md"Note that `for` loops *do not return anything* (or, rather, they return `nothing`)"

# ╔═╡ 21638b14-ebcc-11ea-1761-bbd2f4306a96
new_phil

# ╔═╡ 70cb0e36-e936-11ea-3ade-49fde77cb696
md"""## Element-wise operations: "Broadcasting"

- Julia provides powerful technology for operating element by element: **broadcasting**

- Adding "`.`" applies an operation element by element
"""

# ╔═╡ b3ea975e-e936-11ea-067d-81339575a3cb
begin
	new_phil2 = copy(new_phil)
	new_phil2[100:200, 1:100] .= RGB(0, 1, 0)
	new_phil2
end

# ╔═╡ 918a0762-e93b-11ea-1115-71dbfdb03f27
md"##"

# ╔═╡ daabe66c-e937-11ea-3bc3-d77f2bce406c
new_phil2

# ╔═╡ 095ced62-e938-11ea-1169-939dc7136fd0
md"## Modifying the whole image at once

- We can use the same trick to modify the whole image at once

- Let's **redify** the image

- We define a **function** that turns a colour into just its red component
"

# ╔═╡ 31f3605a-e938-11ea-3a6d-29a185bbee31
function redify(c)
	return RGB(c.r, 0, 0)
end

# ╔═╡ 2744a556-e94f-11ea-2434-d53c24e59285
begin
	color = RGB(0.9, 0.7, 0.2)

	[color, redify(color)]
end

# ╔═╡ 98412a36-e93b-11ea-1954-f1c105c6ed4a
md"##"

# ╔═╡ 3c32efde-e938-11ea-1ae4-5d88290f5311
redify.(philip)

# ╔═╡ 4b26e4e6-e938-11ea-2635-6d4fc15e13b7
md"## Transforming an image

- The main goal of this week will be to transfrom images in more interesting ways

- First let's **decimate** poor Phil
"

# ╔═╡ c12e0928-e93b-11ea-0922-2b590a99ee89
md"##"

# ╔═╡ ff5dc538-e938-11ea-058f-693d6b016640
md"## Experiments come alive with interaction

- We start to get a feel for things when we can **experiment**!
"

# ╔═╡ fa24f4a8-e93b-11ea-06bd-25c9672166d6
md"##"

# ╔═╡ 15ce202e-e939-11ea-2387-93be0ec4cf1f
@bind repeat_count Slider(1:20, show_value=true)

# ╔═╡ bf2167a4-e93d-11ea-03b2-cdd24b459ba9
md"## Summary

- Images are readily-accessible data about the world

- We want to process them to extract information

- Relatively simple mathematical operations can transform images in useful ways
"

# ╔═╡ 58184d88-e939-11ea-2fc8-73b3476ebe92
expand(image, ratio=5) = kron(image, ones(ratio, ratio))

# ╔═╡ 2dd09f16-e93a-11ea-2cdc-13f558e3391d
extract_red(c) = c.r

# ╔═╡ df1b7996-e93b-11ea-1a3a-81b4ec520679
decimate(image, ratio=5) = image[1:ratio:end, 1:ratio:end]

# ╔═╡ 41fa85c0-e939-11ea-1ad8-79805a2083bb
poor_phil = decimate(head, 5)

# ╔═╡ cd5721d0-ede6-11ea-0918-1992c69bccc6
repeat(poor_phil, repeat_count, repeat_count)

# ╔═╡ 9709e500-ee94-11ea-0297-d3f01ebcdd14
md"## Experimenting -- Daniel Fryer"

# ╔═╡ 88caeec0-ee95-11ea-22dd-1f5a5a41d3cd
philip

# ╔═╡ 10ddf1a0-ee95-11ea-0b63-3bf6ef1386f1
function swapify(c)
	return RGB(c.r*c.b, c.g*c.b, c.b*c.b)
end

# ╔═╡ b2bc0b20-ee94-11ea-0e97-61df7e3d45f0
swapify.(philip)

# ╔═╡ f9a4d3b0-ee98-11ea-3b58-3331ca84ec7a

# ╔═╡ f678c110-ee98-11ea-0ca4-cf4d5b21ad4b
new_philip = copy(philip)

# ╔═╡ c963f680-ee99-11ea-0b05-97c1e6dbf39a
function npi(n)
  trunc(Int,rem(10.0^(n-1)*pi,1)*10)
end

# ╔═╡ b8daeea0-ec79-11ea-34b5-3f13e8a56a42
md"# Appendix"

# ╔═╡ bf1bb2c8-ec79-11ea-0671-3ffb34828f3c
md"## Package environment"

# ╔═╡ 69e3aa82-e93c-11ea-23fe-c1103d989cba
md"## Camera input"

# ╔═╡ 739c3bb6-e93c-11ea-127b-efb6a8ab9379
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

# ╔═╡ 9529bc40-e93c-11ea-2587-3186e0978476
@bind raw_camera_data camera_input(;max_size=2000)

# ╔═╡ 832ebd1a-e93c-11ea-1d18-d784f3184ebe

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

# ╔═╡ 9a843af8-e93c-11ea-311b-1bc6d5b58492
grant = decimate(process_raw_camera_data(raw_camera_data), 2)

# ╔═╡ 6aa73286-ede7-11ea-232b-63e052222ecd
[
	grant             grant[:,end:-1:1]
	grant[end:-1:1,:] grant[end:-1:1,end:-1:1]
]

# ╔═╡ Cell order:
# ╠═a50b5f48-e8d5-11ea-1f05-a3741b5d15ba
# ╠═8a6fed4c-e94b-11ea-1113-d56f56fb293b
# ╠═dc53f316-e8c8-11ea-150f-1374dbce114a
# ╟─c3f43d66-e94b-11ea-02bd-23cfeb878ff1
# ╟─c6c77738-e94b-11ea-22f5-1dce3dbcc3ca
# ╟─cf80793a-e94b-11ea-0120-f7913ae06f22
# ╟─d1638d96-e94b-11ea-2ff4-910e399f864d
# ╟─0117246a-e94c-11ea-1a76-c981ce8e725d
# ╟─27060098-e8c9-11ea-2fe0-03b39b1ddc32
# ╟─4fc58814-e94b-11ea-339b-cb714a63f9b6
# ╟─f067d3b8-e8c8-11ea-20cb-474709ffa99a
# ╠═37c1d012-ebc9-11ea-2dfe-8b86bb78f283
# ╟─a0a97214-e8d2-11ea-0f46-0bfaf016ab6d
# ╟─1697a756-e93d-11ea-0b6e-c9c78d527993
# ╟─af28faca-ebb7-11ea-130d-0f94bf9bd836
# ╠═9529bc40-e93c-11ea-2587-3186e0978476
# ╟─ee1d1596-e94a-11ea-0fb4-cd05f62471d3
# ╠═6aa73286-ede7-11ea-232b-63e052222ecd
# ╠═9a843af8-e93c-11ea-311b-1bc6d5b58492
# ╟─8ab9a978-e8c9-11ea-2476-f1ef4ba1b619
# ╟─38c54bfc-e8cb-11ea-3d52-0f02452f8ba1
# ╟─983f8270-e8c9-11ea-29d2-adeccb5a7ffc
# ╟─2fcaef88-e8ca-11ea-23f7-29c48580f43c
# ╟─7636c4b0-e8d1-11ea-2051-757a850a9d30
# ╟─bca22176-e8ca-11ea-2004-ebeb103116b5
# ╟─0ad91f1e-e8d2-11ea-2c18-93f66c906a8b
# ╠═de373816-ec79-11ea-2772-ebdca52246ac
# ╠═552129ae-ebca-11ea-1fa1-3f9fa00a2601
# ╠═54c1ba3c-e8d2-11ea-3564-bdaca8563738
# ╠═6e0fefb6-e8d4-11ea-1f9b-e7a3db40df39
# ╠═9c359212-ec79-11ea-2d7e-0124dad5f127
# ╠═7703b032-ebca-11ea-3074-0b80a077078e
# ╠═7eff3522-ebca-11ea-1a65-59e66a4e72ab
# ╠═c9cd6c04-ebca-11ea-0990-5fa19ff7ed97
# ╟─0d873d9c-e93b-11ea-2425-1bd79677fb97
# ╠═6b09354a-ebb9-11ea-2d5a-3b75c5ae7aa9
# ╟─2d6c434e-e93b-11ea-2678-3b9db4975089
# ╠═2b14e93e-e93b-11ea-25f1-5f565f80e778
# ╟─0bdc6058-e8d5-11ea-1889-3f706cea7a1f
# ╠═e61db924-ebca-11ea-2f79-f9f1c121b7f5
# ╠═ef60fcc4-ebca-11ea-3f69-155afffe8ea8
# ╠═fac550ec-ebca-11ea-337a-dbc16848c617
# ╟─42aa8cfe-e8d5-11ea-3cb9-c365b98e7a8c
# ╠═4eea5710-e8d5-11ea-3978-af66ee2a137e
# ╟─57b3a0c2-e8d5-11ea-15aa-8da4549f849b
# ╠═03a7c0fc-ebba-11ea-1c71-79d750c97b16
# ╟─e6fd68fa-e8d8-11ea-3dc4-274caceda222
# ╠═63a1d282-e8d5-11ea-0bba-b9cdd32a218b
# ╟─fc5e1af0-e8d8-11ea-1077-07216ff96d29
# ╟─c79dd836-e8e8-11ea-029d-57be9899979a
# ╠═ae260168-e932-11ea-38fd-4f2c6f43e21c
# ╠═47d1bc04-ebcb-11ea-3643-d1ba8dea57c8
# ╠═72400458-ebcb-11ea-26b6-678ae1de8e23
# ╟─f57ea7c2-e932-11ea-0d52-4112187bcb38
# ╠═740ed2e2-e933-11ea-236c-f3c3f09d0f8b
# ╟─6128a5ba-e93b-11ea-03f5-f170c7b90b25
# ╠═78eafe4e-e933-11ea-3539-c13feb894ef6
# ╟─bf3f9050-e933-11ea-0df7-e5dcff6bb3ee
# ╟─212e1f12-e934-11ea-2f35-51c7a6c8dff1
# ╠═117a98c0-e936-11ea-3aac-8f66337cea68
# ╟─8004d076-e93b-11ea-29cc-a1bfcc75e87f
# ╠═3ac63296-e936-11ea-2144-f94bdbd60eaf
# ╠═3e3f841a-e936-11ea-0a81-1b95fe0faa83
# ╟─5978db50-e936-11ea-3145-059a51be2281
# ╠═21638b14-ebcc-11ea-1761-bbd2f4306a96
# ╟─70cb0e36-e936-11ea-3ade-49fde77cb696
# ╠═b3ea975e-e936-11ea-067d-81339575a3cb
# ╟─918a0762-e93b-11ea-1115-71dbfdb03f27
# ╠═daabe66c-e937-11ea-3bc3-d77f2bce406c
# ╠═095ced62-e938-11ea-1169-939dc7136fd0
# ╠═31f3605a-e938-11ea-3a6d-29a185bbee31
# ╠═2744a556-e94f-11ea-2434-d53c24e59285
# ╟─98412a36-e93b-11ea-1954-f1c105c6ed4a
# ╠═3c32efde-e938-11ea-1ae4-5d88290f5311
# ╟─4b26e4e6-e938-11ea-2635-6d4fc15e13b7
# ╠═41fa85c0-e939-11ea-1ad8-79805a2083bb
# ╟─c12e0928-e93b-11ea-0922-2b590a99ee89
# ╟─ff5dc538-e938-11ea-058f-693d6b016640
# ╠═fbe11200-e938-11ea-12e9-6125c1b56b25
# ╟─fa24f4a8-e93b-11ea-06bd-25c9672166d6
# ╠═15ce202e-e939-11ea-2387-93be0ec4cf1f
# ╠═cd5721d0-ede6-11ea-0918-1992c69bccc6
# ╟─bf2167a4-e93d-11ea-03b2-cdd24b459ba9
# ╟─5e688928-e939-11ea-0e16-fbc80af390ab
# ╠═58184d88-e939-11ea-2fc8-73b3476ebe92
# ╟─2dd09f16-e93a-11ea-2cdc-13f558e3391d
# ╟─df1b7996-e93b-11ea-1a3a-81b4ec520679
# ╟─9709e500-ee94-11ea-0297-d3f01ebcdd14
# ╠═88caeec0-ee95-11ea-22dd-1f5a5a41d3cd
# ╠═b2bc0b20-ee94-11ea-0e97-61df7e3d45f0
# ╠═10ddf1a0-ee95-11ea-0b63-3bf6ef1386f1
# ╠═f9a4d3b0-ee98-11ea-3b58-3331ca84ec7a
# ╠═f678c110-ee98-11ea-0ca4-cf4d5b21ad4b
# ╠═c963f680-ee99-11ea-0b05-97c1e6dbf39a
# ╟─b8daeea0-ec79-11ea-34b5-3f13e8a56a42
# ╟─bf1bb2c8-ec79-11ea-0671-3ffb34828f3c
# ╟─69e3aa82-e93c-11ea-23fe-c1103d989cba
# ╟─739c3bb6-e93c-11ea-127b-efb6a8ab9379
# ╟─832ebd1a-e93c-11ea-1d18-d784f3184ebe


text:
__precompile__(false)
"""
Submodule Julia Client Echo
Intended for running in its own worker process.
HTTP and WebSockets need to be loaded in the calling context.
LOAD_PATH must include the directory logutils_ws.
See comment at the end of this file for debugging code.
"""
module ws_jce
using ..HTTP
using ..WebSockets
# We want to log to a separate file, and so use our own
# instance of logutils_ws in this process
import logutils_ws: logto, clog, zlog, zflush, clog_notime
const SRCPATH = Base.source_dir() == nothing ? Pkg.dir("WebSockets", "benchmark") : Base.source_dir()
const LOGFILE = "ws_jce.log"

const PORT = 8000
const SERVER = "ws://127.0.0.1:$(PORT)"
const CLOSEAFTER = Base.Dates.Second(30)

"""
Opens a client, echoes with an optional delay, an integer in milliseconds.
Stores time records for received messages and before sending messages.
Specify the delay in milliseconds by sending a message on the websocket:
    send(ws_jce, "delay|15")
Echoes any message except "exit" and "delay".

Delays to reading, in the websocket use situation, would be caused by usefully spent
calculation time between reads. However, they may be interpreted by the underlying protocol
as transmission problems and cause large slowdowns. Hence the interest in testing
with delays. A countermeasure for optimizing speed might be to run a websocket
reading function in a parallel, not asyncronous process, putting messages on an internal queue.

At exit or after CLOSEAFTER, this function sends one message containing two vectors of
timestamps [ns].
"""
function echowithdelay_jce()
    # This will be run in a worker process. Even so, individual console log
    # output will be redirected to process 1 and prefixed
    # with a "From worker 2".
    # It will also be interspersed with process 1 output,
    # sometimes before a line is finished printing.
    # We use :green to distinguish more easily.
    id = "echowithdelay_jce"
    f = open(joinpath(SRCPATH, "logs", LOGFILE), "w")
    try
        logto(f)
        clog(id, :green, "Open client on ", SERVER, "\nclient side handler ", _jce)
        zflush()
        WebSockets.open(_jce, SERVER)
        zlog(id, :green, " Websocket closed, control returned.")
    catch err
        clog(id, :red, err)
        clog_notime.(stacktrace(catch_backtrace())[1:4])
        zflush()
    finally
        clog(id, :green, " Closing log ", LOGFILE)
        zflush()
        logto(Base.DevNullStream())
        close(f)
    end
end
"
Handler for client websocket, defined by echowithdelay_jce
"
function _jce(ws)
    id = "_jce"
    clog(id, :green, ws)
    zflush()
    receivetimes = Vector{Int64}()
    replytime = Int64(0)
    replytimes = Vector{Int64}()
    msg = Vector{UInt8}()
    delay = 0 # Integer milliseconds
    t1 = now() + CLOSEAFTER
    while isopen(ws) && now() < t1
        # read, record receive time
        msg = read(ws)
        ti = time_ns()
        push!(receivetimes, Int64(ti < typemax(Int64) ? ti : 0 ))
        # break out when receiving 'exit'
        length(msg) == 4 && msg == Vector{UInt8}("exit") && break
        # react to delay instruction
        if length(msg) < 16 && msg[1:6] == Vector{UInt8}("delay=")
            delay = Meta.parse(Int, String(msg[7:end]))
            clog(id, :green, " Changing delay to ", delay, " ms")
            zflush()
        end
        sleep(delay / 1000)
        # record send time, echo
        replytime = time_ns()
        write(ws, msg)
        # clean record of instruction message
        if length(msg) > 16 && msg[1:6] != Vector{UInt8}("delay=")
            push!(replytimes, Int64(replytime < typemax(Int64) ? replytime : 0 ))
        elseif msg[1:6] != Vector{UInt8}("delay=")
            push!(replytimes,  Int64(replytime < typemax(Int64) ? replytime : 0 ))
        end
    end
    if length(receivetimes) > 1
        # Don't include receive time for the "exit" message. Reply time was not recorded
        pop!(receivetimes)
        if length(receivetimes) > 1
            # Send one message with serialized receive and reply times.
            buf = IOBuffer()
            serialize(buf, (receivetimes, replytimes))
            if isopen(ws)
                write(ws, take!(buf))
                zlog(id, :green, " Sent serialized receive and reply times.")
                zflush()
            end
        end
    end
    clog(id, :green, " Exit, close websocket.")
    zflush()
    # Exiting this function starts a closing handshake
    nothing
end
end # module

#=
For debugging in a separate terminal:

using WebSockets
const SRCPATH = Base.source_dir() == nothing ? Pkg.dir("WebSockets", "benchmark") :Base.source_dir()
const LOGGINGPATH = realpath(joinpath(SRCPATH, "../logutils/"))
# for finding local modules
SRCPATH ∉ LOAD_PATH && push!(LOAD_PATH, SRCPATH)
LOGGINGPATH ∉ LOAD_PATH && push!(LOAD_PATH, LOGGINGPATH)

import ws_jce.echowithdelay_jce
ws_jce.echowithdelay_jce()
=#

ext: jl


text:
######################################################################
# c_interface.jl
######################################################################

# internal pointer for ccall stuff
mutable struct rtlsdr_dev
end

mutable struct RTLSDRError <: Exception
    message::String
end

# returns a pointer to
function rtlsdr_open(index::Int64=0)
    rd = Array{Ptr{rtlsdr_dev},1}(undef,1)
    ret = ccall( (:rtlsdr_open, "librtlsdr"),
                 Cint,
                 (Ptr{Ptr{rtlsdr_dev}}, UInt32),
                 rd,
                 index
               )

    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error opening device.")); end

    rf = rd[1]

    # resetting buffers is critical to reading bytes
    rtlsdr_reset_buffer(rf)

    return rf
end

function rtlsdr_reset_buffer(rf::Ref{rtlsdr_dev})
    ret = ccall( (:rtlsdr_reset_buffer, "librtlsdr"),
                 Cint,
                 (Ref{rtlsdr_dev},),
                 rf
               )
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error resetting buffer (error code $ret).")); end
end

function rtlsdr_close(rf::Ref{rtlsdr_dev})
    ret = ccall( (:rtlsdr_close, "librtlsdr"), Cint, (Ref{rtlsdr_dev},), rf)
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error closing device (error code $ret).")); end
end

# center frequency
function rtlsdr_set_center_freq(rf::Ref{rtlsdr_dev}, freq)
    ret = ccall( (:rtlsdr_set_center_freq, "librtlsdr"), Cint, (Ref{rtlsdr_dev}, UInt32), rf, Cint(freq))
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error setting center frequency (error code $ret).")); end
end

function rtlsdr_get_center_freq(rf::Ref{rtlsdr_dev})
    ret = ccall( (:rtlsdr_get_center_freq, "librtlsdr"), UInt32, (Ref{rtlsdr_dev},), rf)
    if ret < 0
        throw(RTLSDRError("RTLSDR.jl reports: Error getting center frequency (error code $ret)."));
    end
    return ret
end

# sample rate
function rtlsdr_get_sample_rate(rf::Ref{rtlsdr_dev})
    ret = ccall( (:rtlsdr_get_sample_rate, "librtlsdr"),
                 UInt32,
                 (Ref{rtlsdr_dev},),
                 rf
               )
    if ret < 0
        throw(RTLSDRError("RTLSDR.jl reports: Error getting sample rate (error code $ret)."));
    end
    return ret
end

function rtlsdr_set_sample_rate(rf::Ref{rtlsdr_dev}, sample_rate)
    ret = ccall( (:rtlsdr_set_sample_rate, "librtlsdr"),
                 Cint,
                 (Ref{rtlsdr_dev},UInt32),
                 rf,
                 sample_rate
               )
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error satting sample rate (error code $ret).")); end
end

# gain
function rtlsdr_set_tuner_gain(rf::Ref{rtlsdr_dev}, gain)
    ret = ccall( (:rtlsdr_set_tuner_gain, "librtlsdr"),
                 Cint,
                 (Ref{rtlsdr_dev},Cint),
                 rf,
                 gain
               )
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error setting tuner gain (error code $ret).")); end
end

function rtlsdr_get_tuner_gain(rf::Ref{rtlsdr_dev})
    gain = ccall( (:rtlsdr_get_tuner_gain, "librtlsdr"),
                 Cint,
                 (Ref{rtlsdr_dev},),
                 rf
               )
    return gain[1]
end

function rtlsdr_get_tuner_gains(rf::Ref{rtlsdr_dev})
    buf = Vector{Cint}(undef, 50) # max 50 gain settings. no RTL has this many.
    num_gains = ccall( (:rtlsdr_get_tuner_gains, "librtlsdr"),
                 Cint,
                 (Ref{rtlsdr_dev}, Ref{Cint}),
                 rf,
                 buf
                )
    if num_gains == 0
        throw(RTLSDRError("RTLSDR.jl reports: Error getting gains."))
    end
    avail_gains = Vector{Int32}(undef, num_gains)
    for i = 1:num_gains
        avail_gains[i] = buf[i]
    end
    return avail_gains
end

function rtlsdr_set_tuner_gain_mode(rf::Ref{rtlsdr_dev}, manual)
    ret = ccall( (:rtlsdr_set_tuner_gain_mode, "librtlsdr"),
                 Cint,
                 (Ref{rtlsdr_dev},Cint),
                 rf,
                 manual
               )
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error setting tuner gain mode (error code $ret).")); end
end

function rtlsdr_set_agc_mode(rf::Ref{rtlsdr_dev}, on)
    ret = ccall( (:rtlsdr_set_tuner_gain, "librtlsdr"),
                 Cint,
                 (Ref{rtlsdr_dev},Cint),
                 rf,
                 on
                )
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error setting AGC mode (error code $ret).")); end
end

# reading
function read_bytes(rf::Ref{rtlsdr_dev}, num_bytes)
    buf = Vector{Cuchar}(undef,num_bytes)
    bytes_read = Ref{Cint}(0)
    ret = ccall( (:rtlsdr_read_sync, "librtlsdr"),
                 Int32,
                 (Ref{rtlsdr_dev}, Ref{Cuchar}, Cint, Ref{Cint}),
                 rf,
                 buf,
                 num_bytes,
                 bytes_read
                )
    if ret != 0; throw(RTLSDRError("RTLSDR.jl reports: Error while reading from device (error code $ret).")); end

    return buf
end


text:
# Autogenerated wrapper script for AFL_jll for x86_64-linux-gnu-cxx03
export afl_gcc_exe, afl_qemu_trace_aarch64_exe, afl_qemu_trace_arm_exe, afl_qemu_trace_i386_exe, afl_qemu_trace_x86_64_exe, afl_test_instr_exe

using Glib_jll
using Pixman_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `afl_gcc_exe`
const afl_gcc_exe_splitpath = ["bin", "afl-gcc"]

# This will be filled out by __init__() for all products, as it must be done at runtime
afl_gcc_exe_path = ""

# afl_gcc_exe-specific global declaration
function afl_gcc_exe(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(afl_gcc_exe_path)
    end
end

# Relative path to `afl_qemu_trace_aarch64_exe`
const afl_qemu_trace_aarch64_exe_splitpath = ["bin", "afl-qemu-trace-aarch64"]

# This will be filled out by __init__() for all products, as it must be done at runtime
afl_qemu_trace_aarch64_exe_path = ""

# afl_qemu_trace_aarch64_exe-specific global declaration
function afl_qemu_trace_aarch64_exe(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(afl_qemu_trace_aarch64_exe_path)
    end
end

# Relative path to `afl_qemu_trace_arm_exe`
const afl_qemu_trace_arm_exe_splitpath = ["bin", "afl-qemu-trace-arm"]

# This will be filled out by __init__() for all products, as it must be done at runtime
afl_qemu_trace_arm_exe_path = ""

# afl_qemu_trace_arm_exe-specific global declaration
function afl_qemu_trace_arm_exe(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(afl_qemu_trace_arm_exe_path)
    end
end

# Relative path to `afl_qemu_trace_i386_exe`
const afl_qemu_trace_i386_exe_splitpath = ["bin", "afl-qemu-trace-i386"]

# This will be filled out by __init__() for all products, as it must be done at runtime
afl_qemu_trace_i386_exe_path = ""

# afl_qemu_trace_i386_exe-specific global declaration
function afl_qemu_trace_i386_exe(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(afl_qemu_trace_i386_exe_path)
    end
end

# Relative path to `afl_qemu_trace_x86_64_exe`
const afl_qemu_trace_x86_64_exe_splitpath = ["bin", "afl-qemu-trace-x86_64"]

# This will be filled out by __init__() for all products, as it must be done at runtime
afl_qemu_trace_x86_64_exe_path = ""

# afl_qemu_trace_x86_64_exe-specific global declaration
function afl_qemu_trace_x86_64_exe(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(afl_qemu_trace_x86_64_exe_path)
    end
end

# Relative path to `afl_test_instr_exe`
const afl_test_instr_exe_splitpath = ["bin", "afl-test-instr"]

# This will be filled out by __init__() for all products, as it must be done at runtime
afl_test_instr_exe_path = ""

# afl_test_instr_exe-specific global declaration
function afl_test_instr_exe(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(afl_test_instr_exe_path)
    end
end

# Inform that the wrapper is available for this platform
wrapper_available = true

"""
Open all libraries
"""
function __init__()
    # This either calls `@artifact_str()`, or returns a constant string if we're overridden.
    global artifact_dir = find_artifact_dir()

    global PATH_list, LIBPATH_list
    # Initialize PATH and LIBPATH environment variable listings
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (Glib_jll.PATH_list, Pixman_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (Glib_jll.LIBPATH_list, Pixman_jll.LIBPATH_list,))

    global afl_gcc_exe_path = normpath(joinpath(artifact_dir, afl_gcc_exe_splitpath...))

    push!(PATH_list, dirname(afl_gcc_exe_path))
    global afl_qemu_trace_aarch64_exe_path = normpath(joinpath(artifact_dir, afl_qemu_trace_aarch64_exe_splitpath...))

    push!(PATH_list, dirname(afl_qemu_trace_aarch64_exe_path))
    global afl_qemu_trace_arm_exe_path = normpath(joinpath(artifact_dir, afl_qemu_trace_arm_exe_splitpath...))

    push!(PATH_list, dirname(afl_qemu_trace_arm_exe_path))
    global afl_qemu_trace_i386_exe_path = normpath(joinpath(artifact_dir, afl_qemu_trace_i386_exe_splitpath...))

    push!(PATH_list, dirname(afl_qemu_trace_i386_exe_path))
    global afl_qemu_trace_x86_64_exe_path = normpath(joinpath(artifact_dir, afl_qemu_trace_x86_64_exe_splitpath...))

    push!(PATH_list, dirname(afl_qemu_trace_x86_64_exe_path))
    global afl_test_instr_exe_path = normpath(joinpath(artifact_dir, afl_test_instr_exe_splitpath...))

    push!(PATH_list, dirname(afl_test_instr_exe_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()


text:
const GPS_CA_CODES = extend_front_and_back(read_in_codes(
    joinpath(dirname(pathof(GNSSSignals)), "..", "data", "codes_gps_l1.bin"),
    37,
    1023
))

"""
$(SIGNATURES)

Get codes of type GPSL1 as a Matrix where each column
represents a PRN.
```julia-repl
julia> get_code(GPSL1)
```
"""
function get_codes(::Type{GPSL1})
    @view GPS_CA_CODES[2:end - 2, :]
end

"""
$(SIGNATURES)

Get code length of GNSS system GPSL1.
```julia-repl
julia> get_code_length(GPSL1)
```
"""
@inline function get_code_length(::Type{GPSL1})
    1023
end

"""
$(SIGNATURES)

Get secondary code length of GNSS system GPSL1.
```julia-repl
julia> get_secondary_code_length(GPSL1)
```
"""
@inline function get_secondary_code_length(::Type{GPSL1})
    1
end

"""
$(SIGNATURES)

Get center frequency of GNSS system GPSL1.
```julia-repl
julia> get_center_frequency(GPSL1)
```
"""
@inline function get_center_frequency(::Type{GPSL1})
    1_575_420_000Hz
end

"""
$(SIGNATURES)

Get code frequency of GNSS system GPSL1.
```julia-repl
julia> get_code_frequency(GPSL1)
```
"""
@inline function get_code_frequency(::Type{GPSL1})
    10_230_000Hz
end

"""
$(SIGNATURES)

Get data frequency of GNSS system GPSL1.
```julia-repl
julia> get_data_frequency(GPSL1)
```
"""
@inline function get_data_frequency(::Type{GPSL1})
    50Hz
end

"""
$(SIGNATURES)

Get code of GNSS system GPSL1 at phase `phase` of prn `prn`.
The phase will not be wrapped by the code length. The phase has to smaller
than the code length and must be an integer.
```julia-repl
julia> get_code_unsafe(GPSL1, 10, 1)
```
"""
Base.@propagate_inbounds function get_code_unsafe(
    ::Type{GPSL1},
    phase::Integer,
    prn::Integer
)
    GPS_CA_CODES[2 + phase, prn]
end


text:
#---------------------------------------------------------
# # [Plot basic things](@id plots)
#---------------------------------------------------------

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/2_plot.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/2_plot.ipynb)

#md # !!! note
#md #     All the AIBECS tutorials and how-to guides are available as Jupyter notebooks.
#md #     You can view them with [nbviewer](https://nbviewer.jupyter.org/)or execute them online with [binder](https://mybinder.org/) by clicking on the badges above!
#md #     (Note that binder can be slow to launch and its memory caps can be a problem when running.)

# This guide is organized as follows
# - [Horizontal maps](@ref horizontal-plots)
# - [Vertical slices](@ref vertical-plots)
# - [Depth profiles](@ref profile-plots)

# In this guide we will focus on how-to plot things using AIBECS' built-in recipes for [Plots.jl](https://github.com/JuliaPlots/Plots.jl).
# These recipes are implemented using [RecipesBase.jl](https://github.com/JuliaPlots/RecipesBase.jl), which are explained in [Plots.jl's documentation](https://docs.juliaplots.org/latest/recipes/).

# Throughout we will use the OCIM2 grid and we will create a `dummy` tracer as a function of location to showcase each plot, just for the sake of the examples herein.

using AIBECS, Plots
grd, _ = OCIM2.load()
dummy = cosd.(latvec(grd))

#--------------------------------------------
# ## [Horizontal plots](@id horizontal-plots)
#--------------------------------------------

# ### Horizontal slice

# The most common thing you plot after a simulation of marine tracers is a horizontal slice.
# In this case, you just need to provide the tracer (`dummy` here), the grid object `grd`, and the depth at which you want to plot.

plothorizontalslice(dummy, grd, depth=10)

# You can supply units for the depth at which you want to see the horizontal slice.

plothorizontalslice(dummy, grd, depth=10u"m")

# And the units should be understood under the hood.

plothorizontalslice(dummy, grd, depth=3u"km")

# If your tracer is supplied with units, those will show in the colorbar label

plothorizontalslice(dummy * u"mol/m^3", grd, depth=10u"m")

# The advantage of Plots.jl recipes like this one is that you can specify other pieces of the plot as you would with built-in functions.
# The advantage of Plots.jl recipes like this one is that you can specify other pieces of the plot as you would with built-in functions.
# For example, you can chose the colormap with the `color` keyword argument.

dummy .*= cosd.(lonvec(grd))
plt = plothorizontalslice(dummy, grd, depth=100, color=:balance)

# And you can finetune attributes after the plot is created.

plot!(plt, xlabel="Lon", ylabel="Lat", colorbar_title="dummy value", title="The pacific as a whole")

#----------------------------------------
# ## [Vertical plots](@id vertical-plots)
#----------------------------------------

# Exploring the vertical distribution of tracers is important after all.

# ### Zonal slices

# You must specify the longitude

dummy = cosd.(latvec(grd))
dummy .+= sqrt.(depthvec(grd)) / 30
plotmeridionalslice(dummy, grd, lon=330)

# ### Zonal averages

# #### Global zonal average

zonalaverage(dummy, grd)

# If you want a zonal average over a specific region, you can just mask it out

# #### Basin zonal average

# This is experimental at this stage and relies on [OceanBasins.jl](https://github.com/briochemc/OceanBasins.jl).
# You can create basin masks using this package with

using OceanBasins
OCEANS = oceanpolygons()
basins = sum(i*isbasin(latvec(grd), lonvec(grd), OCEANS) for (i,isbasin) in enumerate([isindian2, ispacific2, isatlantic2, isantarctic]))
plothorizontalslice(basins, grd, depth=0, seriestype=:heatmap, color=:lightrainbow)

# and you can mask a specific region with the `mask` keyword argument

mPAC = ispacific(latvec(grd), lonvec(grd), OCEANS)
plotzonalaverage(dummy, grd, mask=mPAC)

# ### Meridional slices

# Just as you should expect at this stage, you can plot a meridional slice with

plotmeridionalslice(dummy, grd, lon=-30)

#----------------------------------------------------
# ## [Depth profiles](@id profile-plots)
#----------------------------------------------------

# Sometimes you want a profile at a given station or location

plotdepthprofile(dummy, grd, lonlat=(-30,30))


text:
function writesdpasparse{T<:Real}(io::IO, sdp::SparseSDP{T})
    if !normalized(sdp)
        sdp, = normalize(sdp)
    end

    if !ismaximizationproblem(sdp)
        sdp.obj = -one(T)*sdp.obj
    end

    bsizes = blocksizes(sdp)

    println(io, ncons(sdp))
    println(io, length(bsizes))
    for i = 1:length(bsizes)
        print(io, bsizes[i])
        print(io, ' ')
    end
    println(io)

    rhs = rhsdense(sdp)
    for r in rhs[1:end-1]
        print(io, r)
        print(io, ' ')
    end
    if length(rhs) > 0
        println(io, rhs[end])
    end

    for (bi, block) in obj(sdp)
        for ((i, j), v) in block
            if i <= j || !haskey(block.entries, (j, i))
                print(io, '0')
                print(io, ' ')
                print(io, bi)
                print(io, ' ')
                print(io, min(i, j))
                print(io, ' ')
                print(io, max(i, j))
                print(io, ' ')
                print(io, (block[i,j] + block[j,i])/2)
                println(io)
            end
        end
    end

    for (ci, bm) in cons(sdp)
        for (bi, block) in bm
            for ((i, j), v) in block
                if i <= j || !haskey(block.entries, (j, i))
                    print(io, ci)
                    print(io, ' ')
                    print(io, bi)
                    print(io, ' ')
                    print(io, min(i, j))
                    print(io, ' ')
                    print(io, max(i, j))
                    print(io, ' ')
                    print(io, (block[i,j] + block[j,i])/2)
                    println(io)
                end
            end
        end
    end
end

function readsdpasparse(io::IO)
    sdp = NormalizedSparseSDP()

    s = 1
    for l in readlines(io)
        if !beginswith(l, '#')
            if s == 4
                t = split(l)
                for i = 1:length(t)
                     setrhs!(sdp, i, parsefloat(t[i]))
                end
            elseif s > 4
                t = split(l)
                ci = parseint(t[1])
                bi = parseint(t[2])
                i = parseint(t[3])
                j = parseint(t[4])
                v = parsefloat(t[5])
                if ci == 0
                    setobj!(sdp, bi, i, j, v)
                else
                    setcon!(sdp, ci, bi, i, j, i == j ? v : 2v)
                end
            end
            s += 1
        end
    end

    sdp
end


text:
using Test, MadNLP, MadNLPTests, MINLPTests
import MathOptInterface
import AmplNLReader: AmplModel
import SparseArrays: sparse

@testset "MadNLP test" begin
    @testset "Matrix tools" begin
        include("matrix_test.jl")
    end

    @testset "MOI interface" begin
        include("MOI_interface_test.jl")
    end

    # this is temporarily commented out due to the incompatibility between NLPModels v0.17.2 and AmplNLReader v0.11.0
    # @testset "NLPModels interface" begin
    #     include("nlpmodels_test.jl")
    # end

    @testset "MadNLP test" begin
        include("madnlp_test.jl")
        include("madnlp_dense.jl")
    end

    @testset "MINLP test" begin
        include("minlp_test.jl")
    end
end # @testset


text:
module PokerGUI

end # module


text:
# This is a port of some of the Farmer Toda (2016) MATLAB code
# http://onlinelibrary.wiley.com/doi/10.3982/QE737/abstract

# __precompile__()
#
# module markovProcesses

using NLopt
using Distributions
using Optim
using ProgressMeter
using StatsBase

export VAR_process, VAR_states_transition, sparsify_transition_matrix

# -------------------- types --------------------------

"Includes the states & transition for discretized VAR process"
struct VAR_states_transition{T<:AbstractFloat}
  X::Array{T}
  P::Array{T,2}
end

"VAR process object that we will discretize"
struct VAR_process{T<:AbstractFloat}
  b::Array{T}
  B::Array{T,2}
  Ψ::Array{T,2}
  M::Int
  A::Array{T,2}
  C::Array{T,2}
  Σ::Array{T,2}
  μ::Array{T,2}
end

# -------------------- cosntructors --------------------------

"Wrapper constructs markov Transition matrix & discrete states from VAR_process"
function VAR_states_transition(vpp::VAR_process{T}, states_per_dim::Int=4, num_moments::Int=2, nσ::T=0.0) where {T<:AbstractFloat}
  yy = collect(yspace(vpp, states_per_dim, nσ))
  SS = nprod(yy, vpp.M)
  DD = makeD(SS)

  P, JN, Λ, nM, er = makeTransitionMatrix(yy, SS, DD, vpp, num_moments)
  transitions = probabilityProd(P)
  StatsBase.countmap(nM)

  states = vpp.C*DD .+ vpp.μ

  return VAR_states_transition(states, transitions)
end

"Constructor for univariate AR process"
function VAR_process(b::T, B::T, Ψ::T) where {T<:AbstractFloat}
  o2 = ones(T, 1, 1)
  o1 = ones(T, 1)

  M = 1.0
  C = sqrt(Ψ)   * o2;
  A = B         * o2;
  μ =   b/(1.0-B)   * o2;
  Σ = 1.0/(1.0-B^2) * o2;

  return VAR_process(b*o1, B*o2, Ψ*o2, M, A, C, Σ, μ)
end

"Constructor for VAR process"
function VAR_process(b::Array{T,1}, B::Array{T,2}, Ψ::Array{T,2}) where {T<:AbstractFloat}
  size(B) == size(Ψ)     || error("B and Ψ must be the same size")
  length(b) == size(B,1) || error("b must be same length as B's width")

  M = size(B,1)

  b = reshape(b,M,1)
  tildeC = chol(Ψ)'
  tildeA = tildeC \ B * tildeC
  tildeΣ = reshape( (eye(T,M^2) - kron(tildeA,tildeA)) \ vec(eye(T,M)), M, M)

  # Objective is to make unconditional varances of y_k's equal
  function f(Vvec::Array{T,1}, grad::Array{T,1}) where {T}
    V = reshape(Vvec,M,M)
    return norm( diag(V'*tildeΣ*V) .- trace(tildeΣ)/M )
  end

  # Constraint is that U is orthognoal (U'U = I)
  function c(res::Array{T,1}, Vvec::Array{T,1}, grad::Array{T,2}) where {T}
    V = reshape(Vvec,M,M)
    res[:] = vec(V'*V - eye(M))
  end

  # find best U matrix such that y has same unconditional variance
  opt = Opt(:LN_COBYLA, M^2)
  min_objective!(opt, f)
  equality_constraint!(opt, c, 1e-12*ones(M^2))
  xtol_rel!(opt, 1e-12)
  (fObjOpt, xOpt, flag) = NLopt.optimize(opt, vec(eye(M)))

  # make transformation matrices
  U = reshape(xOpt,M,M)
  C = tildeC*U
  A = U'*tildeA*U
  μ = (eye(M) - B)\b
  Σ = U'*tildeΣ*U

  return VAR_process(b,B,Ψ,M,A,C,Σ,μ)
end

# -------------------- helper functions --------------------------

"Linear space for y given a VAR process and length Nm"
function yspace(vp::VAR_process{T}, Nm::Int=5, nσ::T=0.0) where {T<:AbstractFloat}
  σmin = sqrt(minimum(eigvals(vp.Σ)))
  if nσ <= zero(T)
    nσ = sqrt(convert(T, Nm) - one(T))
  end
  return range(-σmin*nσ, stop=σmin*nσ, length=Nm)
end

"Make Cartesian Product given state space"
function nprod(x::Array{T,1}, n::Int) where {T<:Real}
    if n<1
      return error("n must be >= 1")
    end
    if n==1
      return x
    elseif n==2
      return Base.Iterators.Prod2(x, nprod(x,n-1))
    else
      return Base.Iterators.Prod(x, nprod(x,n-1))
    end
end

"Make state space matrix from vector"
function makeD(S::Vector{T}) where {T<:AbstractFloat}
  return reshape(S, (1, length(S) ) )
end

"Make state space matrix from product iterator"
function makeD(S::Base.Iterators.ProductIterator)
  T = eltype(eltype(S))
  D = zeros(T, length(size(S)), length(S) )
  for (j,s) in enumerate(S)
    D[:,j] = collect(s)
  end
  return D
end

"Moments of scaled normal distribution"
function Tbar(L::Int, δ_m::AbstractFloat)
   out = zeros(eltype(δ_m), L)
   nz = filter(iseven, 1:L)
   out[nz] = [2^(-k/2) * factorial(k) / factorial(k/2) * (1/δ_m)^k for k in nz]
  return(out)
end

"Given 3-d matrix of (independent) probability distributions get the product of them"
function probabilityProd(x::Array{T,3}) where {T<:AbstractFloat}
  J =size(x,1)
  XX = Array(T, J, J)

  function myprod2(y::Array{T,2}, m::Int)
    Mm = size(y,1)-m+1
    if m==1
      return y[Mm,:]
    elseif m==2
      return Base.Iterators.Prod2(y[Mm,:], myprod2(y, m-1))
    else
      return Base.Iterators.Prod(y[Mm,:], myprod2(y, m-1))
    end
  end

  for j in 1:J
    XX[j,:] = map(prod, myprod2(x[j,:,:], size(x,2)))
  end
  return XX
end

# -------------------- big function that makes transition matrix --------------------------

# return a PDF
normpdf(x::Real) = pdf(Normal(), x)
normpdf(x::Vector{<:Real}) = normpdf.(x)

# objective
f(x::Vector{T}, q::Vector{T}, ΔT::Matrix{T}) where {T} = dot(q,exp(x'*ΔT[1:l,:]))

# gradient
function g!(grad::Vector{T}, x::Vector{T}, q::Vector{T}, ΔT::Matrix{T}) where {T}
  grad .= zero(T)
  for i in 1:length(q)
    grad .+= q[i] * exp( dot(x, ΔT[1:l,i]) ) * ΔT[1:l,i]
  end
end

# hessian
function h!(hess::Matrix{T}, x::Vector{T}, q::Vector{T}, ΔT::Matrix{T}) where {T}
  hess .= zero(T)
  for i in 1:length(q)
    hess .+= q[i] * exp( dot(x, ΔT[1:l,i]) ) * (ΔT[1:l,i]*ΔT[1:l,i]')
  end
end

"""
    makeTransitionMatrix(y::Vector{T}, S::Union{Vector{T}, Base.AbstractProdIterator}, D::Matrix{T}, vp::VAR_process, L::Int=2, κ::T=1e-8)

Make transition matrix for discretized VAR process with `M` dimensions. The discretization
allows for each dimension to take on `Nm = length(y)` possible values. This means that there
are `Nm^M` possible states `j`. We try to match up to `L` moments.

Each row `i` of the `Nm^M` x `Nm^M` Markov transition matrix is formed from `P::Array{T}(Nm^M, M, Nm)`
by taking the Kronecker product of the `M` columns of matrix `P[i, :, :]`.

# Arguments
* `y::Vector` provides a 'generic' grid of length which is transformed into the state space
* `S::Union{Vector{T}, Base.AbstractProdIterator}`, the primitive of the state-space
* `D::Matrix{T}` is a `M` by `Nm^M` grid of states
* `vp::VAR_process`
* `L::Int=2` the number of moments we try to match
* `κ::T=1e-8` we set the initial guess for the distribution as the max of the pdf over the
    discrete states & `κ`

# Returns
* `P::Array{T}(Nm^M, M, Nm)` matrix is probability, given we are in one state `i ∈ Nm^M`
    states that the `m ∈ M` component is equal to `j ∈ Nm`
* `JN::Matrix{T}(Nm^M, M)` is set of `min` for the objective function
* `Λ::Matrix{T}(Nm^M, M, L)` is set of `argmin`s
* `numMoments::Matrix{Int}(Nm^M, M)` number of moments successfully matched
* `approxErr::Matrix{T}(Nm^M, M, L)` approsimation error for each moment
"""
function makeTransitionMatrix(y::Vector{T}, S::PT, D::Array{T}, vp::VAR_process{T}, L::Int=2, κ::T=1e-8) where {T<:AbstractFloat, PT<:Union{AbstractVector{T}, Base.AbstractProdIterator}}

  if L >= length(y)
    warn("Using fewer moments. Gave ", L, " using ", length(y)-1)
    L = length(y)-1
  end

  # Initialize elements that will be returned
  Λ          = zeros(T  , length(S), vp.M, L)
  P          = zeros(T  , length(S), vp.M, length(y))
  JN         = zeros(T  , length(S), vp.M)
  approxErr  = zeros(T  , length(S), vp.M, L)
  numMoments = zeros(Int, length(S), vp.M)

  # preallocate these, which will be updated each iteration
  ΔT   = Array(T,L,length(y))
  dev  = similar(y)
  q    = similar(y)

  δ = y[end]  # a scaling factor

  # Conditional on each state (j), examine each of the m state vs
  prog = Progress(length(S), 5)

  for j in 1:length(S)
    for m in 1:vp.M
      dev .= (y - dot(vp.A[m,:], D[:,j]))
      q .= max(normpdf.(dev), κ)

      for l in L:-1:1
        ΔT[1:l,:] = (dev ./ δ )'.^(1:l) .- Tbar(l, δ)

        # closures
        f_cl = f(x::Vector) = f(x, q, ΔT)
        g_cl!(grad::Vector, x::Vector) = g!(grad, x, q, ΔT)
        h_cl!(hess::Matrix, x::Vector) = h!(hess, x, q, ΔT)

        # optimize to match moments
        try
          res = Optim.optimize(f_cl, g_cl!, h_cl!, ones(T,l))
          λ = Optim.minimizer(res)
          J_candidate = Optim.minimum(res)
          grad = zeros(T,l)
          g!(λ,grad)

          # if we like the results, update and break
          if ( norm( grad ./ J_candidate ) < 1e-5 ) & all(isfinite.(grad)) & all(isfinite.(λ)) & (J_candidate > 0.0)
            JN[j,m] = J_candidate
            Λ[j,m,1:l] .= λ
            for k in 1:length(y)
              P[j,m,k] = q[k] * exp( dot(λ, ΔT[1:l,k]) ) / J_candidate
            end
            approxErr[j,m,1:l] .= grad ./ J_candidate
            numMoments[j,m] = l
            break
          end # if statment
        catch
        end

      end   # loop over moment number of conditions (l=L:1)
    end     # loop over number of states (m = 1:M)
    ProgressMeter.next!(prog)
  end       # loop over state space (j=1:J)

  return P, JN, Λ, numMoments, approxErr
end

# -------------------- make striped versions --------------------------

"""
    striped_bool(m::Integer, n::Integer, ur::UnitRange)::Matrix{Bool}(m,n)

Returns striped array where trues are on diagonals `ur`
"""
function striped_bool(m::Integer, n::Integer, ur::UnitRange)
  out   = fill(false, m, n)
  trues = fill(true, m, n)
  for i in ur
    out .+= diagm(diag(trues, i), i)
  end
  return out
end

"""
    stripe_and_sparsify_transition_matrix{T<:AbstractFloat}(P::Matrix{T}, minp::T, which_diags::UnitRange)::SparseMatrixCSC{T}

Construct sparse transition matrix from `P` where elements not on diagonals `which_diags`
or less than `minp` are zeroed out. (Ensures that rows sum to 1)
"""
function sparsify_transition_matrix(P::Matrix{T}, minp::T, which_diags::UnitRange) where {T<:AbstractFloat}

  sb = striped_bool(size(P, 1), size(P, 2), which_diags)
  P_big = P .> minp
  P_new = P .* sb .* P_big
  P_new_rowsums = sum(P_new, dims=2)
  for i in 1:size(P_new, 1)
    P_new[i, :] .= P_new[i, :] ./ P_new_rowsums[i]
  end

  return sparse(P_new)

end

function sparsify_transition_matrix(P::Matrix{T}, minp::T)::SparseMatrixCSC{T} where T<:AbstractFloat

  P_big = P .> minp
  P_new = P .* P_big
  P_new_rowsums = sum(P_new, dims=2)
  for i in 1:size(P_new, 1)
    P_new[i, :] .= P_new[i, :] ./ P_new_rowsums[i]
  end

  return sparse(P_new)

end

# # End module
# end


text:
module Saw

import GreyDecision.GreyNumbers.GreyNumber
import GreyDecision.Utility: unitize, rowsums, weightize, reverseminmax

function saw(decisionMat::Array{GreyNumber{T}, 2}, weights::Array{Float64,1}, fns::Array{Function,1};  v::Float64=0.5) where {T <: Real}
    n, p = size(decisionMat)

    normalizedDecisionMat = similar(decisionMat)

    w = unitize(weights)

    colminmax = zeros(GreyNumber{Float64}, p)

    for i in 1:p
        colminmax[i] = decisionMat[:, i] |> fns[i]
        if fns[i] == maximum
            normalizedDecisionMat[:, i] = map(x -> x / colminmax[i], decisionMat[:, i])
        elseif fns[i] == minimum
            normalizedDecisionMat[:, i] = map(x -> colminmax[i] / x, decisionMat[:, i])
        else
            @error fns[i]
            error("Function not found")
        end
    end

    scores = weightize(normalizedDecisionMat, w) |> rowsums

    rankings = scores |> sortperm |> reverse

    bestIndex = rankings |> first

    return Dict(
        :scores => scores,
        :rankings => rankings,
        :bestindex => bestIndex
    )
end

end # End of module saw


text:
# This implements viscous creep laws and routines to compute with them
#
# Note that various simple creep laws are defined in this file;
# more complex ones (such as DislocationCreep) are in separate files
# in the same directory
#
# In case you want to add new creep laws, have a look at how the ones
# here are implemented. Please add tests as well!

export  DislocationCreep,
        SetDislocationCreep

# Dislocation Creep ------------------------------------------------
"""
    DislocationCreep(n = 1.0NoUnits, r = 0.00.0NoUnits, A = 1.5MPa/s, E = 476.0kJ/mol, V = 6e-6m^3/mol, apparatus = "AxialCompression" )

Defines the flow law parameter of a dislocation creep law

The (isotropic) dislocation creep law, as used by experimtalists, is given by
```math
     \\dot{\\gamma} = A \\sigma_\\mathrm{d}^n f_\\mathrm{H2O}^r \\exp(-\\frac{E+PV}{RT})
```

where ``n`` is the power law exponent,
``r`` is the exponent of fugacity dependence,
``A`` is a pre-exponential factor [MPa^(n+r)] (if manually defined, n and r must be either pre-defined or substituted),
``E`` is the activation energy [kJ/mol], ``V`` is the activation volume [m^3/mol]. ``\\dot{\\gamma}`` is the ordinary strain rate [1/s],
and ``\\sigma_\\mathrm{d}`` is the differential stress which are converted into second invariants using the apparatus type that can be
either "AxialCompression", "SimpleShear" or "Invariant".
If the flow law paramters are already given as a function of second invariants, choose apparatus = "Invariant"

```julia-repl
julia> x2      =   DislocationCreep(n=3)
DislocationCreep: n=3, r=0.0, A=1.5 MPa^-3 s^-1, E=476.0 kJ mol^-1, V=6.0e-6 m^3 mol^-1, Apparatus=AxialCompression
```
"""
@with_kw_noshow mutable struct DislocationCreep <: AbstractCreepLaw
    equation::LaTeXString   =   L"\tau_{ij} = 2 \eta  \dot{\varepsilon}_{ij}" # TO BE UPDATED
    n::GeoUnit        = 1.0NoUnits         # power-law exponent
    r::GeoUnit        = 0.0NoUnits         # exponent of water-fugacity dependence
    A::GeoUnit        = 1.5MPa^(-n-r)/s    # pre-exponential factor
    E::GeoUnit        = 476.0kJ/mol        # activation energy
    V::GeoUnit        = 6e-6m^3/mol        # activation volume
    R::GeoUnit        = 8.314J/mol/K       # Universal gas constant
    Apparatus         = "AxialCompression" # type of experimental apparatus, either AxialCompression, SimpleShear or Invariant
    Comment::String   = ""                 # Some remarks you want to add about this creep law implementation
    BibTex_Reference  = ""                 # BibTeX reference
end

# Calculation routines for linear viscous rheologies
# All inputs must be non-dimensionalized (or converted to consitent units) GeoUnits
function ComputeCreepLaw_EpsII(TauII, a::DislocationCreep, p::CreepLawVariables)
    @unpack n           = a
    @unpack r           = a
    @unpack A           = a
    @unpack E           = a
    @unpack V           = a
    @unpack R           = a
    @unpack Apparatus   = a
    @unpack P           = p
    @unpack T           = p
    @unpack f           = p
    if Apparatus == "AxialCompression"
        FT = sqrt(3.0)NoUnits               # relation between differential stress recorded by apparatus and TauII
        FE = 2.0/sqrt(3.0)NoUnits            # relation between gamma recorded by apparatus and EpsII
    elseif Apparatus == "SimpleShear"
        FT = 2.0NoUnits                      # it is assumed that the flow law parameters were derived as a function of differential stress, not the shear stress. Must be modidified if it is not the case
        FE = 2.0NoUnits
    elseif Apparatus == "Invariant"
        FT = 1.0NoUnits
        FE = 1.0NoUnits
    end
    return A.val*(TauII.val*FT)^n.val*f.val^r.val*exp(-(E.val+P.val*V.val)/(R.val*T.val))/FE;
end

# EpsII .= A.*(TauII.*FT).^n.*f.^r.*exp.(-(E.+P.*V)./(R.*T))./FE; Once we have a
# All inputs must be non-dimensionalized (or converted to consitent units) GeoUnits
function ComputeCreepLaw_TauII(EpsII, a::DislocationCreep, p::CreepLawVariables)
    @unpack n           = a
    @unpack r           = a
    @unpack A           = a
    @unpack E           = a
    @unpack V           = a
    @unpack R           = a
    @unpack Apparatus   = a
    @unpack P           = p
    @unpack T           = p
    @unpack f           = p
    if Apparatus == "AxialCompression"
        FT = sqrt(3.0)               # relation between differential stress recorded by apparatus and TauII
        FE = 2.0/sqrt(3.0)           # relation between gamma recorded by apparatus and EpsII
    elseif Apparatus == "SimpleShear"
        FT = 2.0                     # it is assumed that the flow law parameters were derived as a function of differential stress, not the shear stress. Must be modidified if it is not the case
        FE = 2.0
    elseif Apparatus == "Invariant"
        FT = 1.0
        FE = 1.0
    end
    return A.val^(-1/n.val)*(EpsII.val*FE)^(1/n.val)*f.val^(-r.val/n.val)*exp((E.val+P.val*V.val)/(n.val*R.val*T.val))/FT;
end

# Print info
function show(io::IO, g::DislocationCreep)
    print(io, "DislocationCreep: n=$(g.n.val), r=$(g.r.val), A=$(g.A.val), E=$(g.E.val), V=$(g.V.val), Apparatus=$(g.Apparatus)" )
end
#-------------------------------------------------------------------------

# Add pre-defined creep laws
"""
    SetDislocationCreep["Name of Dislocation Creep"]

This is a dictionary with pre-defined creep laws
"""
SetDislocationCreep = Dict([

# Olivine rheology
("Dry Olivine | Hirth & Kohlstedt (2003)",
# after Hirth, G. & Kohlstedt (2003), D. Rheology of the upper mantle and the mantle wedge: A view from the experimentalists.
# Inside the subduction Factory 83?105. Table 1, "dry dislocation" parameters
    DislocationCreep(
        n = 3.5NoUnits,
        A = 1.1e5MPa^(-3.5)/s,
        E = 530kJ/mol,
        V = 15e-6m^3/mol,
        Apparatus =   "AxialCompression",
        r = 0NoUnits,
        Comment = "Still to be verified with the original publication (BK). Values checked, plots are not reproduced (DK).",
        BibTex_Reference = parse_bibtex("""
            @incollection{eiler_rheology_2003,
            address = {Washington, D. C.},
            title = {Rheology of the upper mantle and the mantle wedge: {A} view from the experimentalists},
            volume = {138},
            isbn = {978-0-87590-997-4},
            shorttitle = {Rheology of the upper mantle and the mantle wedge},
            url = {http://www.agu.org/books/gm/v138/138GM06/138GM06.shtml},
            language = {en},
            urldate = {2019-10-09},
            booktitle = {Geophysical {Monograph} {Series}},
            publisher = {American Geophysical Union},
            author = {Hirth, Greg and Kohlstedt, David},
            editor = {Eiler, John},
            year = {2003},
            doi = {10.1029/138GM06},
            pages = {83--105},
            }
        """)
    )
)

# Olivine rheology
("Wet Olivine | Hirth & Kohlstedt (2003)",
    # After Hirth, G. & Kohlstedt (2003), D. Rheology of the upper mantle and the mantle wedge: A view from the experimentalists.
    #   Inside the subduction Factory 83?105. Table 1, "wet dislocation" parameters
    #  Note that this assumes C_OH=1000
    DislocationCreep(
        n = 3.5NoUnits,
        A = 90MPa^(-3.5)/s,
        E = 480kJ/mol,
        V = 11e-6m^3/mol,
        r   = 1.2NoUnits,
        Apparatus =   "AxialCompression",
        Comment = "Still to be verified with the original publication (BK). Values checked, plots are not reproduced (DK).",
        BibTex_Reference = parse_bibtex("""
            @incollection{HirthKohlstedt_OlivineRheology_2003,
            address = {Washington, D. C.},
            title = {Rheology of the upper mantle and the mantle wedge: {A} view from the experimentalists},
            volume = {138},
            isbn = {978-0-87590-997-4},
            shorttitle = {Rheology of the upper mantle and the mantle wedge},
            url = {http://www.agu.org/books/gm/v138/138GM06/138GM06.shtml},
            language = {en},
            urldate = {2019-10-09},
            booktitle = {Geophysical {Monograph} {Series}},
            publisher = {American Geophysical Union},
            author = {Hirth, Greg and Kohlstedt, David},
            editor = {Eiler, John},
            year = {2003},
            doi = {10.1029/138GM06},
            pages = {83--105},
            }
        """);
    )
)

]); # end of setting pre-defined creep laws


text:
@testset "365.water-and-jug-problem.jl" begin
    @test can_measure_water(3, 5, 4)
    @test !can_measure_water(2, 6, 5)
    @test can_measure_water(11, 3, 13)
    @test !can_measure_water(11, 3, 20)
end


text:
######################
# 1: The Julia type for ToricVarieties
######################

struct toric_variety
           GapToricVariety::GapObj
end
export toric_variety

######################
# 2: Generic constructors
######################

function toric_variety( rays::Vector{Vector{Int}}, cones::Vector{Vector{Int}} )
    # construct the toric variety in GAP
    gap_rays = GapObj( rays, recursive = true )
    gap_cones = GapObj( cones, recursive = true )
    fan = GAP.Globals.Fan( gap_rays, gap_cones )
    variety = GAP.Globals.ToricVariety( fan )

    # wrap it into a struct and return
    return toric_variety( variety )
end
export toric_variety

######################
# 3: Special constructors
######################

function projective_space( x )
    # construct the projective space in gap
    variety = GAP.Globals.ProjectiveSpace( x )

    # wrap it and return
    return toric_variety(  variety )
end
export projective_space

######################
# 4: Properties
######################

function is_normal_variety( v::toric_variety )
    return true
end
export is_normal_variety

function is_affine( v::toric_variety )
    return GAP.Globals.IsAffine(v.GapToricVariety)::Bool
end
export is_affine

function is_projective( v::toric_variety )
    return GAP.Globals.IsProjective( v.GapToricVariety )::Bool
end
export is_projective

function is_smooth( v::toric_variety )
    return GAP.Globals.IsSmooth( v.GapToricVariety )::Bool
end
export is_smooth

function is_complete( v::toric_variety )
    return GAP.Globals.IsComplete( v.GapToricVariety )::Bool
end
export is_complete

function has_torusfactor( v::toric_variety )
    return GAP.Globals.HasTorusfactor( v.GapToricVariety )::Bool
end
export has_torusfactor

function has_no_torusfactor( v::toric_variety )
    return GAP.Globals.HasNoTorusfactor( v.GapToricVariety )::Bool
end
export has_no_torusfactor

function is_orbifold( v::toric_variety )
    return GAP.Globals.IsOrbifold( v.GapToricVariety )::Bool
end
export is_orbifold

function is_simplicial( v::toric_variety )
    return GAP.Globals.IsSimplicial( v.GapToricVariety )::Bool
end
export is_simplicial

function is_isomorphic_to_projective_space( v::toric_variety )
    return GAP.Globals.IsIsomorphicToProjectiveSpace( v.GapToricVariety )::Bool
end
export is_isomorphic_to_projective_space

function is_direct_product_of_projective_spaces( v::toric_variety )
    return GAP.Globals.IsDirectProductOfPNs( v.GapToricVariety )::Bool
end
export is_direct_product_of_projective_spaces

max_stars_repo_name: fingolfin/JToric

max_issues_repo_name: fingolfin/JToric

max_forks_repo_name: fingolfin/JToric


text:
module test_errors_error

using Test

@testset "error testset" begin

@test true

g() = must_be_an_error
f() = g()
f()

end

@testset "error testset 2" begin

@test false
@test false
@test false
throw(ErrorException("error"))

end

end # module test_errors_error


text:
module GeometryModule
export Geometry, readGeometryXYZ, bohr2angstrom, angstrom2bohr, angstrom2bohr!, computeEnergyInteratomicRepulsion, computeNumberElectrons

using ..BaseModule
using ..AtomModule

type Geometry
	atoms::Array{Atom,1}
end

"""
    readGeometryXYZ(filename::AbstractString; unit::Symbol=:Bohr)

reads an .xyz-Geometry file (acc. to the unofficial standard the coordinates
of which are interpreted as Angstrom) and yields a Geometry object in `unit`.
"""
function readGeometryXYZ(filename::AbstractString; unit::Symbol=:Bohr)
	geo = Geometry([])
	lines = open(readlines,filename)
	for line in lines[3:end]	# line 1 is n. of atoms, line 2 is comment
		columns = split(line)
		elem = Element(columns[1])
		pos = Position(map(parse,columns[2:4])...)	# float is not recommended for use on tuples
		append!(geo.atoms,[Atom(elem,pos)])
	end
	if (unit == :Bohr) angstrom2bohr!(geo) end	# Should be default as Bohr is an atomic unit
	return geo
end

function Geometry(filename::AbstractString)
  if isfile(filename)
    return readGeometryXYZ(filename)
  end
  for ext in ("xyz", "XYZ")
    if isfile("$filename.$ext")
      return readGeometryXYZ("$filename.$ext")
    end
  end
  error("Couldn't find file: $filename")
end

angstrom2bohr(x::Real) = x/0.529177210
bohr2angstrom(x::Real) = x*0.529177210
angstrom2bohr(pos::Position) = (return Position(angstrom2bohr(pos.x),angstrom2bohr(pos.y),angstrom2bohr(pos.z)))
function angstrom2bohr!(geo::Geometry)
  for atom in geo.atoms
    atom.position = angstrom2bohr(atom.position)
  end
end

function computeEnergyInteratomicRepulsion(
  atom1::Atom,
  atom2::Atom)
  q1 = atom1.element.atomicNumber
  q2 = atom2.element.atomicNumber
  A = atom1.position
  B = atom2.position
  return q1*q2/(distance(A,B))
end

function computeEnergyInteratomicRepulsion(geo::Geometry)
  result = 0.
  for A in 1:length(geo.atoms), B in 1:A-1
      result += computeEnergyInteratomicRepulsion(geo.atoms[A],geo.atoms[B])
  end
  return result
end

function computeNumberElectrons(geo::Geometry,charge::Int=0)
  sum(at.element.atomicNumber for at in geo.atoms) - charge
end

end #module


text:
export done, next, isdone, isnext, Continuation, Done, Next

@enum ContinuationState Done Next

struct Continuation{T <: Tuple}
    state::ContinuationState
    data::T
end # Continuation

isdone(continuation::Continuation) = continuation.state == Done
isnext(continuation::Continuation) = continuation.state == Next

done(args...) = Continuation(Done, tuple(args...))
next(args...) = Continuation(Next, tuple(args...))

function Base.show(io::IO, next::Continuation)

    if next.state == Done
        print(io, "done(")
    else
        print(io, "next(")
    end

    comma = false
    for v in next
        if comma
            print(io, ", ")
        end
        comma = true

        show(io, v)
    end

    print(io, ")")
end

Base.iterate(c::Continuation) = Base.iterate(c.data)
Base.iterate(c::Continuation, state) = Base.iterate(c.data, state)
# Base.start(c::Continuation) = Base.start(c.data)
# Base.done(c::Continuation, state) = Base.done(c.data, state)
# Base.next(c::Continuation, state) = Base.next(c.data, state)


text:
# is included from test_simulate.jl
function plot2d(pos, reltime=0.0; zoom=true, front=false, segments=se().segments)
    x = Float64[]
    z = Float64[]
    for i in 1:length(pos)
        if front
            push!(x, pos[i][2])
        else
            push!(x, pos[i][1])
        end
        push!(z, pos[i][3])
    end
    x_max = maximum(x)
    z_max = maximum(z)
    if zoom
        xlabel = "x [m]"
        if front xlabel = "y [m]" end
        plot(x,z, xlabel=xlabel, ylabel="z [m]", legend=false, xlims = (x_max-15.0, x_max+5), ylims = (z_max-15.0, z_max+5))
        annotate!(x_max-10.0, z_max-3.0, "t=$(round(reltime,digits=1)) s")
    else
        plot(x,z, xlabel="x [m]", ylabel="z [m]", legend=false)
        annotate!(x_max-10.0, z_max-3.0, "t=$(round(reltime,digits=1)) s")
    end
    if length(pos) > segments+1
        s=segments
        plot!([x[s+1],x[s+4]],[z[s+1],z[s+4]], legend=false) # S6
        plot!([x[s+2],x[s+5]],[z[s+2],z[s+5]], legend=false) # S8
        plot!([x[s+3],x[s+5]],[z[s+3],z[s+5]], legend=false) # S7
        plot!([x[s+2],x[s+4]],[z[s+2],z[s+4]], legend=false) # S2
        plot!([x[s+1],x[s+5]] ,[z[s+1],z[s+5]],legend=false) # S5
    end
    plot!(x, z, seriestype = :scatter)
end


text:
let pointtypes = (wkbPoint, wkbPoint25D, wkbPointM, wkbPointZM),
    multipointtypes =
        (wkbMultiPoint, wkbMultiPoint25D, wkbMultiPointM, wkbMultiPointZM),
    linetypes =
        (wkbLineString, wkbLineString25D, wkbLineStringM, wkbLineStringZM),
    multilinetypes = (
        wkbMultiLineString,
        wkbMultiLineString25D,
        wkbMultiLineStringM,
        wkbMultiLineStringZM,
    ),
    polygontypes = (wkbPolygon, wkbPolygon25D, wkbPolygonM, wkbPolygonZM),
    multipolygontypes = (
        wkbMultiPolygon,
        wkbMultiPolygon25D,
        wkbMultiPolygonM,
        wkbMultiPolygonZM,
    ),
    collectiontypes = (
        wkbGeometryCollection,
        wkbGeometryCollection25D,
        wkbGeometryCollectionM,
        wkbGeometryCollectionZM,
    )

    function GeoInterface.geotype(g::AbstractGeometry)::Symbol
        gtype = getgeomtype(g)
        return if gtype in pointtypes
            :Point
        elseif gtype in multipointtypes
            :MultiPoint
        elseif gtype in linetypes
            :LineString
        elseif gtype == wkbLinearRing
            :LinearRing
        elseif gtype in multilinetypes
            :MultiLineString
        elseif gtype in polygontypes
            :Polygon
        elseif gtype in multipolygontypes
            :MultiPolygon
        elseif gtype in collectiontypes
            :GeometryCollection
        else
            @warn "unknown geometry type" gtype
            :Unknown
        end
    end

    function GeoInterface.coordinates(g::AbstractGeometry)
        gtype = getgeomtype(g)
        ndim = getcoorddim(g)
        return if gtype in pointtypes
            if ndim == 2
                Float64[getx(g, 0), gety(g, 0)]
            elseif ndim == 3
                Float64[getx(g, 0), gety(g, 0), getz(g, 0)]
            else
                error("getcoorddim($g) returned $ndim: expected 2 or 3")
            end
        elseif gtype in multipointtypes
            Vector{Float64}[
                GeoInterface.coordinates(getgeom(g, i - 1)) for i in 1:ngeom(g)
            ]
        elseif gtype in linetypes || gtype == wkbLinearRing
            Vector{Float64}[
                collect(getpoint(g, i - 1)[1:ndim]) for i in 1:ngeom(g)
            ]
        elseif gtype in multilinetypes || gtype in polygontypes
            Vector{Vector{Float64}}[
                GeoInterface.coordinates(getgeom(g, i - 1)) for i in 1:ngeom(g)
            ]
        elseif gtype in multipolygontypes
            Vector{Vector{Vector{Float64}}}[
                GeoInterface.coordinates(getgeom(g, i - 1)) for i in 1:ngeom(g)
            ]
        end
    end
end


text:
"""
This module is the application's entry point. It is responsible for connecting
to discord and routing messages through handlers.
"""
module Edgar

using Discord

include("./Types.jl")
include("./Actions.jl")
include("./Interactive.jl")

"""
Start the discord bot using the provided bot secret from the developer portal.
The secret can be provided as an argument to `main()` or it can come from the
`DISCORD_BOT_SECRET` environment variable.

Example usage:

    > julia --project=. -e 'using Edgar; Edgar.main()'

"""
function main(secret=ENV["DISCORD_BOT_SECRET"])
    client = Client(secret; presence=(game=(name="λ(x) = x", type=AT_GAME),))

    add_handler!(client, MessageCreate, on_message_create)

    open(client)
    wait(client)
end

"""
This function is invoked any time that a message is created in a server which
uses this bot.
"""
function on_message_create(client::Client, e::MessageCreate)
    if from_self(client, e.message) return end

    action = Interactive.handle(e.message.content, string(e.message.author.id))

    take_action(client, e.message, action)
end

"Evaluates to TRUE when the message was sent by this bot."
function from_self(client::Client, message::Message)::Bool
    bot_user = Discord.me(client)
    result = bot_user.id == message.author.id
    @debug """
    Message '$(message.content)'
    Is $(result ? "from" : "not from") this bot.
    """
    result
end

end


text:
mutable struct MipDiscretSolver
    inst::Instance
    loglevel::Int         # niveau de verbosité
    bestsol::Solution     # meilleure Solution rencontrée
    # Les attributs spécifiques au modèle
    model::Model  # Le modèle MIP
    x         # vecteur des variables d'atterrissage
    b         # vecteur des variables de précédence (before)
    cost      # variable du cout de la solution
    costs     # variables du cout de chaque avion
    # Le constructeur
    function MipDiscretSolver(inst::Instance)
        this=new()
        this.inst = inst
        this.loglevel = Args.get("level")
        # this.bestsol = Solution(inst) # plante dans le proto car timing non dispo !
        this.bestsol = Solution(inst, update=false) # pas besoin de résoudre le timing !
        this.model = new_model(mode=:mip, log_level=this.loglevel)
        return this
    end
end

"""Création du modèle discrétisé (plus lent, mais polyvalent)"""
function buildSeqataModel(sv::MipDiscretSolver)

    # nb d'avion, pour abréger l'écriture
    planes = sv.inst.planes
    n = sv.inst.nb_planes
    model = sv.model

    # Quelques raccourcis car utilisés un peu partout
    lbmin = lb_min(sv.inst)
    ubmax = ub_max(sv.inst)

    ##########################################################
    # Création des variables
    ##########################################################

    # y: indicateur des dates d'atterrissage
    # y=1 si l'avion i atterrit à la date t dans [0,T-1]
    # @variable(model, y[p in planes, 1:ubmax], Bin) ORI pour indice Int et Plane
    @variable(model, y[p in planes, lbmin:ubmax], Bin)

    # variable dérivée : date d'atterrissage effective de l'avion i
    @expression(model, x[p in planes], sum(t*y[p,t] for t in lbmin:ubmax))
    sv.x = x

    # Coût de chaque avion compte tenu de sa date d'atterrissage
    @expression(model, costs[p in planes],
            sum(get_cost(p,t)*y[p,t] for t in p.lb:p.ub))

    # Objectif Minimiser le coût de pénalité total
    @expression(model, total_cost,
                sum(get_cost(p,t)*y[p,t] for p in planes for t in lbmin:ubmax))

    sv.costs = costs
    sv.cost = total_cost

    ##########################################################
    # Définition de l'objectif
    ##########################################################

    @objective(model, Min, total_cost)

    ##########################################################
    # Création des contraintes
    ##########################################################

    # Contrainte sur les bornes de la date d'atterrissage
    # AMPL: C1 {i in 1..n} : E[i] <= x[i] <= L[i];
    @constraint(model, c1[p in planes], p.lb <= x[p] <= p.ub)

    # Chaque avion ne peut avoir qu'une seule date d'atterrissage.
    # AMPL: C2 {i in 1..n} : sum{t in T_SET} y[i,t] = 1;
    @constraint(model, c2[p in planes], sum(y[p,t] for t in lbmin:ubmax) == 1)

    # Deux dates d'atterrissage ne doivent pas être trop rapprochées
    # Principe :
    # - on prend toutes les paires d'avions possibles
    # - pour toutes les dates d'atterissage possibles de p1
    #   - on interdit a p2 d'atterrir dans la fenêtre trop près avant
    #
    for p1 in planes, p2 in planes
        if !(p1.id<p2.id) continue end # évite de traiter le cas symétrique
        # ici p1 peut être avant ou après p2
        for t_i in p1.lb:p1.ub
            t_lb_j = max(p2.lb, t_i - get_sep(sv.inst, p2, p1) + 1) # p2 avant p1
            t_ub_j = min(p2.ub, t_i + get_sep(sv.inst, p1, p2) - 1) # p2 après p1
            for t_j in t_lb_j:t_ub_j
                # Si besoin de nommé les contraintes : c3[p1.id,p2.id,t_i,t_j]
                @constraint(model, y[p1,t_i] + y[p2,t_j] <= 1)
            end
        end
    end
end

# Résoud le problème complet : calcul l'ordre, le timing (dates d'atterrissage)
# et le cout total de la solution optimale au problème SEQATA (coûts linéaires)
# Cette fonction crée le modèle PLNE, le résoud puis met l'objet solution à
# jour.
#
function solve(sv::MipDiscretSolver)
    ln2("BEGIN solve(MipDiscretSolver)")
    ln2( "="^60 )

    lg2("Construction du modèle ($(ms())) ... ")
    buildSeqataModel(sv)
    ln2("fait ($(ms())).")

    lg2("Lancement de la résolution ($(ms())) ... ")
    optimize!(sv.model)
    ln2("fait ($(ms())).")

    lg2("Test de validité de la solution ($(ms())) ... ")
    if  JuMP.termination_status(sv.model) != MOI.OPTIMAL
        print("ERREUR : pas de solution pour :\n    ")
        @show JuMP.termination_status(sv.model)
        println( to_s(sv.bestsol) )
        exit(1)
    end

    lg2("Exploitation des résultats (mise à jour de la solution)($(ms())) ... ")
    # Extraction des valeurs entières des variables d'atterrissage

    # Il reste maintenant à mettre à jour notre objet solution à partir du
    # résultat de la résolution MIP
    for (i, p) in enumerate(sv.inst.planes)
        sv.bestsol.planes[i] = p
        sv.bestsol.x[i] = round(Int, value(sv.x[p]))  # value() retourne Float64 !
        sv.bestsol.costs[i] = value(sv.costs[p])
    end
    sv.bestsol.cost = round(value(sv.cost), digits=Args.args[:cost_precision])

    # On trie juste la solution par date d'atterrissage croissante des avions
    # pour améliorer la présentation de la solution
    sort!(sv.bestsol)
    ln2("fait. ($(ms()))")

    ln2("END solve(MipDiscretSolver)")
end


text:
using TermWin

Profile.init( 1000000, 0.02 )
@profile tshow( TermWin, title="TermWin" )
Profile.print()


text:
"""
    MatHeatDiffModule

Module for linear heat diffusion material models.
"""
module MatHeatDiffModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict

# Type for heat diffusion models of materials.
struct MatHeatDiff
  thermal_conductivity::Array{FFlt, 2};# Thermal conductivity
  specific_heat::FFlt;# Specific heat per unit volume
end

function MatHeatDiff(thermal_conductivity)
    return MatHeatDiff(thermal_conductivity, 0.0)
end

# """
#     bar
#
# Compute.
# """
# function update!(self, ms, gradtheta, Rm; output = nothing, outputRm =nothing)
#
#   kappa = self.property.thermal_conductivity;
#   gradtheta = reshape(context.gradtheta, [],1);# make column vector
#   flux = - kappa * gradtheta;# in local material orientation coordinates
#   if outputRm!= nothing #  output coordinate system supplied?
#     if Rm!=nothing
#       flux = (Rm*flux);# in global coordinate system
#     end
#     flux = outputRm'*flux;# in output coordinate system
#   end
#   if output!= nothing
#     if output=="flux"
#       out = flux;
#     else
#       out = nothing;
#     end
#   else
#     out = flux;
#   end
#   return out, ms
# end

end


text:
#! /usr/bin/julia

# Rosetta Code, Integer overflow

s = subtypes(Signed)
u = subtypes(Unsigned)

println("Integer Type Limits")
for i in 1:length(s)
    println(s[i], " [", typemin(s[i]), ", ", typemax(s[i]), "]")
    println(u[i], " [", typemin(u[i]), ", ", typemax(u[i]), "]")
end

println()
println("Add to typemax")
for t in s
    over = typemax(t) + 1
    println(t, " => ", over, " (", typeof(over), ")")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
#=
all_tests:
- Julia version:
- Author: Berk
- Date: 2020-07-01
Note: All tests should be run from a julia REPR within the OCTHaGOn folder, using:
      julia --project=.
      include("test/all.jl")

      To see coverage, run with:
      julia --project=. --code-coverage=tracefile-%p.info --code-coverage=user
      include("test/all.jl")
=#

include("load.jl");

@testset "OCTHaGOn" begin
    include(string(OCTHaGOn.PROJECT_ROOT, "/test/src.jl"))

    include(string(OCTHaGOn.PROJECT_ROOT, "/test/imports.jl"))

#     include(string(OCTHaGOn.PROJECT_ROOT, "/test/cbf.jl"))

    include(string(OCTHaGOn.PROJECT_ROOT, "/test/algorithms.jl"))

#     include(string(OCTHaGOn.PROJECT_ROOT, "/test/lse.jl"))

end

# Other tests to try later

# include("/test/test_transonic.jl");


text:
module DEACgaussiantest
using SpecialFunctions
using JLD
using NPZ

export heaviside, genGaussianDSF_HS, genGaussianDSF2_HS, genGaussianDSF, genGaussianISF, genGaussianMoment0,
genGaussianMoment1, genGaussianMoment3, genGaussianMomentn1

#See https://journals.aps.org/prb/pdf/10.1103/PhysRevB.82.174510
heaviside(x::AbstractFloat) = ifelse(x < 0, zero(x), ifelse(x > 0, one(x), oftype(x,0.5)))

function genGaussianDSF_HS(omega,p,mu,alpha)
    heaviside.(omega).*(p*exp.(-((omega.-mu).^2)/2/(alpha^2))/sqrt(2*pi)/alpha)
end

function genGaussianDSF2_HS(o1,o2,p,mu,alpha)
    (heaviside.(o1).*(p*exp.(-((o1.-mu).^2)/2/(alpha^2))/sqrt(2*pi)/alpha) + heaviside.(o2).*(p*exp.(-((o2.-mu).^2)/2/(alpha^2))/sqrt(2*pi)/alpha))/2
end

function genGaussianDSF(omega,p,mu,alpha)
    p .* exp.(-((omega.-mu).^2) ./ 2 ./ (alpha^2)) ./ sqrt(2*pi) ./ alpha
end

function genGaussianISF(tau,p,mu,alpha)
    p .* exp.(tau .* (((alpha^2) .* tau) .- (2*mu)) ./ 2)
end

function genGaussianMoment0(p,mu,alpha)
    p
end

function genGaussianMoment1(p,mu,alpha)
    mu*p
end

function genGaussianMoment3(p,mu,alpha)
    mu * ( (3*(alpha^2)) + mu^2 ) * p
end

function genGaussianMomentn1(p,mu,alpha)
    sqrt(2) * p * dawson(mu/alpha/sqrt(2))/alpha
end

function getMoment(moment_function::Function,p::Float64,mu::Float64,alpha::Float64)
    moment_function(p,mu,alpha)
end

function getMoment(moment_function::Function,p::Array{Float64,1},mu::Float64,alpha::Float64)
    m = 0.0;
    for pp in p
        m += moment_function(pp,mu,alpha);
    end
    m
end

function getMoment(moment_function::Function,p::Float64,mu::Array{Float64},alpha::Float64)
    m = 0.0;
    for u in mu
        m += moment_function(p,u,alpha);
    end
    m
end

function getMoment(moment_function::Function,p::Float64,mu::Float64,alpha::Array{Float64})
    m = 0.0;
    for a in alpha
        m += moment_function(p,mu,a);
    end
    m
end

function getMoment(moment_function::Function,p::Array{Float64},mu::Array{Float64},alpha::Float64)
    size(p) == size(mu) || throw(DimensionMismatch("size of p not equal to size of mu"))
    m = 0.0;
    @inbounds for i = 1:size(p,1)
        pp = p[i]
        u = mu[i]
        m += moment_function(pp,u,alpha);
    end
    m
end

function getMoment(moment_function::Function,p::Array{Float64},mu::Float64,alpha::Array{Float64})
    size(p) == size(alpha) || throw(DimensionMismatch("size of p not equal to size of mu"))
    m = 0.0;
    @inbounds for i = 1:size(p,1)
        pp = p[i]
        a = alpha[i]
        m += moment_function(pp,mu,a);
    end
    m
end

function getMoment(moment_function::Function,p::Float64,mu::Array{Float64},alpha::Array{Float64})
    size(mu) == size(alpha) || throw(DimensionMismatch("size of mu not equal to size of alpha"))
    m = 0.0;
    @inbounds for i = 1:size(mu,1)
        u = mu[i]
        a = alpha[i]
        m += moment_function(p,u,a);
    end
    m
end

function getMoment(moment_function::Function,p::Array{Float64},mu::Array{Float64},alpha::Array{Float64})
    size(p) == size(mu) || throw(DimensionMismatch("size of p not equal to size of mu"))
    size(p) == size(alpha) || throw(DimensionMismatch("size of p not equal to size of alpha"))
    m = 0.0;
    @inbounds for i = 1:size(p,1)
        pp = p[i]
        u = mu[i]
        a = alpha[i]
        m += moment_function(pp,u,a);
    end
    m
end

function printMoments(p,mu,alpha)
    m = getMoment(genGaussianMomentn1,p,mu,alpha);
    println("inverse first moment: $m")
    m = getMoment(genGaussianMoment0,p,mu,alpha);
    println("zeroeth moment: $m")
    m = getMoment(genGaussianMoment1,p,mu,alpha);
    println("first moment: $m")
    m = getMoment(genGaussianMoment3,p,mu,alpha);
    println("third moment: $m")
end

function getISF(tau::Array{Float64},p,mu,alpha)
    getSpectrum(genGaussianISF,tau,p,mu,alpha)
end

function getDSF(omega::Array{Float64},p,mu,alpha)
    getSpectrum(genGaussianDSF,omega,p,mu,alpha)
end

function getSpectrum(spectrum_function::Function,x::Array{Float64},p::Float64,mu::Float64,alpha::Float64)
    spectrum_function(x,p,mu,alpha)
end

function getSpectrum(spectrum_function::Function,x::Array{Float64},p::Array{Float64,1},mu::Float64,alpha::Float64)
    s = zeros(size(x,1));
    for pp in p
        s += spectrum_function(x,pp,mu,alpha);
    end
    s
end

function getSpectrum(spectrum_function::Function,x::Array{Float64},p::Float64,mu::Array{Float64},alpha::Float64)
    s = zeros(size(x,1));
    for u in mu
        s += spectrum_function(x,p,u,alpha);
    end
    s
end

function getSpectrum(spectrum_function::Function,x::Array{Float64},p::Float64,mu::Float64,alpha::Array{Float64})
    s = zeros(size(x,1));
    for a in alpha
        s += spectrum_function(x,p,mu,a);
    end
    s
end

function getSpectrum(spectrum_function::Function,x::Array{Float64},p::Array{Float64},mu::Array{Float64},alpha::Float64)
    size(p) == size(mu) || throw(DimensionMismatch("size of p not equal to size of mu"))
    s = zeros(size(x,1));
    @inbounds for i = 1:size(p,1)
        pp = p[i]
        u = mu[i]
        s += spectrum_function(x,pp,u,alpha);
    end
    s
end

function getSpectrum(spectrum_function::Function,x::Array{Float64},p::Array{Float64},mu::Float64,alpha::Array{Float64})
    size(p) == size(alpha) || throw(DimensionMismatch("size of p not equal to size of mu"))
    s = zeros(size(x,1));
    @inbounds for i = 1:size(p,1)
        pp = p[i]
        a = alpha[i]
        s += spectrum_function(x,pp,mu,a);
    end
    s
end

function getSpectrum(spectrum_function::Function,x::Array{Float64},p::Float64,mu::Array{Float64},alpha::Array{Float64})
    size(mu) == size(alpha) || throw(DimensionMismatch("size of mu not equal to size of alpha"))
    s = zeros(size(x,1));
    @inbounds for i = 1:size(mu,1)
        u = mu[i]
        a = alpha[i]
        s += spectrum_function(x,p,u,a);
    end
    s
end

function getSpectrum(spectrum_function::Function,x::Array{Float64},p::Array{Float64},mu::Array{Float64},alpha::Array{Float64})
    size(p) == size(mu) || throw(DimensionMismatch("size of p not equal to size of mu"))
    size(p) == size(alpha) || throw(DimensionMismatch("size of p not equal to size of alpha"))
    s = zeros(size(x,1));
    @inbounds for i = 1:size(p,1)
        pp = p[i]
        u = mu[i]
        a = alpha[i]
        s += spectrum_function(x,pp,u,a);
    end
    s
end

function get_tau(tau_0::Float64,tau_max::Float64,tau_steps::Int64)
    collect(tau_0:tau_max/(tau_steps - 1):tau_max);
end

function saveDataFile(tau_0::Float64,tau_max::Float64,tau_steps::Int64,p,mu,alpha,isf_error; filename::String="qmcdata.jld", print_moments::Bool=false)
    tau = get_tau(tau_0,tau_max,tau_steps);
    saveDataFile(tau,p,mu,alpha,isf_error,filename=filename,print_moments=print_moments);
end

function saveDataFile(tau::Array{Float64,1},p,mu,alpha,isf_error; filename::String="qmcdata.jld", print_moments::Bool=false)
    isf = getISF(tau,p,mu,alpha);
    if print_moments
        printMoments(p,mu,alpha);
    end
    saveDataFile(tau,isf,isf_error,filename=filename);
end

function saveDataFile(tau::Array{Float64,1},isf::Array{Float64,1},isf_error::Float64; filename::String="qmcdata.jld")
    size(tau) == size(isf) || throw(DimensionMismatch("size of tau not equal to size of isf"))
    isf_error_array = fill(isf_error,size(tau,1));
    saveDataFile(tau,isf,isf_error_array,filename=filename);
end

function saveDataFile(tau::Array{Float64,1},isf::Array{Float64,1},isf_error::Array{Float64,1}; filename="qmcdata.jld")
    size(tau) == size(isf) || throw(DimensionMismatch("size of tau not equal to size of isf"))
    size(tau) == size(isf_error) || throw(DimensionMismatch("size of tau not equal to size of isf"))

    # Get filename extension and save using correct method
    filename_base, filename_extension = splitext(filename);
    if filename_extension == ".jld"
        save(filename,
             "tau",tau,
             "isf",isf,
             "isf_error",isf_error);
    elseif filename_extension == ".npz"
        npz_dict = Dict("tau" => tau,
                        "isf" => isf,
                        "isf_error" => isf_error);
        npzwrite(filename, npz_dict);
    else
        throw(ArgumentError("unrecoginized file extension '$(filename_extension)', please choose from [.jld, .npz]"))
    end
    nothing
end
end


text:
module BridgeSDEInference

using Bridge, StaticArrays, Distributions
using Statistics, Random, LinearAlgebra
using ForwardDiff
using ForwardDiff: value

# fitzHughNagumo.jl
export FitzhughDiffusion, FitzhughDiffusionAux
export regularToAlter, alterToRegular, regularToConjug, conjugToRegular, display

# types.jl
export ImproperPrior, NoChangePt, SimpleChangePt

# mcmc.jl
export mcmc, PartObs, FPT, FPTInfo, ConjugateUpdt, MetropolisHastingsUpdt

# ODE solvers:
export Ralston3, RK4, Tsit5, Vern7

# random_walk.jl
export RandomWalk

# priors.jl
export Priors

# blocking_schedule.jl
export NoBlocking, ChequeredBlocking

# starting_pt.jl
export KnownStartingPt, GsnStartingPt

include("fitzHughNagumo.jl")
include("fitzHughNagumo_conjugateUpdt.jl")

include("types.jl")
include("ralston3.jl")
include("rk4.jl")
include("tsit5.jl")
include("vern7.jl")

include("priors.jl")

include("guid_prop_bridge.jl")
include("random_walk.jl")
include("blocking_schedule.jl")
include("starting_pt.jl")
include("mcmc.jl")
include("path_to_wiener.jl")

end


text:
struct Variable
    data::Array
end


text:
"""
`Maker.FileTarget`

The type created by `file()`. Fields expected to be accessed publicly
include:

- `name::UTF8String`
- `dependencies::Vector{UTF8String}`
- `description::UTF8String`

"""
type FileTarget <: AbstractTarget
    name::UTF8String
    dependencies::Vector{UTF8String}
    description::UTF8String
    action::Function
    timestamp::DateTime
    funhash::UInt64
    isstale::Bool
end

"""
```julia
resolvedependency(s::AbstractString)
```

Return the target registered under name `s`. If no target is registered and a
file of name `s` exists, return a new `FileTarget` for that file.
"""
function resolvedependency(s::AbstractString)
    if haskey(tasks(), utf8(s))
        tasks()[utf8(s)]
    elseif isfile(s)
        FileTarget(s, UTF8String[], "", () -> nothing,
                   Dates.unix2datetime(mtime(s)), 0, false)
    else
        error("no rule for target '$s'")
    end
end

function isstale(t::FileTarget)
    if !isfile(t.name) || t.isstale
        return true
    end
    ds = dependencies(t)
    # length(ds) > 0 && @show maximum([timestamp(d) for d in ds])
    # length(ds) > 0 && @show timestamp(t)
    # length(ds) > 0 && @show timestamp(t) < maximum([timestamp(d) for d in ds])
    isempty(ds) ? false : timestamp(t) < maximum([timestamp(d) for d in ds])
end

timestamp(t::FileTarget) = max(t.timestamp, Dates.unix2datetime(mtime(t.name)))

originaltimestamp(::Type{FileTarget}, name::AbstractString) =
    isfile(name) ? Dates.unix2datetime(mtime(name)) : DateTime()

"""
```julia
file(action::Function, name::AbstractString, dependencies=[])
file(name::AbstractString, dependencies=[])
```
Define and register a `file` target for Maker.jl.

- `action` is the function that operates when the target is
  executed.

- `name` refers to the name of the task or target.

- `dependencies` refers to the name (`AbstractString`) or names
  (`Vector{AbstractString}`) of targets that need to be satisfied
  for this target.

Targets are registered globally.

`file` targets use the name of the file as the name of the target.
File targets use timestamps to determine when targets need to be
updated. File paths are relative to the current working directory.

If the `action` or `dependencies` of a target are redefined, the
target will be marked as stale, and the action will be updated
at the next target check.

See also `make`, `directory`,, `task`, and `variable`. `file`
registers a `FileTarget` type.
"""
file(action::Function, name::AbstractString, dependencies=UTF8String[]) =
    target(FileTarget, name, action, dependencies)
file(name::AbstractString, dependencies=UTF8String[]) =
    target(FileTarget, name, () -> nothing, dependencies)


text:
# The SteinMaternRadialKernel is given by
#
# k(x, y) = (1 + gamma * r) * exp{-gamma * r}
#
# where r = ||x-y||. This is a norm for the order
# 2 Sobolev space H_2. This is the Matern covariance
# associated with v = 3/2.
#
# See http://www.gaussianprocess.org/gpml/chapters/RW4.pdf
# or http://arxiv.org/pdf/1204.6448.pdf for more details.

mutable struct SteinMaternRadialKernel <: SteinMaternWeightedKernel
    gamma::Float64
end

# default is use gamma sqrt{3}
SteinMaternRadialKernel() = SteinMaternRadialKernel(sqrt(3.0))

function Q(ker::SteinMaternRadialKernel, x::Array{Float64,1}, y::Array{Float64,1})
    1.0
end

function gradxQ(ker::SteinMaternRadialKernel, x::Array{Float64,1}, y::Array{Float64,1})
    d = length(x)
    zeros(d)
end

function gradxyQ(ker::SteinMaternRadialKernel, x::Array{Float64,1}, y::Array{Float64,1})
    0.0
end


text:
"""
Materialize an `AggregateVector` into a `DataFrame`. Currently only eager
materialization is allowed, which requires two steps:

1. Materialize the source node into a `DataFrame`.
2. Use `DataFrames.combine` to aggregate this `DataFrame`.
"""
function materialize(node::AggregateVector)
    DataFrames.combine(
        materialize(node.source),
        map(
            e -> [e.input_columns...] => e.vector_form => e.alias,
            node.aggregates,
        )...,
    )
end


text:
const TEST_NAME = "ARS1x128"

include("common.jl")

using RandomNumbers.Random123

r = ARS1x(123)

test_all(r, 100_000_000)


text:
module FixedPointAcceleration
    using LinearAlgebra: cond, pinv
    using GLM: fit, LinearModel
    using Dates: now
    include("0_Structs_And_Enums.jl")
    include("1_MainFunctions.jl")
    export fixed_point, fixed_point_new_input # Main functionality.
    export FunctionEvaluationResult, FixedPointResults # Structs.
    export put_together_without_jumps, execute_function_safely # Auxillery functions. Mainly exported for testing.
end


text:
# Make plot of positive definiteness of Gaussian approximation:
using PyPlot
using ForwardDiff

#function posdef()
ntime = 10000
time = collect(linspace(1e-2,1000,ntime))
#nc = 3
#amp = [2.7518074963750734,-0.3546256365636662,-1.3968344706386124]
#b = [1.5312499407621147,4.753504292338284,1.7967431129812566]
#omega = [0.34148888178746056,2.042571639590954,1.7641684073510928]
nc = 4
#amp = [2.953823582534888,-0.04961614737314587,-1.6900719772140012,-0.21413545794774125]
#b = [1.589510614489265,3.4729345651560872,1.8698492189074047,6.363451205766626]
#omega = [0.3329119766576696,4.5999414581901386,1.7115890219170917,2.205939216106191]

amp = [2.953823582534889, -0.04961614737314587, -1.6900719772140012, -0.21413545794774125]
b =[1.589510614489265, 3.4729345651560872, 1.8698492189074047, 6.363451205766626]
omega =[0.3329119766576696, 4.5999414581901386,1.7115890219170917, 2.205939216106191]

function func(x)
fout=0.0
for k=1:nc
  fout +=amp[k]*exp(-b[k]*sqrt(x))*cos(omega[k]*sqrt(x))
end
#fout=exp(-x*0.5)
return fout
end

y = zeros(ntime)
dy = zeros(ntime)
d2y = zeros(ntime)
for i=1:ntime
  y[i]=func(time[i])
  dy[i]=ForwardDiff.derivative(func,time[i])
  d2y[i]=ForwardDiff.derivative(z -> ForwardDiff.derivative(func,z),time[i])
end
plot(time,y)
plot(time,dy)
plot(time,d2y)

#return
#end

max_stars_repo_name: dfm/ess

max_issues_repo_name: dfm/ess

max_forks_repo_name: dfm/ess