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

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

# ╔═╡ fbde6364-4f30-11eb-1ece-712293996c04
using Colors: RGB

# ╔═╡ 486457d8-4f37-11eb-306c-57d650508136
using Images

# ╔═╡ 69dc67fa-4cff-11eb-331e-25ffdced4323
let
	using Plots

	blue = "#8DC0FF"
	red = "#FFAEA6"

	t = collect(0:10:3040)
	ϵ₁ = randn(length(t))*15     # noise on Dolphin IQ
	ϵ₂ = randn(length(t))*20     # noise on Human IQ

	Y₁ = dolphinsIQ = t/12 + ϵ₁
	Y₂ = humanIQ = t/20 + ϵ₂

	sc = scatter(t,Y₁; label="Dolphins", color=blue,
	  ylabel="IQ (-)", xlabel ="Time (year BC)", legend=:topleft)
	sc = scatter!(sc, t,Y₂; label="Humans", color=red)

	λₛ = 0.01

qb3 = QuestionBlock(;
	title=md"**Question 5: Ridge Regression**",
	description = md"""
Ridge regression can be seen as an extension of ordinary least squares regression,

$\beta X =b\, ,$

where a matrix $\beta$ is sought which minimizes the sum of squared residuals between the model and the observations,

$SSE(\beta) = (y - \beta X)^T (y - \beta X)$

In some cases it is adviceable to add a regularisation term to this objective function,

$SSE(\beta) = (y - \beta X)^T (y - \beta X) + \lambda \left\lVert X \right\rVert^2_2 \, ,$

this is known as ridge regression. The matrix $\beta$ that minimises the objective function can be computed analytically.

$\beta = \left(X^T X + \lambda I \right)^{-1}X^T y$

Let us look at an example. We found some data on the evolution of human and dolphin intelligence.

```julia
using Plots

blue = "#8DC0FF"
red = "#FFAEA6"

t = collect(0:10:3040)
ϵ₁ = randn(length(t))*15     # noise on Dolphin IQ
ϵ₂ = randn(length(t))*20     # noise on Human IQ

Y₁ = dolphinsIQ = t/12 + ϵ₁
Y₂ = humanIQ = t/20 + ϵ₂

sc = scatter(t,Y₁; label="Dolphins", color=blue,
  ylabel="IQ (-)", xlabel ="Time (year BC)", legend=:topleft)
sc = scatter!(sc, t,Y₂; label="Humans", color=red)
```

$(sc)

> "For instance, on the planet Earth, man had always assumed that he was more intelligent than dolphins because he had achieved so much - the wheel, New York, wars and so on - whilst all the dolphins had ever  done was muck about in the water having a good time. But conversely, the dolphins had always believed that they were far more intelligent than man - for precisely the same reasons."
>
> ~ *Hitchhikers guide to the galaxy* ~

Plot the trend of human vs. dolphin intelligence by implementing the analytical solution for ridge regression. For this, you need the uniform scaling operator `I`, found in the `LinearAlgebra` package. Use a lambda of 0.01.

""",
	questions = [Question(;description=md"",validators=Bool[], status=md"")]
)
end

# ╔═╡ 7308bc54-e6cd-11ea-0eab-83f7535edf25
# edit the code below to set your name and UGent username

student = (name = "Jimmy Janssen", email = "Jimmy.Janssen@UGent.be");

# press the ▶ button in the bottom right of this cell to run your edits
# or use Shift+Enter

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

# ╔═╡ cdff6730-e785-11ea-2546-4969521b33a7
begin
	using Pkg; Pkg.activate("..")
	using DSJulia;
	using PlutoUI
	tracker = ProgressTracker(student.name, student.email);
	md"""

	Submission by: **_$(student.name)_**
	"""
end

# ╔═╡ a2181260-e6cd-11ea-2a69-8d9d31d1ef0e
md"""
# Notebook 2: Collections
"""

# ╔═╡ 2222fe0c-4c1d-11eb-1e63-f1dbc90a813c
fyi(md"""In programming, a collection is a class used to represent a set of similar data type items as a single unit. These unit classes are used for grouping and managing related objects. A collection has an underlying data structure that is used for efficient data manipulation and storage.

	[Techopedia.com](https://www.techopedia.com/definition/25317/collection)""")

# ╔═╡ 44542690-4c1d-11eb-2eea-49f28ed7fd90
md"""
## 1. Arrays

Let us start with `Array`'s. It is very similar to lists in Python, though it can have more than one dimension. An `Array` is defined as follows,

"""

# ╔═╡ 8f8c7b44-4c1d-11eb-0cd8-3bb82c75c086
A = []            # empty array

# ╔═╡ b0420f36-5a71-11eb-01f1-f16b115f5895
md"Note that the element type is `Any`, since the compiler has no idea what we would store in `A`!"

# ╔═╡ a0c22de6-4c1d-11eb-34a2-aff57cfd22a1
X = [1, 3, -5, 7] # array of integers

# ╔═╡ a81d1f22-4c1d-11eb-1b76-2929f30565bf
md"### Indexing and slicing"

# ╔═╡ bb79bf28-4c1d-11eb-35bf-379ac0cd16b6
md"Let's start by eating the frog. Julia uses [1-based indexing](https://docs.google.com/document/d/11ZKaR0a6hvc6xmYLfmslAAPnkVRSZFGz5GZYRNmxmsQ/edit)..."

# ╔═╡ eaf59a7e-4c1d-11eb-2db3-fd7f995db3e4
names = ["Arthur", "Ford", "Zaphod", "Marvin", "Trillian", "Eddie"]

# ╔═╡ efa60284-4c1d-11eb-1c08-09993363e4a8
names[0]        # this does not work, sorry Pythonista's

# ╔═╡ efa6f180-4c1d-11eb-1ab0-3d1ca0b4bc57
names[1]        # hell yeah!

# ╔═╡ efb6718c-4c1d-11eb-0dff-e55e6a676e39
names[end]      # last element

# ╔═╡ efb69e46-4c1d-11eb-1ce7-ed428db8ff44
names[end-1]    # second to last element

# ╔═╡ 0a02b730-4c1e-11eb-1e8a-872dcfc8ab81
md"Slicing arrays is intuitive,"

# ╔═╡ 15035dee-4c1e-11eb-123f-a961fdd48445
names[3:6]

# ╔═╡ 1960de28-4c1e-11eb-1c84-ffe0cbaac940
md"and slicing with assignment too."

# ╔═╡ 4db9d648-4c1e-11eb-1063-e78c78ef5c4b
names[end-1:end] = ["Slartibartfast", "The Whale and the Bowl of Petunias"]

# ╔═╡ 4fb1a53e-4c1e-11eb-381b-1f86a5ed97a1
names

# ╔═╡ 619fd3bc-5a72-11eb-22f7-49ff5486e0fa
md"What the first index again? And how to access the last element? If you forgot, just use `first` and `last`:"

# ╔═╡ 7e460c3e-5a72-11eb-0c52-cf9583c70759
first(names), last(names)

# ╔═╡ 354fa70c-5a74-11eb-31fc-ad21a845d3b0
fyi(md"Pluto might not execute the lines of code in order. As you are changing `names` later on, the results might be a bit unexpected.")

# ╔═╡ 56f5f21e-4c1e-11eb-004e-f19aa9029b01
md"### Types and arrays"

# ╔═╡ 79e0f212-4c1e-11eb-0d64-87308d762180
md"Julia arrays can be of mixed type."

# ╔═╡ 86dcfb26-4c1e-11eb-347f-ffbd8b396f09
Y = [42, "Universe", []]

# ╔═╡ 8a0e3e36-4c1e-11eb-0ec0-d19fdc3c89d8
md"The type of the array changes depending on the elements that make up the array. With `Any` being the most general type. Working with such general objects is not efficient though. In practice, your code will become as slow as Python."

# ╔═╡ 90d3dc80-4c1e-11eb-2a11-3fe581f0b5f7
typeof(Y)

# ╔═╡ 996ba666-4c1e-11eb-3c5c-4bf8673de6bc
X

# ╔═╡ 90d369a8-4c1e-11eb-3c16-c5fb02bdf3bb
typeof(X)

# ╔═╡ f247382c-4c1e-11eb-229c-efe48b7a4d7f
md"When the elements of the arrays are mixed, the type is promoted to the closest common ancestor. For `Y` this is `Any`.  But an array of an integer and a float becomes an..."

# ╔═╡ fd06f130-4c1e-11eb-37cf-03af9372ae45
B = [1.1, 1]

# ╔═╡ 00237b9a-4c1f-11eb-3b19-73c0b8e4cbed
typeof(B)

# ╔═╡ 0425e980-4c1f-11eb-2477-e35a924b8018
eltype(B)   # gives the type of the elements

# ╔═╡ 069cccc4-4c1f-11eb-39b3-b94136c1b468
md"... array of floats."

# ╔═╡ 0fd48034-4c1f-11eb-06a9-0d7353b2a0d6
md"""
Julia allows the flexibility of having mixed types, though this will hinder performance, as the compiler can no longer optimize for the type. If you process an array of `Any`'s, your code will be as slow as Python.

To create an array of a particular type, just use `Type[]`.

"""

# ╔═╡ 2262e4fc-4c1f-11eb-07b8-0b9732b93d86
Float64[1, 2, 3]

# ╔═╡ 3d0107c4-4c1f-11eb-1b5b-ed954348d0aa
md"""### Initialisation
Arrays can be initialised in all the classic, very Pythonesque ways.

"""

# ╔═╡ ad79420c-64a2-11eb-0ab3-4dfce430f6c3
begin

qb_ai = QuestionBlock(;
	title=md"**Assignment: array intialisation**",
	description = md"""

	Test the following statements.

	```julia
	C = []  # empty

	zeros(5)      # row vector of 5 zeroes

	ones(3,3)     # 3×3 matrix of 1's, will be discussed later on

	fill(0.5, 10) # in case you want to fill a matrix with a specific value

	rand(2)       # row vector of 2 random floats [0,1]

	randn(2)      # same but normally-distributed random numbers

	rand(Bool, 10)  # specify the type, many packages overload rand for other stuff

	```
	""")

end

# ╔═╡ 4b3317da-4c1f-11eb-19d5-03570c4d65df

# ╔═╡ 503c9da0-4c1f-11eb-292a-db7b8ce9f458

# ╔═╡ 503d455c-4c1f-11eb-3af2-8f200db1fd30

# ╔═╡ 504d8aca-4c1f-11eb-3600-d77038b0f2bc

# ╔═╡ 505874a8-4c1f-11eb-1132-3bbba81ae1db

# ╔═╡ 5071c430-4c1f-11eb-226b-634abae6082f

# ╔═╡ c2ccb916-5a72-11eb-16d9-15283727d6cf

# ╔═╡ 52a8a6ec-4c1f-11eb-386c-a99ef05b41b0
md"Often it is better to provide a specific type for initialization. For numeric elements `Float64` is the default."

# ╔═╡ 507254b0-4c1f-11eb-2b2c-8bc88e58e0b3
zeros(Int8, 5)

# ╔═╡ 27b18866-4c22-11eb-22da-656ca8a4c01d
md"### Comprehensions and list-like operations"

# ╔═╡ 387c3e70-4c22-11eb-37e4-bb6c36600074
md"Comprehensions are a concise and powerful way to construct arrays and are very loved by the Python community."

# ╔═╡ 3da23eb8-4c22-11eb-1ec4-c350d615322f
Z = [1, 2, 3, 4, 5, 6, 8, 9, 8, 6, 5, 4, 3, 2, 1]

# ╔═╡ 3f5436d2-4c22-11eb-342a-35b7a29ef146
dZ = [Z[i-1] - 2Z[i] + Z[i+1] for i=2:length(Z)-1] # central difference

# ╔═╡ 63fea2f6-4c22-11eb-0802-37fd7653cdb5
md"""General $N$-dimensional arrays can be constructed using the following syntax:

```julia
[ F(x,y,...) for x=rx, y=ry, ... ]
```

Note that this is similar to using set notation. For example:

```julia
[i * j for i in 1:4, j in 1:5]
```

"""

# ╔═╡ 10267fee-5a73-11eb-2947-279f6be1a3fe
md"We can also use logic to specify which elements to include."

# ╔═╡ 27051614-5a73-11eb-1d22-35ec8ebc1fd8
[i^2 for i in 1:20 if isodd(i)]  # squares of the odd integers below 20

# ╔═╡ a5f17ccc-4c22-11eb-2cb8-7b130e1e811f
md"Arrays behave like a stack. So pushing, appending and popping are valid operations. Elements can be added to the back of the array, (*Pushing and appending*)"

# ╔═╡ f1b481f4-4c22-11eb-39b7-39ffdd5bbccc
push!(names, "Eddie") # add a single element

# ╔═╡ f61fc4b0-4c22-11eb-30b3-154ed1aa43bd
append!(names, ["Sam", "Gerard"]) # add an array

# ╔═╡ fcbeda22-4c22-11eb-2d35-a356b98bbc46
md""" "Eddie" was appended as the final element of the Array along with "Sam" and "Gerard". Remember, a "!" is used to indicate an in-place function. `pop()` is used to return and remove the final element of an array"""

# ╔═╡ 3516a722-4c23-11eb-3ee7-fb8d582c8ce0
md"Removing the last element or first element (*popping*)"

# ╔═╡ 08bb725e-4c23-11eb-3338-03370f49dd11
pop!(names)

# ╔═╡ 0cfc84ca-4c23-11eb-124b-5397430fd203
names

# ╔═╡ 4fbecdfe-4c23-11eb-0da7-5945a49c3a2a
popfirst!(names)

# ╔═╡ 562b751e-4c23-11eb-2b8f-73f710bf3520
names

# ╔═╡ 5ed7284a-4c23-11eb-1451-0ff763f52bc7
md"## 2. Matrices"

# ╔═╡ 0186eab2-4c24-11eb-0ff6-d7f8af343647
md"Let's add a dimension and go to 2D Arrays, matrices. It is all quite straightforward,"

# ╔═╡ 097e96e8-4c24-11eb-24c4-31f4d23d3238
P = [0 1 1; 2 3 5; 8 13 21; 34 55 89]

# ╔═╡ 028be3bc-661f-11eb-251e-73abd3abb9fe
size(P)

# ╔═╡ 08e5589e-661f-11eb-262a-dd917f77f56b
size(P, 1)  # size first dimention

# ╔═╡ 0b9bad58-4c24-11eb-26a8-1d04d7b2be61
P[3,2]  # indexing

# ╔═╡ 0b9d0bf8-4c24-11eb-2beb-0763c66e6a20
P[1,:]  # slicing

# ╔═╡ b9a9a730-5a73-11eb-0d17-7bf1aa935697
md"You  will likely not be surprised that most functions just work for two (or more) dimensions:"

# ╔═╡ d4405d8c-5a73-11eb-21f7-37c4d7ac537b
ones(4, 5)

# ╔═╡ da5edf54-5a73-11eb-26ff-2f6af8adceed
fill(0.5, 8, 2)

# ╔═╡ def334c0-5a73-11eb-3006-e7437155cdef
rand(Bool, 3, 3)

# ╔═╡ 25d53e6c-4c24-11eb-02a2-a71d4b2a7974
md"It is important to know that arrays and other collections are copied by reference."

# ╔═╡ 6191c72e-4c24-11eb-21bb-a59e880a3573
@terminal let
	println(); println()
    P = [0 1 1; 2 3 5; 8 13 21; 34 55 89]
	@show R = P
	@show P
	@show P[1, 1] = 42
	println()
	@show R
end

# ╔═╡ 07220d0a-4c4a-11eb-0ae3-298cf03a0bf6
md"`deepcopy()` can be used to make a fully dereferenced object."

# ╔═╡ 26f6f852-4c4a-11eb-3a5c-e7d788713ab8
md"""### Concatenation
Arrays can be constructed and also concatenated using the following functions,
"""

# ╔═╡ 54b81ed8-4c4a-11eb-1d47-99d5823f2ab1
I = [0 1 1; 2 3 5; 8 13 21; 34 55 89]

# ╔═╡ 56fd77e2-4c4a-11eb-1ab1-4793cd9b220c
W = rand(4,3)

# ╔═╡ 56ff198a-4c4a-11eb-1604-8f08c9cf868c
cat(I, W, dims=2)                # concatenation along a specified dimension

# ╔═╡ 57226e44-4c4a-11eb-26fd-fbd6f993bb72
cat(I, W, dims=1) == vcat(I, W) == [I; W]   # vertical concatenation

# ╔═╡ 57327ab4-4c4a-11eb-219f-f70dd02f170c
cat(I, W, dims=2) == hcat(I, W) == [I W]   # horizontal concatenation

# ╔═╡ 969bc7a0-4c4a-11eb-3db6-892f68020468
md"Note that `;` is an operator to use `vcat`, e.g."

# ╔═╡ a8eba678-4c4a-11eb-2866-1135e65bc4df
[zeros(2, 2) ones(2, 1); ones(1, 3)]

# ╔═╡ afd0b9a6-4c4a-11eb-270b-133ddc3e753b
md"This simplified syntax can lead to strange behaviour. Explain the following difference."

# ╔═╡ b4de01c4-4c4a-11eb-3b0b-275ec9ddf5bf
[1  2-3]

# ╔═╡ b76b87e0-4c4a-11eb-3b21-f1365960fdd0
[1 2 -3]

# ╔═╡ be174110-4c4a-11eb-05cf-17e2527dfad8
md"Sometimes, `vcat` and `hcat` are better used to make the code unambiguous..."

# ╔═╡ cad05c2a-4c4a-11eb-0c89-13cddd2aa35f
md"By default, the `*` operator is used for matrix multiplication"

# ╔═╡ d43ae2ee-4c4a-11eb-281a-b353dc1de640
E = [2 4 3; 3 1 5]

# ╔═╡ d5f77750-4c4a-11eb-0624-11ca5cc2a84e
R = [ 3 10; 4 1 ;7 1]

# ╔═╡ d5fa6276-4c4a-11eb-068e-b114e33e5d8f
E * R

# ╔═╡ 45e86286-4c4b-11eb-1516-0140dc69ab58
md"### Element-wise operations"

# ╔═╡ ee779c1a-4c4a-11eb-1894-d743aeff7f44
md"""This is the Julian way since functions act on the objects, and element-wise operations are done with "dot" operations. For every function or binary operation like `^` there is a "dot" operation `.^` to perform element-by-element exponentiation on arrays."""

# ╔═╡ 0697987c-4c4b-11eb-3052-df54b72dec52
T = [10 10 10; 20 20 20]

# ╔═╡ 14bceb32-4c4b-11eb-06b2-5190f7ebb9c2
T.^2

# ╔═╡ 1ace7720-4c4b-11eb-0338-37e7a7227a68
md"""Under the hood, Julia is looping over the elements of `Y`. So a sequence of dot-operations is fused into a single loop."""

# ╔═╡ 2146ac4c-4c4b-11eb-288f-edb3eacff0eb
T.^2 .+ cos.(T)

# ╔═╡ 28f5c018-4c4b-11eb-3530-8b592f2abeda
md"""
Did you notice that dot-operations are also applicable to functions, even user-defined functions? As programmers, we are by lazy by definition and all these dots are a lot of work. The `@.` macro does this for us.
"""

# ╔═╡ 32351fb8-4c4b-11eb-058b-5bb348e8dfb7
T.^2 .+ cos.(T) == @. T^2 + cos(T)

# ╔═╡ be557eda-64a3-11eb-1562-35ad48531ebd

# ╔═╡ eed4faca-4c1f-11eb-3e6c-b342b48080eb
md""" ### Intermezzo: Colors.jl and Images.jl

As has been mentioned before, everything has a type. We also know that functions can behave differently for each type. With this in mind, let us look at two interesting packages. *Colors.jl* and *Images.jl*
"""

# ╔═╡ 42254aa6-4f37-11eb-001b-f78d5383e36f
RGB(0.5,0.2,0.1)

# ╔═╡ 0977a54e-4f31-11eb-148c-1d44be4f6853
fyi(md"*Colors.jl* provides a wide array of functions for dealing with color. This includes conversion between colorspaces, measuring distance between colors, simulating color blindness, parsing colors, and generating color scales for graphics.")

# ╔═╡ 4c21ed0e-4f37-11eb-0a90-3120e1ee7936
url = "https://i.imgur.com/BJWoNPg.jpg"

# ╔═╡ 5686c59e-4f37-11eb-21d5-47bdbcf75805
download(url, "bluebird.jpg") # download to a local file

# ╔═╡ a5dc2904-4f37-11eb-24c1-d7837c8bd487
img = load("bluebird.jpg")

# ╔═╡ cf3ec91e-4f37-11eb-0706-e9f726532654
typeof(img)

# ╔═╡ d7d1b940-4f37-11eb-3f74-9dbe02774e54
eltype(img)

# ╔═╡ f90725be-4f37-11eb-1905-e99a65ad3e07
size(img)

# ╔═╡ dd781738-4f37-11eb-09fb-b7f3390a49b2
md"So an image is basically a two-dimensional array of Colors. Which means it can be processed just like any other array."

# ╔═╡ ab395348-4f39-11eb-0a3c-1d8af3b6442e
sqr_img = img[1:1500, 201:1700]

# ╔═╡ 7acaf008-4f3a-11eb-17c2-dfaa1e9f2ce7
md"Because of this type system, a lot of interesting feature work out of the box."

# ╔═╡ 8ce0ab98-4f3a-11eb-37b2-dd7dda63ad5f
md"""
☼
$(@bind brightness html"<input type=range min=0.0 max=300.0>")
☾
"""

# ╔═╡ ac62d6e0-5a74-11eb-1538-09d157738257
brightness

# ╔═╡ d73eba40-4f3a-11eb-0aa8-617fc22d5ca3
img[1:1500, 201:1700] ./ (brightness/100)

# ╔═╡ 21db9766-64a4-11eb-3ec1-4956431e7a09

# ╔═╡ 5064c592-4c4b-11eb-0dee-5186caf2b1f6
md"### Higher dimensional arrays"

# ╔═╡ 598980b8-4c4b-11eb-0c5b-b7064b189e97
md"Matrices can be generalized to multiple dimensions."

# ╔═╡ 5fcfb5dc-4c4b-11eb-0be6-e7f66ea1839e
H = rand(3, 3, 3)

# ╔═╡ bf6b9fc4-5a74-11eb-2676-bda580c65877
md"That is all there is to see about matrices, feel free to create arrays of any dimension."

# ╔═╡ 6e7d5a94-4c4b-11eb-3e2d-353177d6bca5
md"### Ranges"

# ╔═╡ 79129c92-4c4b-11eb-28f6-633aedabd990
md"The colon operator `:` can be used to construct unit ranges, e.g., from 1 to 20:"

# ╔═╡ 7f693b34-4c4b-11eb-134d-af855593f45e
ur = 1:20

# ╔═╡ 83939d80-4c4b-11eb-3dc8-1b559c61a43b
md"Or by increasing in steps:"

# ╔═╡ 884b0d4a-4c4b-11eb-1d39-5f67ba9e92fe
str = 1:3:20

# ╔═╡ 9483861e-4c4b-11eb-156b-2501ef2c54d0
md"Similar to the `range` function in Python, the object that is created is not an array, but an iterator. This is actually the term used in Python. Julia has many different types and structs, which behave a particular way. Types of `UnitRange` only store the beginning and end value (and stepsize in the case of `StepRange`). But functions are overloaded such that it acts as an array. This can really improve the execution speed since the conversion from `Range` to explicit array is only performed where it is necessary and avoids copying large matrices from function to function."

# ╔═╡ 9fd1be0a-4c4b-11eb-299b-f7f0d8797f71
@terminal let
	for i in ur
	  println(i)
	end
end

# ╔═╡ a2104d08-4c4b-11eb-0ccc-b588e99a2057
str[3]

# ╔═╡ a214405c-4c4b-11eb-118b-b7ec852f9257
length(str)

# ╔═╡ a806a9dc-4c4b-11eb-1101-6f75be7a610c
md"All values can be obtained using `collect`:"

# ╔═╡ af514422-4c4b-11eb-1cfe-cba6029ec52f
collect(str)

# ╔═╡ b3035158-4c4b-11eb-1b8d-1fc4070fa132
md"Such implicit objects can be processed much smarter than naive structures. Compare! (You might need to run the lines of code again as not to measure compile time."

# ╔═╡ c4575438-4c4b-11eb-1da2-97acca3f3e99
@elapsed sum((i for i in 1:100_000_000))

# ╔═╡ 0a26ad38-5a75-11eb-0424-4d26aa905de5
md"Note that we use round brackets to create our own iterator, where every element is processed one by one by `sum`."

# ╔═╡ c64e66c8-4c4b-11eb-1686-8fa64d8b2505
@elapsed sum(1:100_000_000)

# ╔═╡ 3a6c466a-5a75-11eb-07e2-ffbf9ec3ffe4
sum(1:100_000_000) == sum((i for i in 1:100_000_000))

# ╔═╡ 03d82c7a-4c58-11eb-0071-bb9ea16bfbb3
md"`StepRange` and `UnitRange` also work with floats."

# ╔═╡ 0bec2d28-4c58-11eb-0a51-95bf50bbfd79
0:0.1:10

# ╔═╡ 6aa73154-661f-11eb-2b88-578eb2dd2ec2
(0:1:100) / 10  # equivalent

# ╔═╡ ae6064e6-64a4-11eb-24b5-0b0b848aa2d6

# ╔═╡ 0fd08728-4c58-11eb-1b71-c9710d398fab
md"## 3. Other collections"

# ╔═╡ 2c6097f4-4c58-11eb-0807-d5d8cbfbd62c
md"Some of the other collections include tuples, dictionaries, and others."

# ╔═╡ 9505a4d4-4c58-11eb-1e2e-0d080437fa23
md"### Tuples"

# ╔═╡ 23dac6cc-4c58-11eb-2c66-f1f79db08536
tupleware = ("tuple", "ware") # tuples

# ╔═╡ 4f80f0a8-4c58-11eb-3679-c186c61c5a14
fyi(md"A tuple is an array with a fixed size. It is not possible to perform operations that change the size of a tuple.")

# ╔═╡ 6d496b44-4c58-11eb-33b6-5b4d6315b6ea
pop!(tupleware)

# ╔═╡ 56e8f6b4-5a75-11eb-3eeb-ffec491be69c
fyi(md"In contrast to arrays however, the types at positions should not be the same, since the compiler will create a new type for every combination!")

# ╔═╡ 8245e46e-5a75-11eb-2d0a-27ef6a1f2492
mixedtuple = (9, "nine")

# ╔═╡ 942b88b4-5a75-11eb-3e7b-4534bf4a7b12
typeof(mixedtuple)

# ╔═╡ 7bc7bdf4-4c58-11eb-1fd8-376ac6da5ab2
md"Indexing and slicing is the same as arrays,"

# ╔═╡ 74d97654-4c58-11eb-344b-8d6df24323d5
tupleware[end]

# ╔═╡ 9bb1e83a-5a75-11eb-0fcc-59f6cc50bf6a
md"One can also define named tuples where the fields have names."

# ╔═╡ bcd5696a-5a75-11eb-0ec1-f116216aa682
namedtuple = (trainer="Ash", pokemon="Pikachu")

# ╔═╡ d74ff724-5a75-11eb-39a1-a963cd64948b
namedtuple[:trainer]

# ╔═╡ dd8978a4-5a75-11eb-287e-03e4272b6f2c
namedtuple[1]  # but indexing still works

# ╔═╡ 9cecd9a6-4c58-11eb-22dc-33cd2559d815
md"### Dictionaries"

# ╔═╡ c6d236da-4c58-11eb-2714-3f5c43583a3d
fyi(md"A dictionary is a collection that stores a set of values with their corresponding keys internally for faster data retrieval. The operation of finding the value associated with a key is called lookup or indexing. [Techopedia.com](https://www.techopedia.com/definition/10263/dictionary-c)")

# ╔═╡ 3253ab74-4c58-11eb-178e-83ea8aba9c8f
scores = Dict("humans" => 2, "vogons" => 1) # dictionaries

# ╔═╡ 32593936-4c58-11eb-174c-0bb20d93dde5
scores["humans"]

# ╔═╡ fcfb511c-5a75-11eb-2181-33d147ab1806
scores["mice"] = 3  # adding a key

# ╔═╡ 1882840a-5a76-11eb-3392-81c2915487f5
scores

# ╔═╡ 3fc787d6-5a76-11eb-06e9-5378d27ce011
delete!(scores, "humans")  # removing a key, earth was destroyed

# ╔═╡ ebb09172-4c58-11eb-1cc9-91193c57677d
md"## 4. Exercises"

# ╔═╡ ee9069e2-63a7-11eb-12b9-97ae270506f4
hint(md"Remember, `.` is not only used for decimals...")

# ╔═╡ 3de1f1aa-58bd-11eb-2ffc-0de292b13840
function riemannsum(f, a, b; n=100)
  return missing
end

# ╔═╡ 5f47cdf0-58be-11eb-1bca-a3d0941b9bea
begin
	integral1 = missing #...
	integral2 = missing #...
end;

# ╔═╡ 3aa37510-58bb-11eb-2ecb-37ce4428269c
begin

	q1 = Question(
			validators = @safe[
				norm(Solutions.riemannsum(sin, 0, 2pi) - riemannsum(sin, 0, 2pi)) ≤ abs(0.1 * Solutions.riemannsum(sin, 0, 2pi)),

				(Solutions.riemannsum(x->x*sin(x), 0, 2pi) - riemannsum(x->x*sin(x), 0, 2pi)) ≤ abs(0.1 * Solutions.riemannsum(x->x*sin(x), 0, 2pi)),

				(Solutions.riemannsum(x->(sqrt(1-x^2)), 0, 1, n=1000) - riemannsum(x->(sqrt(1-x^2)), 0, 1, n=1000)) ≤ abs(0.1 * Solutions.riemannsum(x->(sqrt(1-x^2)), 0, 1, n=1000))
			]
		)

	q2 = Question(
			validators = @safe[
				integral1 == Solutions.riemannsum(x->x*sin(x), 0, 2pi)
			],
			description = md" **Integral 1:**  $\int_0^{2\pi} x\,\sin(x)\,dx$ (n=100)")

	q3 = Question(
			validators = @safe[
			integral2 == Solutions.riemannsum(x->(sqrt(1-x^2)), 0, 1, n=1000)
			],
			description = md" **Integral 2:**  $\int_0^1 \sqrt{1 - x^2}\,dx$ (n=1000)")

qb1 = QuestionBlock(;
	title=md"**Question 1: do you still remember how to integrate?**",
	description = md"""
Integrating for dummies. Compute the Riemann sum **without** making use of a for-loop.

Riemann approximates the integration of a function in the interval [a, b],

$$\int_a^b f(x)\, dx \approx \sum_{i=1}^n f(x_i) \,\Delta x$$

which is the sum of the function $f(x)$ evaluated over an array of x-values in the interval [a,b] multiplied by the $\Delta x$ which is,

$$\Delta x = \cfrac{(b-a)}{n}$$

Complete the function `riemannsum(f, a, b,; n=100)` where the arguments are the function to integrate (f) the boundaries of the interval a, b and the number of bins with a default value of 100, n.

	""",
	questions = [q1, q2, q3]
)
	validate(qb1, tracker)
end

# ╔═╡ c1e377c4-64a4-11eb-3e7f-b163cb465057

# ╔═╡ 5619fd6c-4cfe-11eb-1512-e1800b6c7df9
function mydet(A)
	size(A,1) != size(A,2) && throw(DimensionMismatch)
	return missing
end

# ╔═╡ b1a00da4-4cfe-11eb-0aff-69099e40d28f
let
	using LinearAlgebra
	M₁ = [1 2; 3 4]
	M₂ = rand(10, 10)

q2 = Question(;
	description=md"""
Write a function `mydet` to compute the determinant of a 2x2 square matrix. Remember, for a $2 \times 2$ matrix, the determinant is computed as

${\displaystyle|A|={\begin{vmatrix}a&b\\c&d\end{vmatrix}}=ad-bc.}$
""",
	validators = @safe[det(M₁) == mydet(M₁)]
)

q3 = QuestionOptional{Hard}(;
	description=md"""For larger matrices, there is a recursive way of computing the determinant based on the minors, i.e. the determinants of the submatrices. See [http://mathworld.wolfram.com/Determinant.html](http://mathworld.wolfram.com/Determinant.html).

Update `mydet` to compute the determinant of a general square matrix.
""",
	validators = @safe[det(M₁) ≈ mydet(M₁), det(M₂) ≈ mydet(M₂)]
)

qb2 = QuestionBlock(;
	title=md"**Question 2: determinining the determinant**",
	description = md"""
	""",
	questions = [q2, q3]
)
	validate(qb2, tracker)
end

# ╔═╡ e5293248-64a4-11eb-0d30-53a15bec0d01

# ╔═╡ cb20fffe-58cf-11eb-1b65-49699f2d3699
function estimatepi(n)
	missing
end

# ╔═╡ c6e16d7a-58cf-11eb-32a4-3372939066e3
begin
q91 = Question(
		  validators = @safe[
			abs(estimatepi(100) - π) < 1.0,
			abs(estimatepi(100000) - π) < 1e-2
		  ],
		)

q92 = QuestionOptional{Easy}(
		validators = @safe[],
		description = md"Did you use a for loop? If so, try to do this without an explicit for-loop")

qb90 = QuestionBlock(;
	title=md"**Question 3: it is pi 'o clock**",
	description = md"""
	Estimate pi through Monte Carlo sampling. Do this by simulating throwing `n` pebbles in the [-1, 1] x [-1, 1] square and track the fraction that land in the unit circle. Complete the function `estimatepi` below.
	""",
	questions = [q91, q92],
	hints=[
		hint(md"""
			[Check this image](http://www.pythonlikeyoumeanit.com/_images/circle_square_small.png)
			 """),
		hint(md"""
			Because each throw falls randomly within the square, you realize that the probability of a dart landing within the circle is given by the ratio of the circle’s area to the square’s area:

			$$P_{circle} = \frac{Area_{circle}}{Area_{square}} = \frac{\pi r^2}{(2r)^2}$$

			Furthermore, we can interpret Pcircle as being approximated by the fraction of darts thrown that land in the circle. Thus, we find:

			$$\frac{N_{circle}}{N_{total}} \approx \frac{\pi r^2}{(2r)^2} = \frac{\pi}{4}$$

			where $N_{total}$ is the total number of darts thrown, and $N_{circle}$ is the number of darts that land within the circle. Thus simply by keeping tally of where the darts land, you can begin to estimate the value of π!

			[source:](http://www.pythonlikeyoumeanit.com/Module3_IntroducingNumpy/Problems/Approximating_pi.html)  pythonlikeyoumeanit.com
			 """),

		]
)
	validate(qb90, tracker)
end

# ╔═╡ cee388d2-58cf-11eb-3b88-971b4b85e957
function estimatepi2(n)
	missing
end

# ╔═╡ 41b19e20-4d0f-11eb-1c3c-572cc5243d99

# ╔═╡ 04aff640-58bb-11eb-1bb6-69ad9fc32314
md"## 5. Extra exercises"

# ╔═╡ 75d14674-58ba-11eb-3868-172fc00a0eb8
function markdowntable(table, header)
	missing
end

# ╔═╡ 0c91ce30-58b9-11eb-3617-4d87682831dd
begin

	q71 = Question(
			validators = @safe[
			  markdowntable(ones(8, 4), ["A", "B", "C", "D"])==
			  Solutions.markdowntable(ones(8, 4), ["A", "B", "C", "D"]),
			  markdowntable(ones(20, 2), ["Foo", "Bar"])==
			  Solutions.markdowntable(ones(20, 2), ["Foo", "Bar"]),
			]
		  )
	q72 = QuestionOptional{Easy}(
			description = md"Make sure the user does not put in invalid table specifications by asserting the number of columns in the content equal to the number of elements in the header."
		  )

	qb70 = QuestionBlock(;
		title=md"**Question 4: markdown tables**",
		description = md"""
	Markdown is a lightweight markup language that you can use to add formatting elements to plaintext text documents. It is also the markup language used in this notebook. Markdown is really easy to learn (see the example below). The problem with markdown is that tables generation is a tedious process... Write a small julia package (read function) that generates a markdown table that takes a an array of strings for the header and a n-by-m array of table values. Complete `markdowntable()` below. The function should both return a string of the markdown table and should automatically copies this to the clipboard using the `clipboard()` function. Just for completion you should end your table with a newline (\n).

		```MD
			# Header 1
			## Header 2
			### Header 3

			**This text is bold** and *this is italic*.

			* This
			* is
			* a
			* list.

			| This | is | a | table |
			| :--|:--|:--|:--|
			| 5 | 10 | 10 | 3 |
			| 9 | 3 | 1 | 5 |
			| 8 | 4 | 7 | 6 |
		```
	# Header 1
	## Header 2
	### Header 3

	**This text is bold** and *this is italic*.

	* This
	* is
	* a
	* list.

	| This | is | a | table |
	| :--|:--|:--|:--|
	| 5 | 10 | 10 | 3 |
	| 9 | 3 | 1 | 5 |
	| 8 | 4 | 7 | 6 |

	-----------------

		""",
		questions = [q71, q72],
		hints = [
			hint(md""" The `join` and `repeat`-functions might come in handy """),
			hint(md""" The `@assert` macro should get you close to solving the second part."""),
		]

	)

	validate(qb70, tracker)
end

# ╔═╡ 9f1a2834-4d0f-11eb-3c3e-b7ff55f65dd3
begin
	t = collect(0:10:3040)
	ϵ₁ = randn(length(t))*15     # noise on Dolphin IQ
	ϵ₂ = randn(length(t))*20     # noise on Human IQ

	Y₁ = dolphinsIQ = t/12 + ϵ₁
	Y₂ = humanIQ = t/20 + ϵ₂
end;

# ╔═╡ 85fb018e-4c1d-11eb-2519-a5abe100748e
begin
	β₁ = missing    # replace with the correct way to compute β
	β₂ = missing    # replace with the correct way to compute β

	Yₚ₁ = β₁.*t      # Dolphin IQ
	Yₚ₂ = β₂.*t      # Human IQ
end;

# ╔═╡ 00121c4e-64a5-11eb-2993-61c695c4e6a1

# ╔═╡ a8837ec2-5a4b-11eb-2930-55e48850b7db
vandermonde(α, n) = missing

# ╔═╡ b56686ec-4cfa-11eb-2b14-a5d49a137cc5
let

q1 = Question(;
	description=md"""
Write a one-liner function `vandermonde` to generate this matrix. This function takes as a vector `α` and `n`, the number of powers to compute.
""",
	validators = @safe[vandermonde(1:20, 5) == Solutions.vandermonde(1:20, 5)]
)

qb1 = QuestionBlock(;
	title=md"**Question 6: Vandermonde matrix**",
	description = md"""Write a function to generate an $n \times m$ [Vandermonde matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix) for a given vector $\alpha=[\alpha_1,\alpha_2,\ldots,\alpha_m]^T$. This matrix is defined as follows

${\displaystyle V={\begin{bmatrix}1&\alpha _{1}&\alpha _{1}^{2}&\dots &\alpha _{1}^{n-1}\\1&\alpha _{2}&\alpha _{2}^{2}&\dots &\alpha _{2}^{n-1}\\1&\alpha _{3}&\alpha _{3}^{2}&\dots &\alpha _{3}^{n-1}\\\vdots &\vdots &\vdots &\ddots &\vdots \\1&\alpha _{m}&\alpha _{m}^{2}&\dots &\alpha _{m}^{n-1}\end{bmatrix}},}$

or

$V = [\alpha_i^{j-1}] .$

""",
	questions = [q1]
)
	validate(qb1)
end

# ╔═╡ 16ec4ee4-64a5-11eb-26f3-15313b8b5acb

# ╔═╡ a0026d1c-bcdb-4a2e-b1a1-11c5235a4956
md"""## Answers:
If you would like to take a look at the answers, you can do so by checking the boxo the question you would like to see. The function will be shown just below the question you want to look at.

| Question | Show solution |
|-----|:---------:|
| Question 1 | $(@bind answ_q1 CheckBox()) |
| Question 2 | $(@bind answ_q2 CheckBox()) |
| Question 3 | $(@bind answ_q3 CheckBox()) |
| Question 4 | $(@bind answ_q4 CheckBox()) |
| Question 5 | $(@bind answ_q5 CheckBox()) |
| Question 6 | $(@bind answ_q6 CheckBox()) |
"""

# ╔═╡ 6d43af49-f127-4ffe-ba97-0f04fb792efb
if answ_q1 == true
	md"""
	```Julia
	function riemannsum(f, a, b; n=100)
	  dx = (b - a) / n
	  return sum(f.(a:dx:b)) * dx
	end
	```
	"""
end

# ╔═╡ eae5611f-913c-48e4-bc1b-86d33908ac46
if answ_q2 == true
	md"""
	```Julia
	function mydet(A)
		size(A,1) != size(A,2) && throw(DimensionMismatch)
		return A[1,1]*A[2,2]-[1,2]*A[2,1]
	end
	```
	"""
end

# ╔═╡ be4c9bbe-0ec1-4477-b422-dcf308cd6f5b
if answ_q3 == true
	md"""
	```Julia
	function estimatepi(n)
	  hits = 0
	  for i in 1:n
		  x = rand()
		  y = rand()
		  if x^2 + y^2 ≤ 1.0
			  hits += 1
		  end
	  end
	  return 4hits/n
	end

	estimatepi2(n) = 4count(x-> x ≤ 1.0, sum(rand(n, 2).^2, dims=2)) / n
	```
	"""
end

# ╔═╡ cf2b6aa4-b0e8-471a-95fc-646bbddb989a
if answ_q4 == true
	md"""
	```Julia
	function markdowntable(table, header)
	  n, m = size(table)
	  table_string = ""
	  # add header
	  table_string *= "| " * join(header, " | ") * " |\n"
	  # horizontal lines for separating the columns
	  table_string *= "| " * repeat(":--|", m) * "\n"
	  for i in 1:n
		  table_string *= "| " * join(table[i,:], " | ") * " |\n"
	  end
	  #clipboard(table_string)
	  return table_string
	end
	```
	"""
end

# ╔═╡ 65a45d44-7942-4808-bed0-0369077c5edb
if answ_q5 == true
	md"""
	```Julia
	β₁ = inv(transpose(t)*t + UniformScaling(0.01)) * transpose(t)* Y₁
	β₂ = inv(transpose(t)*t + UniformScaling(0.01)) * transpose(t)* Y₂

	Yₚ₁ = β₁.*t      # Dolphin IQ
	Yₚ₂ = β₂.*t      # Human IQ
	```
	"""
end

# ╔═╡ 80159cd5-2e46-4375-b899-a44e45a948da
if answ_q6 == true
	md"""
	```Julia
	vandermonde(α, n) = [αᵢ^j for αᵢ in α, j in 0:n-1]
	```
	"""
end

# ╔═╡ 2e7973b6-4d0f-11eb-107c-cdaf349428c0
md""" ## 5. References

- [Julia Documentation](https://juliadocs.github.io/Julia-Cheat-Sheet/)
- [Introduction to Julia UCI data science initiative](http://ucidatascienceinitiative.github.io/IntroToJulia/)
- [Month of Julia](https://github.com/DataWookie/MonthOfJulia)
- [Why I love Julia, Next Journal](https://nextjournal.com/kolia/why-i-love-julia)

"""

# ╔═╡ Cell order:
# ╟─cdff6730-e785-11ea-2546-4969521b33a7
# ╠═7308bc54-e6cd-11ea-0eab-83f7535edf25
# ╟─a2181260-e6cd-11ea-2a69-8d9d31d1ef0e
# ╟─2222fe0c-4c1d-11eb-1e63-f1dbc90a813c
# ╟─44542690-4c1d-11eb-2eea-49f28ed7fd90
# ╠═8f8c7b44-4c1d-11eb-0cd8-3bb82c75c086
# ╟─b0420f36-5a71-11eb-01f1-f16b115f5895
# ╠═a0c22de6-4c1d-11eb-34a2-aff57cfd22a1
# ╟─a81d1f22-4c1d-11eb-1b76-2929f30565bf
# ╟─bb79bf28-4c1d-11eb-35bf-379ac0cd16b6
# ╠═eaf59a7e-4c1d-11eb-2db3-fd7f995db3e4
# ╠═efa60284-4c1d-11eb-1c08-09993363e4a8
# ╠═efa6f180-4c1d-11eb-1ab0-3d1ca0b4bc57
# ╠═efb6718c-4c1d-11eb-0dff-e55e6a676e39
# ╠═efb69e46-4c1d-11eb-1ce7-ed428db8ff44
# ╟─0a02b730-4c1e-11eb-1e8a-872dcfc8ab81
# ╠═15035dee-4c1e-11eb-123f-a961fdd48445
# ╟─1960de28-4c1e-11eb-1c84-ffe0cbaac940
# ╠═4db9d648-4c1e-11eb-1063-e78c78ef5c4b
# ╠═4fb1a53e-4c1e-11eb-381b-1f86a5ed97a1
# ╟─619fd3bc-5a72-11eb-22f7-49ff5486e0fa
# ╠═7e460c3e-5a72-11eb-0c52-cf9583c70759
# ╟─354fa70c-5a74-11eb-31fc-ad21a845d3b0
# ╟─56f5f21e-4c1e-11eb-004e-f19aa9029b01
# ╟─79e0f212-4c1e-11eb-0d64-87308d762180
# ╠═86dcfb26-4c1e-11eb-347f-ffbd8b396f09
# ╟─8a0e3e36-4c1e-11eb-0ec0-d19fdc3c89d8
# ╠═90d3dc80-4c1e-11eb-2a11-3fe581f0b5f7
# ╠═996ba666-4c1e-11eb-3c5c-4bf8673de6bc
# ╠═90d369a8-4c1e-11eb-3c16-c5fb02bdf3bb
# ╟─f247382c-4c1e-11eb-229c-efe48b7a4d7f
# ╠═fd06f130-4c1e-11eb-37cf-03af9372ae45
# ╠═00237b9a-4c1f-11eb-3b19-73c0b8e4cbed
# ╠═0425e980-4c1f-11eb-2477-e35a924b8018
# ╟─069cccc4-4c1f-11eb-39b3-b94136c1b468
# ╟─0fd48034-4c1f-11eb-06a9-0d7353b2a0d6
# ╠═2262e4fc-4c1f-11eb-07b8-0b9732b93d86
# ╟─3d0107c4-4c1f-11eb-1b5b-ed954348d0aa
# ╟─ad79420c-64a2-11eb-0ab3-4dfce430f6c3
# ╠═4b3317da-4c1f-11eb-19d5-03570c4d65df
# ╠═503c9da0-4c1f-11eb-292a-db7b8ce9f458
# ╠═503d455c-4c1f-11eb-3af2-8f200db1fd30
# ╠═504d8aca-4c1f-11eb-3600-d77038b0f2bc
# ╠═505874a8-4c1f-11eb-1132-3bbba81ae1db
# ╠═5071c430-4c1f-11eb-226b-634abae6082f
# ╠═c2ccb916-5a72-11eb-16d9-15283727d6cf
# ╟─52a8a6ec-4c1f-11eb-386c-a99ef05b41b0
# ╠═507254b0-4c1f-11eb-2b2c-8bc88e58e0b3
# ╟─27b18866-4c22-11eb-22da-656ca8a4c01d
# ╟─387c3e70-4c22-11eb-37e4-bb6c36600074
# ╠═3da23eb8-4c22-11eb-1ec4-c350d615322f
# ╠═3f5436d2-4c22-11eb-342a-35b7a29ef146
# ╟─63fea2f6-4c22-11eb-0802-37fd7653cdb5
# ╟─10267fee-5a73-11eb-2947-279f6be1a3fe
# ╠═27051614-5a73-11eb-1d22-35ec8ebc1fd8
# ╟─a5f17ccc-4c22-11eb-2cb8-7b130e1e811f
# ╠═f1b481f4-4c22-11eb-39b7-39ffdd5bbccc
# ╠═f61fc4b0-4c22-11eb-30b3-154ed1aa43bd
# ╟─fcbeda22-4c22-11eb-2d35-a356b98bbc46
# ╟─3516a722-4c23-11eb-3ee7-fb8d582c8ce0
# ╠═08bb725e-4c23-11eb-3338-03370f49dd11
# ╠═0cfc84ca-4c23-11eb-124b-5397430fd203
# ╠═4fbecdfe-4c23-11eb-0da7-5945a49c3a2a
# ╠═562b751e-4c23-11eb-2b8f-73f710bf3520
# ╟─5ed7284a-4c23-11eb-1451-0ff763f52bc7
# ╟─0186eab2-4c24-11eb-0ff6-d7f8af343647
# ╠═097e96e8-4c24-11eb-24c4-31f4d23d3238
# ╠═028be3bc-661f-11eb-251e-73abd3abb9fe
# ╠═08e5589e-661f-11eb-262a-dd917f77f56b
# ╠═0b9bad58-4c24-11eb-26a8-1d04d7b2be61
# ╠═0b9d0bf8-4c24-11eb-2beb-0763c66e6a20
# ╟─b9a9a730-5a73-11eb-0d17-7bf1aa935697
# ╠═d4405d8c-5a73-11eb-21f7-37c4d7ac537b
# ╠═da5edf54-5a73-11eb-26ff-2f6af8adceed
# ╠═def334c0-5a73-11eb-3006-e7437155cdef
# ╟─25d53e6c-4c24-11eb-02a2-a71d4b2a7974
# ╠═6191c72e-4c24-11eb-21bb-a59e880a3573
# ╟─07220d0a-4c4a-11eb-0ae3-298cf03a0bf6
# ╟─26f6f852-4c4a-11eb-3a5c-e7d788713ab8
# ╠═54b81ed8-4c4a-11eb-1d47-99d5823f2ab1
# ╠═56fd77e2-4c4a-11eb-1ab1-4793cd9b220c
# ╠═56ff198a-4c4a-11eb-1604-8f08c9cf868c
# ╠═57226e44-4c4a-11eb-26fd-fbd6f993bb72
# ╠═57327ab4-4c4a-11eb-219f-f70dd02f170c
# ╟─969bc7a0-4c4a-11eb-3db6-892f68020468
# ╠═a8eba678-4c4a-11eb-2866-1135e65bc4df
# ╟─afd0b9a6-4c4a-11eb-270b-133ddc3e753b
# ╠═b4de01c4-4c4a-11eb-3b0b-275ec9ddf5bf
# ╠═b76b87e0-4c4a-11eb-3b21-f1365960fdd0
# ╟─be174110-4c4a-11eb-05cf-17e2527dfad8
# ╟─cad05c2a-4c4a-11eb-0c89-13cddd2aa35f
# ╠═d43ae2ee-4c4a-11eb-281a-b353dc1de640
# ╠═d5f77750-4c4a-11eb-0624-11ca5cc2a84e
# ╠═d5fa6276-4c4a-11eb-068e-b114e33e5d8f
# ╟─45e86286-4c4b-11eb-1516-0140dc69ab58
# ╟─ee779c1a-4c4a-11eb-1894-d743aeff7f44
# ╠═0697987c-4c4b-11eb-3052-df54b72dec52
# ╠═14bceb32-4c4b-11eb-06b2-5190f7ebb9c2
# ╟─1ace7720-4c4b-11eb-0338-37e7a7227a68
# ╠═2146ac4c-4c4b-11eb-288f-edb3eacff0eb
# ╟─28f5c018-4c4b-11eb-3530-8b592f2abeda
# ╠═32351fb8-4c4b-11eb-058b-5bb348e8dfb7
# ╟─be557eda-64a3-11eb-1562-35ad48531ebd
# ╟─eed4faca-4c1f-11eb-3e6c-b342b48080eb
# ╠═fbde6364-4f30-11eb-1ece-712293996c04
# ╠═42254aa6-4f37-11eb-001b-f78d5383e36f
# ╟─0977a54e-4f31-11eb-148c-1d44be4f6853
# ╠═486457d8-4f37-11eb-306c-57d650508136
# ╠═4c21ed0e-4f37-11eb-0a90-3120e1ee7936
# ╠═5686c59e-4f37-11eb-21d5-47bdbcf75805
# ╠═a5dc2904-4f37-11eb-24c1-d7837c8bd487
# ╠═cf3ec91e-4f37-11eb-0706-e9f726532654
# ╠═d7d1b940-4f37-11eb-3f74-9dbe02774e54
# ╠═f90725be-4f37-11eb-1905-e99a65ad3e07
# ╟─dd781738-4f37-11eb-09fb-b7f3390a49b2
# ╠═ab395348-4f39-11eb-0a3c-1d8af3b6442e
# ╟─7acaf008-4f3a-11eb-17c2-dfaa1e9f2ce7
# ╟─8ce0ab98-4f3a-11eb-37b2-dd7dda63ad5f
# ╠═ac62d6e0-5a74-11eb-1538-09d157738257
# ╠═d73eba40-4f3a-11eb-0aa8-617fc22d5ca3
# ╟─21db9766-64a4-11eb-3ec1-4956431e7a09
# ╟─5064c592-4c4b-11eb-0dee-5186caf2b1f6
# ╟─598980b8-4c4b-11eb-0c5b-b7064b189e97
# ╠═5fcfb5dc-4c4b-11eb-0be6-e7f66ea1839e
# ╟─bf6b9fc4-5a74-11eb-2676-bda580c65877
# ╟─6e7d5a94-4c4b-11eb-3e2d-353177d6bca5
# ╟─79129c92-4c4b-11eb-28f6-633aedabd990
# ╠═7f693b34-4c4b-11eb-134d-af855593f45e
# ╟─83939d80-4c4b-11eb-3dc8-1b559c61a43b
# ╠═884b0d4a-4c4b-11eb-1d39-5f67ba9e92fe
# ╟─9483861e-4c4b-11eb-156b-2501ef2c54d0
# ╠═9fd1be0a-4c4b-11eb-299b-f7f0d8797f71
# ╠═a2104d08-4c4b-11eb-0ccc-b588e99a2057
# ╠═a214405c-4c4b-11eb-118b-b7ec852f9257
# ╟─a806a9dc-4c4b-11eb-1101-6f75be7a610c
# ╠═af514422-4c4b-11eb-1cfe-cba6029ec52f
# ╟─b3035158-4c4b-11eb-1b8d-1fc4070fa132
# ╠═c4575438-4c4b-11eb-1da2-97acca3f3e99
# ╟─0a26ad38-5a75-11eb-0424-4d26aa905de5
# ╠═c64e66c8-4c4b-11eb-1686-8fa64d8b2505
# ╠═3a6c466a-5a75-11eb-07e2-ffbf9ec3ffe4
# ╟─03d82c7a-4c58-11eb-0071-bb9ea16bfbb3
# ╠═0bec2d28-4c58-11eb-0a51-95bf50bbfd79
# ╠═6aa73154-661f-11eb-2b88-578eb2dd2ec2
# ╟─ae6064e6-64a4-11eb-24b5-0b0b848aa2d6
# ╟─0fd08728-4c58-11eb-1b71-c9710d398fab
# ╟─2c6097f4-4c58-11eb-0807-d5d8cbfbd62c
# ╟─9505a4d4-4c58-11eb-1e2e-0d080437fa23
# ╠═23dac6cc-4c58-11eb-2c66-f1f79db08536
# ╟─4f80f0a8-4c58-11eb-3679-c186c61c5a14
# ╠═6d496b44-4c58-11eb-33b6-5b4d6315b6ea
# ╟─56e8f6b4-5a75-11eb-3eeb-ffec491be69c
# ╠═8245e46e-5a75-11eb-2d0a-27ef6a1f2492
# ╠═942b88b4-5a75-11eb-3e7b-4534bf4a7b12
# ╟─7bc7bdf4-4c58-11eb-1fd8-376ac6da5ab2
# ╠═74d97654-4c58-11eb-344b-8d6df24323d5
# ╟─9bb1e83a-5a75-11eb-0fcc-59f6cc50bf6a
# ╠═bcd5696a-5a75-11eb-0ec1-f116216aa682
# ╠═d74ff724-5a75-11eb-39a1-a963cd64948b
# ╠═dd8978a4-5a75-11eb-287e-03e4272b6f2c
# ╟─9cecd9a6-4c58-11eb-22dc-33cd2559d815
# ╟─c6d236da-4c58-11eb-2714-3f5c43583a3d
# ╠═3253ab74-4c58-11eb-178e-83ea8aba9c8f
# ╠═32593936-4c58-11eb-174c-0bb20d93dde5
# ╠═fcfb511c-5a75-11eb-2181-33d147ab1806
# ╠═1882840a-5a76-11eb-3392-81c2915487f5
# ╠═3fc787d6-5a76-11eb-06e9-5378d27ce011
# ╟─ebb09172-4c58-11eb-1cc9-91193c57677d
# ╟─3aa37510-58bb-11eb-2ecb-37ce4428269c
# ╟─ee9069e2-63a7-11eb-12b9-97ae270506f4
# ╠═3de1f1aa-58bd-11eb-2ffc-0de292b13840
# ╠═5f47cdf0-58be-11eb-1bca-a3d0941b9bea
# ╟─6d43af49-f127-4ffe-ba97-0f04fb792efb
# ╟─c1e377c4-64a4-11eb-3e7f-b163cb465057
# ╟─b1a00da4-4cfe-11eb-0aff-69099e40d28f
# ╠═5619fd6c-4cfe-11eb-1512-e1800b6c7df9
# ╟─eae5611f-913c-48e4-bc1b-86d33908ac46
# ╟─e5293248-64a4-11eb-0d30-53a15bec0d01
# ╟─c6e16d7a-58cf-11eb-32a4-3372939066e3
# ╠═cb20fffe-58cf-11eb-1b65-49699f2d3699
# ╠═cee388d2-58cf-11eb-3b88-971b4b85e957
# ╟─be4c9bbe-0ec1-4477-b422-dcf308cd6f5b
# ╟─41b19e20-4d0f-11eb-1c3c-572cc5243d99
# ╟─04aff640-58bb-11eb-1bb6-69ad9fc32314
# ╟─0c91ce30-58b9-11eb-3617-4d87682831dd
# ╠═75d14674-58ba-11eb-3868-172fc00a0eb8
# ╟─cf2b6aa4-b0e8-471a-95fc-646bbddb989a
# ╟─69dc67fa-4cff-11eb-331e-25ffdced4323
# ╠═9f1a2834-4d0f-11eb-3c3e-b7ff55f65dd3
# ╠═85fb018e-4c1d-11eb-2519-a5abe100748e
# ╟─65a45d44-7942-4808-bed0-0369077c5edb
# ╟─00121c4e-64a5-11eb-2993-61c695c4e6a1
# ╟─b56686ec-4cfa-11eb-2b14-a5d49a137cc5
# ╠═a8837ec2-5a4b-11eb-2930-55e48850b7db
# ╟─80159cd5-2e46-4375-b899-a44e45a948da
# ╟─16ec4ee4-64a5-11eb-26f3-15313b8b5acb
# ╟─a0026d1c-bcdb-4a2e-b1a1-11c5235a4956
# ╟─2e7973b6-4d0f-11eb-107c-cdaf349428c0


text:
using LinearAlgebra

# integrate sinh^(n)(ax) from 0 to r
function integrate_sinh(n; r=1.0, a=1.0)
    if n==0
        return r
    elseif n==1
        return (cosh(a*r)-1)/a
    else
        return (sinh(a*r)^(n-1))*cosh(a*r)/(a*n) - (n-1)/n * integrate_sinh(n-2,r=r, a=a)
    end
end

function hyp_radial_density(r, d; curvature=-1.0, radius=1.0)
#    k = 1/(curvature^2)
    k = 1.0
    hyp_r = (sinh(r/k))^(d-1)/integrate_sinh(d-1,r=radius, a=radius/k)
    return hyp_r
end

function euc_radial_density(r, d; radius=1.0)
    return d*(r^(d-1))/radius^d
end

# rejection sampling n=numofpts points from density dens, where the argument lies between 0 and maxval

function rejection_sampling(dens::Function, maxval,numofpts=1)
    max_val = dens(maxval)
    iter = 1
    rands = Array{Float64,1}(undef, numofpts)
    while iter <= numofpts
        x = rand()*maxval
        val = dens(x)
        u = rand()
        if u*max_val < val
            rands[iter] = x
            iter+=1
        end
    end
    return rands
end

function sample_hyp_rad(d, numofpts=1; curvature=-1.0, radius=1.0)
    rands = rejection_sampling(x->hyp_radial_density(x,d,curvature=curvature, radius=radius), radius,numofpts)

    # ...radius within the Poincare ball
    euc_rands = map(x->tanh(x/2.0),rands)
    return euc_rands
end

function sample_euc_rad(d, numofpts=1; radius=1.0)
    rands = rejection_sampling(x->euc_radial_density(x,d,radius=radius), radius, numofpts)
    return rands
end

function sample_sph(d, numofpts=1; curvature=1.0)
    rands = []
    i=0
    while i<=numofpts
        vec = randn(d+1)
        if vec[d+1]>0
            push!(rands, normalize(vec))
            i+=1
        end
    end
    return rands
end

function sample_sphere(d, numofpts=1)
    vecs = randn(d, numofpts)
    rands = []
    for i=1:numofpts
        push!(rands, normalize(vecs[:,i]))
    end
    return rands
end

function sample_hyp(d, numofpts=1; radius=1.0, curvature=-1)
    sphere_pts = sample_sphere(d,numofpts)
    radii = sample_hyp_rad(d, numofpts, radius=radius, curvature=curvature)
    ball_pts = [radii[i]*sphere_pts[i] for i=1:numofpts]
    return ball_pts
end

function sample_euc(d, numofpts=1; radius=1.0)
    sphere_pts = sample_sphere(d,numofpts)
    radii = sample_euc_rad(d, numofpts, radius=radius)
    ball_pts = [radii[i]*sphere_pts[i] for i=1:numofpts]
    return ball_pts
end

function sample_ball(d, numofpts=1; radius=1.0, curvature=0.0)
    sphere_pts = sample_sphere(d,numofpts)
    if curvature < 0
        radii = sample_hyp_rad(d, numofpts, radius=radius, curvature=curvature)
        ball_pts = [radii[i]*sphere_pts[i] for i=1:numofpts]
    elseif curvature == 0.0
        radii = sample_euc_rad(d, numofpts, radius=radius)
        ball_pts = [radii[i]*sphere_pts[i] for i=1:numofpts]
    elseif curvature > 0
        ball_pts = sample_sph(d, numofpts, curvature=curvature)
    end
end

function hyp_distance(pts; curvature=-1.0)
    distances = zeros(length(pts), length(pts))
    for i=1:length(pts)
        for j=1:i-1
            nx = 1-(norm(pts[i]))^2
            ny = 1-(norm(pts[j]))^2
            delta = 2 * norm(pts[i]-pts[j])^2/(nx*ny)
            distances[i,j] = acosh(1+delta)
            distances[j,i] = distances[i,j]
        end
    end
    return distances
end

function euc_distance(pts)
    distances = zeros(length(pts), length(pts))
    for i=1:length(pts)
        for j=1:i-1
            distances[i,j] = norm(pts[i]-pts[j])
            distances[j,i] = distances[i,j]
        end
    end
    return distances
end

function sph_distance(pts; curvature=1.0)
    distances = zeros(length(pts), length(pts))
    for i=1:length(pts)
        for j=1:i-1
            distances[i,j] = acos(dot(pts[i],pts[j]))
            distances[j,i] = distances[i,j]
        end
    end
    return distances
end

function distance_matrix(pts; curvature=0.0)
    if curvature < 0
        return hyp_distance(pts, curvature=curvature)
    elseif curvature == 0
        return euc_distance(pts)
    elseif curvature > 0
        return sph_distance(pts, curvature=curvature)
    end
end

function to_density(matr)
    dens_matr = zeros(size(matr,1), size(matr,2))
    n = size(matr)[1]
    all_entries = sort(setdiff(unique(matr), 0.0))
    total = binomial(n,2)
    for i=1:n
        for j=i+1:n
            dens_matr[i,j] = (findfirst(x->x==matr[i,j], all_entries))/total
            dens_matr[j,i] = dens_matr[i,j]
        end
    end
    return dens_matr
end

function to_density!(matr)
    matr = to_density(matr)
    return matr
end

ext: jl


text:
using JuMP
using Gurobi
using JLD2
using Suppressor
using Printf
using XLSX

max_stars_repo_name: nimazzi/PlanModelDet

max_issues_repo_name: nimazzi/PlanModelDet

max_forks_repo_name: nimazzi/PlanModelDet


text:
# define functions used
# give type , market_type, action_type, t_shock

println("market_type a ", sim_set["market_type"])
println("action_type a ", sim_set["action_type"])

if econ_set["prod_trans_type"] == 1
    q_draw_productivity_gen(productivity, product_trans, econ_set) = q_draw_productivity_norm(productivity, product_trans, econ_set)
else
    q_draw_productivity_gen(productivity, product_trans, econ_set) = q_draw_productivity_discrete(productivity, product_trans, econ_set)
end

# choose utility function
q_utility_gen(consume) = q_utility_log(consume)

if sim_set["type"] == 0
    if sim_set["action_type"] == 0  # simple number
        q_action_gen_one(input, para, scale, shock_p) = q_action_simple_one(input, para, scale, shock_p)
        n_decide_para = 1
    end
    if sim_set["action_type"] == 1 # linear
        q_action_gen_one(input, para, scale, shock_p) = q_action_lin_one_three(input, para, scale, shock_p)
        n_decide_para = 3
        println("n_decide_para ", n_decide_para)
    end
    if sim_set["action_type"] == 2 # polynomial inputs
        q_action_gen_one(input, para, scale, shock_p) = q_action_po_one_six(input, para, scale, shock_p)
        n_decide_para = 6
        println("n_decide_para ", n_decide_para)
    end
    if sim_set["action_type"] == 3 # neural net inputs
        q_action_gen_one(input, para, scale, shock_p) = q_action_nn_one_six(input, para, scale, shock_p)
        n_decide_para = 6
        println("n_decide_para ", n_decide_para)
    end
    if sim_set["action_type"] == 4  # linear shock dummy only
        q_action_gen_one(input, para, scale, shock_p) = q_action_lin_one_four(input, para, scale, shock_p)
        n_decide_para = 4
        println("n_decide_para ", n_decide_para)
    end
    if sim_set["action_type"] == 5 # linear inputs plus shock dummy
        q_action_gen_one(input, para, scale, shock_p) = q_action_lin_one_six(input, para, scale, shock_p)
        n_decide_para = 6
        println("n_decide_para ", n_decide_para)
    end
    if sim_set["action_type"] == 6  # poly inputs linear shock dummy
        q_action_gen_one(input, para, scale, shock_p) = q_action_po_one_nine(input, para, scale, shock_p)
        n_decide_para = 9
        println("n_decide_para ", n_decide_para)
    end

    # market type does not matter for rep agent model
    q_market_gen(supply, demand, post_price, price_para) = q_market_clear(supply, demand, post_price, price_para)
    price_para = [0, 0, 0]  # does not matter for rep agent model

    q_gen_period(state, s_para, econ_set, decide_set, account_a, shock_p) = q_rep_period(state, s_para, econ_set, decide_set, account_a, shock_p)
end

#


text:
using SymbolicMarkov
using SymbolicDiff
using SparseMatrix
using NMarkov
using LinearAlgebra
using Test

# include("test_prob.jl")
# include("test_ctmc.jl")
include("test_markov.jl")
# include("test_mexp.jl")
# include("test_ctmctran.jl")


text:
# ---
# title: 1347. Minimum Number of Steps to Make Two Strings Anagram
# id: problem1347
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/minimum-number-of-steps-to-make-two-strings-anagram/description/>
# hidden: true
# ---
#
# Given two equal-size strings `s` and `t`. In one step you can choose **any
# character** of `t` and replace it with **another character**.
#
# Return _the minimum number of steps_ to make `t` an anagram of `s`.
#
# An  **Anagram**  of a string is a string that contains the same characters
# with a different (or the same) ordering.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "bab", t = "aba"
#     Output: 1
#     Explanation: Replace the first 'a' in t with b, t = "bba" which is anagram of s.
#
#
# **Example 2:**
#
#
#
#     Input: s = "leetcode", t = "practice"
#     Output: 5
#     Explanation: Replace 'p', 'r', 'a', 'i' and 'c' from t with proper characters to make t anagram of s.
#
#
# **Example 3:**
#
#
#
#     Input: s = "anagram", t = "mangaar"
#     Output: 0
#     Explanation: "anagram" and "mangaar" are anagrams.
#
#
# **Example 4:**
#
#
#
#     Input: s = "xxyyzz", t = "xxyyzz"
#     Output: 0
#
#
# **Example 5:**
#
#
#
#     Input: s = "friend", t = "family"
#     Output: 4
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 50000`
#   * `s.length == t.length`
#   * `s` and `t` contain lower-case English letters only.
#
#
## @lc code=start
using LeetCode

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


text:
function joInvertibleMatrix(A::AbstractMatrix{F};DDT::DataType=F,RDT::DataType=promote_type(F,DDT)) where {F<:Number}
    (m,n) = size(A)
    m==n || throw(ArgumentError("A must be square"))
    Fact = lu(A);
    (L,U,p,q,R) = Fact.L, Fact.U, Fact.p, Fact.q, Fact.Rs
    P = joPermutation(p,DDT=DDT)
    Q = joPermutation(q,DDT=DDT)
    Ut = U'
    Lt = L'
    forw_div = v->Q'*(U\(L\(P*(R.*v))))
    adj_div = v->conj(R).*(P'*(Lt\(Ut\(Q*v))))

    return joLinearFunction_A(n,n,
                              v->A*v,
                              v->A'*v,
                              forw_div,
                              adj_div,
                              DDT,RDT,name="joInvertibleMatrix",fMVok=true,iMVok=true)
end


text:
using Base.isless
# Let's define the less than operator for any two vectors that have the same type:
# This does lexicographic comparison, we use it on vectors of chars in this task.
function Base.isless(t1, t2)
    for (a, b) in zip(t1, t2) # zip only to the shorter length
        if !isequal(a, b)
            return isless(a, b)
        end
    end
    return length(t1) < length(t2)
end

# The sort function of Julia doesn't work on strings, so we write one:
# This returns a sorted vector of the chars of the given string
sortchars(s::AbstractString) = sort(collect(Char, s))

# Custom comparator function for sorting the loaded wordlist
sortanagr(s1::AbstractString, s2::AbstractString) =
    if length(s1) != length(s2) length(s1) < length(s2) else sortchars(s1) < sortchars(s2) end

# Tests if two strings are deranged anagrams, returns a bool:
# in our case s2 is never longer than s1
function deranged(s1::AbstractString, s2::AbstractString)
    # Tests for derangement first
    for (a, b) in zip(s1, s2)
        if a == b return false end
    end
    # s1 and s2 are deranged, but are they anagrams at all?
    return sortchars(s1) == sortchars(s2)
end

# Task starts here, we load the wordlist line by line, strip eol char, and sort the wordlist
# in a way that ensures that longer words come first and anagrams go next to each other
words = sort(open(readlines, "./data/unixdict.txt"), rev = true, lt = sortanagr)

# Now we just look for deranged anagrams in the neighbouring words of the sorted wordlist
for i in 1:length(words)-1
    if deranged(words[i], words[i+1])
        # The first match is guaranteed to be the longest due to the custom sorting
        println("The longest deranged anagrams are $(words[i]) and $(words[i+1])")
        break
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using JuliaPBF.StructPBF

function GenParticles(inAnsDataStruct::AnalysisDataStruct)
    npart = inAnsDataStruct.nMaxParticles
    ar_ParticleData = Vector{ParticleData}(undef,npart)
    for ii in 1:npart
        ar_ParticleData[ii] = ParticleData(
                                inAnsDataStruct.nMaxParticles+1,  # gridID :: Int64
                                Vec2(),                           # vel    :: Vec2
                                Vec2(),                           # pos    :: Vec2
                                Vec2(),                           # temppos:: Vec2
                                0.0,                              # Lambda :: Float64
                                0.0,                              # mass   :: Float64
                                0.0)                              # phase  :: Int64
    end
    FillInitialBox(inAnsDataStruct, ar_ParticleData)
    return ar_ParticleData
end

function GenGrids(inAnsDataStruct::AnalysisDataStruct)
    start_point = inAnsDataStruct.analysisBox.staPoint
    end_point = inAnsDataStruct.analysisBox.endPoint
    grid_size = 2.0 * inAnsDataStruct.kernelRadius
    n_Grid_x = Int(ceil((end_point.x - start_point.x)/grid_size))
    n_Grid_y = Int(ceil((end_point.y - start_point.y)/grid_size))

    ini_nNearGrid      = Vector{Int64}(undef,n_Grid_x*n_Grid_y)
    ini_nGridParticles = Vector{Int64}(undef,n_Grid_x*n_Grid_y)
    ini_NearGridIndex  = Vector{Int64}(undef,9*(n_Grid_x-2)*(n_Grid_y-2) + 6*(2*(n_Grid_x-2)+2*(n_Grid_y-2)) + 4*4)
    nong = 0
    for iy in 1:n_Grid_y
        for ix in 1:n_Grid_x
            # 아랫줄
            if iy > 1
                if ix > 1
                    nong+=1
                    ini_NearGridIndex[nong] = ix-1 + (iy-2)*n_Grid_x
                end
                nong+=1
                ini_NearGridIndex[nong] = ix + (iy-2)*n_Grid_x
                if ix < n_Grid_x
                    nong+=1
                    ini_NearGridIndex[nong] = ix+1 + (iy-2)*n_Grid_x
                end
            end
            # 자기줄
            if ix > 1
                nong+=1
                ini_NearGridIndex[nong] = ix-1 + (iy-1)*n_Grid_x
            end
            nong+=1
            ini_NearGridIndex[nong] = ix + (iy-1)*n_Grid_x
            if ix < n_Grid_x
                nong+=1
                ini_NearGridIndex[nong] = ix+1 + (iy-1)*n_Grid_x
            end
            # 윗줄
            if iy < n_Grid_y
                if ix > 1
                    nong+=1
                    ini_NearGridIndex[nong] = ix-1 + iy*n_Grid_x
                end
                nong+=1
                ini_NearGridIndex[nong] = ix + iy*n_Grid_x
                if ix < n_Grid_x
                    nong+=1
                    ini_NearGridIndex[nong] = ix+1 + iy*n_Grid_x
                end
            end
            ii = ix + (iy-1)*n_Grid_x
            ini_nNearGrid[ii] = nong
            ini_nGridParticles[ii] = 0
        end
    end

    return GridData(ini_nNearGrid,ini_NearGridIndex,ini_nGridParticles)
end

function GenPhases(inAnsDataStruct::AnalysisDataStruct)
    ## CREATE BOUNDARY LINES ##
    # Boundary line is a object which shape never
    # change, but only 회전/병진. at this time, it
    # won't be added since defining the method how
    # treat the boundary condition.
    # line collision? dummy particles?
    return []
end

function GenBoundaries(inAnsDataStruct::AnalysisDataStruct)
    ## CREATE BOUNDARY LINES ##
    # Boundary line is a object which shape never
    # change, but only 회전/병진. at this time, it
    # won't be added since defining the method how
    # treat the boundary condition.
    # line collision? dummy particles?
    return []
end

# private

function FillInitialBox(inAnsDataStruct::AnalysisDataStruct, arParticleData::Vector{ParticleData})
    ID = 1
    for initBox_i in 1:length(inAnsDataStruct.a_initialBox)
        initBox_phaseID = inAnsDataStruct.a_initialBox[initBox_i].phaseID
        initBox_box     = inAnsDataStruct.a_initialBox[initBox_i].box
        initBox_vel     = inAnsDataStruct.a_initialBox[initBox_i].vel
        initBox_type    = inAnsDataStruct.a_phases[initBox_phaseID].type
        if initBox_type == "Fluid"
            initBox_name      = inAnsDataStruct.a_phases[initBox_phaseID].Fluid.name
            initBox_size      = inAnsDataStruct.a_phases[initBox_phaseID].Fluid.size
            initBox_density   = inAnsDataStruct.a_phases[initBox_phaseID].Fluid.density
            initBox_viscosity = inAnsDataStruct.a_phases[initBox_phaseID].Fluid.viscosity
            # initBox_mass      = initBox_density * (π*(initBox_size*0.5)^2.0) # Circle
            # initBox_mass      = initBox_density * (initBox_size^2.0) # Box
            initBox_mass      = initBox_density * ((initBox_size^2.0)*0.85 + (π*(initBox_size*0.5)^2.0)*0.15)
        end
        boxDX = initBox_box.endPoint.x - initBox_box.staPoint.x
        boxDY = initBox_box.endPoint.y - initBox_box.staPoint.y
        for ix in 1:div(boxDX+initBox_size*0.5,initBox_size)
            for iy in 1:div(boxDY+initBox_size*0.5,initBox_size)
                xx = initBox_box.staPoint.x + (ix-0.5)*initBox_size
                yy = initBox_box.staPoint.y + (iy-0.5)*initBox_size
                arParticleData[ID] = ParticleData(
                                        0,              # gridID :: Int64
                                        initBox_vel,    # vel    :: Vec2
                                        Vec2(xx,yy),    # pos    :: Vec2
                                        Vec2(xx,yy),    # temppos:: Vec2
                                        0.0,            # lambda :: Float64
                                        initBox_mass,   # mass   :: Float64
                                        initBox_phaseID)# phase  :: Int64
                ID += 1
            end
        end
    end
end

max_stars_repo_name: Hipgineer/JuliaPBF

max_issues_repo_name: Hipgineer/JuliaPBF

max_forks_repo_name: Hipgineer/JuliaPBF


text:
import PhysicsEngine

@testset "Simulation" begin
  sim = PhysicsEngine.newSimulation()

  # Time test
  initialTime = sim.time
  PhysicsEngine.simulationStep!(sim)
  finalTime = sim.time
  @test finalTime == initialTime + 1

  # Test add particle
  particle = PhysicsEngine.Electron([1, 1, 1], [2, 2, 2])
  @test (particle in sim.objects) == false
  PhysicsEngine.addParticle!(sim, particle)
  @test (particle in sim.objects) == true

  # Remove particle
  PhysicsEngine.removeParticle!(sim, particle)
  @test (particle in sim.objects) == false
end

max_stars_repo_name: jonathanbuchanan/PhysicsEngine

max_issues_repo_name: jonathanbuchanan/PhysicsEngine

max_forks_repo_name: jonathanbuchanan/PhysicsEngine


text:
const CEND      = "\33[0m"
const CBOLD     = "\33[1m"
const CITALIC   = "\33[3m"
const CURL      = "\33[4m"
const CBLINK    = "\33[5m"
const CBLINK2   = "\33[6m"
const CSELECTED = "\33[7m"

const CBLACK  = "\33[30m"
const CRED    = "\33[31m"
const CGREEN  = "\33[32m"
const CYELLOW = "\33[33m"
const CBLUE   = "\33[34m"
const CVIOLET = "\33[35m"
const CBEIGE  = "\33[36m"
const CWHITE  = "\33[37m"

const CBLACKBG  = "\33[40m"
const CREDBG    = "\33[41m"
const CGREENBG  = "\33[42m"
const CYELLOWBG = "\33[43m"
const CBLUEBG   = "\33[44m"
const CVIOLETBG = "\33[45m"
const CBEIGEBG  = "\33[46m"
const CWHITEBG  = "\33[47m"

const CGREY    = "\33[90m"
const CRED2    = "\33[91m"
const CGREEN2  = "\33[92m"
const CYELLOW2 = "\33[93m"
const CBLUE2   = "\33[94m"
const CVIOLET2 = "\33[95m"
const CBEIGE2  = "\33[96m"
const CWHITE2  = "\33[97m"

const CGREYBG    = "\33[100m"
const CREDBG2    = "\33[101m"
const CGREENBG2  = "\33[102m"
const CYELLOWBG2 = "\33[103m"
const CBLUEBG2   = "\33[104m"
const CVIOLETBG2 = "\33[105m"
const CBEIGEBG2  = "\33[106m"
const CWHITEBG2  = "\33[107m"

max_stars_repo_name: GabrielTrettel/tetris

max_issues_repo_name: GabrielTrettel/Portran

max_forks_repo_name: GabrielTrettel/Portran


text:
"""
A matrix of colors is interpreted as an image
"""
_default(::Node{Array{RGBA{N0f8}, 2}}, ::Style{:default}, ::Dict{Symbol,Any})

function _default(main::MatTypes{T}, ::Style, data::Dict) where T <: Colorant
    @gen_defaults! data begin
        spatialorder = "yx"
    end
    if !(spatialorder in ("xy", "yx"))
        error("Spatial order only accepts \"xy\" or \"yx\" as a value. Found: $spatialorder")
    end
    ranges = get(data, :ranges) do
        const_lift(main, spatialorder) do m, s
            (0:size(m, s == "xy" ? 1 : 2), 0:size(m, s == "xy" ? 2 : 1))
        end
    end
    delete!(data, :ranges)
    @gen_defaults! data begin
        image = main => (Texture, "image, can be a Texture or Array of colors")
        position_x = nothing => Texture
        position_y = nothing => Texture
        primitive = const_lift(ranges) do r
            x, y = minimum(r[1]), minimum(r[2])
            xmax, ymax = maximum(r[1]), maximum(r[2])
            return Rect2f(x, y, xmax - x, ymax - y)
        end => to_uvmesh
        preferred_camera = :orthographic_pixel
        fxaa = false
        shader = GLVisualizeShader("fragment_output.frag", "image.vert", "texture.frag",
            view = Dict("uv_swizzle" => "o_uv.$(spatialorder)"))
    end
end

function to_uvmesh(geom)
    return NativeMesh(const_lift(GeometryBasics.uv_mesh, geom))
end

function to_plainmesh(geom)
    return NativeMesh(const_lift(GeometryBasics.triangle_mesh, geom))
end

"""
A matrix of Intensities will result in a contourf kind of plot
"""
function gl_heatmap(main::MatTypes{T}, data::Dict) where T <: AbstractFloat
    @gen_defaults! data begin
        intensity = main => Texture
        primitive = Rect2(0f0,0f0,1f0,1f0) => native_triangle_mesh
        nan_color = RGBAf(1, 0, 0, 1)
        highclip = RGBAf(0, 0, 0, 0)
        lowclip = RGBAf(0, 0, 0, 0)
        color_map = nothing => Texture
        color_norm = nothing
        stroke_width::Float32 = 0.0f0
        levels::Float32 = 0f0
        stroke_color = RGBA{Float32}(0,0,0,0)
        shader = GLVisualizeShader("fragment_output.frag", "heatmap.vert", "intensity.frag")
        fxaa = false
    end
    return data
end

#Volumes
const VolumeElTypes = Union{Gray, AbstractFloat}

const default_style = Style{:default}()

using .GLAbstraction: StandardPrerender

struct VolumePrerender
    sp::StandardPrerender
end
VolumePrerender(a, b) = VolumePrerender(StandardPrerender(a, b))

function (x::VolumePrerender)()
    x.sp()
    glEnable(GL_CULL_FACE)
    glCullFace(GL_FRONT)
end

function _default(main::VolumeTypes{T}, s::Style, data::Dict) where T <: VolumeElTypes
    @gen_defaults! data begin
        volumedata = main => Texture
        hull = Rect3f(Vec3f(0), Vec3f(1)) => to_plainmesh
        model = Mat4f(I)
        modelinv = const_lift(inv, model)
        color_map = default(Vector{RGBA}, s) => Texture
        color_norm = color_map == nothing ? nothing : const_lift(extrema2f0, main)
        color = color_map == nothing ? default(RGBA, s) : nothing

        algorithm = MaximumIntensityProjection
        absorption = 1f0
        isovalue = 0.5f0
        isorange = 0.01f0
        shader = GLVisualizeShader("fragment_output.frag", "util.vert", "volume.vert", "volume.frag")
        prerender = VolumePrerender(data[:transparency], data[:overdraw])
        postrender = () -> glDisable(GL_CULL_FACE)
    end
    return data
end

function _default(main::VolumeTypes{T}, s::Style, data::Dict) where T <: RGBA
    @gen_defaults! data begin
        volumedata = main => Texture
        hull = Rect3f(Vec3f(0), Vec3f(1)) => to_plainmesh
        model = Mat4f(I)
        modelinv = const_lift(inv, model)
        # These don't do anything but are needed for type specification in the frag shader
        color_map = nothing => Texture
        color_norm = nothing
        color = color_map === nothing ? default(RGBA, s) : nothing

        algorithm = AbsorptionRGBA
        shader = GLVisualizeShader("fragment_output.frag", "util.vert", "volume.vert", "volume.frag")
        prerender = VolumePrerender(data[:transparency], data[:overdraw])
        postrender = () -> glDisable(GL_CULL_FACE)
    end
end


text:
# This file was generated, do not modify it.

import MLJ: schema, std, mean, median, coerce, coerce!, scitype
using DataFrames
using UrlDownload
using PyPlot
ioff() # hide

raw_data = urldownload("https://github.com/tlienart/DataScienceTutorialsData.jl/blob/master/data/wri_global_power_plant_db_be_022020.csv?raw=true")
data = DataFrame(raw_data);

schema(data)

is_active(col) = !occursin(r"source|generation", string(col))
active_cols = [col for col in names(data) if is_active(col)]
select!(data, active_cols);

select!(data, Not([:wepp_id, :url, :owner]))
schema(data)

describe(data)

# The describe() function shows that there are several features with missing values.

capacity = select(data, [:country, :primary_fuel, :capacity_mw]);
first(capacity, 5)

cap_gr = groupby(capacity, [:country, :primary_fuel]);

cap_mean = combine(cap_gr, :capacity_mw => mean)
cap_sum = combine(cap_gr, :capacity_mw => sum)
first(cap_sum, 3)

ctry_selec = r"BEL|FRA|DEU"
tech_selec = r"Solar"

cap_sum_plot = cap_sum[occursin.(ctry_selec, cap_sum.country) .& occursin.(tech_selec, cap_sum.primary_fuel), :]

sort!(cap_sum_plot, :capacity_mw_sum, rev=true)

figure(figsize=(8,6))

plt.bar(cap_sum_plot.country, cap_sum_plot.capacity_mw_sum, width=0.35)
plt.xticks(rotation=90)

savefig(joinpath(@OUTPUT, "D0-processing-g1.svg")) # hide

cap_sum_ctry_gd = groupby(capacity, [:country]);
cap_sum_ctry = combine(cap_sum_ctry_gd, :capacity_mw => sum);

cap_sum = DataFrame(cap_sum);
cap_sum_ctry = DataFrame(cap_sum_ctry);
cap_share = leftjoin(cap_sum, cap_sum_ctry, on = :country, makeunique = true)
cap_share.capacity_mw_share = cap_share.capacity_mw_sum ./ cap_share.capacity_mw_sum_1;

nMissings = length(findall(x -> ismissing(x), data.commissioning_year))

nMissings_share = nMissings/size(data)[1]

typeof(data.commissioning_year)

data_nmiss = dropmissing(data, :commissioning_year);

map!(x -> round(x, digits=0), data_nmiss.commissioning_year, data_nmiss.commissioning_year);

# We can now calculate plant age for each plant (worth remembering that the dataset only contains active plants)

current_year = fill!(Array{Float64}(undef, size(data_nmiss)[1]), 2020);
data_nmiss[:, :plant_age] = current_year - data_nmiss[:, :commissioning_year];

mean_age = mean(skipmissing(data_nmiss.plant_age))
median_age = median(skipmissing(data_nmiss.plant_age))

figure(figsize=(8,6))

plt.hist(data_nmiss.plant_age, color="blue", edgecolor="white", bins=100,
      density=true, alpha=0.5)
plt.axvline(mean_age, label = "Mean", color = "red")
plt.axvline(median_age, label = "Median")

plt.legend()

plt.xlim(0,)

savefig(joinpath(@OUTPUT, "D0-processing-g2.svg")) # hide

age = select(data_nmiss, [:country, :primary_fuel, :plant_age])
age_mean = combine(groupby(age, [:country, :primary_fuel]), :plant_age => mean)

coal_means = age_mean[occursin.(ctry_selec, age_mean.country) .& occursin.(r"Coal", age_mean.primary_fuel), :]
gas_means = age_mean[occursin.(ctry_selec, age_mean.country) .& occursin.(r"Gas", age_mean.primary_fuel), :]

width = 0.35  # the width of the bars

fig, (ax1, ax2) = plt.subplots(1,2)

fig.suptitle("Mean plant age by country and technology")

ax1.bar(coal_means.country, coal_means.plant_age_mean, width, label="Coal")
ax2.bar(gas_means.country, gas_means.plant_age_mean, width, label="Gas")

ax1.set_ylabel("Age")

ax1.set_title("Coal")
ax2.set_title("Gas")

savefig(joinpath(@OUTPUT, "D0-processing-g3.svg")) # hide


text:
export LogScore

"""
    struct LogScore{I} <: Score{I}

A Log score.
"""
struct LogScore{I} <: Score{I}
    logscores :: Dict{I, Float64}
    function LogScore(ls::Dict{I, Float64}) where I
        new{I}(ls)
    end
    function LogScore(vs::Vector{I}, ss::Vector{Float64}) where I
        # must handle non-unique values correctly
        d = Dict{I, Float64}()
        for (v,s) in zip(vs,ss)
            d[v] = logsumexp([get(d, v, -Inf), s])
        end
        new{I}(d)
    end
end

@impl begin
    struct LogScoreGetScore end
    function get_score(sf::LogScore{I}, i::I)::AbstractFloat where {I}
        return exp(get_log_score(sf, i))
    end
end

@impl begin
    struct LogScoreGetLogScore end
    function get_log_score(sf::LogScore{I}, x::I)::AbstractFloat where {I}
        return x in keys(sf.logscores) ? sf.logscores[x] : -Inf
    end
end


text:
#===================================================================================================
Abstract types and common interfaces.
===================================================================================================#
"""
    Evolver

Abstract type of all DE core.
"""
abstract type Evolver end

"""
    Options

Abstract type of all DE options.
"""
abstract type Options end

abstract type Controller end

"""
    DEEvolver <: Evolver

An abstract normal (single-population) DE evolution type
"""
abstract type DEEvolver <: Evolver end

"""
    DEOptions <: Options

Abstract type of options for a DE, particually, a single-population one. See also [`DEEvolver`](@ref).
"""
abstract type DEOptions <: Options end
abstract type DEController <: Controller end

_throw_unimplemented_error() = error("Unimplemented")

# all subtypes of `DEEvolver` should have the following fields
"""
    pop(e::DEEvolver)

Get the population of a DE, which is typically a matrix.
"""
pop(e::DEEvolver) = e.pop

"""
    fitness(e::DEEvolver)

Get the fitness of the population, which is generally a vector for single-objective DE or a matrix
for multi-objective DE.
"""
fitness(e::DEEvolver) = e.fitness

"""
    senses(e::DEEvolver)

Get the senses of optimization. In traditional evolutionary algorithms, the fitness of an individual
is going to be maximized. If minimization is required instead, then set a negative sense.
"""
senses(e::DEEvolver) = e.senses

"""
    bounds(e::DEEvolver)

Get the bound constraints,typically a ``n x 2`` matrix for `n` variables.
"""
bounds(e::DEEvolver) = e.bounds

"""
    options(e::DEEvolver)

Get the options of this DE.
"""
options(e::DEEvolver) = e.options

"""
    nindividuals(e::DEEvolver)

Number of individuals in the population.
"""
function nindividuals(e::DEEvolver)
    return size(pop(e), 2)
end

"""
    ndims(de::DEEvolver)

Number of variables (dimension) of the DE.
"""
Base.ndims(de::DEEvolver) = size(pop(de), 1)

"""
    nobjectives(e::Evolver)

Number of objectives.
"""
function nobjectives(e::DEEvolver)
    f = fitness(e)
    if ndims(f) == 1
        return 1
    end
    return size(f, 1)
end

"""
    best_individual(e::Evolver; all=false)

Get individuals that have the best fitness. If `all` is `true`, then all such individuals are returned;
otherwise, only one of them is fetched.
"""
function best_individual(e::Evolver; all=false)
    _throw_unimplemented_error()
end

"""
    best_fitness(e::Evolver)

Get the best fitness among the individuals of the evolver.
"""
function best_fitness(e::Evolver)
    _throw_unimplemented_error()
end


text:
#=

 This file implements the functions presented in the paper

 * [2] Wexler, A. and Hyland, R. W., "Formulations for the thermodynamic properties of dry air from 173.15 K to 473.15 K, and of saturated moist air from 173.15 K to 372.15 K at pressures to 5 MPa

=#

"""
    ```Bm(Tk,xv)```

First virial coefficient Bm of moist air of saturated vapor eq 2 [2]

 * `Tk` Temperature in K
 * `xv` Molar fraction of water vapor in the moist air
 * Output: Bm in m^3/mol
"""
function Bm(Tk, xv)
    xa = 1 - xv
    xa*xa*Baa(Tk) + 2*xa*xv*Baw(Tk) + xv*xv*Bww(Tk)
end

"""
    ```dBm(Tk,xv)```

Derivative of first virial coefficient Bm of moist air of saturated vapor eq 2 [2]

 * `Tk` Temperature in K
 * `xv` Molar fraction of water vapor in the moist air
 * Output: dBm/dT in m^3/mol/K
"""
function dBm(Tk, xv)
    xa = 1.0 - xv
    xa*xa*dBaa(Tk) + 2*xa*xv*dBaw(Tk) + xv*xv*dBww(Tk)
end

"""
    ```Cm(Tk,xv)```

Second virial coefficient Cm of moist air of saturated vapor eq 3 [2]

 * `Tk` Temperature in K
 * `xv` Molar fraction of water vapor in the moist air
 * Output: Cm in m^6/mol^2
"""
function Cm(Tk, xv)
    xa = 1-xv
    xa*xa*xa*Caaa(Tk) + 3*xa*xa*xv*Caaw(Tk) + 3*xa*xv*xv*Caww(Tk) + xv*xv*xv*Cwww(Tk)
end

"""
    ```dCm(Tk,xv)```

Derivative of second virial coefficient Cm of moist air of saturated vapor eq 3 [2]

 * `Tk` Temperature in K
 * `xv` Molar fraction of water vapor in the moist air
 * Output: dCm/dT in m^6/mol^2/K
"""
function dCm(Tk, xv)
    xa = 1-xv
    xa*xa*xa*dCaaa(Tk) + 3*xa*xa*xv*dCaaw(Tk) + 3*xa*xv*xv*dCaww(Tk) + xv*xv*xv*dCwww(Tk)
end

"""
    ```Baa(Tk)```

Virial coefficient Bₐₐ of dry air eq 10 [2]

 * `Tk` Temperature in K
 * Output: Bₐₐ in m^3/mol

"""
Baa(Tk) = 0.349568e-4 + (1.0/Tk)*(-0.668772e-2 + (1.0/Tk) * (-0.210141e1 + 0.924746e2/Tk))

"""
    ```dBaa(Tk)```

Derivative of Bₐₐ(T). Calculated from equation 10 [2].

 * `Tk` Temperature in K
 * Output: Bₐₐ in m^3/mol/K

"""
dBaa(Tk) = 1.0/(Tk*Tk) * (0.668772e-2 + (1.0/Tk)*(0.420282e1 - 2.774238e2/Tk))

"""
    ```Caa(Tk)```

Second virial coefficient Cₐₐₐ of dry air eq 11 [2].

 * `Tk` Temperature in K
 * Output: Bₐₐ in m^6/mol^2

"""
Caaa(Tk) = (0.125975e-8 + 1.0/Tk * (-0.190905e-6 + 0.632467e-4/Tk))

"""
    ```dCaa(Tk)```

Derivative of second virial coefficient Cₐₐₐ of dry air eq 11 [2].

 * `Tk` Temperature in K
 * Output: Bₐₐ in m^6/mol^2/K

"""
dCaaa(Tk) = 1.0/(Tk*Tk) * (0.190905e-6 - 1.264934e-4/Tk)

"""
    ```Zair(Tk, P, EPS, MAXITER)```

Compressibility factor of dry air.

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * `EPS`: Acceptable error
 * `MAXITER`: Maixmum number of iterations
 * Output: Z (nondimensional)
"""
function Zair(Tk, P, EPS=1e-8, MAXITER=100)
    vm0 = R*Tk/P
    b0 = Baa(Tk) / vm0
    c0 = Caaa(Tk) / (vm0*vm0)
    z = calcz(b0, c0, EPS, MAXITER)
end

"""
    molarvolumeair(Tk)

Molar volume of dry air. Equation 12 of reference [1].
This function requires iteration to compute the volume.

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * `EPS`: Acceptable error
 * `MAXITER`: Maixmum number of iterations
 * Output: molar volume in m^3/mol
"""
function molarvolumeair(Tk, P, EPS=1e-8, MAXITER=100)

    vm0 = R*Tk/P

    b0 = Baa(Tk)/vm0
    c0 = Caaa(Tk)/(vm0*vm0)
    z = calcz(b0, c0, EPS, MAXITER)

    return z*vm0
end

"""
    volumeair

Specific volume of dry air. Uses `molarvolumeair`.

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * `EPS`: Acceptable error
 * `MAXITER`: Maixmum number of iterations
 * Output: molar volume in m^3/kg
"""
volumeair(Tk, P, EPS=1e-8, MAXITER=100) = molarvolumeair(Tk, P, EPS, MAXITER) / Ma

"""
    molarenthalpyair(Tk)

Molar enthalpy of dry air. Equation 13 of reference [2]

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * `EPS`: Acceptable error
 * `MAXITER`: Maixmum number of iterations
 * Output: m^3/mol
"""
function molarenthalpyair(Tk, P, EPS=1e-8, MAXITER=100)
    h1 = -0.79078691e4 + Tk*(0.28709015e2 +
                             Tk*(0.26431805e-2 +
                                 Tk*(-0.10405863e-4 +
                                     Tk*(0.18660410e-7 - 0.97843331e-11*Tk))))
    va = molarvolumeair(Tk, P, EPS, MAXITER)

    h2 = R*Tk/va * ( (Baa(Tk) - Tk*dBaa(Tk)) + (Caaa(Tk) - 0.5*Tk*dCaaa(Tk))/va )

    return (h1 + h2)

end

"""
    enthalpyair(Tk)

Specific enthalpy of dry air. Equation 13 of reference [2]

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * `EPS`: Acceptable error
 * `MAXITER`: Maixmum number of iterations
 * Output: m^3/kg
"""
function enthalpyair(Tk, P, EPS=1e-8, MAXITER=100)
    return molarenthalpyair(Tk, P, EPS, MAXITER) / Ma

end

const ℓ =  (-0.16175159e3, 0.52863609e-2, -0.15608795e-4,
            0.24880547e-7, -0.12230416e-10, 0.28709015e2)

"""
    molarentropyair(Tk)

Specific entropy of dry air. Equation 14 of reference [2]

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * `EPS`: Acceptable error
 * `MAXITER`: Maixmum number of iterations
 * Output: m^3/(kg.K)
"""
function molarentropyair(Tk, P, EPS=1e-8, MAXITER=100)

    s1 = @polyeval(Tk, ℓ, 5) + ℓ[6]*log(Tk) - R*log(P/101325.0)
    va = molarvolumeair(Tk, P, EPS, MAXITER)
    s2 = R*log(P*va/R/Tk) - R/va * ( (Baa(Tk) + Tk*dBaa(Tk)) + 0.5/va * (Caaa(Tk) + Tk*dCaaa(Tk)))

    return (s1 + s2)
end

"""
    entropyair(Tk)

Specific entropy of dry air. Equation 14 of reference [2]

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * `EPS`: Acceptable error
 * `MAXITER`: Maixmum number of iterations
 * Output: m^3/(kg.K)
"""
function entropyair(Tk, P, EPS=1e-8, MAXITER=100)

    return molarentropyair(Tk, P, EPS, MAXITER) / Ma
end

"""
    ```Baw(Tk)```

Cross virial coefficient Baw of saturated vapor eq 15 [2]

 * `Tk` Temperature in K
 * Output: Baw in m^3/mol
"""
Baw(Tk) = 0.32366097e-4 + 1.0/Tk * (-0.141138e-1 + 1.0/Tk * (-0.1244535e1 - 0.2348789e4/(Tk*Tk)))

"""
    ```Baw(Tk)```

Derivative of cross virial coefficient Bww of saturated vapor eq 15 [2]

 * `Tk` Temperature in K
 * Output: dBaw/dT in m^3/mol/K
"""
dBaw(Tk) = 1.0/(Tk*Tk) * (0.141138e-1 + 1.0/Tk * (0.248907e1 + 0.93951568e4 /(Tk*Tk)))

"""
    ```Caaw(Tk)```

Cross virial coefficient Caaw of saturated vapor eq 16 [2]

 * `Tk` Temperature in K
 * Output: Caaw in m^6/mol^2
"""
Caaw(Tk) = 0.482737e-9 + 1.0/Tk * (0.105678e-6 + 1.0/Tk *
                                   (-0.656394e-4 + 1.0/Tk*
                                    (0.294442e-1 - 0.319317e1/Tk)))

"""
    ```dCaaw(Tk)```

Derivative of cross virial coefficient dCaaw/dT of saturated vapor eq 16 [2]

 * `Tk` Temperature in K
 * Output: dCaaw/dT in m^6/mol^2/K
"""
dCaaw(Tk) = 1.0/(Tk*Tk) * (-0.105678e-6 + 1.0/Tk*(1.312788e-4 + 1.0/Tk*(-0.883326e-1 + 1.277268e1/Tk)))

"""
    ```Caww(Tk)```

Cross virial coefficient Caaw of saturated vapor eq 17 [2]

 * `Tk` Temperature in K
 * Output: Caww in m^6/mol^2
"""
Caww(Tk) = -1e-6 * exp( -0.10728876e2 + 1.0/Tk * (0.347802e4 + 1.0/Tk*(-0.383383e6 +  0.33406e8/Tk)))

"""
    ```dCaww(Tk)```

Derivative of cross virial coefficient dCaaw/dT of saturated vapor eq 17 [2]

 * `Tk` Temperature in K
 * Output: dCaww/dT in m^6/mol^2/K
"""
dCaww(Tk) = 1.0/(Tk*Tk) * (-0.347802e4 + 1.0/Tk * (2*0.383383e6 - 3*0.33406e8/Tk)) * Caww(Tk)

"""

    ```e_factor(Tk,P)```

Enhancement factor of moist air.
Uses an iterative procedure to calculate the Enhancement factor defined
in equation 18 [2]

"""
function efactor2(Tk, P, EPS=1e-8, MAXITER=200)
    f = 1.0
    err = 0.0
    for i = 1:MAXITER
        xas = (P-f*Pws(Tk)) / P
        fnovo = exp(lnf(Tk,P,xas))
        err = abs(fnovo-f)

        if err < EPS
            if fnovo < 1
                fnovo = 1.0
            end
            println(i)
            return fnovo
        end

        f = fnovo
    end
    if fnovo < 1.0
        fnovo = 1.0
    end
    throw(ConvergenceError("Enhancement factor calculation did not converge!", fnovo, MAXITER, err))
end

function efactor(Tk, P, relax=1.0, EPS=1e-8, MAXITER=200)

    if Tk < 273.16
        vc = volumeice(Tk) * Mv
        p = Pws_s(Tk)
        k = 0.0
    else
        vc = volumewater(Tk) * Mv
        p = Pws_l(Tk)
        k = henryk(Tk)
    end

    κ = kappa_f(Tk)

    RT = R * Tk

    baa = Baa(Tk)
    bww = Bww(Tk)
    baw = Baw(Tk)
    caaa = Caaa(Tk)
    caaw = Caaw(Tk)
    caww = Caww(Tk)
    cwww = Cwww(Tk)
    fun = f -> (log(f) - lnf(f, P, p, k, κ, vc, RT, baa, baw, bww, caaa, caaw, caww, cwww))

    ef = 1e-8
    f = 1.0
    for i = 1:MAXITER
        g = fun(f)
        dg = (fun(f+ef) - g) / ef
        df = -g/dg

        if abs(df) < EPS
            f = f + df
            if f < 1
                f = 1.0
            end
            return f
        end

        f = f + relax*df
    end
    throw(ConvergenceError("Enhancement factor calculation did not converge!", f, MAXITER, abs(df)))
end

function lnf(f, P, p, k, κ, vc, RT, baa, baw, bww, caaa, caaw, caww, cwww)
    P2 = P*P
    p2 = p*p
    xas = (P - f*p) / P
    xas2 = xas*xas
    RT2 = RT*RT

    t1 = vc/RT * ( (1 + κ*p)*(P-p) - 0.5*κ*(P2 - p2) )
    t2 = log(1.0 - k*xas*P) + (xas2*P/RT)*baa - (2*xas2*P/RT)*baw
    t3 = -(P-p-xas2*P)/RT*bww + xas2*xas*P*P/RT2 * caaa
    t4 = 3*xas2*(1-2xas)*P2/(2*RT2) * caaw - (3xas2*(1-xas)*P2)/RT2 * caww
    t5 = -( (1+2xas)*(1-xas)^2*P2 - p2 )/(2*RT2) * cwww - (xas2*(1-3xas)*(1-xas)*P2)/RT2 * baa*bww
    t6 = -(2xas2*xas*(2-3xas)*P2) / RT2 * baa*baw + ( 6xas2*(1-xas)^2*P2 )/RT2 * bww*baw
    t7 = -3xas2*xas2*P2/(2RT2)*baa*baa - ( 2xas2*(1-xas)*(1-3xas)*P2 )/RT2 * baw*baw
    t8 = -( p2 - (1+3xas)*(1-xas)^3*P2 ) / (2*RT2) * bww*bww

   return t1+t2+t3+t4+t5+t6+t7+t8

end

"""
    ```lnf(Tk, P, xas)```

Auxiliary function used to calculate the enhancement factor.
Actually, this function implements the RHS of eq. 18 of [2].
"""
function lnf2(Tk, P, xas)
    if Tk < 273.16
        vc = volumeice(Tk) * Mv
        p = Pws_s(Tk)
        k = 0.0
    else
        vc = volumewater(Tk) * Mv
        p = Pws_l(Tk)
        k = henryk(Tk)
    end

    κ = kappa_f(Tk)

    RT = R*Tk
    RT2 = RT*RT
    xas2 = xas*xas
    p2 = p*p
    P2 = P*P

    baa = Baa(Tk)
    bww = Bww(Tk)
    baw = Baw(Tk)
    caaa = Caaa(Tk)
    caaw = Caaw(Tk)
    caww = Caww(Tk)
    cwww = Cwww(Tk)

    t1 = vc/RT * ( (1 + κ*p)*(P-p) - 0.5*κ*(P2 - p2) )
    t2 = log(1.0 - k*xas*P) + (xas2*P/RT)*baa - (2*xas2*P/RT)*baw
    t3 = -(P-p-xas2*P)/RT*bww + xas2*xas*P*P/RT2 * caaa
    t4 = 3*xas2*(1-2xas)*P2/(2*RT2) * caaw - (3xas2*(1-xas)*P2)/RT2 * caww
    t5 = -( (1+2xas)*(1-xas)^2*P2 - p2 )/(2*RT2) * cwww - (xas2*(1-3xas)*(1-xas)*P2)/RT2 * baa*bww
    t6 = -(2xas2*xas*(2-3xas)*P2) / RT2 * baa*baw + ( 6xas2*(1-xas)^2*P2 )/RT2 * bww*baw
    t7 = -3xas2*xas2*P2/(2RT2)*baa*baa - ( 2xas2*(1-xas)*(1-3xas)*P2 )/RT2 * baw*baw
    t8 = -( p2 - (1+3xas)*(1-xas)^3*P2 ) / (2*RT2) * bww*bww

   return t1+t2+t3+t4+t5+t6+t7+t8

end

"""
    molarfracmoist_sat(Tk, P)

Molar fraction of water vapor in saturated moist air

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * Output: molar fraction

"""
molarfracmoist_sat(Tk, P) = efactor(Tk, P) * Pws(Tk) / P

"""
    ```Bww(Tk)```

Virial coefficient Bww of saturated vapor eq 19 [2]

 * `Tk` Temperature in K
 * Output: Bww in m^3/mol
"""
Bww(Tk) = R * Tk * Blin(Tk)

"""
    ```dBww(Tk)```

Derivative of virial coefficient dBww/dT of saturated vapor eq 19 [2]

 * `Tk` Temperature in K
 * Output: dBww/dT in m^3/mol/K
"""
dBww(Tk) = R * (Tk * dBlin(Tk) + Blin(Tk))

"""
    ```Cwww(Tk)```

Second virial coefficient Cwww of saturated vapor eq 20 [2]

 * `Tk` Temperature in K
 * Output: Cwww in m^3/mol
"""
Cwww(Tk) = R*R*Tk*Tk * (Clin(Tk) + Blin(Tk)^2)

"""
    ```dCwww(Tk)```

Derivative of second virial coefficient dCww/dT of saturated vapor eq 20 [2]

 * `Tk` Temperature in K
 * Output: dCwww/dT in m^3/mol/K
"""
dCwww(Tk) = R*Tk*R*Tk * (dClin(Tk) + 2*Blin(Tk) * dBlin(Tk)) + (2*R*R*Tk) * (Clin(Tk) + Blin(Tk)^2)

"""
    kappa_l(Tk)

Isothermal compressibility of saturated liquid water. Equation 21, ref. [2].
Range of applicability 0°C - 150°C but can be extended up to 200°C
More details in reference [4].
 * `Tk` Temperature in K
 * Output: κ in 1/Pa
"""
function kappa_l(Tk)
    Tc = Tk - 273.15

    #@check_range Tc 0 200
    if Tc < 100.0
        k = (50.88496 + Tc*(0.6163813 +
                            Tc*(1.459187e-3 +
                                Tc*(20.08438e-6 + Tc*(-58.47727e-9 + Tc*0.4104110e-9))))) /
                                (1.0 + 0.1967348e-1*Tc)

    else
        k = (50.884917 + Tc*(0.62590623 +
                            Tc*(1.3848668e-3 +
                                Tc*(21.603427e-6 + Tc*(-0.72087667e-7 + Tc*0.46545054e-9))))) /
                                (1.0 + 0.1967348e-1*Tc)
    end

    k*1e-11
end

"""
    kappa_s(Tk)

Isothermal compressibility of ice. Equation 22, ref. [2].

 * `Tk` Temperature in K
 * Output: κ in 1/Pa
"""
kappa_s(Tk) = (8.875 + 0.0165*Tk)*1e-11

kappa_f(Tk) =
    if Tk < 273.16
        kappa_s(Tk)
    else
        kappa_l(Tk)
    end

"""
    henryk_O2(Tk)

Henry's coefficient for solubility of O2 in water.
Implements equation 23 from ref. [2]. See reference [3]
Valid for temperatures between 273.15 K and 500 K

 * `Tk` Temperature in K
 * Output: k in atm/mol * 10^(-4)
"""
function henryk_O2(Tk)

    τ = 1000.0/Tk
    α = -0.0005943
    β = -0.1470
    γ = -0.05120
    δ = -0.1076
    ɛ = 0.8447
    a2 = α;
    a1 = γ*τ + δ
    a0 = β*τ*τ + ɛ*τ - 1.0

    10^( (-a1-sqrt(a1*a1 - 4.0*a2*a0)) / (2.0*a2))
end

"""
    henryk_N2(Tk)

Henry's coefficient for solubility of N2 in water.
Implements equation 23 from ref. [2]. See reference [3]
Valid for temperatures between 273.15 K and 500 K

 * `Tk` Temperature in K
 * Output: k in atm/mol * 10^(-4)
"""
function henryk_N2(Tk)

    τ = 1000.0/Tk
    α = -0.1021
    β = -0.1482
    γ = -0.019
    δ = -0.03741
    ɛ = 0.851
    a2 = α;
    a1 = γ*τ + δ
    a0 = β*τ*τ + ɛ*τ - 1.0

    10^( (-a1-sqrt(a1*a1 - 4.0*a2*a0)) / (2.0*a2))
end

"""
    henryk(Tk)

Henry's coefficient for solubility of air in water.
Implements equations 24 and 25 from ref. [2]. See reference [3]
Valid for temperatures between 273.15 K and 500 K

 * `Tk` Temperature in K
 * Output: k in atm/mol fract * 10^(-4)
"""
function henryk(Tk)

    kO2 = henryk_O2(Tk)
    kN2 = henryk_N2(Tk)
    xO2 = 0.22
    xN2 = 0.78

    k = 1.0 / (xO2/kO2 + xN2/kN2)

    1e-4/k * 1.0 /101325.0
end

"""
    ```molarvolumemoist```

Molar volume of moist air given the molar fraction of water vapor.
Assumes a real gas.

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * `xv`  Molar fraction of water vapor
 * Output: molar volume in m^3/mol
"""
molarvolumemoist(Tk, P, xv, EPS=1e-8, MAXITER=100) = Zmoist(Tk, P, xv, EPS, MAXITER) * R * Tk / P

"""
    ```volumemoist```

Specific volume of moist air given the molar fraction of water vapor.
Assumes a real gas. The specific volume is defined as the volume occupied
by a mixture of dry air and water vapor per kg of dry air.

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * `xv`  Molar fraction of water vapor
 * Output: molar volume in m^3/kg of dry air.
"""
volumemoist(Tk, P, xv, EPS=1e-8, MAXITER=100) = molarvolumemoist(Tk,P,xv,EPS,MAXITER) / (Ma*(1-xv))

"""
    ```Zmoist(Tk, P, xv)```

Compressibility factor of moist air.

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * `xv` molar fraction of water vapor.
 * `EPS`: Acceptable error
 * `MAXITER`: Maixmum number of iterations
 * Output: Z (nondimensional)
"""
function Zmoist(Tk, P, xv, EPS=1e-8, MAXITER=500)

    vm0 = R*Tk/P
    b0 = Bm(Tk, xv)/vm0
    c0 = Cm(Tk, xv)/(vm0*vm0)

    return calcz(b0, c0, EPS, MAXITER, 1.0)
end

"Coefficients a_i to calculate enthalpy of moist air ref[2]"
const a = (0.63290874e1, 0.28709015e2, 0.26431805e-2,
           -0.10405863e-4, 0.18660410e-7, -0.97843331e-11)
"Coefficients d_i to calculate enthalpy of moist air ref[2]"
const d = (-0.5008e-2, 0.32491829e2, 0.65576345e-2,
           -0.26442147e-4, 0.51751789e-7, -0.31541624e-10)

"""
    molarenthalpymoist(Tk, P, xv, EPS=1e-8, MAXITER=100)

Molar enthalphy of moist air defined as enthalpy dry air.

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * `xv` molar fraction of water vapor.
 * `EPS`: Acceptable error
 * `MAXITER`: Maixmum number of iterations
 * Output: J/mol

"""
function molarenthalpymoist(Tk, P, xv, EPS=1e-8, MAXITER=100)
    xa = 1.0-xv
    h1 = xa*(@polyeval(Tk, a, 6) -7914.1982)
    h2 = xv*(@polyeval(Tk, d, 6) + 35994.17)
    vm = molarvolumemoist(Tk, P, xv, EPS, MAXITER)
    h3 = R*Tk/vm * (Bm(Tk, xv) - Tk*dBm(Tk,xv)
                    + 1/vm*( Cm(Tk,xv) - 0.5*Tk*dCm(Tk,xv) ) )
    return (h1 + h2 + h3)
end

"""
    enthalpymoist(Tk, P, xv, EPS=1e-8, MAXITER=100)

Specific enthalphy of moist air defined as enthalpy per kg of
dry air.

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * `xv` molar fraction of water vapor.
 * `EPS`: Acceptable error
 * `MAXITER`: Maixmum number of iterations
 * Output: J/kg of dry air

"""
function enthalpymoist(Tk, P, xv, EPS=1e-8, MAXITER=100)
    hm = molarenthalpymoist(Tk, P, xv, EPS, MAXITER)
    return hm / ((1-xv)*Ma)
end

"Coefficients g_i to calculate enthalpy of moist air ref[2]"
const g = (0.34373874e2, 0.52863609e-2, -0.15608795e-4,
           0.24880547e-7, -0.12230416e-10, 0.28709015e2)
"Coefficients k_i to calculate enthalpy of moist air ref[2]"
const k = (0.2196603e1, 0.19743819e-1, -0.70128225e-4,
           0.14866252e-6, -0.14524437e-9, 0.55663583e-13,
           0.32284652e2)

"""
    molarentropymoist(Tk, P, xv, EPS=1e-8, MAXITER=100)

Molar entropy of moist air defined as entropy per kg of
dry air.

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * `xv` molar fraction of water vapor.
 * `EPS`: Acceptable error
 * `MAXITER`: Maixmum number of iterations
 * Output: J/mol/K of dry air

"""
function molarentropymoist(Tk, P, xv, EPS=1e-8, MAXITER=100)
    xa = 1.0-xv
    h1 =  @polyeval(Tk, g, 5) + g[6]*log(Tk) - 196.125465
    h2 =  @polyeval(Tk, k, 6) + k[7]*log(Tk) - 63.31449

    z = Zmoist(Tk, P, xv, EPS, MAXITER)
    vm = z*R*Tk/P

    h3 = -R * log(P / 101325.0) + xa*R*log(z/xa)
    if xv > 1e-8
        h3 = h3 + xv*R*log(z/xv)
    end

    h4 = -R/vm * (  (Bm(Tk, xv) + Tk*dBm(Tk,xv)) + 0.5/vm*(Cm(Tk, xv) + Tk*dCm(Tk, xv)) )

    return (xa*h1 + xv*h2 + h3 + h4)
end

"""
    entropymoist(Tk, P, xv, EPS=1e-8, MAXITER=100)

Specific entropy of moist air defined as entropy per kg of
dry air.

 * `Tk` Temperature in K
 * `P` Pressure in Pa
 * `xv` molar fraction of water vapor.
 * `EPS`: Acceptable error
 * `MAXITER`: Maixmum number of iterations
 * Output: J/kg/K of dry air

"""
function entropymoist(Tk, P, xv, EPS=1e-8, MAXITER=100)

    return molarentropymoist(Tk, P, xv, EPS, MAXITER) / ( (1-xv)*Ma )

end


text:
# ==============================================================================
# Five-dimensional system
#
# Example from Althoff's thesis, Section 3.2.3
#
#
# system type: continuous LTI system
# state dimension: 5
# input dimension: 1
#
# This five-dimensional system is taken from Example 4.1 in [1].
#
# [1] Colas Le Guernic: Reachability analysis of hybrid systems with linear
#     continuous dynamics. PhD thesis, 2009.
# ==============================================================================

function linear5D_homog()
    A = Float64[-1 -4   0   0   0;
                 4 -1   0   0   0;
                 0  0  -3   1   0;
                 0  0  -1  -3   0;
                 0  0   0   0  -2]
    X0 = Hyperrectangle(low=fill(0.9, 5), high=fill(1.1, 5))
    prob = @ivp(x' = Ax, x(0) ∈ X0)
    tspan = (0.0, 5.0)
    return prob, tspan
end

function linear5D()
    # system matrix
    D = [-1.0 -4.0  0.0  0.0  0.0;
          4.0 -1.0  0.0  0.0  0.0;
          0.0  0.0 -3.0  1.0  0.0;
          0.0  0.0 -1.0 -3.0  0.0;
          0.0  0.0  0.0  0.0 -2.0]
    P = [ 0.6 -0.1  0.1  0.7 -0.2;
         -0.5  0.7 -0.1 -0.8  0.0;
          0.9 -0.5  0.3 -0.6  0.1;
          0.5 -0.7  0.5  0.6  0.3;
          0.8  0.7  0.6 -0.3  0.2]
    A = P * D * inv(P)
    X = Universe(5) # state domain
    U = Ball2(zeros(5), 0.01) # input domain

    X0 = BallInf([1.0, 0.0, 0.0, 0.0, 0.0], 0.1)
    prob = @ivp(x' = Ax + u, x ∈ X, u ∈ U, x(0) ∈ X0) # continuous LTI system
    tspan = (0.0, 5.0)
    return prob, tspan
end


text:
# rounding.jl

function Base.floor(qdt::QDate, p::Dates.Year)
    value(p) < 1 && throw(DomainError())
    years = year(qdt)
    return QDate(years - mod(years, value(p)))
end

function Base.floor(qdt::QDate, p::Dates.Month)
    value(p) < 1 && throw(DomainError())
    y, m, l = yearmonthleap(qdt)
    qdt0 = QDate(y, m, l, 1)
    if value(p) == 1
        return qdt0
    end
    months = round(Int, (value(qdt0) - FIRST_VALUE) / _DAYS_IN_MONTH_F)
    return qdt0 - Dates.Month(mod(months, value(p)))
end

function Base.floor(qdt::QDate, p::Day)
    value(p) < 1 && throw(DomainError())
    days = value(qdt) - FIRST_VALUE
    return QDate(UTD(days - mod(days, value(p)) + FIRST_VALUE))
end


text:
using DICOM

# todo


text:
function copy(m::MRG32k3a)
    MRG32k3a(copy(m.Cg), copy(m.Bg), copy(m.Ig))
end


text:
# HSA error handling utilities

export HSAError

struct HSAError <: Exception
    code::Status
end

"""
    description(err::HSAError)

Gets the string description of an error code.
"""
function description(err::HSAError)
    str_ref = Ref{Ptr{Int8}}()
    HSA.status_string(err.code, str_ref)
    unsafe_string(reinterpret(Cstring, str_ref[]))
end

function Base.showerror(io::IO, err::HSAError)
    print(io, "HSA error (code #$(Int(err.code)), $(description(err)))")
end

Base.show(io::IO, err::HSAError) = print(io, "HSAError($(Int(err.code)), $(description(err)))")

"""
    check(result)

Check the return code of an HSA call, and throw an `HSAError` exception
if it is not the success code.
"""
function check(result)
    if result != HSA.STATUS_SUCCESS
        throw(HSAError(result))
    end
end


text:
module TestErrors

using Distributed: RemoteChannel
using Folds
using FoldsPreview
using Test
using Transducers

using ..Utils: @test_error, foreach_executor, with_racer, print_error

function test_early_error_in_body(ex)
    th = IntervalThrottle(Inf)
    local reach_body::Bool = false
    err = @test_error with_preview(identity, ex, th) do previewer
        reach_body = true
        error("ERROR FROM BODY")
    end
    @test reach_body
    @test any(err) do e
        occursin("ERROR FROM BODY", sprint(showerror, e))
    end
end

function test_early_error_in_body()
    foreach_executor(test_early_error_in_body)
end

function test_early_error_in_on_preview(ex)
    th = IntervalThrottle(0)
    err = nothing
    ok = with_racer(:reach_on_preview) do win
        function on_preview(_)
            win()
            error("ERROR FROM ON_PREVIEW")
        end
        err = @test_error with_preview(on_preview, ex, th) do previewer
            Folds.collect(previewer(1:10), ex)
        end
    end

    @test err isa CompositeException
    if !ok
        print_error(err)
    end

    @test any(err) do e
        occursin("ERROR FROM ON_PREVIEW", sprint(showerror, e))
    end
end

function test_early_error_in_on_preview()
    foreach_executor(test_early_error_in_on_preview)
end

function test_late_error_in_on_preview(ex)
    counter = Ref(0)
    barrier = if ex isa DistributedEx
        RemoteChannel(() -> Channel{Nothing}())
    else
        Channel{Nothing}()
    end
    err = nothing
    ok = with_racer(:iter1, :iter2, :iter3) do win1, win2, win3
        function on_preview(_)
            i = counter[] += 1
            @debug "`on_preview`: i = $i"
            if i == 1
                win1()
            elseif i == 2
                win2()
            elseif i == 3
                win3()
            else
                put!(barrier, nothing)
                error("ERROR FROM ON_PREVIEW")
            end
        end
        err = @test_error with_preview(on_preview, ex, 0) do previewer
            itr = (
                begin
                    # The first iteration is required for starting `on_preview`.
                    # So, block on the second:
                    x == 2 && take!(barrier)
                    x
                end
                for x in 1:10
            )
            Folds.collect(previewer(itr), ex)
        end
    end

    @test err isa CompositeException
    if !ok
        print_error(err)
    end

    @test any(err) do e
        occursin("ERROR FROM ON_PREVIEW", sprint(showerror, e))
    end
end

function test_late_error_in_on_preview()
    foreach_executor(test_late_error_in_on_preview)
end

end  # module


text:
module Sample

using ..Space, ..Util
import Random: AbstractRNG
import Random

export defrandalg, condomegasample, randsample, omegarandsample

"Default sampling algorithm"
function defrandalg end

function omegarandsample end
function omegarandsample1 end

"""
`n` random samples from `x` such that `yₓ(ω)` is true where `yₓ` are conditions of `x
using `alg` algorithm

```
x = 1 ~ StdNormal()
randsample(x |ᶜ x >ₚ 2.0, 3; alg = RejectionSample)
'''
"""
randsample(rng::AbstractRNG, x, n; alg = defrandalg(x), ΩT = defΩ(), kwargs...) =
  randsample(rng, ΩT, x, n, alg; kwargs...)

"""
Sample `n` `ω::AbstractΩ` according to function `logenergy(ω)` using inference algorithm `alg`

# Arguments
`rng` : random number generator
`logenergy` :

Returns
Vector of `n` samples
"""
omegarandsample(rng::AbstractRNG, logenergy, n; ΩT = defΩ(), alg, kwargs...) =
  omegarandsample(rng, ΩT, logenergy, n, alg; kwargs...)

omegarandsample(logenergy, n; ΩT = defΩ(), alg, kwargs...) =
  omegarandsample(Random.GLOBAL_RNG, logenergy, n; ΩT, alg, kwargs...)

# Convenience methods

"`randsample` using GLOBAL_RNG"
randsample(x, n::Integer; kwargs...) =
  randsample(Random.GLOBAL_RNG, x, n; kwargs...)

randsample(x; kwargs...) =
  first(randsample(x, 1, kwargs...))

randsample(rng::AbstractRNG, xs::Tuple, args...; kwargs...) =
  randsample(rng, ω -> mapf(ω, xs), args...; kwargs...)

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect