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

using Markdown
using InteractiveUtils

# ╔═╡ 065a1e67-6b63-43df-9d6d-303af08d8434
begin
	using Pkg
	Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
    using PlutoUI
    PlutoUI.TableOfContents()
end

# ╔═╡ 034f7cc2-eacd-11eb-29e2-810a573295cf
using LinearAlgebra

# ╔═╡ 034f7c90-eacd-11eb-1b40-8722209689fa
using DataFrames # load the DataFrames package

# ╔═╡ 034f7d1a-eacd-11eb-1d48-591daf597cd6
using Random

# ╔═╡ 034f7d3a-eacd-11eb-1144-65f9acfcf24c
begin
	using Plots
	plotly()
end

# ╔═╡ 034f7d58-eacd-11eb-3189-bff16ec27739
using Distributions

# ╔═╡ 034f7d76-eacd-11eb-24cc-fd8d53d736c2
using Statistics

# ╔═╡ d6e5fe02-21a5-486c-a237-878be1d95439
begin
    using MLCourse
    MLCourse.list_notebooks(@__FILE__)
end

# ╔═╡ f3508747-da29-47c9-a98e-22ea15caaf2f
Markdown.parse("Hi and welcome to an introduction to the Julia programming language.

$(haskey(ENV, "html_export") ? "This page was created with an interactive [Pluto notebook](https://plutojl.org/). If you want to run it locally on your machine, follow the instructions [here](https://github.com/jbrea/MLCourse)." :
"This is an interactive Pluto notebook. You can create new cells by clicking any \"+\" above or below existing cells. In these cells you can write code and run it by clicking the play button or [Shift] + [Enter] on your keyboard (or [Ctrl] + [Enter]).
The output gets displayed above the cell. Have a look at [this website](https://plutojl.org/) if you want to learn more about Pluto notebooks.

To get help, please open the Live docs at the bottom right of this page and click on the code you want to get help for. For example, click on the ÷ symbol in the 8th cell below.")

If you want to learn more about Julia visit [julialang.org](https://julialang.org).
In the following sections you find links to different chapters of the [manual](https://docs.julialang.org/en/v1/). There is also this [cheat sheet](https://juliadocs.github.io/Julia-Cheat-Sheet/).")

# ╔═╡ d92c0729-c7ab-42cc-b713-30f00e237833
md"# Mathematical Operations

Here is the [link](https://docs.julialang.org/en/v1/manual/mathematical-operations/) to the respective section in the Julia manual."

# ╔═╡ b8630ee9-d2c0-4cd0-b3f7-77f66e997a80

# ╔═╡ 034f7b46-eacd-11eb-1aa2-f380d08206cc
1 + 2

# ╔═╡ 034f7b64-eacd-11eb-1618-6de6def66405
2 * 3

# ╔═╡ 034f7b6e-eacd-11eb-1733-650ad4da2730
3 ^ 2

# ╔═╡ 034f7b78-eacd-11eb-003a-5b5801c19d80
3 == 9 / 3 # test for equality

# ╔═╡ 034f7b6e-eacd-11eb-184a-d504fab31681
4 / 3 # compute the fraction as a floating point number

# ╔═╡ fb3677a9-1610-49dd-96e4-3bd7badc58db
4 // 3 # leave the fraction as a fraction

# ╔═╡ 884f568b-079c-47d6-a597-34b67925f195
4//3 * 3//2 * 1//7 # simplify fractions

# ╔═╡ 034f7b84-eacd-11eb-2836-11cebedb05b9
5 ÷ 2 # Write this symbol with \div [TAB]. Try it in a new cell!

# ╔═╡ 52efdece-3c21-43d4-a826-3b96e3a23e60
md"# Variables and Simple Functions

See also [here](https://docs.julialang.org/en/v1/manual/variables/)."

# ╔═╡ 034f7ba0-eacd-11eb-3fc1-9dca0e072755
x = 2

# ╔═╡ 034f7ba0-eacd-11eb-23ae-cd327563ec60
y = 4.5

# ╔═╡ 034f7ba0-eacd-11eb-3347-49deb13761d9
z = x * y

# ╔═╡ 121b3dbb-f86c-44b8-8b0f-f6b4d0b9992f
md"In the following cells we apply simple built-in functions to the variables defined above. If you change the values of these variables, the output of the cells below will also change, thanks to the reactivity of this notebook. Try it by assigning another value to `x` in the cell above."

# ╔═╡ 034f7baa-eacd-11eb-33f2-5f619e4a7832
sin(x)

# ╔═╡ 034f7bb2-eacd-11eb-006b-7b935e35af66
exp(x)

# ╔═╡ 034f7bb2-eacd-11eb-214b-cbfc83965c26
log(x)

# ╔═╡ ec1dd4c4-d61d-4fdb-b576-3fe3b059eb3a
typeof(2)

# ╔═╡ 1a3387ee-b59d-49d7-9cca-2a8bd04eaaad
typeof(2.0)

# ╔═╡ 034f7b96-eacd-11eb-0c81-ad33210de120
typeof(4//3)

# ╔═╡ 034f7b96-eacd-11eb-2ee1-cd8e6813a726
typeof(2f0)

# ╔═╡ 0f4e3f45-71dd-41b2-a0b3-c55a04ace61f
md"In the following cells you see 3 different ways to define custom functions. See also [here](https://docs.julialang.org/en/v1/manual/functions/), and, in particular, [here](https://docs.julialang.org/en/v1/manual/functions/#man-anonymous-functions) for anonymous functions."

# ╔═╡ 034f7bbe-eacd-11eb-2e95-a93ad2d5e867
myfunc(x, y) = 2 * x + 3 - y # creates a new function called `myfunc`

# ╔═╡ 034f7bbe-eacd-11eb-2c36-d3949367ce5a
myfunc(x, y)

# ╔═╡ e0ec7ba3-4858-4572-994d-953caa98080b
anonfunc = x -> x^2 # creates an anonymous function called `anonfunc`

# ╔═╡ 04da39a1-e3bc-4086-b555-6cea2cf1689d
anonfunc(4)

# ╔═╡ 034f7bc8-eacd-11eb-0396-5b743b505b2a
function longfunc(x, y)
    tmp = x + y
    tmp ^ 2 - 0.25 # the result of the last line is the output
end

# ╔═╡ 034f7bc8-eacd-11eb-31f9-c119d1dfb8a5
longfunc(x, y)

# ╔═╡ ce3c7687-c7a9-4b5d-9a55-51514ed90fd9
md"We will sometimes use the pipe operator to chain functions."

# ╔═╡ 034f7bd2-eacd-11eb-28f3-b707012e1db0
x |> sin |> exp |> log10

# ╔═╡ 034f7bd2-eacd-11eb-1789-7380a89d5339
log10(exp(sin(x)))

# ╔═╡ 034f7cd6-eacd-11eb-0f8c-11179c3ab480
md" # Strings and Markdown

See also [here](https://docs.julialang.org/en/v1/manual/strings/)
and [here](https://docs.julialang.org/en/v1/stdlib/Markdown/).
"

# ╔═╡ 034f7cd6-eacd-11eb-0f8c-11179c3ab48f
s1 = "Hi! "

# ╔═╡ 034f7ce0-eacd-11eb-3c49-e39b8d3e917f
s2 = "Welcome to Machine Learning"

# ╔═╡ 034f7ce0-eacd-11eb-18a1-bd2d1849b9bc
s1 * s2

# ╔═╡ 034f7cec-eacd-11eb-3262-bdc0f2e13a3a
s1 ^ 3 * s2

# ╔═╡ 034f7cec-eacd-11eb-27da-9bb861643df2
split(s2)

# ╔═╡ 832ce046-ab34-4623-9d02-793dc5748208
md"Sometimes we want to paste something computed into a string. This is called [interpolation](https://docs.julialang.org/en/v1/manual/strings/#string-interpolation)."

# ╔═╡ 91b0de0c-221c-468e-9077-5233aedeb4ca
answer = 42

# ╔═╡ 0c7ae3ec-a40f-4e3e-85d8-3f79f091e3b4
"The answer is $answer." # we interpolate with the dollar sign

# ╔═╡ 732ae082-b90f-4108-9872-36077fb2c54c
md"Whenever you see just text here, it is actually the output of a markdown cell.
Markdown cells can be written in two ways: either as a markdown string `md\"this is a markdown string\"` (note the trailing `md`) or as `Markdown.parse(\"this is a markdown string\")`.
You can look at the code of this cell by toggling its visibility with the \"eye\"
button on the left. For you to take notes, it may be useful to know some of the
markdown features. We look at some in the cell below.
"

# ╔═╡ 6bfb60b4-964e-46ec-adc2-ec9e9a0f8158
md"#### This is a fourth level header

[This is a link](https://docs.julialang.org/en/v1/)

- here we have a list
- of multiple items
   1. with different levels.
   1. The number you put here doesn't matter.
   0. The markdown interpreter knows how to count.

Below you see a Julia code block:

```julia
my_fancy_function(x) = x^x # here I define my function
my_fancy_function(3) == 27
```

and here is the same thing in python:

```python
def my_fancy_function(x): # here I define my function
    return x**x

my_fancy_function(3) == 27
```

This is a formula ``E = mc^2``.

And below we have some nicely colored text boxes.

!!! note

    Sometimes I just note random stuff.

!!! tip

    It is useful to know markdown.

!!! warning \"Don't forget the exercises!\"

    No pain, no gain!

!!! danger \"Machine Learning Can Be Addictive\"

    Rumours say that in the past there were students that liked their first
    Machine Learning course so much that they decided to pursue a career in
    Machine Learning.
"

# ╔═╡ 0efb1f1a-8051-4843-aa68-1f4ad8618f92
md"# Vectors, Matrices and Arrays

See also [here](https://docs.julialang.org/en/v1/manual/arrays/)."

# ╔═╡ 034f7bdc-eacd-11eb-1999-9b62efdb4128
v = [1, 2, 3] # this is a vector

# ╔═╡ 034f7cba-eacd-11eb-2f1a-836ccc185dd9
push!(v, 12)

# ╔═╡ 034f7cba-eacd-11eb-3901-2bc3d0ff0e39
pop!(v)

# ╔═╡ 50d83a01-f4c5-45c1-9c81-ed71d8313642
md"If you wonder about the ! in the function names: have a look [here](https://docs.julialang.org/en/v1/manual/style-guide/#bang-convention)."

# ╔═╡ 034f7bfa-eacd-11eb-1839-371285e2c6ed
w = [1 2 3] # this is a matrix with 1 row and 3 columns

# ╔═╡ 034f7be4-eacd-11eb-1616-2df92f244e73
v[1] # access to the first element of the vector

# ╔═╡ 034f7bf0-eacd-11eb-0b3f-69f172657051
v[2] = 4 # assign value 4 to the second element of the vector

# ╔═╡ 034f7bf0-eacd-11eb-012e-997132cd41c6
v # here we check that the second element is indeed 4 now

# ╔═╡ 034f7bfa-eacd-11eb-2bf8-c795453a8959
m = [1 2 3
     4 5 6
     7 8 9] # this is a 3 x 3 matrix.

# ╔═╡ 034f7c04-eacd-11eb-3448-832bd75c3356
m[3, 2] # accessing the element in the 3 row, second column

# ╔═╡ 9c34e4ef-4eba-459f-a7ff-b20708155e65
m' # transpose the matrix

# ╔═╡ e576bf44-9cc1-4053-b40f-b82c8c39c89d
md"Any function can be applied element-wise to arrays using the [dot syntax](https://docs.julialang.org/en/v1/manual/functions/#man-vectorized)."

# ╔═╡ 034f7c0e-eacd-11eb-114f-9dbad3e74dd6
sin.(v)

# ╔═╡ 034f7c0e-eacd-11eb-1b73-27539a68f26d
exp.(m)

# ╔═╡ 034f7c16-eacd-11eb-0d0b-518b28558fda
v .+ 3

# ╔═╡ 034f7c16-eacd-11eb-2cb0-8b1198bfcb78
myfunc.(v, v)

# ╔═╡ 034f7ccc-eacd-11eb-1613-cbdf2386b60b
m * v  # matrix vector multiplication

# ╔═╡ 034f7ccc-eacd-11eb-2ab4-599e0ed86304
v' * v # inner product

# ╔═╡ 034f7cd6-eacd-11eb-37bb-f1930f2f2772
eigvals(m) # compute the Eigen values

# ╔═╡ c997a07f-bc89-4c29-9945-62046b6889d2
md"To concatenate vectors and matrices we have the following syntax (see also [here](https://docs.julialang.org/en/v1/manual/arrays/#man-array-concatenation))."

# ╔═╡ 07fd5baf-ae0b-4e74-97a6-5a6e14644122
[[1, 2, 3]; [77, 88, 99]] # the semicolon concatenates the two vectors

# ╔═╡ 11422ca4-c096-4167-aa53-de312615348e
md"The same can be achieved with"

# ╔═╡ 0ace8c81-e481-4323-9b20-116de412330a
vcat([1, 2, 3], [77, 88, 99])

# ╔═╡ 5d09ccbf-5a51-41f7-a7ef-84cf2cc5fa1f
[[1, 2, 3] [77, 88, 99]] # alternatively you can use hcat([1, 2, 3], [77, 88, 99])

# ╔═╡ 9fd6e8c1-7b84-4d76-abef-9e3da4a56dc7
md"Sometimes it is easier to construct a vector or matrix with [comprehension](https://docs.julialang.org/en/v1/manual/arrays/#man-comprehensions)."

# ╔═╡ e33c0eea-4f5c-4fc8-800d-fa16335eedd0
[i^2 for i in [3, 4, 9, 7]]

# ╔═╡ 315e5a3b-dc28-4c3b-8fc2-3442cb7590be
[i^2 for i in [3, 4, 9, 7] if isodd(i)]

# ╔═╡ 45deab46-9a8c-4d22-9460-980a3900396f
["$i, $j" for i in 1:5, j in 1:5] # see string interpolation above for the $

# ╔═╡ b91eda3c-95cb-4372-9c6b-03fa6a6b02b2
md"# Ranges

See also [here](https://docs.julialang.org/en/v1/base/math/#Base.::)."

# ╔═╡ 034f7c16-eacd-11eb-3877-6ba06a933698
collect(1:10) # 1:10 is a range; the function `collect` turns it into a vector

# ╔═╡ 034f7c22-eacd-11eb-3ab7-b30813d11e94
r1 = 1:.1:10 # a range from 1 to 10 with step-size 0.1

# ╔═╡ 034f7c2c-eacd-11eb-1aae-97b030af71e9
r3 = 10:-.1:1 # a range from 10 to 1 with step-size -0.1

# ╔═╡ 034f7c2c-eacd-11eb-1f56-052874c5c6ab
collect(r3)

# ╔═╡ 034f7c36-eacd-11eb-0188-cd307154555c
exp.(r1) # apply exp to the range

# ╔═╡ 034f7c40-eacd-11eb-26d9-956b65235a10
v[2:3] # you can use ranges to access multiple elements of a vector or matrix

# ╔═╡ 4237f5df-208d-44c4-90d7-db0efbe7b710
v[[1, 3]] # a vector of indices can also be used

# ╔═╡ 034f7c40-eacd-11eb-1beb-bb9e9bb325db
v[2:end]

# ╔═╡ 034f7c4a-eacd-11eb-0de4-537297318c85
v[begin:2]

# ╔═╡ 034f7c4a-eacd-11eb-03ed-bfc7b532501b
m[:, 2] # the colon : alone is a short-hand for begin:end

# ╔═╡ 034f7c54-eacd-11eb-2918-ef8f2befddc1
m[1:2, 2:3]

# ╔═╡ 034f7c54-eacd-11eb-2fe9-5330d3cd3cd5
function is_larger_than_10(x)
	result = ""
    for i in 1:10 # ranges are useful to construct for loops
        if x < i
            return result * " is smaller than 10"
        else
            result *= "|"
        end
    end
    result * "... is larger than 10"
end

# ╔═╡ 034f7c5e-eacd-11eb-24c6-d7e7213a7a5e
is_larger_than_10(4)

# ╔═╡ 034f7c5e-eacd-11eb-1ff7-55f09fb796fe
is_larger_than_10(78)

# ╔═╡ d9ad1d6d-a3b5-4a83-bda9-05ee43ded9d6
md"# Tuples

See also [here](https://docs.julialang.org/en/v1/manual/functions/#Tuples)."

# ╔═╡ 034f7c68-eacd-11eb-1e11-d77379c677c8
t = (2, 3, 4., "bla")

# ╔═╡ 034f7c68-eacd-11eb-12e6-51022c54f6ed
typeof(t)

# ╔═╡ 034f7c72-eacd-11eb-28a5-d7d6f3b9f75a
t[1]

# ╔═╡ 034f7c72-eacd-11eb-1cc2-91eaa2de5e25
t[4]

# ╔═╡ 034f7c7c-eacd-11eb-1aaa-297c951dd8cd
nt = (year = 1789,
	  event = "French Revolution",
	  slogan = "liberté, égalité, fraternité")

# ╔═╡ 034f7c88-eacd-11eb-28f2-c1c7e05279c6
typeof(nt)

# ╔═╡ 034f7c88-eacd-11eb-14fe-b3d49ff91622
nt.year

# ╔═╡ 034f7c90-eacd-11eb-0d1c-5dbd7b62d55c
nt.slogan

# ╔═╡ df8588cb-400b-4423-9766-e6cac2c9717a
md"# DataFrames

DataFrames are tables with named columns.
We will constantly use DataFrames to organize data.
See also [the documentation of the DataFrames.jl package](https://dataframes.juliadata.org/stable/)."

# ╔═╡ 034f7c9a-eacd-11eb-2e85-29c5bffbc39a
df = DataFrame(year = [1789, 1863],
               event = ["French Revolution", "Foundation of Red Cross"],
               slogan = ["liberté, égalité, fraternité", missing])

# ╔═╡ 034f7cae-eacd-11eb-39dc-ab39e97cc716
push!(df, [1789,
           "George Washington becomes first president of the USA",
           "Deeds, not Words."]) # append a new row to the already created dataframe

# ╔═╡ 034f7c9a-eacd-11eb-3aee-0bab7ab5d502
df.year

# ╔═╡ 034f7ca4-eacd-11eb-2f3b-c1dd9be6ad84
df[:, 1] # we can also access columns like for matrices

# ╔═╡ 034f7ca4-eacd-11eb-054e-1f930c99d304
dropmissing(df) # removes all rows with missing entries.

# ╔═╡ 034f7cae-eacd-11eb-38c0-57df4ba09992
df[1, :] # show the first row

# ╔═╡ 034f7cc2-eacd-11eb-0f0f-ff16d1ac454b
df[df.year .== 1789, [:year, :slogan]] # select all rows where year == 1789 and columns year and slogan

# ╔═╡ 034f7cf4-eacd-11eb-302f-55b164df3b90
md"# Random Numbers

See also [here](https://docs.julialang.org/en/v1/stdlib/Random/).
"

# ╔═╡ 034f7cf4-eacd-11eb-302f-55b164df3b9f
rand() # a sample from the uniform distribution over the interval [0, 1).

# ╔═╡ 034f7cf4-eacd-11eb-003c-133f34cca508
rand() # different number than in the cell above.

# ╔═╡ 034f7cfe-eacd-11eb-2514-4b6fdbc42e8f
rand(10)

# ╔═╡ 034f7d08-eacd-11eb-2694-51c1baa997e6
rand(3, 3)

# ╔═╡ 034f7d12-eacd-11eb-1673-c94969427bbf
randn() # a sample from the normal distribution with mean 0 and variance 1

# ╔═╡ 034f7d12-eacd-11eb-261c-a5de2f06f420
rand((:bla, "bli", 3, 1.2))

# ╔═╡ 034f7d1a-eacd-11eb-1e34-21d5eef62940
Random.seed!(123) # sets the seed of the random number generator to 123

# ╔═╡ 034f7d26-eacd-11eb-066c-538e5b614536
rand()

# ╔═╡ 034f7d30-eacd-11eb-081f-970c4daee184
Random.seed!(123) # resets the seed

# ╔═╡ 034f7d30-eacd-11eb-0b09-7b3b9ac47570
rand() # same result as when `rand()` was called last time.

# ╔═╡ 034f7d3a-eacd-11eb-3fbe-6dc4ab1efd50
md"# Plotting

See also the [documentation of the Plots.jl package](http://docs.juliaplots.org/latest/)
"

# ╔═╡ 034f7d3a-eacd-11eb-3fbe-6dc4ab1efd56
plot(1:4, rand(4))

# ╔═╡ 034f7d3a-eacd-11eb-1ef3-4f89b8e5aee9
scatter(rand(100), rand(100))

# ╔═╡ 034f7d44-eacd-11eb-0ab2-2f0d2caed9d6
with_terminal(plotattr) # get some help on plotting attributes

# ╔═╡ 034f7d4c-eacd-11eb-3960-d10d1c16075d
with_terminal(plotattr, :Series)

# ╔═╡ 034f7d4c-eacd-11eb-3841-97ca37e34c9f
with_terminal(plotattr, :Series, "label")

# ╔═╡ 034f7d58-eacd-11eb-123f-9b563a7ae27e
scatter(rand(100), rand(100), label = "my data", xlabel = "X1", ylabel = "X2")

# ╔═╡ 4ca1cd98-a24b-491f-b9ac-28a7daf96d50
md"To compose figures with multiple elements, you can use the `!` version of the plotting functions"

# ╔═╡ 65e73d27-3709-4894-87ab-432343d314c6
begin # we use begin to write a cell with multiple lines
	scatter(1:50, sqrt.(1:50) .+ randn(50), label = "data points")
	plot!(sqrt, color = :red, label = "square root function")
end

# ╔═╡ 9f69c2ab-2940-40ce-b164-07bfa6cc1697
md"Instead of `begin`-`end`-blocks we will also use `let`-`end`-blocks. Variables computed in `let`-blocks are local to that cell, whereas variables defined in `begin`-blocks are global and can be accessed by other cells. (With the `;` at the end of the cell we suppress the output)."

# ╔═╡ e7ed9e83-7f17-4d9a-a49b-e87bc0747035
begin
	a_global_variable = 17
	another_global_variable = "Hi"
end;

# ╔═╡ 0d66592a-d098-4ab5-9323-503db8fc73f6
a_global_variable * 3

# ╔═╡ a9a3f29e-4931-4464-9bc5-af6598815f24
"$another_global_variable there!"

# ╔═╡ 6e65247f-b285-49e5-a16e-31b4541225bf
let
	a_local_variable = 17
	another_local_variable = "Hi"
end;

# ╔═╡ fcd77620-f18a-499b-87c4-fffaa876e59f
a_local_variable * 3

# ╔═╡ 034f7d62-eacd-11eb-3840-51c0622eb8a0
md"# Distributions

See also the [documentation of the Distributions.jl package](https://juliastats.org/Distributions.jl/stable/)
"

# ╔═╡ 034f7d62-eacd-11eb-3840-51c0622eb8a9
d1 = Normal(3, 2.5)

# ╔═╡ 034f7d6c-eacd-11eb-0eb6-277bf88c13d3
rand(d1) # draw a sample from this normal distribution

# ╔═╡ 034f7d6c-eacd-11eb-2792-d773a2c77cb6
rand(d1, 10) # draw 10 samples from this distriution

# ╔═╡ 034f7d6c-eacd-11eb-0dc0-3bd2b5a24d38
mean(rand(d1, 10^5))

# ╔═╡ 034f7d76-eacd-11eb-141d-0f90eef73f08
std(rand(d1, 10^5)) # standard deviation

# ╔═╡ 034f7d7e-eacd-11eb-0945-2dda6be01349
d2 = Bernoulli(.7)

# ╔═╡ 034f7d7e-eacd-11eb-2a1b-933e1231c220
rand(d2, 10)

# ╔═╡ 034f7d8a-eacd-11eb-1dc0-d99d481ab6fa
md"# Other Packages

Here are some other popular Julia packages. They are not relevant for this course, but maybe you want to explore them at some point.
* [Symbolics](https://github.com/JuliaSymbolics/Symbolics.jl)
* [DifferentialEquations](https://github.com/SciML/DifferentialEquations.jl)
* [Optim](https://github.com/JuliaNLSolvers/Optim.jl)
* [Images](https://github.com/JuliaImages/Images.jl)
* [Turing](https://github.com/TuringLang/Turing.jl)
"

# ╔═╡ 4a03cfae-9876-4cf0-a498-d750853191cb
md"""# Exercises

1. Create a data frame with 3 columns named A, B and C.
   1. Column A contains 5 random numbers sampled from a Bernoulli distribution with rate 0.3, column B contains 5 random numbers from the uniform distribution over the interval [0, 1), and column C contains 5 samples from the set `(:hip, :hop)`.
   2. Create a vector whose i'th element contains the sum of the i'th entries of columns A and B of the data frame created in 1.
   3. Select all rows with `:hop` in column C and display the resulting data frame.
1.
   1. Use comprehension (see "Vectors, Matrices, Arrays") to create a vector where with all numbers of the form ``x^y`` with ``x=1, \ldots, 10``, ``y = 2, \ldots, 7`` and ``y > x``.
   2. Compute the sum of the square root of this numbers.
2.
   1. Write a function that returns the smallest entry of a vector (without using the built-in function `minimum`, `argmin` or `findmin`).
   2. Test your function on a vector of 10 randomly sampled integers in the range 1 to 100.
3.
   1. Plot the `cos` function on the interval 0 to 4π. Hint: type `\pi + [Tab]` to enter the symbol π. To learn how to place custom tick labels on the x-axis, type `xticks` in a cell and open the \"Live docs\" at the bottom-right.
   2. Add a scatter plot with 100 points whose ``x`` and ``y`` coordinates are randomly sampled from the interval ``[0, 1)`` on top of the figure with the cosine.
"""

# ╔═╡ 0314376e-ff8c-4ad0-8a4b-f94f04f31f2c
MLCourse.footer()

# ╔═╡ Cell order:
# ╟─f3508747-da29-47c9-a98e-22ea15caaf2f
# ╟─d92c0729-c7ab-42cc-b713-30f00e237833
# ╠═b8630ee9-d2c0-4cd0-b3f7-77f66e997a80
# ╠═034f7b46-eacd-11eb-1aa2-f380d08206cc
# ╠═034f7b64-eacd-11eb-1618-6de6def66405
# ╠═034f7b6e-eacd-11eb-1733-650ad4da2730
# ╠═034f7b78-eacd-11eb-003a-5b5801c19d80
# ╠═034f7b6e-eacd-11eb-184a-d504fab31681
# ╠═fb3677a9-1610-49dd-96e4-3bd7badc58db
# ╠═884f568b-079c-47d6-a597-34b67925f195
# ╠═034f7b84-eacd-11eb-2836-11cebedb05b9
# ╟─52efdece-3c21-43d4-a826-3b96e3a23e60
# ╠═034f7ba0-eacd-11eb-3fc1-9dca0e072755
# ╠═034f7ba0-eacd-11eb-23ae-cd327563ec60
# ╠═034f7ba0-eacd-11eb-3347-49deb13761d9
# ╟─121b3dbb-f86c-44b8-8b0f-f6b4d0b9992f
# ╠═034f7baa-eacd-11eb-33f2-5f619e4a7832
# ╠═034f7bb2-eacd-11eb-006b-7b935e35af66
# ╠═034f7bb2-eacd-11eb-214b-cbfc83965c26
# ╠═ec1dd4c4-d61d-4fdb-b576-3fe3b059eb3a
# ╠═1a3387ee-b59d-49d7-9cca-2a8bd04eaaad
# ╠═034f7b96-eacd-11eb-0c81-ad33210de120
# ╠═034f7b96-eacd-11eb-2ee1-cd8e6813a726
# ╟─0f4e3f45-71dd-41b2-a0b3-c55a04ace61f
# ╠═034f7bbe-eacd-11eb-2e95-a93ad2d5e867
# ╠═034f7bbe-eacd-11eb-2c36-d3949367ce5a
# ╠═e0ec7ba3-4858-4572-994d-953caa98080b
# ╠═04da39a1-e3bc-4086-b555-6cea2cf1689d
# ╠═034f7bc8-eacd-11eb-0396-5b743b505b2a
# ╠═034f7bc8-eacd-11eb-31f9-c119d1dfb8a5
# ╟─ce3c7687-c7a9-4b5d-9a55-51514ed90fd9
# ╠═034f7bd2-eacd-11eb-28f3-b707012e1db0
# ╠═034f7bd2-eacd-11eb-1789-7380a89d5339
# ╟─034f7cd6-eacd-11eb-0f8c-11179c3ab480
# ╠═034f7cd6-eacd-11eb-0f8c-11179c3ab48f
# ╠═034f7ce0-eacd-11eb-3c49-e39b8d3e917f
# ╠═034f7ce0-eacd-11eb-18a1-bd2d1849b9bc
# ╠═034f7cec-eacd-11eb-3262-bdc0f2e13a3a
# ╠═034f7cec-eacd-11eb-27da-9bb861643df2
# ╟─832ce046-ab34-4623-9d02-793dc5748208
# ╠═91b0de0c-221c-468e-9077-5233aedeb4ca
# ╠═0c7ae3ec-a40f-4e3e-85d8-3f79f091e3b4
# ╟─732ae082-b90f-4108-9872-36077fb2c54c
# ╠═6bfb60b4-964e-46ec-adc2-ec9e9a0f8158
# ╟─0efb1f1a-8051-4843-aa68-1f4ad8618f92
# ╠═034f7bdc-eacd-11eb-1999-9b62efdb4128
# ╠═034f7cba-eacd-11eb-2f1a-836ccc185dd9
# ╠═034f7cba-eacd-11eb-3901-2bc3d0ff0e39
# ╟─50d83a01-f4c5-45c1-9c81-ed71d8313642
# ╠═034f7bfa-eacd-11eb-1839-371285e2c6ed
# ╠═034f7be4-eacd-11eb-1616-2df92f244e73
# ╠═034f7bf0-eacd-11eb-0b3f-69f172657051
# ╠═034f7bf0-eacd-11eb-012e-997132cd41c6
# ╠═034f7bfa-eacd-11eb-2bf8-c795453a8959
# ╠═034f7c04-eacd-11eb-3448-832bd75c3356
# ╠═9c34e4ef-4eba-459f-a7ff-b20708155e65
# ╟─e576bf44-9cc1-4053-b40f-b82c8c39c89d
# ╠═034f7c0e-eacd-11eb-114f-9dbad3e74dd6
# ╠═034f7c0e-eacd-11eb-1b73-27539a68f26d
# ╠═034f7c16-eacd-11eb-0d0b-518b28558fda
# ╠═034f7c16-eacd-11eb-2cb0-8b1198bfcb78
# ╠═034f7ccc-eacd-11eb-1613-cbdf2386b60b
# ╠═034f7ccc-eacd-11eb-2ab4-599e0ed86304
# ╠═034f7cc2-eacd-11eb-29e2-810a573295cf
# ╠═034f7cd6-eacd-11eb-37bb-f1930f2f2772
# ╟─c997a07f-bc89-4c29-9945-62046b6889d2
# ╠═07fd5baf-ae0b-4e74-97a6-5a6e14644122
# ╟─11422ca4-c096-4167-aa53-de312615348e
# ╠═0ace8c81-e481-4323-9b20-116de412330a
# ╠═5d09ccbf-5a51-41f7-a7ef-84cf2cc5fa1f
# ╟─9fd6e8c1-7b84-4d76-abef-9e3da4a56dc7
# ╠═e33c0eea-4f5c-4fc8-800d-fa16335eedd0
# ╠═315e5a3b-dc28-4c3b-8fc2-3442cb7590be
# ╠═45deab46-9a8c-4d22-9460-980a3900396f
# ╟─b91eda3c-95cb-4372-9c6b-03fa6a6b02b2
# ╠═034f7c16-eacd-11eb-3877-6ba06a933698
# ╠═034f7c22-eacd-11eb-3ab7-b30813d11e94
# ╠═034f7c2c-eacd-11eb-1aae-97b030af71e9
# ╠═034f7c2c-eacd-11eb-1f56-052874c5c6ab
# ╠═034f7c36-eacd-11eb-0188-cd307154555c
# ╠═034f7c40-eacd-11eb-26d9-956b65235a10
# ╠═4237f5df-208d-44c4-90d7-db0efbe7b710
# ╠═034f7c40-eacd-11eb-1beb-bb9e9bb325db
# ╠═034f7c4a-eacd-11eb-0de4-537297318c85
# ╠═034f7c4a-eacd-11eb-03ed-bfc7b532501b
# ╠═034f7c54-eacd-11eb-2918-ef8f2befddc1
# ╠═034f7c54-eacd-11eb-2fe9-5330d3cd3cd5
# ╠═034f7c5e-eacd-11eb-24c6-d7e7213a7a5e
# ╠═034f7c5e-eacd-11eb-1ff7-55f09fb796fe
# ╟─d9ad1d6d-a3b5-4a83-bda9-05ee43ded9d6
# ╠═034f7c68-eacd-11eb-1e11-d77379c677c8
# ╠═034f7c68-eacd-11eb-12e6-51022c54f6ed
# ╠═034f7c72-eacd-11eb-28a5-d7d6f3b9f75a
# ╠═034f7c72-eacd-11eb-1cc2-91eaa2de5e25
# ╠═034f7c7c-eacd-11eb-1aaa-297c951dd8cd
# ╠═034f7c88-eacd-11eb-28f2-c1c7e05279c6
# ╠═034f7c88-eacd-11eb-14fe-b3d49ff91622
# ╠═034f7c90-eacd-11eb-0d1c-5dbd7b62d55c
# ╟─df8588cb-400b-4423-9766-e6cac2c9717a
# ╠═034f7c90-eacd-11eb-1b40-8722209689fa
# ╠═034f7c9a-eacd-11eb-2e85-29c5bffbc39a
# ╠═034f7cae-eacd-11eb-39dc-ab39e97cc716
# ╠═034f7c9a-eacd-11eb-3aee-0bab7ab5d502
# ╠═034f7ca4-eacd-11eb-2f3b-c1dd9be6ad84
# ╠═034f7ca4-eacd-11eb-054e-1f930c99d304
# ╠═034f7cae-eacd-11eb-38c0-57df4ba09992
# ╠═034f7cc2-eacd-11eb-0f0f-ff16d1ac454b
# ╟─034f7cf4-eacd-11eb-302f-55b164df3b90
# ╠═034f7cf4-eacd-11eb-302f-55b164df3b9f
# ╠═034f7cf4-eacd-11eb-003c-133f34cca508
# ╠═034f7cfe-eacd-11eb-2514-4b6fdbc42e8f
# ╠═034f7d08-eacd-11eb-2694-51c1baa997e6
# ╠═034f7d12-eacd-11eb-1673-c94969427bbf
# ╠═034f7d12-eacd-11eb-261c-a5de2f06f420
# ╠═034f7d1a-eacd-11eb-1d48-591daf597cd6
# ╠═034f7d1a-eacd-11eb-1e34-21d5eef62940
# ╠═034f7d26-eacd-11eb-066c-538e5b614536
# ╠═034f7d30-eacd-11eb-081f-970c4daee184
# ╠═034f7d30-eacd-11eb-0b09-7b3b9ac47570
# ╟─034f7d3a-eacd-11eb-3fbe-6dc4ab1efd50
# ╠═034f7d3a-eacd-11eb-1144-65f9acfcf24c
# ╠═034f7d3a-eacd-11eb-3fbe-6dc4ab1efd56
# ╠═034f7d3a-eacd-11eb-1ef3-4f89b8e5aee9
# ╠═034f7d44-eacd-11eb-0ab2-2f0d2caed9d6
# ╠═034f7d4c-eacd-11eb-3960-d10d1c16075d
# ╠═034f7d4c-eacd-11eb-3841-97ca37e34c9f
# ╠═034f7d58-eacd-11eb-123f-9b563a7ae27e
# ╟─4ca1cd98-a24b-491f-b9ac-28a7daf96d50
# ╠═65e73d27-3709-4894-87ab-432343d314c6
# ╟─9f69c2ab-2940-40ce-b164-07bfa6cc1697
# ╠═e7ed9e83-7f17-4d9a-a49b-e87bc0747035
# ╠═0d66592a-d098-4ab5-9323-503db8fc73f6
# ╠═a9a3f29e-4931-4464-9bc5-af6598815f24
# ╠═6e65247f-b285-49e5-a16e-31b4541225bf
# ╠═fcd77620-f18a-499b-87c4-fffaa876e59f
# ╟─034f7d62-eacd-11eb-3840-51c0622eb8a0
# ╠═034f7d58-eacd-11eb-3189-bff16ec27739
# ╠═034f7d62-eacd-11eb-3840-51c0622eb8a9
# ╠═034f7d6c-eacd-11eb-0eb6-277bf88c13d3
# ╠═034f7d6c-eacd-11eb-2792-d773a2c77cb6
# ╠═034f7d6c-eacd-11eb-0dc0-3bd2b5a24d38
# ╠═034f7d76-eacd-11eb-24cc-fd8d53d736c2
# ╠═034f7d76-eacd-11eb-141d-0f90eef73f08
# ╠═034f7d7e-eacd-11eb-0945-2dda6be01349
# ╠═034f7d7e-eacd-11eb-2a1b-933e1231c220
# ╟─034f7d8a-eacd-11eb-1dc0-d99d481ab6fa
# ╟─4a03cfae-9876-4cf0-a498-d750853191cb
# ╟─d6e5fe02-21a5-486c-a237-878be1d95439
# ╟─065a1e67-6b63-43df-9d6d-303af08d8434
# ╟─0314376e-ff8c-4ad0-8a4b-f94f04f31f2c

max_stars_repo_name: ablaom/MLCourse

max_issues_repo_name: ablaom/MLCourse

max_forks_repo_name: ablaom/MLCourse


text:
"""
**Main simulation loop**

    simulate(parameters, biomass; start::Int64=0, stop::Int64=500, use::Symbol=:stiff)

This function takes two mandatory arguments:

- `p` is a `Dict` as returned by `make_parameters`
- `biomass` is an `Array{Float64, 1}` with the initial biomasses of every species

Internally, the function will check that the length of `biomass` matches
with the size of the network.

In addition, the function takes three optional arguments:

- `start` (defaults to 0), the initial time
- `stop` (defaults to 500), the final time
- `use` (defaults to `:stiff`), a hint to select the solver

The integration method is, by default, `:stiff`, and can be changed to
`:nonstiff`. This is because internally, this function used the
`DifferentialEquations` package to pick the most appropriate algorithm.

The `simulate` function returns a `Dict{Symbol, Any}`, with three
top-level keys:

- `:p`, the parameters that were given as input
- `:t`, the timesteps
- `:B`, an `Array{Float64, 2}` with the biomasses

The array of biomasses has one row for each timestep, and one column for
each species.
"""
function simulate(parameters, biomass; n_concentration::Vector{Float64}=rand(Float64, 2).*10, start::Int64=0, stop::Int64=500, use::Symbol=:nonstiff, cb_interp_points::Int64=100, extinction_threshold::Float64=1e-6)
  @assert stop > start
  @assert length(biomass) == size(parameters[:A],1)
  @assert length(n_concentration) == 2
  if parameters[:productivity] == :nutrients
      biomass = vcat(biomass, n_concentration)
  end

  @assert use ∈ vec([:stiff :nonstiff])
  alg = use == :stiff ? Rodas4(autodiff=false) : Tsit5()

  S = size(parameters[:A], 1)

  # Pre-allocate the timeseries matrix
  tspan = (float(start), float(stop))
  t_keep = collect(start:0.25:stop)

  # Perform the actual integration
  prob = ODEProblem(dBdt, biomass, tspan, parameters)

  ϵ = []

  function species_under_extinction_threshold(u, t, integrator)
    workingbm = deepcopy(u)
    sort!(ϵ)
    deleteat!(workingbm, unique(ϵ))
    cond = any(x -> x < extinction_threshold, workingbm) ? -0.0 : 1.0
    return cond
  end

  function remove_species!(integrator)
    u = integrator.u
    idϵ = findall(x -> x < extinction_threshold, u)
    for e in idϵ
        if !(e ∈ ϵ)
            u[e] = 0.0
            append!(ϵ,e)
        end
    end
    sort!(ϵ)
    nothing
  end

  function remove_species_and_update!(integrator)
    remove_species!(integrator)
    if parameters[:productivity] == :nutrients
      workingbm = deepcopy(integrator.u[1:end-2])
    else
      workingbm = deepcopy(integrator.u)
    end
    parameters = update_rewiring_parameters(parameters, workingbm, integrator.t)
  end

  function remove_target_and_update!(u, t, integrator)
    remove_species!(integrator)
    if parameters[:productivity] == :nutrients
      workingbm = deepcopy(integrator.u[1:end-2])
    else
      workingbm = deepcopy(integrator.u)
    end
    parameters = BioEnergeticFoodWebs.update_rewiring_parameters(parameters, workingbm, integrator.t)
  end

  cb = species_under_extinction_threshold
  affect_function = remove_species_and_update!
  if parameters[:rewire_method] == :ADBM
    if parameters[:adbm_trigger] == :interval
      Δt = parameters[:adbm_interval]
      cb1 = PeriodicCallback(affect_function, Δt)
    else
      cb1 = ContinuousCallback(cb, affect_function, interp_points = cb_interp_points)
    end
  else
    cb1 = ContinuousCallback(cb, affect_function, interp_points = cb_interp_points)
  end

  is_any_extinct = any(biomass .< extinction_threshold)
  if is_any_extinct
    cb2 = FunctionCallingCallback(remove_target_and_update!, funcat = [1.0])
    extinction_callback = CallbackSet(cb1, cb2)
  else
    extinction_callback = cb1
  end

  sol = solve(prob, alg, callback = extinction_callback, saveat=t_keep, dense=false, save_timeseries=false, force_dtmin=false)

  B = hcat(sol.u...)'

  if parameters[:productivity] == :nutrients
      output = Dict{Symbol,Any}(
      :p => parameters,
      :t => sol.t,
      :B => B[:,1:S],
      :C => B[:,S+1:end]
      )
  else
      output = Dict{Symbol,Any}(
      :p => parameters,
      :t => sol.t,
      :B => B)
  end

  return output

end

ext: jl


text:
#=
throw:
- Julia version: 1.0
- Author: singularitti
- Date: 2019-07-18
=#
export throw_which_in,
    throw_which_occursin,
    throw_which_match

function throw_which_in(xs, each)
    for x in xs
        if x in each
            return x
        else
            continue
        end  # if-else
    end  # for
    return nothing
end  # function throw_which_in

function throw_which_occursin(xs, str::AbstractString)
    for x in xs
        if occursin(x, str)
            return x
        else
            continue
        end  # if-else
    end  # for
    return nothing
end  # function throw_which_occursin

function throw_which_match(rs, line)
    for r in rs
        m = match(r, line)
        if m ≢ nothing
            return r
        else
            continue
        end  # if-else
    end  # for
    return nothing
end  # function throw_which_match


text:
using Test, SimplePolynomials, Polynomials, Mods, Multisets

a = SimplePolynomial(1, 1, 1)
b = SimplePolynomial(1, 2, 3)
c = SimplePolynomial(0, 1, -3, 4, 1)
x = getx()
p = (1 + x + x^3 - x^2)
q = (1 + 3x) * (1 - x) + 5

@testset "Polynomial Constructors" begin
    @test p == SimplePolynomial(1, 1, -1, 1)
    @test p == SimplePolynomial([1, 1, -1, 1])
    @test c == SimplePolynomial(c)
    @test c == SimplePolynomial(c.data)
    @test x + Mod{13}(2) == SimplePolynomial([Mod{13}(2), Mod{13}(1)])
end

@testset "Rational Function Constructors" begin
    f = p / q
    @test SimpleRationalFunction(p, q) == f
    @test f(1) == p(1) // q(1)
    f = inv(p)
    @test f(1) * p(1) == 1
end

@testset "String" begin
    p = getx()+1
    @test string(p) == "1 + x"
end

@testset "Polynomial Arithmetic" begin
    @test a + a == 2a
    @test (a + a + a) / 3 == a
    @test a - a == 0
    @test a * a == SimplePolynomial(1, 2, 3, 2, 1)
    @test a * a * a == a^3
    @test eltype(b) == BigInt
end

@testset "Rational Function Arithmetic" begin
    f = p / q
    g = a / b
    @test f * g == (p * a) / (b * q)
    @test f + g == (p * b + a * q) // (b * q)
    @test f - g == (-g) + f
    @test f / g == (p * b) / (a * q)
    @test inv(f) == q / p
end

@testset "Complex numbers" begin
    z = a + im
    @test z(-5) == a(-5) + im
    f = (im + x) / (2 * im * x)
    @test 2f(1) == 1 - im
    z = (x-im)*(x+im)
    @test z == x^2 + 1
end

@testset "Mod numbers" begin
    z = SimplePolynomial(Mod{13}.([1,2,3]))
    @test z(1) == Mod{13}(6)
    f = inv(z)
    @test f(1) == inv(Mod{13}(6))
end

@testset "Polynomial GCD" begin
    p = (x-1)^2 * (x-2)
    q = (x-1) * (x-3)
    @test gcd(p,q) == x-1

    p = (x-1)^3 * (x-3)
    @test gcd(p,p') == (x-1)^2

    c = Mod{13}.(coeffs(p))
    p = SimplePolynomial(c)
    @test gcd(p,p') == (x-1)^2
end

# @testset "Size change" begin
#     # B = big(b)
#     @test B == b
#     @test eltype(B) == BigInt
#     @test eltype(small(B)) == Int

#     r = a/b
#     @test eltype(numerator(r)) == Rational{BigInt}
#     # r = big(r)
#     # @test eltype(denominator(r)) == Rational{BigInt}
# end

@testset "PolyCalculus" begin
    p = 3x // 2 - 5x^2 + (2 - im) * x^3
    @test derivative(integral(p)) == p
end

@testset "Mod Polys" begin
    p = SimplePolynomial(Mod{13}.(0:4))
    @test p(0) == Mod{13}(0)
    @test 2p == p + p
    @test 3 * (p / 3) + 1 == Mod{13}(1) + p
end

@testset "Rational Basics" begin
    f = a / b
    g = b / a
    @test f * g == 1
    @test 2f == (a + a) / b
    @test (f - g)^2 == f^2 - 2f * g + g^2
    @test (f / g)^2 == (f^2) // (g^2)
end

@testset "Rational derivative" begin
    f = inv(1 + x)
    g = derivative(f)
    @test g == -1 // (1 + x)^2
end

@testset "Rational Mod" begin
    x = getx(Mod{13})
    p = Mod{13}(1) + x
    q = Mod{13}(2) + 3x
    f = p / q
    @test f(0) == Mod{13}(1 // 2)
end

@testset "Functionize" begin
    x = getx()
    p = 3 + 5x + 8x^2
    P = make_function(p)
    @test p(10) == P(10)

    f = p / (1-x)
    F = make_function(f)
    @test F(10) == f(10)
    A = [ 2 3 ; 0 -1]
    p = -2 - x + x^2
    @test all(p(A).==0)
end

@testset "Rational Roots" begin
    x = getx()
    p = (x^2+1)*(2x-1)^2*(3x+5)
    R = rational_roots(p)
    @test R == Multiset([1//2, -5//3, 1//2])
end

@testset "Numerical Roots" begin
    x = getx()
    p = (x^2+1)*(2x-1)^2*(3x+5)*(x^2+2)
    R = roots(p)
    Z = p.(R)
    @test sum(abs(z) for z in Z) < 10^-10

    R = newton_roots(p)
    Z = p.(R)
    @test sum(abs(z) for z in Z) < 10^-50

end

nothing


text:
function _convex_combination(Q::AbstractMatrix, κ::Real)::AbstractMatrix
    """
    src: https://arxiv.org/pdf/1310.3396.pdf
    ---------
    Parameters:
    Q: Covariance matrix

    """
    Q_prim = κ * Q + (1-κ) * I(size(Q)[1])
    return Q_prim
end


text:
using ConvexTests, SDPA

using MathOptInterface, Dualization
const MOI = MathOptInterface

function opt(mode)
    () -> begin
        opt = SDPA.Optimizer(Mode=mode)
        MOI.set(opt, MOI.Silent(), true)
        opt
    end
end

exclusions = [r"mip", r"exp", r"sdp_Complex_Semidefinite_constraint"]

@info "Starting SDPA tests with `Mode=PARAMETER_DEFAULT`"
do_tests("SDPA", opt(PARAMETER_DEFAULT); exclude = exclusions, variant="`PARAMETER_DEFAULT`", description="Tests with SDPA via SDPA.jl.", last = false, first = true)

@info "Starting SDPA tests with `Mode=PARAMETER_UNSTABLE_BUT_FAST`"
do_tests("SDPA", opt(PARAMETER_UNSTABLE_BUT_FAST); exclude = exclusions, variant="`PARAMETER_UNSTABLE_BUT_FAST`", first = false, last = false)

@info "Starting SDPA tests with `Mode=PARAMETER_STABLE_BUT_SLOW`"
do_tests("SDPA", opt(PARAMETER_STABLE_BUT_SLOW); exclude = exclusions, variant="`PARAMETER_STABLE_BUT_SLOW`", first = false, last = false)

function dual_opt(mode)
    () -> DualOptimizer(SDPA.Optimizer(Mode=mode))
end

@info "Starting SDPA tests with `Mode=PARAMETER_DEFAULT` (dualized)"

do_tests("SDPA", dual_opt(PARAMETER_DEFAULT); exclude = exclusions, variant="`PARAMETER_DEFAULT` (dualized)", description="Tests with SDPA via SDPA.jl using Dualization.jl.", first = false, last = false)

@info "Starting SDPA tests with `Mode=PARAMETER_UNSTABLE_BUT_FAST` (dualized)"

do_tests("SDPA", dual_opt(PARAMETER_UNSTABLE_BUT_FAST); exclude = exclusions, variant="`PARAMETER_UNSTABLE_BUT_FAST` (dualized)", first = false, last = false)

@info "Starting SDPA tests with `Mode=PARAMETER_STABLE_BUT_SLOW` (dualized)"
do_tests("SDPA", dual_opt(PARAMETER_STABLE_BUT_SLOW); exclude = exclusions, variant="`PARAMETER_STABLE_BUT_SLOW` (dualized)", first = false, last = true)


text:
# NOTE missing inf or nan NOT allowed
# NOTE one at a time, syntax for multiple would be confusing
@with_kw mutable struct Scatter3D{T} <: Drawing3D
    data       ::T  # data container
    linestyle  ::LineStyle   = LineStyle()   # line style (color, width, ...)
    markerstyle::MarkerStyle = MarkerStyle() # marker style (color, size, ...)
    label      ::String      = ""            # plot labels (to go in the legend)
end

mutable struct Surface <: Drawing3D

end


text:
"""
Fichier principal
"""

const MAX_CAPAVEHICLE = 20
const SPEED_VEHICLE = 20

const SafeFloat64 = Union{Missing, Float64}

@enum Etat S R E Non

include("structs.jl")
include("dataManager.jl")
include("branchAndBound.jl")
include("tools.jl")

function formul(lat1::Float64, long1::Float64, lat2::Float64, long2::Float64)
    r = 6371008
    distAng = acos( sin(lat1 * pi/180) * sin(lat2 * pi/180) + cos(lat1 * pi/180) * cos(lat2 * pi/180) * cos((long1 - long2)*pi/180) )

    return r * distAng
end

function calculDistLatLong(stations::Vector{Station})

    nbStations = length(stations)
    distStations = Array{Float64, 2}(undef, nbStations, nbStations)

	for iter = 1:nbStations
        distStations[iter, iter] = 0.
    end

    for iter = 1:(nbStations-1)
		dist = formul(stations[iter].latitude, stations[iter].longitude, stations[iter+1].latitude, stations[iter+1].longitude)
        distStations[iter, iter+1] = dist
		distStations[iter+1, iter] = dist
    end

	for iterRow = 1:nbStations
		for iterCol = (iterRow+2):nbStations
			distStations[iterRow, iterCol] = 0.
			distStations[iterCol, iterRow] = 0.
			for iter = 1:(iterCol-iterRow)
				distStations[iterRow, iterCol] += distStations[iter, iter+1]
				distStations[iterCol, iterRow] += distStations[iter, iter+1]
			end
		end
	end
    return distStations
end

function calcLimit(elt::Elt, w::Float64, epsilon::Float64, tripDate::Array{Float64, 2}, indReq::Int, distStation::Array{Float64, 2}, dictStation::Dict{String, Station}, requests::Vector{Request})
	if elt.state == S
		return w - tripDate[indReq, 1] + tripDate[elt.idReq, 1]
	elseif elt.state == E && tripDate[elt.idReq, 2] != -1.
		return (1+epsilon)*distStation[dictStation[requests[elt.idReq].departureStation].id, dictStation[requests[elt.idReq].arrivalStation].id] + requests[indReq].t - tripDate[elt.idReq, 2]
	elseif elt.state == E
		return -1
	else
		error("La tu essayes de calculer la limite d'un element R")
	end
end

function calcCMin(matTime::Array{Float64, 2}, iter::Int)
	mini = Inf
	taille = size(matTime)[1]
	for ind = 1:taille
		if ind != iter+1
			if mini > matTime[iter+1, ind]
				mini = matTime[iter+1, ind]
			end
		end
	end
	return mini
end

function simulate(nameSimul::String, nameFirstStation::String, w::Float64, epsilon::Float64; verbose)
	# On récupère toutes les données relatives au problème
    allStation = getStations("stations.dat")
    distStation = calculDistLatLong(allStation) / SPEED_VEHICLE

    # Landreau, Cousteau, Perray, ZI1Garde, Lycee, Carquefou-Gare et Carquefou-Centre
    indActualStation = [1, 2, 3 ,4 ,7, 9, 12]

    actualStation = allStation[indActualStation]									# Les stations utilisées
    dictStation = Dict(station.name => station for station in allStation) 	# Dictionnaire pour avoir les stations depuis leur nom
    actualDistance = distStation[indActualStation, indActualStation]				# La distance entre ces stations
    requests, nbShuttles = parseRequests(nameSimul)									# On récupère la liste des requètes
    nbRequests = length(requests)													# Calcul du nombre de requete

    tripDate = Array{Float64, 2}(undef, nbRequests, 3)								# Matrice qui retiendra le temps de reception d'une demande, de la
    for iter = 1:nbRequests*3    													# récupération et de l'amenage du client
    	tripDate[iter] = -1.
    end
    etatToInt = Dict{Etat, Int}(R => 1, S => 2, E => 3)

    actualTime = requests[1].t
    tripDate[1, 1] = actualTime														# Temps réel de la navette
    iterTrip = 2																	# Indice du premier élement non effectué

    #= Pour une navette pour le moment =#
    iterFinTrip = 3
    trip = Vector{Elt}(undef, 3*nbRequests) 										# Liste du trip pour une navette
    trip[1] = Elt(1, R, 1, false, false, 0., "r1", 0.)								# On fixe la première demande au début du trip de la première navette
    trip[2] = Elt(2, S, 1, true, false, 0., "s1", 0.)
    trip[3] = Elt(3, E, 1, false, true, 0., "e1", 0.)
    #First request assigned to the first shuttle

    tripDate[1] = requests[1].t
    indTimeStation = dictStation[nameFirstStation].id
    stopStation = (-1, -1)
    trajTime = 0

    indReq = 2
    for indReq in 2:nbRequests
    	if trip[iterTrip].state == S
    		indNextStation = dictStation[requests[trip[iterTrip].idReq].departureStation].id
    	elseif trip[iterTrip].state == E
    		indNextStation = dictStation[requests[trip[iterTrip].idReq].arrivalStation].id
    	else
    		error("Tu essayes d'aller sur un sommet R ?!")
    	end

    	if length(indTimeStation) == 1
			if  indTimeStation > indNextStation
				trajTime = distStation[indTimeStation, (indTimeStation-1)]
				stopStation = (indTimeStation, indTimeStation-1)
			elseif indTimeStation < indNextStation
				trajTime = distStation[indTimeStation, (indTimeStation+1)]
				stopStation = (indTimeStation, indTimeStation+1)
			else
				trajTime = 0
			end
    	end

    	while actualTime + trajTime < requests[indReq].t && iterTrip <= iterFinTrip

    		if length(indTimeStation) > 1
    			actualTime += trajTime
    			indTimeStation = indTimeStation[1]

    			if  indTimeStation > indNextStation
					trajTime = distStation[indTimeStation, (indTimeStation-1)]
					stopStation = (indTimeStation, indTimeStation-1)
				elseif indTimeStation < indNextStation
					trajTime = distStation[indTimeStation, (indTimeStation+1)]
					stopStation = (indTimeStation, indTimeStation+1)
				else
					trajTime = 0
				end

    		end

			while actualTime + trajTime < requests[indReq].t && iterTrip <= iterFinTrip && indTimeStation != indNextStation

				actualTime += trajTime

				if indTimeStation > indNextStation
					indTimeStation -= 1
				else
					indTimeStation += 1
				end

				if indTimeStation > indNextStation
					trajTime = distStation[indTimeStation, (indTimeStation-1)]
					stopStation = (indTimeStation, indTimeStation-1)
				elseif indTimeStation < indNextStation
					trajTime = distStation[indTimeStation, (indTimeStation+1)]
					stopStation = (indTimeStation, indTimeStation+1)
				else
					trajTime = 0
				end

			end

			if indTimeStation == indNextStation
				tripDate[trip[iterTrip].idReq, etatToInt[trip[iterTrip].state]] = actualTime

				iterTrip += 1

				if iterTrip <= iterFinTrip
					if trip[iterTrip].state == S
						indNextStation = dictStation[requests[trip[iterTrip].idReq].departureStation].id
					elseif trip[iterTrip].state == E
						indNextStation = dictStation[requests[trip[iterTrip].idReq].arrivalStation].id
					else
						error("Tu essayes d'aller sur un sommet R ?!")
					end

					if  indTimeStation > indNextStation
						trajTime = distStation[indTimeStation, (indTimeStation-1)]
						stopStation = (indTimeStation, indTimeStation-1)
					elseif indTimeStation < indNextStation
						trajTime = distStation[indTimeStation, (indTimeStation+1)]
						stopStation = (indTimeStation, indTimeStation+1)
					else
						trajTime = 0
					end
				end
			end
		end

		if iterTrip > iterFinTrip
			actualTime = requests[indReq].t
			trajTime = 0
			stopStation = (indTimeStation, indTimeStation)
			iterTrip = iterFinTrip
		end
		tripDate[indReq, 1] = requests[indReq].t

		iterFinTripBis = iterFinTrip
		tripBis = trip[1:end]

		iterFinTrip += 1
		trip[iterFinTrip] = Elt(iterFinTrip, S, indReq, true, false, 0., "s$indReq", 0.)
		iterFinTrip += 1
		trip[iterFinTrip] = Elt(iterFinTrip, E, indReq, false, true, 0., "e$indReq", 0.)

		listEltTri = sort(trip[iterTrip:iterFinTrip], lt = myIsLess)

		listIndEltPasFini = [dictStation[nameStat].id for nameStat in [elt.state == S ? requests[elt.idReq].departureStation : requests[elt.idReq].arrivalStation for elt in listEltTri]]

		matTime = Array{Float64, 2}(undef, iterFinTrip-iterTrip+2, iterFinTrip-iterTrip+2)
		firstDest = Vector{Int}(undef, iterFinTrip-iterTrip+1)
		timeNextStatR = Vector{Float64}(undef, iterFinTrip-iterTrip+1)
		listeVal1 = Vector{Float64}(undef, iterFinTrip-iterTrip+1)
		listeVal2 = Vector{Float64}(undef, iterFinTrip-iterTrip+1)
		listeChoix = Vector{String}(undef, iterFinTrip-iterTrip+1)
		ecartTemps1 = 0

		matTime[2:end, 2:end] = distStation[listIndEltPasFini, listIndEltPasFini]

		if length(indTimeStation) == 1
			ecartTemps1 = (requests[indReq].t - actualTime)
			for iter =  2:(iterFinTrip-iterTrip+2)
				iterStat = listIndEltPasFini[iter-1]
				val1 = (requests[indReq].t - actualTime) + distStation[stopStation[1], iterStat]
				val2 = trajTime + distStation[stopStation[2], iterStat] - (requests[indReq].t - actualTime)

				listeVal1[iter-1] = val1
				listeVal2[iter-1] = val2
				if val1 < val2
					matTime[iter, 1] = val1
					matTime[1, iter] = val1
					firstDest[iter-1] = stopStation[1]
					timeNextStatR[iter-1] = (requests[indReq].t - actualTime)
					listeChoix[iter-1] = "val1"
				else
					matTime[iter, 1] = val2
					matTime[1, iter] = val2
					firstDest[iter-1] = stopStation[2]
					timeNextStatR[iter-1] = trajTime - (requests[indReq].t - actualTime)
					listeChoix[iter-1] = "val2"
				end
			end
		else
			listeVal1Prec, listeVal2Prec, listeChoixPrec, firstDestPrec = indTimeStation[2:end]
			diffTemps = requests[indReq].t - actualTime
			for iter =  2:(length(listeChoixPrec)+1)
				iterStat = listIndEltPasFini[iter-1]

				val1, val2 = 0, 0

				if listeChoixPrec[iter-1] == "val1"
					val1 = listeVal1Prec[iter-1] - diffTemps
					val2 = listeVal2Prec[iter-1] + diffTemps
				elseif listeChoixPrec[iter-1] == "val2"
					val1 = listeVal1Prec[iter-1] + diffTemps
					val2 = listeVal2Prec[iter-1] - diffTemps
				else
					error("Nom non accapté dans la liste listeChoix !")
				end

				listeVal1[iter-1] = val1
				listeVal2[iter-1] = val2
				if val1 < val2
					matTime[iter, 1] = val1
					matTime[1, iter] = val1
					timeNextStatR[iter-1] = val1 - distStation[firstDestPrec[iter-1], iterStat]
					firstDest[iter-1] = stopStation[1]
					listeChoix[iter-1] = "val1"
				else
					matTime[iter, 1] = val2
					matTime[1, iter] = val2
					timeNextStatR[iter-1] = val2 - distStation[firstDestPrec[iter-1], iterStat]
					firstDest[iter-1] = stopStation[2]
					listeChoix[iter-1] = "val2"
				end
			end
			for iter = (length(listeChoixPrec)+2):(iterFinTrip-iterTrip+2)
				iterStat = listIndEltPasFini[iter-1]
				val1 = trajTime - (requests[indReq].t - actualTime) + distStation[stopStation[1], iterStat]
				val2 = distStation[stopStation[1], stopStation[2]] - trajTime + (requests[indReq].t - actualTime) + distStation[stopStation[2], iterStat]

				listeVal1[iter-1] = val1
				listeVal2[iter-1] = val2
				if val1 < val2
					matTime[iter, 1] = val1
					matTime[1, iter] = val1
					firstDest[iter-1] = stopStation[1]
					timeNextStatR[iter-1] = (requests[indReq].t - actualTime)
					listeChoix[iter-1] = "val1"
				else
					matTime[iter, 1] = val2
					matTime[1, iter] = val2
					firstDest[iter-1] = stopStation[2]
					timeNextStatR[iter-1] = trajTime - (requests[indReq].t - actualTime)
					listeChoix[iter-1] = "val2"
				end
			end
		end

		matTime[1, 1] = 0.

		listElt = Vector{Elt}(undef, iterFinTrip-iterTrip+1)
		for iter = 1:(iterFinTrip-iterTrip+1)
			elt = listEltTri[iter]
			listElt[iter] = Elt(iter+1, elt.state, elt.idReq, elt.isSource, (elt.state == E && tripDate[elt.idReq, 2] == -1), calcLimit(elt, w, epsilon, tripDate, indReq, distStation, dictStation, requests), elt.name, calcCMin(matTime, iter))
		end

		listNextElt, lastElt = branchAndBound(listElt, matTime, epsilon, verbose=verbose)

		if (listNextElt, lastElt) != (false, false)

			trajTime = timeNextStatR[listNextElt[2].id-1]
			indTimeStation = (firstDest[listNextElt[2].id-1], listeVal1, listeVal2, listeChoix, firstDest)

			trip[iterTrip] = Elt(iterFinTrip, R, indReq, false, false, 0., "r$indReq", 0.)
			for iter = 1:(iterFinTrip-iterTrip)
				elt = listNextElt[iter+1]
				trip[iter + iterTrip] = Elt(iter+ iterTrip - 1, elt.state, elt.idReq, elt.isSource, false, 0., elt.name, 0.)
			end
			iterFinTrip += 1
			trip[iterFinTrip] = Elt(iterFinTrip+1, lastElt.state, lastElt.idReq, lastElt.isSource, false, 0., lastElt.name, 0.)

			iterTrip += 1
		else
			println("Impossible d'accepter le client $indReq")
			iterFinTrip = iterFinTripBis
			trip = tripBis[1:end]
		end
		actualTime = requests[indReq].t
		# Recup Next Station
		# Se mettre au bon temps grace à matTime
		# Et on recommence en les ajoutant au Trip
	end

	while iterTrip <= iterFinTrip
		if length(indTimeStation) != 1
			actualTime += trajTime

			indTimeStation = indTimeStation[1]
		end
		elt = trip[iterTrip]
		nextStat = dictStation[elt.state == S ? requests[elt.idReq].departureStation : requests[elt.idReq].arrivalStation].id
		actualTime += distStation[indTimeStation, nextStat]

		tripDate[trip[iterTrip].idReq, etatToInt[trip[iterTrip].state]] = actualTime
		iterTrip += 1
		indTimeStation = nextStat
	end

	return trip, tripDate
end

# fonction principale
function main(nameSimul::String = "simulation2.dat", nameFirstStation::String = "Gare", w::Float64 = 15*60., epsilon::Float64 = 0.5; verbose = false)
	trip, tripDate = simulate(nameSimul,nameFirstStation,w,epsilon,verbose=verbose)

	print("Trajet de la navette : ")
	for i in 1:length(trip)
		if (trip[i].state == R)
			println("")
		end
		print(" -> ")
		print(trip[i].name)
		print(" (",round(tripDate[i],digits=0),")")
	end
end


text:
@testset "ϵ-support" begin
    a, b = ϵsupport(inv, 1)
    @test a ≈ -1 atol = 1e-2
    @test b ≈ +1 atol = 1e-2
end

@testset "Quad Rule" begin
    q = adaptive_G_quad(identity)
    @test q[1] ≈ 0 atol = 1e-14
    q = adaptive_G_quad(identity; a = 0)
    @test q[1] ≈ 1/2 atol = 1e-14
end

max_stars_repo_name: EmileParolin/STrAW

max_issues_repo_name: EmileParolin/STrAW

max_forks_repo_name: EmileParolin/STrAW


text:
const _blacklisted_replacements = Set([ :*, :/, :+, :-, :\ ])
const _irrationals = Set([:π, :ℯ, :γ, :φ, :catalan])

function should_replace(x)
    MacroTools.isexpr(x, :call) || return false
    (x.args[1] ∉ _blacklisted_replacements) || return false

    # HACK: CUDa can't handle `Irrational` so we gotta do something about that.
    # If uniary function and argument is `Irrational`, don't replace
    if length(x.args) == 2 && x.args[2] ∈ _irrationals
        return false
    end

    return true
end

function replace_device_all(ex)
    MacroTools.postwalk(ex) do x
        x = if should_replace(x)
            :(CUDA.cufunc($(x.args[1]))($(x.args[2:end]...)))
        else
            x
        end

        x = CUDA._cuint(x)
        x = CUDA._cupowliteral(x)
        x
    end
end

function tocuname(f)
    # If we have expression of form `Mod1.Mod2.f` then we define `cufunc(Mod1.Mod2.f) = cuf`
    # not `cuMod1.Mod2.f`.
    if MacroTools.isexpr(f, :.)
        Symbol(:cu, f.args[end].value)
    else
        Symbol(:cu, f)
    end
end

function make_cufunc_function(def)
    f = def[:name]
    def[:name] = tocuname(f) # Expr(:., Symbol(mod), QuoteNode(tocuname(f)))
    def[:body] = replace_device_all(def[:body])
    return :($(esc(MacroTools.combinedef(def))))
end

function make_cufunc_register(fname)
    cuname = tocuname(fname)
    return :(@inline CUDA.cufunc(::typeof($(esc(fname)))) = $(esc(cuname)))
end

function make_cufunc_diffrule(mod, f, nargs)
    # TODO: should default be `Base` or `@__MODULE__`?
    m, fname = MacroTools.isexpr(f, :.) ? (Symbol(f.args[end - 1]), f.args[end].value) : (:Base, f)

    # Return early if no diff-rule present
    DiffRules.hasdiffrule(m, fname, nargs) || return nothing

    # TODO: support binary also; currently the issue is that the rule isn't generated
    # properly due to the use of `map` I believe.
    nargs == 1 || throw(ArgumentError("`make_cufunc_diffrule` currently only supports nargs=1"))

    # Get the name
    cuname = tocuname(Expr(:., m, QuoteNode(fname)))

    # Create placeholders for arguments used to define new diffrule
    args = ntuple(i -> gensym(Symbol(:x, i)), nargs)
    # Get the diffrule with `args` as the representations of the arguments
    drule = DiffRules.diffrule(m, fname, args...)

    drule_interpolated = if nargs > 1
        map(drule) do dc
            MacroTools.postwalk(Meta.quot(dc)) do e
                e in args ? Expr(:$, e) : e
            end
        end
    else
        MacroTools.postwalk(Meta.quot(drule)) do e
            e in args ? Expr(:$, e) : e
        end
    end

    # Replace all methods with `cufunc`
    drule_cu = if nargs > 1
        map(replace_device_all, drule_interpolated)
    else
        replace_device_all(drule_interpolated)
    end

    # On the RHS of a DiffRule, we need to make the arguments interpolated, e.g. x -> $x
    rhs = drule_cu

    # Determine which method to use for evaluating the diff-rule created
    def_func = if nargs == 1
        :(ForwardDiff.unary_dual_definition)
    # elseif nargs == 2
    #     :(ForwardDiff.binary_dual_definition)
    else
        throw(ArgumentError("this one you gotta do yourself buddy"))
    end

    return quote
        $(DiffRules).@define_diffrule $(Symbol(mod)).$(cuname)($(esc.(args)...)) = $(esc(rhs))
        eval($(def_func)($(QuoteNode(Symbol(mod))), $(QuoteNode(cuname))))
    end
end

macro cufunc(ex)
    def = MacroTools.splitdef(ex)
    f = def[:name]
    nargs = length(def[:args])
    cufunc_function = make_cufunc_function(def)
    cufunc_register = make_cufunc_register(f)
    # `__module__` corresponds the module of the calling scope, rather than
    # `@__MODULE__` which corresponds to the scope where the code is read, i.e. here.
    cufunc_diffrule = make_cufunc_diffrule(__module__, f, nargs)
    quote
        $cufunc_function
        $cufunc_register
        $cufunc_diffrule
    end
end

"""
    @cufunc_function function f(args...) ... end

Converts the function into a CUDA-broadcastable function by replacing all
functions calls in the body with `cufunc` versions.
"""
macro cufunc_function(ex)
    def = MacroTools.splitdef(ex)
    return make_cufunc_function(def)
end

"""
    @cufunc_register Module.f
    @cufunc_register Module.f(args...)

Registers `Module.f` to its corresponding `cuf` via `CUDA.cufunc`.

Useful in cases where you have multiple different definitions of `Module.f`.
"""
macro cufunc_register(ex)
    MacroTools.@capture(ex, (f_(args_) | f_))
    return make_cufunc_register(f)
end

"""
    @cufunc_diffrule Module.f(args...)
    @cufunc_diffrule TargetModule Module.f(args...)

If it exists, diffrule from `DiffRules` is extracted and `cufunc`ed.

`TargetModule` defines the module for the diffrule to be registered,
as is required by DiffRules.jl.
"""
macro cufunc_diffrule(ex)
    MacroTools.@capture(ex, f_(args__)) || return nothing
    # `__module__` corresponds the module of the calling scope, rather than
    # `@__MODULE__` which corresponds to the scope where the code is read, i.e. here.
    return make_cufunc_diffrule(__module__, f, length(args))
end

macro cufunc_diffrule(mod, ex)
    MacroTools.@capture(ex, f_(args__)) || return nothing
    return make_cufunc_diffrule(mod, f, length(args))
end

"""
    @cufuncf function f(args...) ... end

Same as [CUDAExtensions.@cufunc](@ref) but also includes the function `f`.

This is useful if you're defining a method for the first time and also
want it to be CUDA-broadcastable, rather than making an existing method
CUDA-broadcastable.
"""
macro cufuncf(ex)
    newex = quote
        $(Base).@__doc__ $ex
        CUDAExtensions.@cufunc($ex)
    end
    return esc(newex)
end


text:
function GenDataSize(m,n,r)

## This script creates an sxs data matrix consisting of a discrezation of 4
## 2-D Gaussian bumps. The resulting matrix is rank 4 (= number of bumps)

sc = 1

# Create regular grid
st = 1/n
X = 0:st:1-st
X = sc*X

sy = 1/m
Y = 0:sy:1-sy
Y = sc*Y

c = 50/sc^2

Reg = zeros(m,n)

# Bump 1
for j = 1:r

a = round(rand(1),2)
b = round(rand(1),2)

    for k = 1:m
        for l = 1:n
            Reg[k,l] += dot(exp((-c*([sc*a - X[l]].^2))[1]),exp((-c*([sc*b - Y[k]].^2))[1]))
        end
    end
end

return Reg
end
export GenDataSize


text:
export PARSDMM_multi_level

"""
Computes the projection onto an intersection of sets using a coarse-to-fine grid
multilevel strategy. This wraps around PARSDMM.jl and includes interpolation to finer grids
as well as coarsening of the input model
"""
function PARSDMM_multi_level(
                            m                 :: Vector{TF},
                            TD_OP_levels,
                            AtA_levels,
                            P_sub_levels      ::Vector{Vector{Any}},
                            set_Prop_levels   ::Vector{Any},
                            comp_grid_levels  ::Vector{Any},
                            options,
                            x_ini=zeros(TF,prod(comp_grid_levels[end].n)) ::Vector{TF},
                            l_ini=[],
                            y_ini=[]
                            ) where {TF<:Real}

n_levels = length(TD_OP_levels)

m_levels         = Vector{Vector{TF}}(undef,n_levels) #allocate space for the model we will project at each level
m_levels[1]      = m

#interpolations to finer grids use nearest neighbour interpolation
#from the "Interpolations" julia package. Change (BSpline(Constant())) for other
#types of interpolation
rho_orig = deepcopy(options.rho_ini)

#detect 2D or 3D
if length(comp_grid_levels[1].n)==3 && comp_grid_levels[1].n[3]>1
  dim3 = true #indicate we are working in 3D
else
  dim3 = false
end

#coarsen original model to new grids
for i = 2:n_levels
  itp_m       = interpolate(reshape(m,comp_grid_levels[1].n),BSpline(Constant()))
  if dim3
    m_level   = itp_m(range(1,stop=comp_grid_levels[1].n[1],length=comp_grid_levels[i].n[1]), range(1,stop=comp_grid_levels[1].n[2],length=comp_grid_levels[i].n[2]), range(1,stop=comp_grid_levels[1].n[3],length=comp_grid_levels[i].n[3]))
  else
    m_level   = itp_m(range(1,stop=comp_grid_levels[1].n[1],length=comp_grid_levels[i].n[1]), range(1,stop=comp_grid_levels[1].n[2],length=comp_grid_levels[i].n[2]))
  end
  m_levels[i] = vec(m_level)
end

#start at coarsest grid (zero ini guess for x, y and l)
i=n_levels
#(x,log_PARSDMM,l,y) = compute_projection_intersection_PARSDMM(m_levels[i],AtA_levels[i],TD_OP_levels[i],set_Prop_levels[i],P_sub_levels[i],comp_grid_levels[i],options,FL)
options.zero_ini_guess = true

#solve coarsest level more accurately
(x,log_PARSDMM,l,y) = PARSDMM(m_levels[i],AtA_levels[i],TD_OP_levels[i],set_Prop_levels[i],P_sub_levels[i],comp_grid_levels[i],options,x_ini,l_ini,y_ini)
options.rho_ini = log_PARSDMM.rho[end,:]
for i=(n_levels-1):-1:1

  # resample x, l & y to a finer grid
  itp_x_level = interpolate(reshape(x,comp_grid_levels[i+1].n), BSpline(Constant()))
  if dim3
    x = itp_x_level(range(1,stop=comp_grid_levels[i+1].n[1],length=comp_grid_levels[i].n[1]), range(1,stop=comp_grid_levels[i+1].n[2],length=comp_grid_levels[i].n[2]), range(1,stop=comp_grid_levels[i+1].n[3],length=comp_grid_levels[i].n[3]))
  else
    x = itp_x_level(range(1,stop=comp_grid_levels[i+1].n[1],length=comp_grid_levels[i].n[1]), range(1,stop=comp_grid_levels[i+1].n[2],length=comp_grid_levels[i].n[2]))
  end
  x = vec(x)

  #interpolate y and l
  if options.parallel==true #for now, gather -> interpolate -> distribute, this is inefficient and need to be updated
    l = convert(Vector{Vector{TF}},l)
    y = convert(Vector{Vector{TF}},y)
  end
  (l,y)=interpolate_y_l(l,y,set_Prop_levels,comp_grid_levels,dim3,i)
  if options.parallel==true #for now, gather -> interpolate -> distribute, this is inefficient and need to be updated
    l = distribute(l)
    y = distribute(y)
  end

  # solve on a finer grid
  options.zero_ini_guess = false
  (x,log_PARSDMM,l,y)    = PARSDMM(m_levels[i],AtA_levels[i],TD_OP_levels[i],set_Prop_levels[i],P_sub_levels[i],comp_grid_levels[i],options,x,l,y)
  options.rho_ini        = log_PARSDMM.rho[end,:]

end #end for loop over levels

options.rho_ini = rho_orig
return x, log_PARSDMM, l, y
end #end multi_level_PARSDMM

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
type SAGAinitsd <: StepData
	d
	functions
	y
	I
	getStep::Function
	stepString::AbstractString
	function SAGAinitsd(numVars::Int64,dt::DataType, stepString::AbstractString, numDp::Int64)
		new(zeros(numVars),Array(Function, numDp), Array(dt, numDp),0,  SAGAinitComputation, stepString)
	end
end

function SAGAinit(numVars::Int64, dt::DataType, numDP::Int64)
	SAGAinitsd(numVars,dt, "SAGAinit",numDP)
end

function SAGAinitComputation(x, k, gnum, sd::SAGAinitsd, problem::Problem; outputLevel =0)

	outputLevel >0 && println("In SAGAinitComputation")
	i = rand(1:problem.numTP)

	step = sd.d

	# step 1
	try
		outputLevel >0 && println("Trying to do full SAGA update")
		oldG = sd.y[i]
		step=sd.d/sd.I - oldG
		sd.d -= oldG

		(f,sampleG) = (sd.functions[i])(x)
		sd.y[i] = sampleG
		gnum += 1
		sd.d += sampleG
		step = step + sampleG

	catch
		outputLevel >0 && println("Falling back on SG")
		sd.functions[i] = problem.getSampleFunctionAt(i)
		outputLevel >0 && println("Added a function")
		sd.I+=1

		(f,sampleG) = (sd.functions[i])(x)
		outputLevel >0 && println("Computed new gradient")
		sd.y[i] = sampleG
		outputLevel >0 && println("Stored in y")
		gnum += 1
		sd.d += sampleG

		step = sampleG
	end

	outputLevel >0 && println("Finishing SAGAinitComputation")

	(step, gnum)
end


text:
# to be run on all nodes

function sayhello()
   println("hi I am worker number $(myid()), I live on $(gethostname())")
end
flush(stdout)

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
# ---
# title: 1375. Bulb Switcher III
# id: problem1375
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/bulb-switcher-iii/description/>
# hidden: true
# ---
#
# There is a room with `n` bulbs, numbered from `1` to `n`, arranged in a row
# from left to right. Initially, all the bulbs are turned off.
#
# At moment _k_ (for _k_ from `0` to `n - 1`), we turn on the `light[k]` bulb. A
# bulb **change  color to blue** only if it is on and all the previous bulbs (to
# the left) are turned on too.
#
# Return the number of moments in which **all  turned on** bulbs  **are blue.**
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/02/29/sample_2_1725.png)
#
#
#
#     Input: light = [2,1,3,5,4]
#     Output: 3
#     Explanation: All bulbs turned on, are blue at the moment 1, 2 and 4.
#
#
# **Example 2:**
#
#
#
#     Input: light = [3,2,4,1,5]
#     Output: 2
#     Explanation: All bulbs turned on, are blue at the moment 3, and 4 (index-0).
#
#
# **Example 3:**
#
#
#
#     Input: light = [4,1,2,3]
#     Output: 1
#     Explanation: All bulbs turned on, are blue at the moment 3 (index-0).
#     Bulb 4th changes to blue at the moment 3.
#
#
# **Example 4:**
#
#
#
#     Input: light = [2,1,4,3,6,5]
#     Output: 3
#
#
# **Example 5:**
#
#
#
#     Input: light = [1,2,3,4,5,6]
#     Output: 6
#
#
#
#
# **Constraints:**
#
#   * `n == light.length`
#   * `1 <= n <= 5 * 10^4`
#   * `light` is a permutation of  `[1, 2, ..., n]`
#
#
## @lc code=start
using LeetCode

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


text:
using Gadfly

set_default_plot_size(5.6inch, 5inch)

parse_colorant = Gadfly.parse_colorant
cs = ["red", "green", "blue"]
shps = ["circle", "square", "triangle"]
szs = ["1mm", "2mm", "3mm"]

plot(x=rand(100), y=rand(100), color=rand(cs, 100),
    shape=rand(shps, 100), size=rand(szs, 100),
    Scale.color_discrete(c->parse_colorant.(cs), levels=cs),
    Scale.size_discrete2(s->[1mm, 2mm, 3mm], levels=szs),
    Scale.shape_discrete(levels=shps),
    Theme(key_swatch_color="gray", key_swatch_size=2mm,
      point_shapes=[Shape.circle,Shape.square,Shape.utriangle])
)


text:
using Neo4j, DataFrames
using Test

@testset "Module imports" begin
    @test (@isdefined Neo4j) == true
    @test typeof(Neo4j) == Module
end

# defaults for testing
global username = "neo4j"
global passwd = "neo5j"

global graph = nothing
global conn = nothing
@testset "Creating a connection to localhost" begin
    try
      global graph = getgraph()
    catch
      @info "[TEST] Anonymous connection failed! Creating a Neo4j connection to localhost:7474 with neo4j:$(passwd) credentials..."
      #Trying with security.
      global conn = Neo4j.Connection("localhost"; user=username, password=passwd);
      global graph = getgraph(conn);
    end
    global conn = graph.connection;
end

@testset "Checking version of connected graph = Neo4j $(ascii(graph.version))..." begin
    # Have to account for newer Neo4j! Using version text - ref: http://docs.julialang.org/en/release-0.4/manual/strings/
    @test VersionNumber(graph.version) > v"2.0.0" # convert(VersionNumber, graph.version) >  v"2.0.0"
    # Check that
    @test graph.node == "http://localhost:7474/db/data/node"
end

global barenode = nothing
global propnode = nothing
@testset "Nodes: CRUD, properties, and labels..." begin
    global barenode = Neo4j.createnode(graph)
    @test barenode.self == "http://localhost:7474/db/data/node/$(barenode.id)"

    global propnode = Neo4j.createnode(graph, Dict{AbstractString,Any}("a" => "A", "b" => 1))  #UTF8String
    @test propnode.data["a"] == "A"
    @test propnode.data["b"] == 1

    global gotnode = getnode(graph, propnode.id)
    @test gotnode.id == propnode.id
    @test gotnode.data["a"] == "A"
    @test gotnode.data["b"] == 1

    setnodeproperty(barenode, "a", "A")
    global barenode = getnode(barenode)
    @test barenode.data["a"] == "A"

    global props = getnodeproperties(propnode)
    @test props["a"] == "A"
    @test props["b"] == 1
    @test length(props) == 2

    updatenodeproperties(barenode, Dict{AbstractString,Any}("a" => 1, "b" => "A"))  #UTF8String
    global barenode = getnode(barenode)
    @test barenode.data["a"] == 1
    @test barenode.data["b"] == "A"

    deletenodeproperties(barenode)
    global barenode = getnode(barenode)
    @test length(barenode.data) == 0

    deletenodeproperty(propnode, "b")
    global propnode = getnode(propnode)
    @test length(propnode.data) == 1
    @test propnode.data["a"] == "A"

    addnodelabel(barenode, "A")
    global barenode = getnode(barenode)
    @test getnodelabels(barenode) == ["A"]

    addnodelabels(barenode, ["B", "C"])
    global barenode = getnode(barenode)
    global labels = getnodelabels(barenode)
    @test "A" in labels
    @test "B" in labels
    @test "C" in labels
    @test length(labels) == 3

    updatenodelabels(barenode, ["D", "E", "F"])
    global barenode = getnode(barenode)
    global labels = getnodelabels(barenode)
    @test "D" in labels
    @test "E" in labels
    @test "F" in labels
    @test length(labels) == 3

    deletenodelabel(barenode, "D")
    global barenode = getnode(barenode)
    global labels = getnodelabels(barenode)
    @test "E" in labels
    @test "F" in labels
    @test length(labels) == 2

    global nodes = getnodesforlabel(graph, "E")
    @test length(nodes) > 0
    @test barenode.id in [n.id for n = nodes]

    global labels = getlabels(graph)
end

@testset "Relationships: CRUD, neighbors" begin
    global rel1 = createrel(barenode, propnode, "test"; props=Dict{AbstractString,Any}("a" => "A", "b" => 1)); #UTF8String
    global rel1alt = getrel(graph, rel1.id);
    @test rel1.reltype == "TEST"
    @test rel1.data["a"] == "A"
    @test rel1.data["b"] == 1
    @test rel1.id == rel1alt.id

    global endnode = Neo4j.createnode(graph, Dict{AbstractString,Any}("a" => "A", "b" => 1))  # UTF8String
    global rel2 = createrel(propnode, endnode, "test"; props=Dict{AbstractString,Any}("a" => "A", "b" => 1));  # UTF8String
    @test length(Neo4j.getrels(endnode)) == 1
    @test length(Neo4j.getrels(propnode)) == 2
    @test length(Neo4j.getrels(barenode)) == 1
    @test length(Neo4j.getrels(endnode, incoming=true, outgoing=false)) == 1
    @test length(Neo4j.getrels(endnode, incoming=false, outgoing=true)) == 0
    @test length(Neo4j.getrels(propnode, incoming=true, outgoing=false)) == 1
    @test length(Neo4j.getrels(propnode, incoming=false, outgoing=true)) == 1

    global neighbors = Neo4j.getneighbors(propnode)
    @test length(neighbors) == 2
    global neighbors = Neo4j.getneighbors(propnode, incoming=true, outgoing=false)
    @test length(neighbors) == 1
    @test neighbors[1].metadata["id"] == barenode.metadata["id"]
    global neighbors = Neo4j.getneighbors(propnode, incoming=false, outgoing=true)
    @test length(neighbors) == 1
    @test neighbors[1].metadata["id"] == endnode.metadata["id"]

    global rel1prop = getrelproperties(rel1);
    @test rel1prop["a"] == "A"
    @test rel1prop["b"] == 1
    @test length(rel1prop) == 2
    @test getrelproperty(rel1, "a") == "A"
    @test getrelproperty(rel1, "b") == 1

    deleterel(rel1)
    deleterel(rel2)
    @test_throws ErrorException getrel(graph, rel1.id)
    @test_throws ErrorException getrel(graph, rel2.id)
end

@testset "Nodes: Deleting nodes (cleaning up)" begin
    deletenode(graph, barenode.id)
    deletenode(graph, propnode.id)
    @test_throws ErrorException getnode(graph, barenode.id)
    @test_throws ErrorException getnode(graph, propnode.id)
end

# --- New transaction code from Glesica source ---
function createnode(txn, name, age; submit=false)
  q = "CREATE (n:Neo4jjl) SET n.name = {name}, n.age = {age}"
  txn(q, "name" => name, "age" => age; submit=submit)
end

@testset "Transactions" begin
    global loadtx = transaction(conn)

    @test length(loadtx.statements) == 0

    createnode(loadtx, "John Doe", 20)

    @test length(loadtx.statements) == 1

    createnode(loadtx, "Jane Doe", 20)

    @test length(loadtx.statements) == 2

    global query = "MATCH (n:Neo4jjl) WHERE n.age = {age} RETURN n.name";
    global people = loadtx(query, "age" => 20; submit=true)

    @test length(loadtx.statements) == 0
    @test length(people.results) == 3
    @test length(people.errors) == 0

    global matchresult = people.results[3]
    @test matchresult["columns"][1] == "n.name"
    @test "John Doe" in [row["row"][1] for row = matchresult["data"]]
    @test "Jane Doe" in [row["row"][1] for row = matchresult["data"]]

    global loadresult = commit(loadtx)

    @test length(loadresult.results) == 0
    @test length(loadresult.errors) == 0

    global query = "MATCH (n:Neo4jjl) WHERE n.age = {age} DELETE n"

    global deletetx = transaction(conn)
    deletetx(query, "age" => 20)

    global deleteresult = commit(deletetx)

    @test length(deleteresult.results) == 1
    @test length(deleteresult.results[1]["columns"]) == 0
    @test length(deleteresult.results[1]["data"]) == 0
    @test length(deleteresult.errors) == 0

    global rolltx = transaction(conn)

    global person = createnode(rolltx, "John Doe", 20; submit=true)

    @test length(rolltx.statements) == 0
    @test length(person.results) == 1
    @test length(person.errors) == 0

    rollback(rolltx)

    global rolltx = transaction(conn)
    global rollresult = rolltx("MATCH (n:Neo4jjl) WHERE n.name = 'John Doe' RETURN n"; submit=true)

    @test length(rollresult.results) == 1
    @test length(rollresult.results[1]["columns"]) == 1
    @test length(rollresult.results[1]["data"]) == 0
    @test length(rollresult.errors) == 0

end

# --- New cypherQuery using transaction/commit endpoint ---

@testset "DataFrames with cypherQuery()" begin

    # Open transaction and create node
    global loadtx = transaction(conn)
    createnode(loadtx, "John Doe", 20; submit=true)
    Neo4j.commit(loadtx)

    global matchresult = cypherQuery(conn,
                      "MATCH (n:Neo4jjl {name: {name}}) RETURN n.name AS Name, n.age AS Age;",
                      "name" => "John Doe")
    @test DataFrames.DataFrame(Name = "John Doe", Age = 20) == matchresult

    # Cleanup
    global deletetx = transaction(conn)
    global query = "MATCH (n:Neo4jjl) WHERE n.age = {age} DELETE n"
    deletetx(query, "age" => 20)
    global deleteresult = commit(deletetx)
end


text:
# ---
# title: 756. Pyramid Transition Matrix
# id: problem756
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Bit Manipulation, Depth-first Search
# link: <https://leetcode.com/problems/pyramid-transition-matrix/description/>
# hidden: true
# ---
#
# We are stacking blocks to form a pyramid. Each block has a color which is a
# one letter string.
#
# We are allowed to place any color block `C` on top of two adjacent blocks of
# colors `A` and `B`, if and only if `ABC` is an allowed triple.
#
# We start with a bottom row of `bottom`, represented as a single string. We
# also start with a list of allowed triples `allowed`. Each allowed triple is
# represented as a string of length 3.
#
# Return true if we can build the pyramid all the way to the top, otherwise
# false.
#
# **Example 1:**
#
#
#
#     Input: bottom = "BCD", allowed = ["BCG", "CDE", "GEA", "FFF"]
#     Output: true
#     Explanation:
#     We can stack the pyramid like this:
#         A
#        / \
#       G   E
#      / \ / \
#     B   C   D
#
#     We are allowed to place G on top of B and C because BCG is an allowed triple.  Similarly, we can place E on top of C and D, then A on top of G and E.
#
#
#
# **Example 2:**
#
#
#
#     Input: bottom = "AABA", allowed = ["AAA", "AAB", "ABA", "ABB", "BAC"]
#     Output: false
#     Explanation:
#     We can't stack the pyramid to the top.
#     Note that there could be allowed triples (A, B, C) and (A, B, D) with C != D.
#
#
#
#
# **Constraints:**
#
#   * `bottom` will be a string with length in range `[2, 8]`.
#   * `allowed` will have length in range `[0, 200]`.
#   * Letters in all strings will be chosen from the set `{'A', 'B', 'C', 'D', 'E', 'F', 'G'}`.
#
#
## @lc code=start
using LeetCode

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


text:
## Getindex default.
## Is this the proper default?
"""
    x[i]

Some SymPy Python objects have index notation provided for them through `__getitem__`. This allows Julia's `getindex` to dispatch to Python's `__getitem__`. The index (indices) must be symbolic. This will use 0-based indexing, as it is a simple pass through to Python.

Examples:

```jldoctest utils
julia> using SymPy

julia> i,j = sympy.symbols("i j", integer=True)
(i, j)

julia> x = sympy.IndexedBase("x")
x

julia> a = sympy.Sum(x[i], (i, 1, j))
  j
 ___
 ╲
  ╲
  ╱   x[i]
 ╱
 ‾‾‾
i = 1
```

"""
function Base.getindex(x::Sym, xs::Sym...)
    if pycall_hasproperty(x, :__getitem__)
        return x.__getitem__(xs)
    elseif pycall_hasproperty(x, :__getslice__)
        return x.__getslice__(xs)
    else

        # no indexing defined, but this won't effect broadcasting
        x
    end
end

## Call
## Call symbolic object with natural syntax
## ex(x=>val)
## how to do from any symbolic object?
(ex::Sym)() = ex

## without specification, variables to substitute for come from ordering of `free_symbols`:
function (ex::Sym)(args...)
    xs = ex.free_symbols
    for (var, val) in zip(xs, args)
        ex = ex.subs(var, Sym(val))
    end
    ex
end

## can use a Dict or pairs to specify:
function (ex::Sym)(x::Dict)
    for (k,v) in x
        ex = ex.subs(k, Sym(v))
    end
    ex
end
function (ex::Sym)(kvs::Pair...)
    for (k,v) in kvs
        ex = ex.subs(k, Sym(v))
    end
    ex
end

##################################################
## subs
##
"""
`subs` is used to subsitute a value in an expression with another
value.
Examples:

```jldoctest subs
julia> using SymPy

julia> x,y = symbols("x,y")
(x, y)

julia> ex = (x-y)*(x+2y)
(x - y)⋅(x + 2⋅y)

julia> subs(ex, (y, y^2))
⎛     2⎞ ⎛       2⎞
⎝x - y ⎠⋅⎝x + 2⋅y ⎠

julia> subs(ex, (x,1), (y,2))
-5

julia> subs(ex, (x,y^3), (y,2))
72

julia> subs(ex, y, 3)
(x - 3)⋅(x + 6)
```

There is a curried form of `subs` to use with the chaining `|>` operator

```jldoctest subs
julia> ex |> subs(x,ℯ)
(ℯ - y)⋅(2⋅y + ℯ)
```
The use of pairs gives a convenient alternative:

```jldoctest subs
julia> subs(ex, x=>1, y=>2)
-5

julia> ex |> subs(x=>1, y=>2)
-5
```

"""
subs(ex::T, y::Tuple{Any, Any}; kwargs...)          where {T <: SymbolicObject} = ex.subs(y[1], Sym(y[2]), kwargs...)
subs(ex::T, y::Tuple{Any, Any}, args...; kwargs...) where {T <: SymbolicObject} = subs(subs(ex, y), args...)
subs(ex::T, y::S, val; kwargs...)                   where {T <: SymbolicObject, S<:SymbolicObject} = subs(ex, (y,val))
subs(ex::T, dict::Dict; kwargs...)                  where {T <: SymbolicObject} = subs(ex, dict...)
subs(ex::T, d::Pair...; kwargs...)                  where {T <: SymbolicObject} = subs(ex, ((p.first, p.second) for p in d)...)
subs(exs::Tuple{T, N}, args...; kwargs...)          where {T <: SymbolicObject, N} = map(u -> subs(u, args...;kwargs...), exs)
subs(x::Number, args...; kwargs...) = x

## curried versions to use with |>
subs(x::SymbolicObject, y; kwargs...) = ex -> subs(ex, x, y; kwargs...)
subs(;kwargs...)                      = ex -> subs(ex; kwargs...)
subs(dict::Dict; kwargs...)           = ex -> subs(ex, dict...; kwargs...)
subs(d::Pair...; kwargs...)           = ex -> subs(ex, [(p.first, p.second) for p in d]...; kwargs...)

##################################################
## doit
##
"""
`doit` evaluates objects that are not evaluated by default.

Examples:

```jldoctest doit
julia> using SymPy

julia> @syms x f()
(x, f)

julia> D = Differential(x)
Differential(x)

julia> df = D(f(x))
d
──(f(x))
dx

julia> dfx = subs(df, (f(x), x^2))
d ⎛ 2⎞
──⎝x ⎠
dx

julia> doit(dfx)
2⋅x
```

Set `deep=true` to apply `doit` recursively to force evaluation of nested expressions:

```jldoctest doit
julia> @syms g()
(g,)

julia> dgfx = g(dfx)
 ⎛d ⎛ 2⎞⎞
g⎜──⎝x ⎠⎟
 ⎝dx    ⎠

julia> doit(dgfx)
 ⎛d ⎛ 2⎞⎞
g⎜──⎝x ⎠⎟
 ⎝dx    ⎠

julia> doit(dgfx, deep=true)
g(2⋅x)
```

There is also a curried form of `doit`:
```jldoctest doit
julia> dfx |> doit
2⋅x

julia> dgfx |> doit(deep=true)
g(2⋅x)
```

"""
doit(ex::T; deep::Bool=false) where {T<:SymbolicObject} = ex.doit(deep=deep)
doit(; deep::Bool=false)      where {T<:SymbolicObject} = ex -> doit(ex, deep=deep)

## simplify(ex::SymbolicObject, ...) is exported
"""
    simplify

SymPy has dozens of functions to perform various kinds of simplification. There is also one general function called `simplify` that attempts to apply all of these functions in an intelligent way to arrive at the simplest form of an expression. (See [Simplification](https://docs.sympy.org/latest/tutorial/simplification.html) for details on `simplify` and other related functionality).

For non-symbolic expressions, `simplify` returns its first argument.
"""
simplify(x,args...;kwargs...) = x

##################################################
##
## access documentation of SymPy
"""
    SymPy.Doc(f::Symbol, [module=sympy])

Return docstring of `f` found within the specified module.

Examples
```
SymPy.Doc(:sin)
SymPy.Doc(:det, sympy.matrices)
## add module to query
SymPy.pyimport_conda("sympy.crypto.crypto", "sympy")
SymPy.Doc(:padded_key, sympy.crypto)
```
"""
struct Doc
    u
    m
end
Doc(u::Union{Symbol, Function}) = Doc(Symbol(u), sympy)
function Base.show(io::IO, u::Doc)
    v = getproperty(u.m, Symbol(u.u)).__doc__
    print(io, v)
end

##################################################
# avoid type piracy. After we call `pytype` mappings, some
# objects are automatically converted and no longer PyObjects
function pycall_hasproperty(x::PyCall.PyObject, k)
    PyCall.hasproperty(x, k) && (getproperty(x,k) != nothing)
end

pycall_hasproperty(x::SymbolicObject, k) = pycall_hasproperty(PyCall.PyObject(x), k)
pycall_hasproperty(x, k) = false

# simple helper for boolean properties
# x.is_finite -> is_(:finite, x)
# e.g.: is_(:FiniteSet, x) = hasproperty && get property
function is_(k::Symbol, x::Sym)::Bool
    key = Symbol("is_$k")
    pycall_hasproperty(x, key) && getproperty(x, key) == Sym(true)
end

_get_members(sm) = PyCall.inspect.getmembers(sm)
function _get_member_functions(sm, exclude=())
    mems = _get_members(sm)
    fns = Dict()
    for (u, v) in mems
        # special cases
        u in exclude && continue

        pycall_hasproperty(v, :deprecated) && continue
        !isa(v, PyCall.PyObject) && continue

        ## we grab functions or FunctionClass objects only
        if PyCall.hasproperty(v, :__class__) &&
            v.__class__.__name__ in ("function", "FunctionClass")
            fns[u] = v
        end
    end
    fns
end

function _is_module(x::Symbol)
    try
        typeof(eval(x)) <:  Module && x ≠ :Main
    catch err
        false
    end
end

function loaded_modules()
    nms = names(Main,imported=true)
    nms = filter(nm -> nm != :SymPy, nms)
    Ms = filter(_is_module,  nms)
    eval.(Ms)
end

# default list of modules to search for namespace collicsions
const base_Ms = (Base, SpecialFunctions, Base.MathConstants,
           LinearAlgebra
           )

# default list of methods to exclude from importing
#
# In addition to issues (such as "C") this should list
#
# * julia methods for which the sympy method is different from julia's generic usage (e.g. `div`)
# * `julia` methods that are clearly issues with the package system, though not in base_Ms (e.g., `plot`, `latex`)
# * sympy methods for which a substitute is used (e.g. `lambdify`)
#
# these are still accesible through dot-call syntax.
#
const base_exclude=("C", "lambdify",
              "latex", "eye", "sympify","symbols", "subs",
              "div", "rem", "log", "sinc",
              "dsolve",
              "ask",
              "plot")

"""
    import_from(module, meths; kwargs...)

Import methods from a python module. Implements functionality  of `from module import function` in  Python.

* `module`: a python module, such as `sympy`
* `meths`: nothing or a tuple of symbols to import. If `nothing`, then all member functions of the module are imported (but not constructors or other objects)
* `Ms`: additional Julia Modules to import from. By default, a few base modules are searched for to avoid namespace collisions.
* `typ`: a symbol indicating variable type for first argument that the new function should be restricted to. For most, the default, `:SymbolicObject` will be appropriate
* `exclude`: when importing all (`meths=nothing`), this can be used to avoid importing some methods by name. The default has a few to avoid.

Examples:

```
import_from(sympy)  # bring in functions from sympy (done `import_sympy`)
import_from(sympy, (:sin, :cos))  # just bring in a few methods
import_from(sympy , (:Wild,), typ=:Any) # Allows `Wild("x")`
#
import PyCall
PyCall.pyimport_conda("sympy.physics.wigner", "sympy")
import_from(sympy.physics.wigner)
```

"""
function import_from(sm, meths=nothing;
                     Ms::Tuple=(),
                     typ::Symbol=:SymbolicObject,
                     exclude::Union{Nothing, NTuple{N,Symbol}}=nothing
                     ) where {N}

    if meths == nothing
        _exclude = isa(exclude, Nothing) ? base_exclude : exclude
        fns = _get_member_functions(sm, _exclude)
    else
        mems = PyCall.inspect.getmembers(sm)
        fns = Dict()
        for (m,p) in mems
            if Symbol(m) in meths
                fns[m] = p
            end
        end
    end
     for (k,v) in fns
         meth = Symbol(k)
         inMs = false
         for M in union(base_Ms, Ms)
             M1 = replace("$M",r"^.*\."=>"")
             if isdefined(M, meth)
                 inMs = true
                 ## @show "import", M, k
##                 println("$M.$k(ex::$typ, args...; kwargs...)=getproperty($(sm.__name__), :$k)(ex, args...; kwargs...)")
                 @eval begin
                     ($M.$meth)(ex::($typ), args...; kwargs...) =
                         getproperty($sm,$k)(ex, Sym.(args)...; kwargs...)
                 end
                 break
             end
        end
        if !inMs
            ## need to export
            ## @show "export", k
##            println("$k(ex::$typ, args...; kwargs...)=getproperty($(sm.__name__), :$k)(ex, args...; kwargs...); export $k")
            @eval begin
                ($meth)(ex::($typ), args...; kwargs...) =
                    getproperty($sm,$k)(ex, args...; kwargs...)
            end
            eval(Expr(:export, meth))
        end
    end

end

"""
    free_symbols(ex)
    free_symbols(ex::Vector{Sym})

Return vector of free symbols of expression or vector of expressions. The results are orderded by
`sortperm(string.(fs))`.

Example:

```jldoctest
julia> using SymPy

julia> @syms x y z a
(x, y, z, a)

julia> free_symbols(2*x + a*y) # [a, x, y]
3-element Vector{Sym}:
 a
 x
 y
```
"""
function free_symbols(ex::Union{T, Vector{T}}) where {T<:SymbolicObject}
    pex = PyObject(ex)
    #fs.__class__.__name__ == "set"
    if PyCall.hasproperty(pex, :free_symbols)
        fs = collect(Sym, pex.free_symbols)
        fs[sortperm(string.(fs))]   # some sorting order to rely on
    else
        Sym[]
    end
end

## Module for Introspection
module Introspection

import SymPy: Sym
import PyCall: PyObject, hasproperty, PyNULL, inspect
export args, func, funcname, class, classname, getmembers

# utilities

"""
    funcname(x)

Return name or ""
"""
function funcname(x::Sym)
    y = PyObject(x)
    if hasproperty(y, :func)
        return y.func.__name__
    else
        return ""
    end
end

"""
   func(x)

Return function head from an expression

[Invariant:](http://docs.sympy.org/dev/tutorial/manipulation.html)

Every well-formed SymPy expression `ex` must either have `length(args(ex)) == 0` or
`func(ex)(args(ex)...) = ex`.
"""
func(ex::Sym) = return ex.func

"""
    args(x)

Return arguments of `x`, as a tuple. (Empty if no `:args` property.)
"""
function args(x::Sym)
    if hasproperty(PyObject(x), :args)
        return x.args
    else
        return ()
    end
end

function class(x::T) where {T <: Union{Sym, PyObject}}
    if hasproperty(PyObject(x), :__class__)
        return x.__class__
    else
        return PyNull()
    end
end

function classname(x::T) where {T <: Union{Sym, PyObject}}
    cls = class(x)
    if cls == PyNULL()
        "NULL"
    else
        cls.__name__
    end
end

function getmembers(x::T) where {T <: Union{Sym, PyObject}}
    Dict(u=>v for (u,v) in inspect.getmembers(x))
end

end


text:
using Pkg
Pkg.activate("./")
using RK

using Plots

function OutputSolution()
    # Output numerical solution in the end
    nz = length(history(tspan[1]))
    if nz == 1
        println("y(", tspan[2], ") = ", sol[2][end])
    else
        for i = 1 : nz
            println("y", i, "(", tspan[2], ") = ", (sol[2]')[end - nz + i])
        end
    end
end

# Speed test
stepsize = 1e-2;

# Example 1
idefun(t,y,z,i) = ((t - 1) * exp(t^2) * i) / (exp(-1) * y[1] - 1)
K(t,s,y)        = y * exp(-s * t)
delays_int(t)   = t - 1
history(t)      = exp(t)

tspan = [1.1 5]

print("Elapsed time:")
@time sol = ide_solve(idefun, K, delays_int, history, tspan, stepsize)

OutputSolution()

# Example 2
idefun(t,y,z,i) = -y.^2 - t * exp(t^2) * z^4 * i
K(t,s,y)        = y * exp(s - s * t)
delays(t,y)     = t / 2
delays_int(t)   = t - 1
history(t)      = exp(-t)

print("Elapsed time:")
@time sol = ide_solve(idefun, K, delays_int, history, tspan, stepsize, delays)

OutputSolution()

# Example 3
idefun(t,y,z,i) = exp(1) - exp(t^2) / (z[1]^2) * (i[1] - exp(-2 * t) * i[2]) * (t - 1)
K(t,s,y)        = [ y * exp(-s * t);
                    y * exp(t * (2 - s))]
delays(t,y)     = t - 1
delays_int(t)   = [ t - 1;
                    t - 2 ]
history(t)      = exp(t)

tspan = [0 5]

print("Elapsed time:")
@time sol = ide_solve(idefun, K, delays_int, history, tspan, stepsize, delays)

OutputSolution()

# Example 4
idefun(t,y,z,i) = -z[1]^((t + 1) / 2) * z[2] * y.^2 * (1 + exp(t^2) * t * i[1]) / exp(0.5)
K(t,s,y)        = y * exp(s - s * t)
delays(t,y)     = [ (log(y[1]))^2 / (t + 1) - 0.5;
                    (t - 1) / 4 ]
delays_int(t)   = t / 2 - 1
history(t)      = exp(-t)

tspan = [0 5]

print("Elapsed time:")
@time sol = ide_solve(idefun, K, delays_int, history, tspan, stepsize, delays)

OutputSolution()

# Example 5
idefun(t,y,z,i) = [ -2.5 * i[1];
                    -15 * i[2];
                    0.001 * exp(t)]
K(t,s,y)        = [ sin(y[1]);
                    sin(y[2])]
delays_int(t)   = [ t - 1;
                    t - 1 ]
history(t)      = [ 1.5;
                    1.5;
                    0.001 * exp(t)]

tspan = [0 10]

print("Elapsed time:")
@time sol = ide_solve(idefun, K, delays_int, history, tspan, stepsize)

OutputSolution()

max_stars_repo_name: vitesempl/IDE-Delay-Runge-Kutta

max_issues_repo_name: vitesempl/IDE-Delay-Runge-Kutta

max_forks_repo_name: vitesempl/IDE-Delay-Runge-Kutta


text:
module Doublets

export Doublet

using ..Elements
using ..Motions

struct Doublet{T <: Number} <: Element
    z::ComplexF64
    S::T
end

Elements.kind(::Doublet) = Singleton
Elements.kind(::Type{Doublet{T}}) where T = Singleton

Elements.position(d::Doublet) = d.z
Elements.streamfunction(z::ComplexF64, d::Doublet) = imag(d.S/(z - d.z)/π)
Elements.complexpotential(z::ComplexF64, d::Doublet) = d.S/(z - d.z)/π

function Motions.induce_velocity(z::ComplexF64, d::Doublet, t)
    conj(-d.S/(z - d.z)^2/π)
end

end


text:
const FILTER_CACHE = Dict{UUID, Tuple{UInt64, Any}}()
const SORTED_CACHE = Dict{UUID, UInt64}()

col_access(row, col) = row[col]

function generate_sorter(params, fixed_col_names)
    lts = []
    for sortspec in params
        let col = sortspec["colId"]
            ind = findfirst(==(col), fixed_col_names)
            if ind === nothing
                continue
            end
            if sortspec["sort"] == "asc"
                push!(lts, (r -> col_access(r, ind), (r1, r2) -> col_access(r1, ind) < col_access(r2, ind)))
            else
                push!(lts, (r -> col_access(r, ind), (r1, r2) -> !(col_access(r1, ind) < col_access(r2, ind))))
            end
        end
    end
    function (row1, row2)
        lt = false
        for (accessor, ltf) in lts
            if ltf(row1, row2)
                return true
            end
            if accessor(row1) != accessor(row2)
                return false
            end
        end
        return lt
    end
end

const filterMapping = Dict{String, Function}(
    "equals" => ==,
    "notEqual" => !=,
    "lessThan" => <,
    "lessThanOrEqual" => <=,
    "greaterThan" => >,
    "greaterThanOrEqual" => >=,
)

function generate_filter(params, col)
    filter_type = params["filterType"]
    if filter_type == "number"
        generate_number_filter(params, col)
    elseif filter_type == "text"
        generate_string_filter(params, col)
    elseif filter_type == "date"
        generate_date_filter(params, col)
    else
        (args...) -> true
    end
end

function generate_number_filter(params, col)
    row -> filterMapping[params["type"]](col_access(row, col), params["filter"])
end

function regex_escape(s::AbstractString)
    res = replace(s, r"([()[\]{}?*+\-|^\$\\.&~#\s=!<>|:])" => s"\\\1")
    replace(res, "\0" => "\\0")
end

function generate_string_filter(params, col)
    filtervalue = regex_escape(params["filter"])
    op = params["type"]
    if op == "equals"
        matcher = Regex("^" * filtervalue * "\$", "i")
        row -> occursin(matcher, string(col_access(row, col)))
    elseif op == "notEqual"
        matcher = Regex("^" * filtervalue * "\$", "i")
        row -> !occursin(matcher, string(col_access(row, col)))
    elseif op == "startsWith"
        matcher = Regex("^" * filtervalue, "i")
        row -> occursin(matcher, string(col_access(row, col)))
    elseif op == "endsWith"
        matcher = Regex(filtervalue * "\$", "i")
        row -> occursin(matcher, string(col_access(row, col)))
    elseif op == "contains"
        matcher = Regex(filtervalue, "i")
        row -> occursin(matcher, string(col_access(row, col)))
    elseif op == "notContains"
        matcher = Regex(filtervalue, "i")
        row -> !occursin(matcher, string(col_access(row, col)))
    end
end

function generate_date_filter(params, col)
    format = dateformat"y-m-d H:M:S"
    dateFrom = Date(params["dateFrom"], format)

    op = params["type"]

    if op == "inRange"
        dateTo = Date(params["dateTo"], format)

        row -> >(col_access(row, col), dateFrom) && <(col_access(row, col), dateTo)
    else
        row -> filterMapping[op](col_access(row, col), dateFrom)
    end
end

function generate_bool(op, cond1, cond2, col)
    let f1 = generate_filter(cond1, col), f2 = generate_filter(cond2, col)
        if op == "AND"
            (args...) -> f1(args...) && f2(args...)
        else
            (args...) -> f1(args...) || f2(args...)
        end
    end
end

function generate_filterer(filters::Dict, fixed_col_names)
    funcs = Function[]
    for (col, filter) in filters
        ind = findfirst(==(col), fixed_col_names)
        if ind === nothing
            continue
        end
        op = get(filter, "operator", "")
        if op in ("AND", "OR")
            push!(funcs, generate_bool(op, filter["condition1"], filter["condition2"], ind))
        else
            push!(funcs, generate_filter(filter, ind))
        end
    end
    return funcs
end

max_stars_repo_name: JaapHaitsma/julia-vscode

max_issues_repo_name: JaapHaitsma/julia-vscode

max_forks_repo_name: JaapHaitsma/julia-vscode


text:
module PoissonDGTests

using Test
using Gridap
import Gridap: ∇
#using LinearAlgebra

#domain = (0,1,0,1)
#partition = (4,4)
#model = CartesianDiscreteModel(domain,partition)
#const h = (domain[2]-domain[1]) / partition[1]

using Gridap.Geometry: DiscreteModelMock
model = DiscreteModelMock()
const h = 1

order = 2
const γ = 10

Ω = Triangulation(model)
Γ = BoundaryTriangulation(model)
Λ = SkeletonTriangulation(model)

degree = order
dΩ = Measure(Ω,degree)
dΓ = Measure(Γ,degree)
dΛ = Measure(Λ,degree)

n_Γ = get_normal_vector(Γ)
n_Λ = get_normal_vector(Λ)

u(x) = x[1]^2 + x[2]
∇u(x) = VectorValue( 2*x[1], one(x[2]) )
Δu(x) = 2
f(x) = - Δu(x)
∇(::typeof(u)) = ∇u

V = TestFESpace(model,ReferenceFE(lagrangian,Float64,order),conformity=:L2)
U = TrialFESpace(V,u)

a(u,v) =
  ∫( ∇(v)⋅∇(u) )*dΩ +
  ∫( (γ/h)*v*u  - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u )*dΓ +
  ∫( (γ/h)*jump(v*n_Λ)⋅jump(u*n_Λ) - jump(v*n_Λ)⋅mean(∇(u)) -  mean(∇(v))⋅jump(u*n_Λ) )*dΛ

l(v) =
  ∫( v*f )*dΩ +
  ∫( (γ/h)*v*u - (n_Γ⋅∇(v))*u )*dΓ

op = AffineFEOperator(a,l,U,V)

uh = solve(op)

e = u - uh

l2(u) = sqrt(sum( ∫( u⊙u )*dΩ ))
h1(u) = sqrt(sum( ∫( u⊙u + ∇(u)⊙∇(u) )*dΩ ))

el2 = l2(e)
eh1 = h1(e)
ul2 = l2(uh)
uh1 = h1(uh)

@test el2/ul2 < 1.e-8
@test eh1/uh1 < 1.e-7

end # module


text:
################################################################################
#                              TEST: Euler Angles
################################################################################

for k = 1:samples
    for rot_seq in rot_seq_array
        # Sample three angles form a uniform distribution [-pi,pi].
        θx = -pi + 2*pi*rand()
        θy = -pi + 2*pi*rand()
        θz = -pi + 2*pi*rand()
        Θ  = EulerAngles(θx, θy, θz, rot_seq)

        δΘ₁ = Θ*inv_rotation(Θ)
        δΘ₂ = inv_rotation(Θ)*Θ

        # Check if the rotation is identity.
        q₁ = angle_to_quat(δΘ₁)
        q₂ = angle_to_quat(δΘ₂)

        @test q₁.q0 ≈ 1 atol=1e-7
        @test q₁.q1 ≈ 0 atol=1e-7
        @test q₁.q2 ≈ 0 atol=1e-7
        @test q₁.q3 ≈ 0 atol=1e-7

        @test q₂.q0 ≈ 1 atol=1e-7
        @test q₂.q1 ≈ 0 atol=1e-7
        @test q₂.q2 ≈ 0 atol=1e-7
        @test q₂.q3 ≈ 0 atol=1e-7

        @test δΘ₁.rot_seq == rot_seq
        @test δΘ₂.rot_seq == Symbol(reverse(String(rot_seq)))
    end
end

# Test promotion
# ==============

@test EulerAngles(1,1,1.0,  :XYZ) === EulerAngles(1.0,1.0,1.0,:XYZ)
@test EulerAngles(1,1,1.0f0,:XYZ) === EulerAngles{Float32}(1,1,1,:XYZ)

# Test printing
# =============

# Color test is not working on AppVeyor. Hence, if the OS is windows, then we
# just skip.
if !Sys.iswindows()
    # With colors.
    expected = """
EulerAngles{Float64}:
\e[32;1m  R(X): \e[0m  1.5708 rad (  90.0000 deg)
\e[33;1m  R(Y): \e[0m  1.0472 rad (  60.0000 deg)
\e[34;1m  R(Z): \e[0m  0.7854 rad (  45.0000 deg)"""

    result = sprint((io,angles)->show(io, MIME"text/plain"(), angles),
                    EulerAngles(π/2, π/3, π/4, :XYZ),
                    context = :color => true)
    @test result == expected
end

# Without colors.
expected = """
EulerAngles{Float64}:
  R(X):   1.5708 rad (  90.0000 deg)
  R(Y):   1.0472 rad (  60.0000 deg)
  R(Z):   0.7854 rad (  45.0000 deg)"""

result = sprint((io,angles)->show(io, MIME"text/plain"(), angles),
                EulerAngles(π/2, π/3, π/4, :XYZ),
                context = :color => false)
@test result == expected

expected = """
EulerAngles{Float64}: R(XYZ):   1.5708   1.0472   0.7854 rad"""

result = sprint(print, EulerAngles(π/2, π/3, π/4, :XYZ))
@test result == expected


text:
using LinearAlgebra

export Scale, factor

"""
    Scale{S <: Union{Number, Val}, D, BT <: AbstractBlock{D}} <: TagBlock{BT, D}

`Scale` a block with scalar. it can be either a `Number` or a compile time `Val`.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> 2 * X
[scale: 2] X

julia> im * Z
[+im] Z

julia> -im * Z
[-im] Z

julia> -Z
[-] Z
```
"""
mutable struct Scale{S<:Union{Number,Val},D,BT<:AbstractBlock{D}} <: TagBlock{BT,D}
    alpha::S
    content::BT
end

content(x::Scale) = x.content
factor(x::Scale{<:Number}) = x.alpha
factor(x::Scale{Val{X}}) where {X} = X

# parameter interface
getiparams(s::Scale{<:Number}) = (factor(s),)
setiparams(s::Scale{<:Number}, alpha::Number) = Scale(alpha, s.content)
setiparams!(s::Scale{T1}, alpha::T2) where {T1<:Number, T2<:Number} = (s.alpha = T1(alpha); s)

Base.copy(x::Scale) = Scale(x.alpha, copy(x.content))
Base.adjoint(x::Scale{<:Number}) = Scale(adjoint(x.alpha), adjoint(content(x)))
Base.adjoint(x::Scale{Val{X}}) where {X} = Scale(Val(adjoint(X)), adjoint(content(x)))

YaoBase.ishermitian(s::Scale) =
    (ishermitian(s |> content) && ishermitian(s |> factor)) || ishermitian(mat(s))
YaoBase.isunitary(s::Scale) =
    (isunitary(s |> content) && isunitary(s |> factor)) || isunitary(mat(s))
YaoBase.isreflexive(s::Scale) =
    (isreflexive(s |> content) && isreflexive(s |> factor)) || isreflexive(mat(s))
YaoBase.iscommute(x::Scale, y::Scale) = iscommute(x |> content, y |> content)
YaoBase.iscommute(x::AbstractBlock, y::Scale) = iscommute(x, y |> content)
YaoBase.iscommute(x::Scale, y::AbstractBlock) = iscommute(x |> content, y)

Base.:(==)(x::Scale, y::Scale) = (factor(x) == factor(y)) && (content(x) == content(y))

chsubblocks(x::Scale, blk::AbstractBlock) = Scale(x.alpha, blk)
cache_key(x::Scale) = (factor(x), cache_key(content(x)))

mat(::Type{T}, x::Scale) where {T} = T(x.alpha) * mat(T, content(x))
mat(::Type{T}, x::Scale{Val{S}}) where {T,S} = T(S) * mat(T, content(x))

function _apply!(r::AbstractArrayReg, x::Scale{S}) where {S}
    _apply!(r, content(x))
    regscale!(r, factor(x))
    return r
end