text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Foundations of machine learning")

# # Minimizing functions - how a computer learns

#-

# In the previous notebooks, we've seen that by changing **parameters** in a function, we could find a "best fit" of that function to some data. We use a **loss function** to quantify the "goodness" of a set of parameters and we look for those parameters that **optimize**, in fact **minimize**, the loss function. If we teach a machine how to minimize the loss function on its own, we say that machine is able to **learn** how to model our data.
#
# In the last notebook, we threw around terms like **derivative**, **gradient**, and **gradient descent** to give you a rough sense of how we minimize a function on a computer. In this notebook, we will step through these concepts more carefully, with the aim of being able to implement them using Julia.

#-

# ## Minimizing a 1D function using calculus

#-

# Let's draw the 1D loss function $L_1$ from a previous notebook again:

σ(x) = 1 / (1 + exp(-x))
f(x, w) = σ(w * x)

x1 = 2
y1 = 0.8

L1(w) = (y1 - f(x1, w))^2

#-

using Plots; gr()

#-

plot(L1, -2, 1.5, xlabel="w", ylabel="L1(w)", leg=false)

# By eye, we can see that the minimum is around $w=0.6$. But how can we get the computer to work this out on its own?

#-

# In the previous notebook, we thought of this function plot as a hill, viewed from the side. We could find the minimum by making the hill sticky, and letting a ball roll down it. The ball will find and settle in the minimum of the function.  Now let's see how to teach a computer to do this.
#
# We need to find the downhill direction along the hill, which is related to its *slope* (how steep it is). Calculus provides us with tools to calculate that slope!

#-

# Namely, the slope of a curve $L_1(w)$ at $w$ is given by its **derivative** $L_1'(w)$; geometrically, this is the slope of the **tangent line** to the curve at that point, i.e. the straight line which touches the curve at that point.

#-

# Calculus provides us with some rules to calculate an analytical formula for the derivative, and we will see later how to apply these rules, albeit indirectly, for machine learning.
# To gain understanding, however, we will see here how to get the computer to help us out by calculating the derivatives numerically instead!

#-

# ## Approximating derivatives

#-

# Let's recall that the derivative $L_1'(w)$ of a function is defined as
#
# $$L_1'(w) \simeq \frac{L_1(w+h) - L_1(w)}{h},$$
#
# for a small step size $h$. (Strictly speaking, we must take the limit when $h$ tends to $0$ to obtain the exact value of the derivative.)

#-

# #### Exercise 1
#
# Write a function to calculate the derivative of a function at a given point. Note that in Julia, we can easily pass functions as arguments to other functions!
# The function should take the function, the point $w$ at which to calculate the derivative, and the value of $h$, which should have a default value of 0.001.
#
# *Note*: A function input argument can have a *default value* if we set the input argument equal to that default value when we define the function. For example,
#
# ```julia
# f(x, a = 3) = a * x^2
# ```
#
# The function `f` will square the input we give it and multiply by `a`. However,  if we choose to call `f(x)` *without* passing it an `a` input, it will assume `a` is `3` and return `3*x^2`.

#-

# #### Solution
#
# We could define the `derivative` function as

function derivative(f, w, h=0.001)
    return ( f(w + h) - f(w) ) / h
end
## Or the short form with:
## derivative(f, w, h=0.001) = ( f(w + h) - f(w) ) / h

# **Tests**

@assert derivative(x->x^2, 0.2) ≈ 0.401    # check default h works
@assert derivative(x->x^2, 0.2, 0.5) ≈ .9  # check we can specify our own h
@assert derivative(x->x^3, 0.2) ≈ 0.120601
@assert derivative(x->x^3, -0.2) ≈ 0.119401

# #### Exercise 2
#
# Write an interactive visualization of the tangent line to the graph of $L_1$, so that we can visualize the tangent at any point on $L_1$. Include the current value of the derivative in the title.
#
# *Hint*: Recall that the straight line through the point $(x_0, y_0)$ with slope $m$ is given by
#
# $$ \frac{y - y_0}{x - x_0} = m,$$
#
# so
#
# $$y = y_0 + m*(x - x_0).$$

#-

# #### Solution

wmin = -2
wmax = 1

w = 0.5 # Try manipulating w between wmin and wmax (-2 to 1)
plot(L1, wmin, wmax, label="L_1")

x0 = w
y0 = L1(w)

scatter!([x0], [y0], label="")

m = derivative(L1, w)

plot!(x->y0 + m*(x - x0), ls=:dash, label="tangent line")

ylims!(-0.1, 0.6)
xlims!(wmin, wmax)

title!("Derivative = $m")
xlabel!("w")
ylabel!("L1(w)")

# #### Exercise 3
#
# What is the value of the derivative (slope of the tangent line) at a minimum? Can this happen anywhere else?

#-

# #### Solution

#-

# The derivative is `0` at a minimum. This happens at any minimum, maximum or "point of inflection" of a smooth function.

#-

# ## Minimization by gradient descent

#-

# The tangent line at a point is a good approximation of the function near that point. Therefore the derivative tells us in which direction, and how fast, a function grows or shrinks when we move a small distance away from that point.
#
# As we saw in the previous notebook, we can think of the function as being a hill, and having a ball that we want to move down the hill. Gravity will automatically pull the ball in the direction *down* the hill; we want to emulate this in the computer!

#-

# #### Exercise 4
#
# When the derivative $L_1'(w)$ is positive, that means that $L_1$ increases from left to right at point $w$.
#
# If the derivative $L_1'(w)$ is positive (> 0), in which direction should we move $w$ to *decrease* $L_1$?

#-

# #### Solution
#
# We should move $w$ left, i.e. decrease it.

#-

# #### Exercise 5
#
# If the derivative $L_1'(w)$ is negative (< 0), in which direction should we move $w$ to *decrease* $L_1$?

#-

# We can use this information to tell the computer which way to take a step. This constitutes the numerical algorithm called **gradient descent**; it is called this since we are descending (moving downwards) along the graph of the function by using information about its gradient (slope).

#-

# #### Exercise 6
#
# Implement gradient descent by following this prescription for an **iterative (repetitive) algorithm**:
#
# 1. Start at an initial guess $w_0$ for $w$.
#
# 2. At each step, calculate the derivative, $L_1'(w_n)$ at the current value of $w_n$ using the function that you created above.
#
# 3. Modify the value of $w$ by a small multiple (for example, $\eta=0.01$) of the value of the derivative you just created, via $w_{n+1} = w_n - \eta L_1'(w_n)$.
#
# For this problem, start with $w_0 = -2.0$. Repeat steps 2 and 3 a total of `2000` times.
#
# Package this code into a function called `gradient_descent` that takes two inputs, a function and a range for values of $w$, and returns the final value of $w$ and $L_1(w)$.
#
# Using `L1` and `-2:0.01:1.5` as your inputs to `gradient_descent`, for what value of $w$ is $L_1$ at a minimum?

#-

# #### Solution

function gradient_descent(f, wrange)
    derivative(f, w, h=0.001) = ( f(w + h) - f(w) ) / h
    w = -2.0
    for i in 1:2000
        fderiv = derivative(f, w)
        step = .01 * fderiv
        w = w - step
    end
    return w, f(w)
end

gradient_descent(L1, -2:0.01:1.5)

# **Tests**

@assert 0.6383 < gradient_descent(L1, -2:0.01:1.5)[1] < 0.6385
@assert 0.00032 < gradient_descent(L1, -2:0.01:1.5)[2] < 0.00033

# #### Exercise 7
#
# Modify your code for gradient descent to return the result once it has found an answer within some *tolerance*, rather than taking a set number of steps. The new prescription for this algorithm is:
#
# 1. Start at an initial guess $w_0$ for $w$.
#
# 2. At each step, calculate the derivative, $L_1'(w_n)$ at the current value of $w_n$, using the function that you created above.
#
# 3. Modify the value of $w$ by a small multiple (for example, $\eta=0.01$) of the value of the derivative you just created, via $w_{n+1} = w_n - \eta L_1'(w_n)$.
#
# 4. Check how different $w_{n+1}$ is from $w_n$. If you're satisfied that $L_1(w_{n+1})$ is minimized, return $w_{n+1}$ and $L_1(w_{n+1})$. Otherwise, go to step (2) and continue.
#
# Edit `gradient_descent` so that it takes three inputs: a function, a range for values of $w$, and a tolerance that tells you how close $w_{n+1}$ must be to $w_n$ before you can stop iterating.
#
# Using `L1`, `-2:0.01:1.5`, and `.000001` as your inputs to `gradient_descent`, for what value of $w$ is $L_1$ at a minimum?

#-

# #### Solution
#
#
# Your code should look something like this
#
# $L_1$ is at a minimum when `gradient_descent(L1, -2:0.01:1.5, .000001)[1]`, or $w$, is around `0.69216`.

function gradient_descent(f, wrange, tol)
    w = -2.0
    while true
        fderiv = derivative(f, w)
        step = .01 * fderiv
        w = w - step
        if abs(step) < tol
           return w, f(w)
        end
    end
end

gradient_descent(L1, -2:0.01:1.5, .000001)

# **Tests**

@assert 0.69215 < gradient_descent(L1, -2:0.01:1.5, .000001)[1] < 0.69217

# After you've completed this step, you will have found a good approximation to the value $w^*$ where $L_1$ reaches its minimum, and the minimum value $L_1(w^*)$, which in this case is (almost) $0$.

#-

# #### Exercise 8
#
# Alter the function `gradient_descent` so that it stores the results `(w, L1(w))` at each step of the algorithm as an array and returns this array. How many steps does the algorithm take for input parameters `(L1, -2:0.01:1.5, .000001)` before terminating? You should count your starting $w_0$ as your first step.

#-

# #### Solution

#-

# Your altered code should look like this, getting the number of steps with `length`:

function gradient_descent(f, wrange, tol)
    w = -2.0
    store_w = [(w, f(w))]
    while true
        fderiv = derivative(f, w)
        step = .01 * fderiv
        w = w - step
        push!(store_w, (w, f(w)))
        if abs(step) < tol
           return store_w
        end
    end
end

length(gradient_descent(L1, -2:0.01:1.5, .000001))

# **Tests**

@assert length(gradient_descent(L1, -2:0.01:1.5, .000001)) == 4255

# #### Exercise 9
#
# Overlay the steps taken in the last exercise with a plot of $L_1(w)$ vs. $w$.
#
# Where does our algorithm take the largest steps? (Where does the ball move fastest down the hill?)
#
# A) Between w = -2:-1

# B) Between w = -1:0

# C) Between w = 0:.6
#
# *Hint*: It may be easier to see what's going on if you only plot, for example, every 15th step.

#-

# #### Solution
#
# B) The following code generates the desired figure:

wsteps = gradient_descent(L1, -2:0.01:1.5, .000001)
fig1 = plot(L1, -2, 1.5, xlabel="w", ylabel="C1(w)", leg=false)
for step in wsteps[1:15:end]
   scatter!([step[1]], [step[2]])
end
fig1

# ## Functions of 2 variables and their derivatives

#-

# So far, we have been looking at the minimizing a loss function $L_1(w)$ that depends on a single parameter, $w$. Now let's turn to the cost function $L_2(w, b)$ from a previous notebook, that is a function of *two* parameters, $w$ and $b$, and try to minimize it.
#
# As we've seen, we get a **surface**, instead of a curve, when we graph $L_2$ as a function of both of its parameters.

#-

# **Exercise 10**
#
# Draw a surface plot of $L_2$, given by
#
# $$L_{2}(w, b) = \sum_i(y_i - g(x_i, w, b))^2$$
#
# using the `surface` function from `Plots.jl`. For this plot, use the values of `xs` and `ys` from notebook 5:
#
# We can get a nice interactive 3D plot by using the Plotly backend of `Plots.jl` by executing the command
#
#     plotly()

xs = [2, -3, -1, 1]
ys = [0.8, 0.3, 0.4, 0.4]

# #### Solution

#-

# The following code will generate the desired graphic:

plotly()
σ(x) = 1 / (1 + exp(-x))
g(x, w, b) = σ(w*x) + b

xs = [2, -3, -1, 1]
ys = [0.8, 0.3, 0.4, 0.4]

L2(w, b) = sum(abs2, ys .- g.(xs, w, b))

surface(-2:0.02:2, -2:0.02:2, L2, alpha=0.6, xlabel = "w", ylabel = "b", zlabel = "L2(w, b)")

# ### Finding the minimum

#-

# We can just about see, by rotating the graph, that $L_2$ has a single minimum. We want to find the values of $w$ and $b$ where this minimum value is reached.

#-

# Following what we did for the function $L_1$ above, we expect that we will need to calculate derivatives of $L_2$. Since the function is more complicated, though, the derivatives are too!
#
# It turns out that the right concept is that of the [**gradient**](https://en.wikipedia.org/wiki/Gradient) of $L_2$, denoted $\nabla L_2(w, b)$. This is a **vector** consisting of $2$ numbers if there are $2$ parameters [or, in general, $n$ numbers if there are $n$ parameters].
#
# The numbers that form the gradient $\nabla L_2(w, b)$ are called the **partial derivatives** of $L_2$ with respect to $w$ and $b$,  written as
#
# $$\frac{\partial L_2}{\partial w} \quad \text{and} \quad \frac{\partial L_2}{\partial b}.$$
#
# Although this notation might look complicated, all it means is that we calculate derivatives just like we did before, except that we fix the value of the other variable.
# For example, to calculate $\frac{\partial L_2}{\partial w}$, the partial derivative of $L_2$ with respect to $w$, we fix the value of $b$ and think of the resulting function as a function of the single variable $w$; then we use the formula for derivatives of functions of a single variable.
#
# [Note that $\frac{\partial L_2}{\partial w}$ is itself a function of $w$ and $b$; we could write $\frac{\partial L_2}{\partial w}(w, b)$.]

#-

# #### Exercise 11
#
# Write functions that will allow you to calculate the partial derivatives of a function $f$ of two variables.
#
# In particular, declare functions called `partial_w` and `partial_b`. Each should take four inputs - a function $f$ of two variables, the first input argument to $f$, the second input argument to $f$, and a step size `h` with default value `0.001`. `partial_w` should return the partial derivative of $f$ with respect to its first input argument and `partial_b` should return the partial derivative of $f$ with respect to its second input argument.

#-

# #### Solution
#
# The partial derivatives are

partial_w(f, w, b, h = 0.001) = ( f(w + h, b) - f(w, b) ) / h
partial_b(f, w, b, h = 0.001) = ( f(w, b + h) - f(w, b) ) / h

# **Tests**

## Test that default parameters work:
@assert partial_w((x, y) -> 3 * x * y, 0.3, 0.3) ≈ 0.9
@assert partial_b((x, y) -> 3 * x * y, 0.3, 0.3) ≈ 0.9
## Test that you can specify your own step sizes:
@assert partial_w((x, y) -> 3 * x^2 * y^3, 2.3, 10.2, 1) ≈ 17828.2944
@assert partial_b((x, y) -> 3 * x^2 * y^3, 2.3, 10.2, 1) ≈ 5454.8364
## Test that x and y are treated separately:
@assert partial_w((x, y) -> 3 * x * y^2, 0.3, 0.2) ≈ 0.12
@assert partial_b((x, y) -> 3 * x * y^2, 0.3, 0.2) ≈ 0.3609
@assert partial_w((x, y) -> 3 * x^2 * y, 0.3, 0.2, 0.1) ≈ 0.42
@assert partial_b((x, y) -> 3 * x^2 * y, 0.3, 0.2, 0.1) ≈ 0.27

# #### Exercise 12
#
# Use `partial_b` from the last exercise to find the partial derivative of $L_2$ with respect to $w$ at b = 0.3, $\frac{\partial L_2}{\partial w}|_{b = 0.3}$ for `w = -2:0.01:1`

#-

# #### Solution

wrange = -2:0.01:1
partialw = partial_w.(L2, wrange, 0.3 )

# **Tests**

@assert isapprox(partialw[3], -0.212573, rtol=2e-6)

# #### Exercise 13
#
# Plot the cross section of the surface of $L_2(w, b)$ at $b = 0.3$. Make this plot interactive to show that the function `partial_w` gives the slope of the tangent to this cross section for any point `w` in the range `-2:0.01:1`.
#
# For what value of $w$ in this range is the slope of the cross section closest to -1?

#-

# #### Solution

#-

# The slope is closest to -1.0 at w = -0.55. You can see this by running the following code:

L2_fixedb = L2.(wrange, 0.3)

w = 0.0 # Try manipulating w to pick a value within `wrange`
plot(wrange, L2_fixedb)
scatter!([w], [L2(w, 0.3)])
slope = partial_w(L2, w, 0.3)
tangent_points = slope .* wrange .- slope * w .+ L2.(w, 0.3)
plot!(wrange, tangent_points)
ylims!(0.3, 2.0)
title!("Slope = $slope @ $(w)")

# ## ***Optional**: Functions with $n$ inputs

#-

# If a function $f$ takes $n$ input arguments, we can write them as $p_1, \ldots, p_n$, where $p_i$ means the "$i$th parameter". In Julia, we can wrap them up into a single **vector**. Now we can calculate the partial derivative $\frac{\partial L_2}{\partial p_i}$ with respect to the $i$th variable.

#-

# #### Exercise 14
#
# For the next exercise, you will need to use the splat command, `...`. You can use this command to "open up" a collection and pass all the elements of that collection as inputs to a function.
#
# For example, say you have an array, `numbers`,
#
# ```julia
# numbers = [4, 3, 2]
# ```
#
# and you want to use `numbers` to create a $4\times3\times3$ randomly populated array via `rand`. `rand(numbers)` will not do what you want. You could index into `numbers` to grab the values you want and pass them to `rand`, as in
#
# ```julia
# rand(numbers[1], numbers[2], numbers[3])
# ```
#
# or you could use a splat:
#
# ```julia
# rand(numbers...)
# ```
#
# Use `...` to pass the contents of `inputs`
#
# ```julia
# inputs = [30, 12, "cats"]
# ```
#
# to the function `dreams`
#
# ```julia
# dreams(i, j, perfect_mammal) = "I wish I had $(i + j) $perfect_mammal."
# ```

#-

# #### Solution

inputs = [30, 12, "cats"]
dreams(i, j, perfect_mammal) = "I wish I had $(i + j) $perfect_mammal."
dreams(inputs...)

# **Tests**

@assert dreams(inputs...) == "I wish I had 42 cats."

# #### Exercise 15:
#
# Write a function, `partial`, to calculate the $i$th partial derivative of a function. This function should have four inputs
#
# * a function, $f$, for which you want to compute the partial derivative
# * an array, *p*, specifying the values of all input arguments to $f$ at the point where you want $\frac{\partial f}{\partial p_i}$ computed
# * the index, $i$, of the variable with respect to which you want to calculate the partial derivative of $f$
# * a step size with default value 0.001
#
# Hint: you will need to `copy` and modify `p` within `partial`.

#-

# #### Solution

function partial(f, p, i, h = 0.001)
    p_next = copy(p)
    p_next[i] += h
    return ( f(p_next...) - f(p...) ) / h
end

# **Tests**

@assert partial((x, y, z) -> x^2 * y * z, [1.0, 2.0, 1.0], 1, 0.0001) ≈ 4.0002
@assert partial((x, y, z) -> x^2 * y * z, [1.0, 2.0, 1.0], 1, 1) ≈ 6

@assert partial((w, x, y, z) -> w * x^2 * y^3 * z^4, [1.1, 2.2, 3.3, 4.4], 1) ≈ 65192.53776074038
@assert partial((w, x, y, z) -> w * x^2 * y^3 * z^4, [1.1, 2.2, 3.3, 4.4], 2) ≈ 65207.35424659506
@assert partial((w, x, y, z) -> w * x^2 * y^3 * z^4, [1.1, 2.2, 3.3, 4.4], 3) ≈ 65212.29507069802
@assert partial((w, x, y, z) -> w * x^2 * y^3 * z^4, [1.1, 2.2, 3.3, 4.4], 4) ≈ 65214.7658571339

# ## Gradient descent in 2 dimensions

#-

# It turns out that the gradient vector of the function $L_2(w, b)$ gives the direction in the plane $(w, b)$ in which the function $L_2$ **increases fastest**.
#
# In other words, if you start at the position $(w_0, b_0)$ and take a small step of length $\eta$ to new values $(w', b')$, the value of the function will change to a new value $L_2(w', b')$. How do we get to the minimum of $L_2(w, b)$ as fast as possible? We want to step in the direction where $L_2$ decreases fastest!
#
# In multivariable calculus courses, it is shown that $L_2$ will *increase* fastest if you take a step **in the direction of the gradient $\nabla L_2(w, b)$**! To decrease $L_2$ the fastest, we should take a step in the *opposite* direction, $-\nabla L_2(w, b)$.
#
# Let's now generalize the gradient descent algorithm that we wrote previously to work our 2-dimensional function.

#-

# #### Exercise 16
#
# Extend your 1D implementation of `gradient_descent` so that it will minimize the function $L_2$.
#
# Requirements:
#
# * Your new method for `gradient_descent` will take four input arguments: the function $f$ for which you seek the minimum, the range of values for $f$'s first input argument that you will consider, the range of values for $f$'s second input argument that you will consider, and a tolerance that will specify the maximum allowable step size, $\sum_i \eta \frac{\partial f}{\partial p_i}$
#
# * Use $\eta = .01$. For example, for a function $f(w, b)$, update $w$ such that $w_{n+1} = w_n - 0.01 * \frac{\partial f}{\partial w_n}$
#
# * Seed `gradient_descent` with the starting coordinates [-2.0, -2.0], i.e. $w_0 = -2.0$ and $b_0 = -2.0$.
#
# * Return all steps (their coordinates) taken during gradient descent and the values of the loss function at these coordinates.
#
# Once you have done this, execute
#
# ```julia
# gradient_descent(L2, -2:0.02:2, -2:0.02:2, .001)
# ```
#
# How many steps were taken by gradient descent?
#
# Hint: Do not count your starting coordinates `[-2.0, -2.0]` as a step.

#-

# #### Solution
#
#
# The new code for `gradient_descent` should look something like

function gradient_descent(f, wrange, brange, tol)
    w, b = -2.0, -2
    store_steps = [[w, b]]
    store_C2_vals = [f(w, b)]
    while true
        # Calculate partial derivatives
        fpartial1 = partial(f, [w, b], 1)
        fpartial2 = partial(f, [w, b], 2)
        # Update parameters
        w = w - .01 * fpartial1
        b = b - .01 * fpartial2
        # Store parameters and resulting function
        push!(store_steps, [w, b])
        push!(store_C2_vals, f(w, b))

        if abs(w + b) < tol
           return store_steps, store_C2_vals
        end
    end
end

# **Tests**

steps, values = gradient_descent(L2, -2:0.02:2, -2:0.02:2, .001 )
@assert length(values) == 396

# #### Exercise 17
#
# Use the `surface` and `scatter!` commands to illustrate the path taken by `gradient_descent` from [-2.0, -2.0] to wherever the algorithm terminates.
#
# Where do the scattered points representing the steps of gradient descent appear the most dense?
#
# A) Near the starting point at [-2.0, -2.0]

# B) Near the point [-1.8, 0] where $C_2$ appears nearly flat

# C) Near the minimum of $C_2$

#-

# #### Solution
#
# B) The following code will plot the function `C2` with an overlay of the steps taken in gradient descent
#
# You can change the frequency with which steps are plotted by changing the iterator in the `for` loop from `1:length(values)` to `1:n:length(values)` for varying values of `n`. This will allow you to see more clearly where steps are densely or sparsely packed.

fig = surface(-2:0.02:2, -2:0.02:2, L2, alpha=0.6, xlabel = "w", ylabel = "b", zlabel = "C2(w, b)")
for i in 1:length(values)
   scatter!([steps[i][1]], [steps[i][2]], [values[i]])
end
fig

# ## What we have learnt

#-

# To recap, in this notebook, we have seen that the computer can calculate approximate derivatives, and that we can use those inside a relatively simple algorithm, gradient descent, to minimize functions of 1 and 2 variables!
#
# Now that we've used gradient descent to **learn** the best values of `b` and `w` for a cost function, we can imagine how we could do our classification of apples vs. bananas by allowing the machine to choose our parameters! To do this, we need to define a loss function for apple vs. banana classification, using images that have already been categorised as apples or bananas.

max_stars_repo_path: Courses/Foundations of machine learning/0900.What-is-learning-quiz.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
module ErrTests

 include("../src/AndersonMoore.jl")
using .AndersonMoore

function noErrors()::Bool

    hh = [0.  0.   0.  0. -1.1  0.  0.  0.  1.  1.  0.  0.;
          0. -0.4  0.  0.  0.   1. -1.  0.  0.  0.  0.  0.;
          0.  0.   0.  0.  0.   0.  1.  0.  0.  0.  0.  0.;
          0.  0.   0. -1.  0.   0.  0.  1.  0.  0.  0.  0.]::Array{Float64, 2}

    e = "AndersonMoore: unique solution."

    (bnow,rtsnow,ia,nexact,nnumeric,lgroots,AMAcode) =
        AndersonMooreAlg(hh, 4, 1, 1, 1.0e-8, 1.0 + 1.0e-8)

    err(AMAcode) == e

end # noErrors()

function tooManyRoots()::Bool

    hh = [0.  0.  0.  0. -1.1  0.  0.  0.  1.  1.  0.  0.;
          0.  4.  0.  0.  0.   1. -1.  0.  0.  0.  0.  0.;
          0.  0.  0.  0.  0.   0.  1.  0.  0.  0.  0.  0.;
          0.  0.  0. -1.  0.   0.  0.  1.  0.  0.  0.  0.]::Array{Float64, 2}

    e = "AndersonMoore: too many big roots."

    (bnow, rtsnow, ia, nexact, nnumeric, lgroots, AMAcode) =
        AndersonMooreAlg(hh, 4, 1, 1, 1.0e-8, 1.0 + 1.0e-8)

    err(AMAcode) == e

end # tooManyRoots()

function tooFewRoots()::Bool

    hh = [0.  0.   0.  0. -0.9  0.  0.  0.  1.  1.  0.  0.;
          0. -0.4  0.  0.  0.   1. -1.  0.  0.  0.  0.  0.;
          0.  0.   0.  0.  0.   0.  1.  0.  0.  0.  0.  0.;
          0.  0.   0. -1.  0.   0.  0.  1.  0.  0.  0.  0.]::Array{Float64, 2}

    e = "AndersonMoore: too few big roots."

    (bnow, rtsnow, ia, nexact, nnumeric, lgroots, AMAcode) =
        AndersonMooreAlg(hh, 4, 1, 1, 1.0e-8, 1.0 + 1.0e-8)

    err(AMAcode) == e

end # tooFewRoots()

function tooManyExactShifts()::Bool

    hh = [0.  0.   0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
          0. -0.4  0.  0.  0.  1. -1.  0.  0.  0.  0.  0.;
          0.  0.   0.  0.  0.  0.  1.  0.  0.  0.  0.  0.;
          0.  0.   0. -1.  0.  0.  0.  1.  0.  0.  0.  0.]::Array{Float64, 2}

    e = "AndersonMoore: too many exact shiftrights."

    (bnow, rtsnow, ia, nexact, nnumeric, lgroots, AMAcode) =
        AndersonMooreAlg(hh, 4, 1, 1, 1.0e-8, 1.0 + 1.0e-8)

    err(AMAcode) == e

end # tooManyExactShifts()

function tooManyNumericShifts()::Bool

    hh = [0.  0.  0.  0.  -2.2  0.  0.  0.  2.  2.  0.  0.;
          0.  0.  0.  0.  -1.1  0.  0.  0.  1.  1.  0.  0.;
          0.  0.  0.  0.  -1.1  0.  1.  0.  1.  1.  0.  0.;
          0.  0.  0. -1.  -1.1  0.  0.  1.  1.  1.  0.  0.]::Array{Float64, 2}

    e = "AndersonMoore: too many numeric shiftrights."

    (bnow, rtsnow, ia, nexact, nnumeric, lgroots, AMAcode) =
        AndersonMooreAlg(hh, 4, 1, 1, 1.0e-8, 1.0 + 1.0e-8)

    err(AMAcode) == e

end # tooManyNumericShifts()

function spurious()::Bool

    e = "AndersonMoore: return code not properly specified"

    err(975) == e

end # spurious

end # module

ext: jl


text:
############################################################
# jo Types #################################################
############################################################

############################################################
# includes
include("joTypes/joMiscTypes.jl")
include("joTypes/joAbstractTypes.jl")
include("joTypes/joAbstractLinearOperator.jl")
    include("joTypes/joAbstractParallelableLinearOperator.jl")
include("joTypes/joAbstractFosterLinearOperator.jl")
    include("joTypes/joAbstractLinearOperatorInplace.jl")
include("joTypes/joPAsetup.jl")
include("joTypes/joAbstractSAparallelToggleOperator.jl")
include("joTypes/joAbstractSAparallelLinearOperator.jl")
include("joTypes/joAbstractDAparallelToggleOperator.jl")
include("joTypes/joAbstractDAparallelLinearOperator.jl")


text:
include("modules.jl")

struct VNet;
    block_one;
    block_one_dw;
    block_two;
    block_two_dw;
    block_three;
    block_three_dw;
    block_four;
    block_four_dw;

    block_five;
    block_five_up;
    block_six;
    block_six_up;
    block_seven;
    block_seven_up;
    block_eight;
    block_eight_up;

    block_nine;
    out_conv;
    out_conv2;
end

function VNet(in_channels::Int=3, n_classes::Int=2, n_filters::Int=16, norm::Bool=false)
    if !norm
    println("No normalization is used.")
    VNet(
        RConvBlock(1, in_channels, n_filters),  # block_one
        DownsamplingRConv3d(n_filters, 2*n_filters, 2),

        RConvBlock(2, 2*n_filters, 2*n_filters),    # block_two
        DownsamplingRConv3d(2*n_filters, 4*n_filters, 2),

        RConvBlock(3, 4*n_filters, 4*n_filters),    # block_three
        DownsamplingRConv3d(4*n_filters, 8*n_filters, 2),

        RConvBlock(3, 8*n_filters, 8*n_filters),    # block_four
        DownsamplingRConv3d(8*n_filters, 16*n_filters, 2),

        RConvBlock(3, 16*n_filters, 16*n_filters),  # block_five
        UpsamplingRConvTranspose3d(16*n_filters, 8*n_filters, 2),

        RConvBlock(3, 8*n_filters, 8*n_filters),  # block_six
        UpsamplingRConvTranspose3d(8*n_filters, 4*n_filters, 2),

        RConvBlock(3, 4*n_filters, 4*n_filters),  # block_seven
        UpsamplingRConvTranspose3d(4*n_filters, 2*n_filters, 2),

        RConvBlock(2, 2*n_filters, 2*n_filters),  # block_eight
        UpsamplingRConvTranspose3d(2*n_filters, n_filters, 2),

        RConvBlock(1, n_filters, n_filters),
        Conv3d(n_filters, n_classes, 1, 0, 1)
        Conv3d(n_filters, n_classes, 1, 0, 1)
    )
    else
    println("Batch normalization is used.")
    VNet(
        BatchConvBlock(1, in_channels, n_filters),  # block_one
        DownsamplingBatchConv3d(n_filters, 2*n_filters, 2),

        BatchConvBlock(2, 2*n_filters, 2*n_filters),    # block_two
        DownsamplingBatchConv3d(2*n_filters, 4*n_filters, 2),

        BatchConvBlock(3, 4*n_filters, 4*n_filters),    # block_three
        DownsamplingBatchConv3d(4*n_filters, 8*n_filters, 2),

        BatchConvBlock(3, 8*n_filters, 8*n_filters),    # block_four
        DownsamplingBatchConv3d(8*n_filters, 16*n_filters, 2),

        BatchConvBlock(3, 16*n_filters, 16*n_filters),  # block_five
        UpsamplingBatchConvTranspose3d(16*n_filters, 8*n_filters, 2),

        BatchConvBlock(3, 8*n_filters, 8*n_filters),  # block_six
        UpsamplingBatchConvTranspose3d(8*n_filters, 4*n_filters, 2),

        BatchConvBlock(3, 4*n_filters, 4*n_filters),  # block_seven
        UpsamplingBatchConvTranspose3d(4*n_filters, 2*n_filters, 2),

        BatchConvBlock(2, 2*n_filters, 2*n_filters),  # block_eight
        UpsamplingBatchConvTranspose3d(2*n_filters, n_filters, 2),

        BatchConvBlock(1, n_filters, n_filters),
        Conv3d(n_filters, n_classes, 1, 0, 1)
        Conv3d(n_filters, n_classes, 1, 0, 1)
    )
    end

end

function (c::VNet)(x)
    # encode
    x1 = c.block_one(x)
    xx = c.block_one_dw(x1)

    x2 = c.block_two(xx)
    xx = c.block_two_dw(x2)

    x3 = c.block_three(xx)
    xx = c.block_three_dw(x3)

    x4 = c.block_four(xx)
    xx = c.block_four_dw(x4)

    x5 = c.block_five(xx)
    x5 = dropout(x5, 0.5)

    # decode
    xx = c.block_five_up(x5)
    xx = xx .+ x4

    x6 = c.block_six(xx)
    xx = c.block_six_up(x6)
    xx = xx .+ x3

    x7 = c.block_seven(xx)
    xx = c.block_seven_up(x7)
    xx = xx .+ x2

    x8 = c.block_eight(xx)
    xx = c.block_eight_up(x8)
    xx = xx .+ x1

    x9 = c.block_nine(xx)
    x9 = dropout(x9, 0.5)

    # output layer
    out = c.out_conv(x9)
    out_tanh = tanh.(x9)
    out_seg = c.out_conv2(x9)
    return out_tanh, out_seg
end

# model = VNet(1, 3, 2)
# x = ones(32, 32, 32, 1, 11);
# y = model(x);
# size(y)

max_stars_repo_name: skoc/julia-project

max_issues_repo_name: skoc/julia-project

max_forks_repo_name: skoc/julia-project


text:
using Test
using ModelingToolkit, OrdinaryDiffEq

# Basic electric components
const t = Sym{ModelingToolkit.Parameter{Real}}(:t)
function Pin(;name)
    @variables v(t) i(t)
    ODESystem(Equation[], t, [v, i], [], name=name, defaults=[v=>1.0, i=>1.0])
end

function Ground(name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    ODESystem(eqs, t, [], [], systems=[g], name=name)
end

function ConstantVoltage(name; V = 1.0)
    val = V
    @named p = Pin()
    @named n = Pin()
    @parameters V
    eqs = [
           V ~ p.v - n.v
           0 ~ p.i + n.i
          ]
    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)
end

function Resistor(name; R = 1.0)
    val = R
    @named p = Pin()
    @named n = Pin()
    @variables v(t)
    @parameters R
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           v ~ p.i * R
          ]
    ODESystem(eqs, t, [v], [R], systems=[p, n], defaults=Dict(R => val), name=name)
end

function Capacitor(name; C = 1.0)
    val = C
    @named p = Pin()
    @named n = Pin()
    @variables v(t)
    @parameters C
    D = Differential(t)
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           D(v) ~ p.i / C
          ]
    ODESystem(eqs, t, [v], [C], systems=[p, n], defaults=Dict(C => val), name=name)
end

R = 1.0
C = 1.0
V = 1.0
resistor = Resistor(:resistor, R=R)
capacitor = Capacitor(:capacitor, C=C)
source = ConstantVoltage(:source, V=V)
ground = Ground(:ground)

function connect(ps...)
    eqs = [
           0 ~ sum(p->p.i, ps) # KCL
          ]
    # KVL
    for i in 1:length(ps)-1
        push!(eqs, ps[i].v ~ ps[i+1].v)
    end

    return eqs
end
rc_eqs = [
          connect(source.p, resistor.p)
          connect(resistor.n, capacitor.p)
          connect(capacitor.n, source.n, ground.g)
         ]

rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground], name=:rc)
sys = structural_simplify(rc_model)
@test !isempty(ModelingToolkit.defaults(sys))
u0 = [
      capacitor.v => 0.0
      capacitor.p.i => 0.0
      resistor.v => 0.0
     ]
prob = ODEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Rodas4())

@test sol[resistor.p.i] == sol[capacitor.p.i]
@test sol[resistor.n.i] == -sol[capacitor.p.i]
@test sol[capacitor.n.i] == -sol[capacitor.p.i]
@test iszero(sol[ground.g.i])
@test iszero(sol[ground.g.v])
@test sol[resistor.v] == sol[source.p.v] - sol[capacitor.p.v]

prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())

@test sol[resistor.p.i] == sol[capacitor.p.i]
@test sol[resistor.n.i] == -sol[capacitor.p.i]
@test sol[capacitor.n.i] == -sol[capacitor.p.i]
@test iszero(sol[ground.g.i])
@test iszero(sol[ground.g.v])
@test sol[resistor.v] == sol[source.p.v] - sol[capacitor.p.v]
#using Plots
#plot(sol)


text:
const TJUMP = big(10)^20

"""
    tseed!([rngs=Random.THREAD_RNGs], seed=Random.make_seed(); jump=big(10)^20)

Produce approximately statistically independent rngs by seeding `rngs[1]` with `seed` and
advancing each the state of the remaining rngs by a multiple of `jump`.
"""
@inline function tseed!(seed = Random.make_seed(); jump = TJUMP)
    tseed!(Random.THREAD_RNGs, seed, jump = jump)
end

@noinline function tseed!(
    rngs::Vector{MersenneTwister},
    seed::Union{Integer,Vector{UInt32}} = Random.make_seed();
    jump::Integer = TJUMP,
)
    nt = Threads.nthreads()
    Threads.threadid() == 1 || argerror("Can only call `tseed!` from main thread")
    length(rngs) == nt || argerror("length(rngs) != Threads.nthreads()")

    rngs[1] = MT = MersenneTwister(seed)
    Threads.@threads for _ = 1:nt
        tid = Threads.threadid()
        tid > 1 && (rngs[tid] = randjump(MT, jump * (tid - 1)))
    end
    return rngs
end

getrange(N::Int, np::Int, i::Int = Threads.threadid()) = splitrange(1:N, np, i)
function getrange(range::UnitRange{Int}, np::Int, i::Int = Threads.threadid())
    np > 0 || argerror("np must be > 0")
    0 < i <= np || argerror("i must be in range (0, np]")

    N = length(range)
    start = first(range)
    each, extras = divrem(N, np)
    if each == 0
        from = start + i - 1
        return from:from
    else
        from = (i - 1) * each + min(extras, i - 1) + start
        to = from + each - 1 + ifelse(i <= extras, 1, 0)
        return from:to
    end
end

splitrange(N::Int, np::Int) = splitrange(1:N, np)
function splitrange(range::UnitRange{Int}, np::Int)
    d, r = divrem(length(range), np)
    np = d == 0 ? r : np
    UnitRange{Int}[getrange(range, np, i) for i = 1:np]
end

"""
    $(SIGNATURES)

Get the current number of threads currently used by the BLAS library.
"""
@inline function nblasthreads()
    ccall((BLAS.@blasfunc(openblas_get_num_threads), Base.libblas_name), LinearAlgebra.BlasInt, ())
end

"""
    @with_blasthreads n expression

Temporarily set the number of threads the BLAS library should use to `n` within the `expression`.

See also: [`nblasthreads`](@ref).

# Examples
```jldoctest
julia> nblasthreads()
6

julia> @with_blasthreads 10 println("Using ", nblasthreads(), " BLAS threads")
Using 10 BLAS threads

julia> nblasthreads()
6
```
"""
macro with_blasthreads(nthreads, expr)
    quote
        local n = nblasthreads()
        BLAS.set_num_threads($(esc(nthreads)))
        $(esc(expr))
        BLAS.set_num_threads(n)
    end
end


text:
config = MOI.Test.Config()
optimize_bridged!(mock) = MOIU.mock_optimize!(mock, [[128.0, 0.0, -256.0, 0.0, 160.0, 0.0, -32.0, 0.0, 1.0]; zeros(90)])
optimize_cached!(mock) = MOIU.mock_optimize!(mock, [zeros(90); [128.0, 0.0, -256.0, 0.0, 160.0, 0.0, -32.0, 0.0, 1.0]])
for mock in [bridged_mock(optimize_bridged!), cached_mock(optimize_cached!)]
    Tests.chebyshev_test(mock, config)
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
module PTBLM

function traindata(dir::String)
    xs = readdata(dir, "ptb.train.txt")
    ys = makeys(xs)
    xs, ys
end

function testdata(dir::String)
    xs = readdata(dir, "ptb.test.txt")
    ys = makeys(xs)
    xs, ys
end

function readdata(dir::String, filename::String)
    mkpath(dir)
    url = "https://raw.githubusercontent.com/tomsercu/lstm/master/data"
    path = joinpath(dir, filename)
    isfile(path) || download("$(url)/$(filename)", path)
    lines = open(readlines, path)
    map(l -> Vector{String}(split(chomp(l))), lines)
end

function makeys(xs::Vector{Vector{String}})
    map(xs) do x
        y = copy(x)
        shift!(y)
        push!(y, "<eos>")
    end
end

end


text:
# Read in the data and load packages

using Wavelets
using DataFrames
using Gadfly

# set_default_plot_size(25cm,25cm)

path = Pkg.dir("Wavelets")
data = readcsv(joinpath(path,"tests","data","testdata.csv"))

# Let us first plot the time series and only concentrate on the first one
(N,k) = size(data)
plot(x=[1:N], y=data[:,1], Geom.line)

# Now do the wavelet decomposition using MODWT on four levels using LA8 wavelet
levels = 4
decomposition = modwt(data[:,1], "la8", 4, "periodic")

# Look what do we have inside the dcompostiion object, should be quite strightforward to guess what means what
names(decomposition)

# Now plot the various wavelet coefficients with the last scaling coefficient and the original time series

# Still dont know how to reorder the things in the picture, if anyone has a clue, please help!!!
df = convert(DataFrame, [reshape(decomposition.W[:,:,1], decomposition.L, decomposition.level) decomposition.V[:,decomposition.level,1] decomposition.original[:,1] [1:decomposition.L]])
names!(df, [[symbol("W$i") for i in 1:decomposition.level], symbol("V$levels"), symbol("Originalseries"), symbol("Time")])
plotdata = stack(df, [1:(decomposition.level+2)], [decomposition.level+3])
plot(plotdata, ygroup="variable", y="value", x="Time", Geom.subplot_grid(Geom.line))

# You can do the similar thing with MRA just use the MRA method
(D, S) = mra(decomposition)
df = convert(DataFrame, [reshape(D[:,:,1], decomposition.L, decomposition.level) S[:,decomposition.level,1] decomposition.original[:,1] [1:decomposition.L]])
names!(df, [[symbol("D$i") for i in 1:decomposition.level], symbol("S$levels"), symbol("Originalseries"), symbol("Time")])
plotdata = stack(df, [1:(decomposition.level+2)], [decomposition.level+3])
plot(plotdata, ygroup="variable", y="value", x="Time", Geom.subplot_grid(Geom.line))


text:
using RegressionAndOtherStories

expend = CSV.read(ros_data("HealthExpenditure", "healthexpenditure.csv"),
    DataFrame)

life = CSV.read(ros_data("HealthExpenditure", "lifeexpectancy.csv"),
    DataFrame)

function fmax(a)
    maximum(filter(x ->!ismissing(x) .&& typeof(x) <: Number, a))
end

country = Vector(expend[4:33, 1])
country[5] = "Czech"
country[20] = "N.Zealand"
country[29] = "UK"
country[30] = "USA"

spending = [fmax(x) for x in [Vector(expend[i, 2:end]) for i in 4:33]]
lifespan = [fmax(x) for x in [Vector(life[i, 2:end]) for i in 4:33]]

df = DataFrame(
        country = country,
        spending = spending,
        lifespan = lifespan
    )

df |> display

CSV.write(ros_datadir("HealthExpenditure", "healthdata.csv"), df)


text:
using Base.Test

push!(LOAD_PATH, "./jl")
using RBT
include("../citydb.jl")

setdb(citydb)

function RBT.query(state, expr; params...)
    params_signature = ([(name, typeof(param)) for (name, param) in params]...)
    q = RBT.prepare(state, expr; params_signature...)
    println(repr(q))
    println(repr(q.pipe))
    r = RBT.execute(q; params...)
    #display(r)
    #println()
    return r
end

@query(6*(3+4))
@query(department)
@query(department.name)
@query(department.employee)
@query(department.employee.name)
@query(employee)
@query(employee.name)
@query(employee.department)
@query(employee.department.name)

@query(count(department))
@query(count(employee))
@query(department.count(employee))
@query(count(department.employee))
@query(max(employee.salary))
@query(max(department.count(employee)))
@query(sum(employee.salary))

@query(department:select(name,count(employee)))
@query(department:select(name,count(employee),max(employee.salary)))
@query(employee:select(name,surname,position,department))

@query(record(count(employee), max(employee.salary), employee))
@query(department.record(name, count(employee), max(employee.salary), employee))

@query(employee:filter(salary>200000))
@query(employee:select(name,surname,position):filter(employee.salary>200000))
@query(department:filter(count(employee)>1000).name)
@query(department:select(name,count(employee:filter(salary>100000))))

@query(department:first)
@query(department:select(name,size => count(employee)):sort(size:desc):first)
@query(department:select(name,size => count(employee)):first(size))
@query(department:last)
@query(department:select(name,size => count(employee)):sort(size):last)
@query(department:select(name,size => count(employee)):last(size:desc))
@query(department:take(5))
@query(department:skip(10):take(5))
@query(department:skip(-5))
@query(department:take(count(department)/2))
@query(department:reverse)

@query(department:select(id,name))
@query(department:get(5))
@query(department:get(-1))
@query(department[5]:select(id,name,count(employee)))
@query(department[5].employee)

@query(department.name:sort)
@query(department:desc:sort:select(id,name))
@query(employee:sort(salary))
@query(employee:sort(salary:desc))
@query(employee:sort(salary:desc,surname:asc,name:asc))

@query(department:select(name,count(employee):as(size)))
@query(department:define(size => count(employee)):filter(size>1000):select(name,size))

@query(employee:filter(salary>managed_by.salary):select(name,surname,position,managed_by))
@query(department[26].employee:select(name,surname,position,connect(managed_by),count(connect(manages))))
@query(department[26].employee:sort_connect(managed_by):select(depth(managed_by),name,surname,position))

@query(employee.position:unique)
@query(employee.position:desc:unique)
@query(department:select(name,count(unique(employee.position)),count(employee)))

@query(employee:group(position))
@query(employee:group(position):select(position,count(employee)))
@query(employee:group(position):select(position, size => count(employee)):sort(size:desc))
@query(
    employee:group(position)
    :define(department => unique(employee.department))
    :filter(count(department)>=5)
    :sort(count(department):desc)
    :select(position, department))
@query(employee:group(name):select(name, size => count(employee)):sort(size:desc))
@query(
    employee:group(name)
    :filter(count(employee)>=10)
    :select(name, max_salary => max(employee.salary))
    :sort(max_salary:desc))
@query(
    employee
    :group(department, salary_bracket => salary/10000*10000 :desc)
    :select(department, salary_bracket, salary_bracket+9999, count(employee)))
@query(
    employee
    :group_cube(department, salary_bracket => salary/10000*10000 :desc)
    :select(department, salary_bracket, salary_bracket+9999, count(employee)))
@query(
    employee
    :partition(department:take(3))
    :select(department.name, count(unique(employee.position)), count(employee), max(employee.salary)))
@query(
    employee
    :partition_cube(department:take(3))
    :select(department.name, count(unique(employee.position)), count(employee), max(employee.salary)))
@query(
    employee
    :define(salary_bracket => salary/60000*60000)
    :partition(
        department:take(3),
        name => ["ANTHONY", "BRIAN"],
        salary_bracket => range(0, 60000, max(employee.salary)))
    :select(
        dept => department.name,
        name,
        low => salary_bracket,
        high => salary_bracket+59999,
        pop_position => employee:group(position):first(count(employee)).position,
        num_pos => count(unique(employee.position)),
        num_empl => count(employee)))

@query(
    mix(a => range(2,1,10), b => range(2,1,10), c => range(2,1,10))
    :filter((a <= b) & (b <= c))
    :select(a, b, c, (a*b)*c))
@query(
    mix(department, department)
    :filter((left.id != right.id) & (left.count(employee)/10 == right.count(employee)/10))
    :select(left.name, left.count(employee), right.name, right.count(employee)))
@query(pack(a => range(1,1,5), z => range(95,1,99)))
@query(pack(department, employee))
@query(pack(size => department.count(employee), position => employee.position))
@query(pack(department, employee).department)
@query(pack(department, employee).employee)
@query(pack(department, employee):select(department.name, employee.name))
@query(employee:take(100):filter(salary>max(unlink(employee).salary)/2))
@query(
    employee[10]
    :define(namesake => link((left.id!=right.id)&(left.name==right.name), employee))
    :select(name, count(namesake), namesake))

@query(department:json)
@query(department:select(name,head => employee:first(salary)):json)
@query(department:json(name,head => employee:first(salary)))
@query(employee:dataframe)
@query(department:select(name,size => count(employee), max_salary => max(employee.salary)):dataframe)
@query(department:dataframe(name,size => count(employee), max_salary => max(employee.salary)))

@query((max(employee.salary) > 100000) & (max(employee.salary) < 300000))

@query(X*X, X=5)
@query(X*Y, X=5, Y=Nullable(4))
@query(X*(Y*Z), X=5, Y=Nullable(4), Z=[2,3,4,5])

q = RBT.@prepare(X*(Y*Z), X=Int, Y=Nullable{Int}, Z=Vector{Int})
q(X=5, Y=4, Z=[2,3,4,5])

@query(employee:filter((position==POSITION) & (name==NAME)), POSITION="POLICE OFFICER", NAME="CHARLES")
@query(department:filter(count(employee)>SIZE):select(name,count(employee)-SIZE), SIZE=1000)
@query(employee:filter(department.name in DEPTS), DEPTS=["POLICE", "FIRE"])

@query(department:select(name, previous_name => first(before).name, next_name => first(after).name))
@query(department:select(name, past_names => before.name))
@query(department:select(1+count(before), name))
@query(department:define(size => count(employee)):select(name, size, total => sum(and_before.size)))
@query(department:define(size => count(employee)):filter(size == max(and_around.size)):select(name, size))
@query(
    employee
    :filter(department.name == DEPT)
    :dataframe(
        name, surname, position, salary, salary_diff => max(and_around(position).salary)-salary),
    DEPT="TREASURER")
@query(department:select(name, employee:filter(salary==max(and_around.salary)):frame))
@query(employee:filter(salary==MAX_SALARY):given(MAX_SALARY => max(employee.salary)))
@query(department:select(name, employee:filter(salary==MAX_SALARY):given(MAX_SALARY => max(employee.salary))))

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
using StaticArrays
using Test

using CLIMA
using CLIMA.Atmos
using CLIMA.PlanetParameters
using CLIMA.MoistThermodynamics
using CLIMA.PlanetParameters
using CLIMA.VariableTemplates

function init_test!(bl, state, aux, (x,y,z), t)
    FT = eltype(state)

    z = FT(z)

    # These constants are those used by Stevens et al. (2005)
    qref       = FT(9.0e-3)
    q_pt_sfc   = PhasePartition(qref)
    Rm_sfc     = FT(gas_constant_air(q_pt_sfc))
    T_sfc      = FT(290.4)
    P_sfc      = FT(MSLP)

    # Specify moisture profiles
    q_liq      = FT(0)
    q_ice      = FT(0)

    θ_liq  = FT(289.0)
    q_tot  = qref

    ugeo = FT(7)
    vgeo = FT(-5.5)
    u, v, w = ugeo, vgeo, FT(0)

    # Pressure
    H     = Rm_sfc * T_sfc / grav
    p     = P_sfc * exp(-z / H)

    # Density, Temperature
    ts    = LiquidIcePotTempSHumEquil_given_pressure(θ_liq, p, q_tot)
    ρ     = air_density(ts)

    e_kin = FT(1/2) * FT((u^2 + v^2 + w^2))
    e_pot = grav * z
    E     = ρ * total_energy(e_kin, e_pot, ts)

    state.ρ               = ρ
    state.ρu              = SVector(ρ*u, ρ*v, ρ*w)
    state.ρe              = E
    state.moisture.ρq_tot = ρ * q_tot

    return nothing
end

function main()
    CLIMA.init()

    FT = Float64

    # DG polynomial order
    N = 4

    # Domain resolution and size
    Δh = FT(40)
    Δv = FT(40)
    resolution = (Δh, Δh, Δv)

    xmax = 320
    ymax = 320
    zmax = 400

    t0 = FT(0)
    timeend = FT(10)

    driver_config = CLIMA.Atmos_LES_Configuration("Driver test", N, resolution,
                                                  xmax, ymax, zmax, init_test!)
    solver_config = CLIMA.setup_solver(t0, timeend, driver_config)

    result = CLIMA.invoke!(solver_config)
end

main()

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
using NeighbourLists
using JuLIP

# at = set_pbc!(bulk(:Si, cubic=true) * 3, true)
# length(at)
# cutoff = 3.1 * rnn(:Si)
# nlist = CellList(positions(at), cutoff, cell(at), pbc(at))
#
# f(r, R) = (exp(-8*(r-1)) - 2 * exp( - 4 * (r - 1))) * 1 / (1+r)
#
# out = zeros(length(at))
# mapreduce_sym!(out, f, pairs(nlist))
#
# @time mapreduce_sym!(out, f, pairs(nlist));
# @time NeighbourLists.tmapreduce_sym!(out, f, pairs(nlist));

macro symm(N, ex)
   if N isa Symbol
      N = eval(N)
   end
   @assert ex.head == :for
   @assert length(ex.args) == 2
   ex_for = ex.args[1]
   ex_body = ex.args[2]
   # iteration symbol
   i = ex_for.args[1]
   # lower and upper bound
   a0 = ex_for.args[2].args[1]
   a1 = ex_for.args[2].args[2]
   # create the for-loop without body
   loopstr = "for $(i)1 = ($a0):(($a1)-$(N-1))"
   for n = 2:N
      loopstr *= ", $i$n = $i$(n-1)+1:(($a1)-$(N-n))"
   end
   loopstr *= "\n $i = CartesianIndex($(i)1"
   for n = 2:N
      loopstr *= ", $i$n"
   end
   loopstr *= ") \n end"
   loopex = Meta.parse(loopstr)
   append!(loopex.args[2].args, ex_body.args)
   # return the expression
   esc(quote
      $loopex
   end)
end

@generated function t(b0, b1, Nval::Val{N}) where N
   quote
      @symm $N for j = b0:b1
         println(j)
      end
   end
end

t(3, 8, Val(3))

M = 4
@symm M for j = 1:6
   println(j)
end

getN(::Val{N}) where N = N::Int
getN(n::Integer) = n
getN(::Type{Val{N}}) where N = N::Int

macro mm(NN)
   @show NN
   @show getN(eval(NN))
end

@mm(4)
@mm(Val(4))

getN(Val{4})

M = 4
@mm 4
@mm M

function m1(::Val{N}) where N
   @show N
   @show typeof(N)
end

m1(Val(4))

@generated function m2(v::T)
   @mm(T)
end

function simplex_lengths(I::SVector{N, TI}, X::AbstractVector{T}) where {
                                          N, TI <: Integer, T <: AbstractFloat}
   N2 = (N*(N-1))÷2
   a = zeros(TI, N2)
   b = zeros(TI, N2)
   s = zeros(T, N2)
   n = 0
   for i = 1:N-1, j = i+1:N
      n += 1
      a[n] = I[i]
      b[n] = I[j]
      s[n] = abs(X[a[n]] - X[b[n]])
   end
   return s, a, b
end

I = SVector(3, 5, 7, 2)
X = rand(10)
simplex_lengths(I, X)

function statalloc(I::SVector{N, TI}) where {N, TI}
   N2 = (N*(N-1)) ÷ 2
   a = zero(MVector{N2, TI})
end

statalloc(I)

function simplex_lengths!(s, a, b, I::SVector{N, TI}, X) where {N, TI <: Integer}
   n = 0
   for i = 1:N-1, j = i+1:N
      n += 1
      a[n] = I[i]
      b[n] = I[j]
      s[n] = norm(X[a[n]] - X[b[n]])
   end
   return SVector(s), SVector(a), SVector(b)
end

function simplex_lengths2!(s, a, b, I::SVector{N, TI}, X) where {N, TI <: Integer}
   n = 0
   for i = 1:N-1, j = i+1:N
      n += 1
      a[n] = I[i]
      b[n] = I[j]
      s[n] = norm(X[a[n]] - X[b[n]])
   end
   return s, a, b
end

function test(slen, M)
   I = SVector(3, 50, 74, 22, 10)
   N = length(I)
   N2 = (N*(N-1)) ÷ 2
   X = rand(SVector{3, Float64}, 100)
   s = zero(MVector{N2, Float64})
   a = zero(MVector{N2, Int})
   b = zero(MVector{N2, Int})
   for m = 1:M
      slen(s, a, b, I, X)
   end
end

@btime test(simplex_lengths!, 10_000)
@btime test(simplex_lengths2!, 10_000)
# @btime simplex_lengths2!($s, $a, $b, $I, $X)


text:
" ProbNumDiffEq.jl-specific implementation of OrdinaryDiffEq.jl's `postamble!`."
function OrdinaryDiffEq.postamble!(integ::OrdinaryDiffEq.ODEIntegrator{<:AbstractEK})
    # OrdinaryDiffEq.jl-related calls:
    OrdinaryDiffEq._postamble!(integ)
    copyat_or_push!(integ.sol.k, integ.saveiter_dense, integ.k)
    pn_solution_endpoint_match_cur_integrator!(integ)

    # Calibrate the solution (if applicable)
    if isstatic(integ.cache.diffusionmodel)
        if integ.cache.diffusionmodel.calibrate
            # The estimated global_diffusion is just a scaling factor
            mle_diffusion = integ.cache.global_diffusion
            calibrate_solution!(integ, mle_diffusion)
        else
            constant_diffusion = integ.cache.default_diffusion
            set_diffusions!(integ.sol, constant_diffusion)
        end
    end

    if integ.alg.smooth
        smooth_solution!(integ)
    end

    @assert (length(integ.sol.u) == length(integ.sol.pu))

    return nothing
end

"""
    calibrate_solution!(integ, mle_diffusion)

Calibrate the solution (`integ.sol`) with the specified `mle_diffusion` by (i) setting the
values in `integ.sol.diffusions` to the `mle_diffusion` (see [`set_diffusions!`](@ref)),
(ii) rescaling all filtering estimates such that they have the correct diffusion, and (iii)
updating the solution estimates in `integ.sol.pu`.
"""
function calibrate_solution!(integ, mle_diffusion)

    # Set all solution diffusions; don't forget the initial diffusion!
    set_diffusions!(integ.sol, mle_diffusion * integ.cache.default_diffusion)

    # Rescale all filtering estimates to have the correct diffusion
    @simd ivdep for s in integ.sol.x_filt
        copy!(s.Σ, apply_diffusion(s.Σ, mle_diffusion))
    end

    # Re-write into the solution estimates
    for (pu, x) in zip(integ.sol.pu, integ.sol.x_filt)
        mul!(pu, integ.cache.SolProj, x)
    end
    # [(su[:] .= pu) for (su, pu) in zip(integ.sol.u, integ.sol.pu.μ)]
end

"""
    set_diffusions!(solution::AbstractProbODESolution, diffusion)

Set the contents of `solution.diffusions` to the provided `diffusion`, overwriting the local
diffusion estimates that are in there. Typically, `diffusion` is either a global quasi-MLE
or the specified initial diffusion value if no calibration is desired.
"""
function set_diffusions!(solution::AbstractProbODESolution, diffusion::Number)
    solution.diffusions .= diffusion
    return nothing
end
function set_diffusions!(solution::AbstractProbODESolution, diffusion::Diagonal)
    @simd ivdep for d in solution.diffusions
        copy!(d, diffusion)
    end
    return nothing
end

"""
    smooth_solution!(integ)

Smoothes the solution saved in `integ.sol`, filling `integ.sol.x_smooth` and updating the
values saved in `integ.sol.pu` and `integ.sol.u`. This function handles the iteration and
preconditioning. The actual smoothing step happens in [`smooth!`](@ref).
"""
function smooth_solution!(integ)
    integ.sol.x_smooth = copy(integ.sol.x_filt)

    @unpack A, Q = integ.cache
    @unpack x_smooth, t, diffusions = integ.sol
    @unpack x_tmp, x_tmp2 = integ.cache
    x = x_smooth

    for i in length(x)-1:-1:1
        dt = t[i+1] - t[i]
        if iszero(dt)
            copy!(x[i], x[i+1])
            continue
        end

        make_preconditioners!(integ.cache, dt)
        P, PI = integ.cache.P, integ.cache.PI

        mul!(x_tmp, P, x[i])
        mul!(x_tmp2, P, x[i+1])
        smooth!(x_tmp, x_tmp2, A, Q, integ, diffusions[i])
        mul!(x[i], PI, x_tmp)

        # Save the smoothed state into the solution
        mul!(integ.sol.pu[i], integ.cache.SolProj, x[i])
        integ.sol.u[i][:] .= integ.sol.pu[i].μ
    end
    integ.sol.interp = set_smooth(integ.sol.interp)
    return nothing
end

"Inspired by `OrdinaryDiffEq.solution_match_cur_integrator!`"
function pn_solution_endpoint_match_cur_integrator!(integ)
    if integ.opts.save_end
        if integ.alg.smooth
            OrdinaryDiffEq.copyat_or_push!(
                integ.sol.x_filt,
                integ.saveiter_dense,
                integ.cache.x,
            )
        end

        OrdinaryDiffEq.copyat_or_push!(
            integ.sol.pu,
            integ.saveiter,
            mul!(integ.cache.pu_tmp, integ.cache.SolProj, integ.cache.x),
        )
    end
end

"Extends `OrdinaryDiffEq._savevalues!` to save ProbNumDiffEq.jl-specific things."
function DiffEqBase.savevalues!(
    integ::OrdinaryDiffEq.ODEIntegrator{<:AbstractEK},
    force_save=false,
    reduce_size=true,
)

    # Do whatever OrdinaryDiffEq would do
    out = OrdinaryDiffEq._savevalues!(integ, force_save, reduce_size)

    # Save our custom stuff that we need for the posterior
    # TODO If we don't want dense output, we might not want to save these!
    # It's not completely clear how to specify that though; They are also needed for sampling.
    if integ.alg.smooth
        OrdinaryDiffEq.copyat_or_push!(integ.sol.x_filt, integ.saveiter, integ.cache.x)
    end
    OrdinaryDiffEq.copyat_or_push!(
        integ.sol.diffusions,
        integ.saveiter,
        integ.cache.local_diffusion,
    )
    if integ.opts.save_everystep
        OrdinaryDiffEq.copyat_or_push!(
            integ.sol.pu,
            integ.saveiter,
            mul!(integ.cache.pu_tmp, integ.cache.SolProj, integ.cache.x),
        )
    end

    return out
end


text:
using HomotopyContinuation.ModelKit

@testset "ModelKit" begin
    include("model_kit/symbolic_test.jl")
    include("model_kit/instructions_test.jl")
    include("model_kit/slp_test.jl")
end


text:
using ISA, Plots

include(raw"colorMaps/cubeYF.jl")
include(raw"colorMaps/viridis.jl")

path = raw"C:\Users\hemad\Desktop\NMSU Cources\EE590\images"

fs = 5000#16_000
t = 0.0:1/fs:3.0
𝐶₀ = AMFMtriplet(t->0.2+0.8cos(11t), t->200.0, 0.0)
𝐶₁ = AMFMtriplet(t->exp(-abs(t/3)), t->100t, 0.1)
𝐶₂ = AMFMtriplet(t->𝒩ᵤ(t; μ=1.5, σ=1.0), t->150 + 125sin(5t), π)
𝐶₃ = AMFMtriplet(t->u(t-1.725)-u(t-2.475), t->50, 0.0)
𝑆1 = compSet([𝐶₀,𝐶₁])
𝑆2 = compSet([𝐶₀,𝐶₁,𝐶₂,𝐶₃])

default(margins=1Plots.mm,bottom_margin =5Plots.mm,
        left_margin=13Plots.mm,dpi=300,zlims=(-2,2))

# 1_ 3D 𝑆1
# 2_ TF
# 3_ 3D
# 4_ TR
# 5_ 3D
# 6_ TURN ON PROJECTION
# 7_ 𝑆2
# 8_ TR
# 9_ 3D
# 10_ TF
# 11_ 3D

# 1_ 3D 𝑆1
plot(𝑆1,camera=(20,80.0),timeaxis=t)
CNT = 0
for i in 0:25
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 26

# 2_ TF

# if 1==0
#     for i ∈ 0:20
#         θ1 = 20.0 - i
#         rotate = i*90.0/20
#         plot(𝑆1,camera=(θ1,80.0),yrotation=rotate,timeaxis=t)
#         savefig(path*raw"\plot"*lpad(i+CNT,5,"0")*".png")
#     end;CNT = CNT + 21
#
#     for i in 0:20
#         savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#     end;CNT = CNT + 21
#
#     for i ∈ 0:10
#         θ2 = 80.0 + i
#         plot(𝑆1,camera=(0.0,θ2),yrotation=90.0,timeaxis=t)
#         # (i<5) ? plot(𝑆,camera=(0.0,θ2),yrotation=90.0) : plot(𝑆,camera=(0.0,90),yrotation=90.0,zguidefont = viridis()[1],showaxis=:xy,z_foreground_color_text=viridis()[1])
#         savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#     end;CNT = CNT + 11
# end
for i ∈ 0:10
    θ1 = 20.0 - i
    θ2 = 80.0 + i
    rotate = i*45.0/10
    plot(𝑆1,camera=(θ1,θ2),yrotation=rotate,timeaxis=t)
    # (i<5) ? plot(𝑆,camera=(0.0,θ2),yrotation=90.0) : plot(𝑆,camera=(0.0,90),yrotation=90.0,zguidefont = viridis()[1],showaxis=:xy,z_foreground_color_text=viridis()[1])
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 11

for i ∈ 0:10
    θ1 = 10.0 - i
    rotate = 45 + i*45.0/10
    plot(𝑆1,camera=(θ1,90.0),yrotation=rotate,timeaxis=t)
    savefig(path*raw"\plot"*lpad(i+CNT,5,"0")*".png")
end;CNT = CNT + 11

plot(𝑆1,camera=(0.0,90),yrotation=90.0,zguidefont = viridis()[1],
     showaxis=:xy,z_foreground_color_text=viridis()[1],
     timeaxis=t)

for i in 0:35
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 36

# 3_ 3D
for i ∈ 0:10
    θ1 = i
    θ2 = 90.0 - i
    rotate = 90 - i*45.0/10
    plot(𝑆1,camera=(θ1,θ2),yrotation=rotate,timeaxis=t)
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 11

for i ∈ 0:10
    θ1 = 10.0 + i
    rotate = 45 - i*45.0/10
    plot(𝑆1,camera=(θ1,80.0),yrotation=rotate,timeaxis=t)
    savefig(path*raw"\plot"*lpad(i+CNT,5,"0")*".png")
end;CNT = CNT + 11

plot(𝑆1,camera=(20,80.0),timeaxis=t)
for i in 0:35
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 36

# 4_ TR
for i in 0:20
    θ1 = 20.0 - i
    θ2 = 80.0 - i
    plot(𝑆1,camera=(θ1,θ2),timeaxis=t)
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 21

for i in 0:60
    θ2 = 60.0 - i
    (i<45) ? plot(𝑆1,camera=(0.0,θ2),timeaxis=t) : plot(
    𝑆1,camera=(0.0,θ2),timeaxis=t,yguidefont = viridis()[1],
    showaxis=:xz,y_foreground_color_text=viridis()[1])
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 61

for i in 0:35
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end; CNT = CNT + 36

# 5_ 3D
for i in 0:20
    plot(𝑆1,camera=(i,i),timeaxis=t,yguidefont = viridis()[1],
    showaxis=:xz,y_foreground_color_text=viridis()[1])
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 21

for i in 0:60
    plot(𝑆1,camera=(20.0,20+i),timeaxis=t)
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 61

for i in 1:3
    plot(𝑆1,camera=(20,80.0),timeaxis=t)
    for i in 0:25
        savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
    end; CNT = CNT + 26

    # 6_ TURN ON PROJECTION
    plot(𝑆1,camera=(20,80.0),timeaxis=t,realProj=true)
    for i in 0:25
        savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
    end; CNT = CNT + 26

    # 7_ 𝑆2
    plot(𝑆2,camera=(20,80.0),timeaxis=t,realProj=true)
    for i in 0:35
        savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
    end; CNT = CNT + 36
end

# 8_ TR
for i in 0:20
    θ1 = 20.0 - i
    θ2 = 80.0 - i
    plot(𝑆2,camera=(θ1,θ2),timeaxis=t,realProj=true)
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 21

for i in 0:60
    θ2 = 60.0 - i
    (i<45) ? plot(𝑆2,camera=(0.0,θ2),timeaxis=t, realProj=true) : plot(
    𝑆2,camera=(0.0,θ2),timeaxis=t,yguidefont = viridis()[1],
    showaxis=:xz,y_foreground_color_text=viridis()[1], realProj=true)
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 61

for i in 0:35
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end; CNT = CNT + 36

# 9_ 3D
for i in 0:20
    plot(𝑆2,camera=(i,i),timeaxis=t,yguidefont = viridis()[1],
    showaxis=:xz,y_foreground_color_text=viridis()[1],
    realProj=true)
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 21

for i in 0:60
    plot(𝑆2,camera=(20.0,20+i),timeaxis=t,realProj=true)
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 61

for i in 0:35
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end; CNT = CNT + 36

#-------------------------------------
# 10_ TF
for i ∈ 0:20
    θ1 = 20.0 - i
    rotate = i*90.0/20
    plot(𝑆2,camera=(θ1,80.0),yrotation=rotate,timeaxis=t,realProj=true)
    savefig(path*raw"\plot"*lpad(i+CNT,5,"0")*".png")
end;CNT = CNT + 21

for i in 0:35
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 36

for i ∈ 0:10
    θ2 = 80.0 + i
    plot(𝑆2,camera=(0.0,θ2),yrotation=90.0,timeaxis=t,realProj=true)
    # (i<5) ? plot(𝑆,camera=(0.0,θ2),yrotation=90.0) : plot(𝑆,camera=(0.0,90),yrotation=90.0,zguidefont = viridis()[1],showaxis=:xy,z_foreground_color_text=viridis()[1])
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 11

plot(𝑆2,camera=(0.0,90),yrotation=90.0,zguidefont = viridis()[1],
     showaxis=:xy,z_foreground_color_text=viridis()[1],
     timeaxis=t,realProj=true)

for i in 0:35
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 36

# 11_ 3D
for i ∈ 0:10
    θ1 = i
    θ2 = 90.0 - i
    rotate = 90 - i*45.0/10
    plot(𝑆2,camera=(θ1,θ2),yrotation=rotate,timeaxis=t,realProj=true)
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 11

for i ∈ 0:10
    θ1 = 10.0 + i
    rotate = 45 - i*45.0/10
    plot(𝑆2,camera=(θ1,80.0),yrotation=rotate,timeaxis=t,realProj=true)
    savefig(path*raw"\plot"*lpad(i+CNT,5,"0")*".png")
end;CNT = CNT + 11

plot(𝑆2,camera=(20,80.0),timeaxis=t,realProj=true)
for i in 0:35
    savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
end;CNT = CNT + 36

# if 1==0
#     for i in 0:140 #TR
#         if i in 0:20
#             plot(𝑆,camera=(0,0),
#             yrotation=90.0,yguidefont = viridis()[1],
#             showaxis=:xz,y_foreground_color_text=viridis()[1])
#             savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#         elseif i in 21:40
#             plot(𝑆,camera=(0,0),realProj=true,
#             yrotation=90.0,yguidefont = viridis()[1],
#             showaxis=:xz,y_foreground_color_text=viridis()[1])
#             savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#         elseif i in 41:60
#             plot(𝑆,camera=(0,0),
#             yrotation=90.0,yguidefont = viridis()[1],
#             showaxis=:xz,y_foreground_color_text=viridis()[1])
#             savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#         elseif i in 61:80
#             plot(𝑆,camera=(0,0),realProj=true,
#             yrotation=90.0,yguidefont = viridis()[1],
#             showaxis=:xz,y_foreground_color_text=viridis()[1])
#             savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#         elseif i in 81:100
#             plot(𝑆,camera=(0,0),
#             yrotation=90.0,yguidefont = viridis()[1],
#             showaxis=:xz,y_foreground_color_text=viridis()[1])
#             savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#         elseif i in 101:120
#             plot(𝑆,camera=(0,0),realProj=true,
#             yrotation=90.0,yguidefont = viridis()[1],
#             showaxis=:xz,y_foreground_color_text=viridis()[1])
#             savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#         elseif i in 121:140
#             plot(𝑆,camera=(0,0),
#             yrotation=90.0,yguidefont = viridis()[1],
#             showaxis=:xz,y_foreground_color_text=viridis()[1])
#             savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#         end
#     end;CNT = CNT + 141
#
#     for i ∈ 0:90
#         marg = i*7/90
#         default(margins=(marg)Plots.mm)
#         (i<75) ? plot(𝑆,camera=(i,0.0)) : plot(𝑆,
#                     camera=(i,0.0),xguidefont = viridis()[1],
#                     showaxis=:yz,x_foreground_color_text=viridis()[1])
#                     savefig(path*raw"\plot"*lpad((CNT + i),5,"0")*".png")
#                 end; CNT = CNT + 91
#
#                 for i in 0:25 #FR
#                     savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#                 end;CNT = CNT + 26
#
#
#                 for i ∈ 0:70
#                     marg = i*6/70
#                     default(margins=(7-marg)Plots.mm)
#                     plot(𝑆,camera=(90-i,i))
#                     savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#                 end;
#                 CNT = CNT + 71
#
#                 for i ∈ 0:10
#                     plot(𝑆,camera=(20,70+i))
#                     savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#                 end;CNT = CNT + 11
#
#                 # plot(𝑆,camera=(20.0,80))
#                 # for i in 0:25
#                 #     savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#                 # end;CNT = CNT + 25
#
#                 for i in 0:100
#                     if i in 0:20
#                         plot(𝑆,camera=(20.0,80))
#                         savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#                     elseif i in 21:40
#                         plot(𝑆,camera=(20.0,80),realProj=true)
#                         savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#                     elseif i in 41:60
#                         plot(𝑆,camera=(20.0,80))
#                         savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#                     elseif i in 61:80
#                         plot(𝑆,camera=(20.0,80),realProj=true)
#                         savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#                     elseif i in 81:100
#                         plot(𝑆,camera=(20.0,80))
#                         savefig(path*raw"\plot"*lpad((i + CNT),5,"0")*".png")
#                     end
#                 end;CNT = CNT + 101
#
# end

#ffmpeg -r 30 -f image2 -s 1920x1080 -i plot%05d.png -vcodec libx264 -crf 25  -pix_fmt yuv420p test.mp4

#issue: the following will plot the same plot
# plot(𝑆,camera=(20.0,80))
# plot(𝑆,camera=(20.5,80))

# and the following too:
# plot(𝑆,camera=(20.6,80))
# plot(𝑆,camera=(21.0,80))


text:
"""
    calculateMassMu!(m, g, m_j, mu, mu_sum, m_sum)

Given the gravitational constant and the mass of each body, calculates sums of masses, standard gravitational parameters of the bodies and jacobian masses

```math
mu_{sum_{i}} = \\sum_{j = 1}^{i} m_{i} \\hspace{0.5cm} mu = G \\cdot m_{i} \\hspace{0.5cm} mu_{sum} = G \\cdot m_{sum_{i}} \\hspace{0.5cm} m' = m_{i} \\frac{m_{sum_{i-1}}}{m_{sum_{i}}}
```

# Args

* `m`: array of N elements, each containing the mass of each body.
* `g`: gravitational constant.
* `m_j`: array of N elements. Each element will contain the jacobi mass ``m'`` of a body, as output.
* `mu`: array of N elements. Each element will contain the standard gravitational parameter of a body, as output.
* `mu_sum`: array of N elements, where the standard gravitational parameter of the bodies from 1 to i will be saved saved in the i-th element, as output.
* `m_sum`: array of N elements, where the element i contains the mass of the bodies from 1 to i, as output.
"""
function calculateMassMu!(m, g, m_j, mu, mu_sum, m_sum)

    m_sum[1] = m[1]
    mu_sum[1] = m[1] * g
    mu[1] = m[1] * g

    for i in 2:length(m)
        mu[i] = m[i] * g
        m_sum[i] = m_sum[i-1] + m[i]
        mu_sum[i] = m_sum[i] * g
        m_j[i] = m[i] * m_sum[i-1] / m_sum[i]
    end
    m_j[1] = m_sum[length(m)]

end

max_stars_repo_name: salanueva/KeplerFlow

max_issues_repo_name: salanueva/KeplerFlow

max_forks_repo_name: salanueva/KeplerFlow


text:
global test_suite_polyhedra

using LazySets: _isbounded_stiemke, _isbounded_unit_dimensions

for N in [Float64, Rational{Int}, Float32]
    # random polyhedron
    rand(HPolyhedron)

    p_univ = HPolyhedron{N}()

    # constructor from matrix and vector
    A = [N(1) N(2); N(-1) N(1)]
    b = [N(1), N(2)]
    p = HPolyhedron(A, b)
    c = p.constraints
    @test c isa Vector{LinearConstraint{N, Vector{N}}}
    @test c[1].a == N[1, 2] && c[1].b == N(1)
    @test c[2].a == N[-1, 1] && c[2].b == N(2)

    # convert back to matrix and vector
    A2, b2 = tosimplehrep(p)
    @test A == A2 && b == b2

    # 2D HPolyhedron
    p = HPolyhedron{N}()
    c1 = LinearConstraint(N[2, 2], N(12))
    c2 = LinearConstraint(N[-3, 3], N(6))
    c3 = LinearConstraint(N[-1, -1], N(0))
    c4 = LinearConstraint(N[2, -4], N(0))
    addconstraint!(p, c3)
    addconstraint!(p, c1)
    addconstraint!(p, c4)
    addconstraint!(p, c2)

    # tosimplehrep for other polytopic types
    A, b = tosimplehrep(BallInf(N[0], N(2)))
    @test (A == hcat(N[1; -1]) || A == hcat(N[-1; 1])) && b == N[2, 2]
    # tosimplehrep from list of constraints
    A, b = tosimplehrep([HalfSpace(N[1], N(2)), HalfSpace(N[-1], N(2))])
    @test (A == hcat(N[1; -1]) || A == hcat(N[-1; 1])) && b == N[2, 2]

    # support vector of polyhedron with no constraints
    @test σ(N[1], p_univ) == N[Inf]

    # universality
    @test !isuniversal(p)
    res, w = isuniversal(p, true)
    @test !res && w ∉ p
    @test isuniversal(p_univ) && isuniversal(p_univ, true) == (true, N[])

    # membership
    @test N[5 / 4, 7 / 4] ∈ p && N[4, 1] ∉ p

    # constrained dimensions
    @test constrained_dimensions(p) == [1, 2]
    @test constrained_dimensions(
        HPolyhedron([LinearConstraint(N[1, 0], N(1))])) == [1]
    @test constrained_dimensions(HPolyhedron()) == Int[]

    # concrete linear map with invertible matrix
    linear_map(N[2 3; 1 2], p)

    # translation
    p2 = translate(p, N[1, 2])
    @test p2 isa HPolyhedron{N} && ispermutation(constraints_list(p2),
        [HalfSpace(N[2, 2], N(18)), HalfSpace(N[-3, 3], N(9)),
         HalfSpace(N[-1, -1], N(-3)), HalfSpace(N[2, -4], N(-6))])

    # constraints iterator
    @test ispermutation(collect(constraints(p)), constraints_list(p))

    # vertices iterator
    @test ispermutation(collect(vertices(p)), vertices_list(p))

    if test_suite_polyhedra
        # conversion to and from Polyhedra's VRep data structure
        cl = constraints_list(HPolyhedron(polyhedron(p)))
        @test length(p.constraints) == length(cl)

        # convert hyperrectangle to a HPolyhedron
        H = Hyperrectangle(N[1, 1], N[2, 2])
        P = convert(HPolyhedron, H)
        @test length(vertices_list(P)) == 4

        # checking for emptiness
        P = HPolyhedron([LinearConstraint(N[1, 0], N(0))])    # x <= 0
        @test !isempty(P)

        # concrete linear map with non-invertible matrix
        if N == Float64
            @test linear_map(N[2 3; 0 0], P) isa HPolyhedron
        end
    end

    if !test_suite_polyhedra
        # concrete linear map of a bounded polyhedron by a non-invertible matrix
        # throws an assertion error, since tovrep(HPolytope(...)) is required
        H = Hyperrectangle(N[1, 1], N[2, 2])
        P = convert(HPolyhedron, H)
        if N != Float32
            # conversion to vrep with Float32 fails
            @test_throws ArgumentError linear_map(N[2 3; 0 0], P, algorithm="vrep")
        end

        if N != Rational{Int} # in floating-point we can use elimination
            lm = linear_map(N[2 3; 0 0], P, algorithm="elimination")
            @test lm isa HPolyhedron{Float64}

            B = N[4e8 2; 0 1]
            P = CartesianProduct(BallInf(N[0.01], N(0.08)), Singleton(N[1.0]))
            lm = linear_map(B, P)
            @test lm isa HPolytope{Float64}
        end
    end
end

# default Float64 constructors
@test HPolyhedron() isa HPolyhedron{Float64}

# tests that only work with Float64 and Float32
for N in [Float64, Float32]
    # normalization
    p1 = HPolyhedron([HalfSpace(N[1e5], N(3e5)), HalfSpace(N[-2e5], N(4e5))])
    p2 = normalize(p1)
    for hs in constraints_list(p2)
        @test norm(hs.a) == N(1)
    end

    # 2D HPolyhedron
    p = HPolyhedron{N}()
    c1 = LinearConstraint(N[2, 2], N(12))
    c2 = LinearConstraint(N[-3, 3], N(6))
    c3 = LinearConstraint(N[-1, -1], N(0))
    c4 = LinearConstraint(N[2, -4], N(0))
    addconstraint!(p, c3)
    addconstraint!(p, c1)
    addconstraint!(p, c4)
    addconstraint!(p, c2)

    p_univ = HPolyhedron{N}()

    # boundedness
    @test !isbounded(p_univ)
    @test isbounded(p)
    @test !isbounded(HPolyhedron([HalfSpace(N[1, 0], N(1))]))

    @test _isbounded_stiemke(p_univ)
    @test _isbounded_stiemke(p)
    @test !_isbounded_stiemke(HPolyhedron([HalfSpace(N[1, 0], N(1))]))

    @test _isbounded_unit_dimensions(p_univ)
    @test _isbounded_unit_dimensions(p)
    @test !_isbounded_unit_dimensions(HPolyhedron([HalfSpace(N[1, 0], N(1))]))
end

# Polyhedra tests that only work with Float64
for N in [Float64]
    p = HPolyhedron{N}()
    c1 = LinearConstraint(N[2, 2], N(12))
    c2 = LinearConstraint(N[-3, 3], N(6))
    c3 = LinearConstraint(N[-1, -1], N(0))
    c4 = LinearConstraint(N[2, -4], N(0))
    addconstraint!(p, c3)
    addconstraint!(p, c1)
    addconstraint!(p, c4)
    addconstraint!(p, c2)

    # support vector
    d = N[1, 0]
    @test σ(d, p) == N[4, 2]
    d = N[0, 1]
    @test σ(d, p) == N[2, 4]
    d = N[-1, 0]
    @test σ(d, p) == N[-1, 1]
    d = N[0, -1]
    @test σ(d, p) == N[0, 0]

    # membership
    @test [Inf, Inf] ∉ p

    # an_element
    P = HPolyhedron([HalfSpace(N[3//50, -4//10], N(1)),
                     HalfSpace(N[-1//50, 1//10], N(-1))])
    @test an_element(P) ∈ P

    # an_element for an unbounded polyhedron
    P = HPolyhedron([HalfSpace(N[-1, 0], N(-1))])
    y = an_element(P)
    # check that all entries are finite
    @test all(!isinf, y)
    # check that the points belong to P
    @test y ∈ P

    # boundedness
    @test isbounded(p)

    if test_suite_polyhedra
        p_unbounded = HPolyhedron([LinearConstraint(N[-1, 0], N(0))])
        p_infeasible = HPolyhedron([LinearConstraint(N[1], N(0)),
                                    LinearConstraint(N[-1], N(-1))])

        # support function/vector
        d = N[1, 0]
        @test σ(d, p_unbounded)[1] == N(Inf)
        @test ρ(d, p_unbounded) == N(Inf)
        @test_throws ErrorException σ(N[1], p_infeasible)

        # isempty
        @test !isempty(p_unbounded)
        @test isempty(p_infeasible)
        P = HPolyhedron([LinearConstraint(N[1, 0], N(0))])    # x <= 0
        addconstraint!(P, LinearConstraint(N[-1, 0], N(-1)))  # x >= 1
        @test isempty(P)

        # intersection
        A = N[0 1; 1 0]
        b = N[1, 1]
        p1 = HPolyhedron(A, b)
        A = N[0 -1; -1 0]
        b = N[0, 0]
        p2 = HPolyhedron(A, b)
        cap = intersection(p1, p2)
        @test ispermutation(constraints_list(cap),
                            constraints_list(BallInf(N[0.5, 0.5], N(0.5))))

        # intersection with polytope
        A = N[1 0;    # x <= 1
              0 1;    # y <= 1
              -1 -1]  # x + y >= 1
        b = N[1, 1, -1]
        p = HPolyhedron(A, b)
        b = BallInf(N[1, 1], N(0.5))
        cap1 = intersection(p, b)
        cap2 = intersection(b, p)
        @test ispermutation(constraints_list(cap1), constraints_list(cap2)) &&
              ispermutation(constraints_list(cap1),
                            constraints_list(BallInf(N[0.75, 0.75], N(0.25))))

        # intersection with half-space
        hs = HalfSpace(N[2, 2], N(-1))
        c1 = constraints_list(intersection(p1, hs))
        c2 = constraints_list(intersection(hs, p1))
        @test length(c1) == 3 && ispermutation(c1, c2)

        # convex hull
        ch = convex_hull(p1, p2)
        @test ch isa HPolyhedron{N} && isempty(constraints_list(ch))

        # Cartesian product
        A = N[1 -1]'
        b = N[1, 0]
        p1 = HPolyhedron(A, b)
        p2 = copy(p1)
        cp = cartesian_product(p1, p2)
        cl = constraints_list(cp)
        @test length(cl) == 4

        # vertices_list
        vertices_list(p1) ≈ [N[1.0], N[0.0]]

        # tovrep from HPolyhedron
        @test tohrep(p1) isa HPolyhedron{N} # test no-op

        # removing a redundant constraint in a 2D polyhedron (see #565)
        A = N[0 1;  # y <= 0
              1 0;  # x <= 0
              2 2]  # 2x + 2y <= 1 (the constraint is implied by the other two)
        b = N[0, 0, 1]
        p1 = HPolyhedron(A, b)
        remove_redundant_constraints!(p1)
        Ar, br = tosimplehrep(p1)
        @test Ar == A[1:2, :] && br == b[1:2]

        # removing redundant constraints from a list of constraints
        constraints = [HalfSpace(N[1], N(1)), HalfSpace(N[1], N(0))]
        constraints2 = remove_redundant_constraints(constraints)
        result = remove_redundant_constraints!(constraints)
        @test result
        @test ispermutation(constraints, constraints2)
        constraints = [HalfSpace(N[1], N(0)), HalfSpace(N[-1], N(-1)), HalfSpace(N[-1], N(-1))]
        result = remove_redundant_constraints!(constraints)
        @test !result

        # checking for empty intersection (also test symmetric calls)
        P = convert(HPolytope, BallInf(zeros(N, 2), N(1)))
        Q = convert(HPolytope, BallInf(ones(N, 2), N(1)))
        R = HPolyhedron([HalfSpace(N[1, 0], N(3)), HalfSpace(N[-1, 0], N(-2))])
        res, w = isdisjoint(P, Q, true)
        @test !isdisjoint(P, Q) && !res && w ∈ P && w ∈ Q
        res, w = isdisjoint(Q, P, true)
        @test !isdisjoint(Q, P) && !res && w ∈ P && w ∈ Q
        res, w = isdisjoint(Q, R, true)
        @test !isdisjoint(Q, R) && !res && w ∈ Q && w ∈ R
        res, w = isdisjoint(R, Q, true)
        @test !isdisjoint(R, Q) && !res && w ∈ Q && w ∈ R
        res, w = isdisjoint(P, R, true)
        @test isdisjoint(P, R) && res && w == N[]
        res, w = isdisjoint(R, P, true)
        @test isdisjoint(R, P) && res && w == N[]

        Punbdd = HPolyhedron([HalfSpace(N[0.68, 1.22], N(-0.76)),
                              HalfSpace(N[-0.75, -0.46], N(0.68))])
        @assert !isbounded(Punbdd)

        Pbdd = HPolyhedron([HalfSpace(N[0.68, 1.22], N(-0.76)),
                            HalfSpace(N[-0.75, -0.46], N(0.68)),
                            HalfSpace(N[1.72, 0.33], N(0.37)),
                            HalfSpace(N[-1.60, -0.41], N(0.67)),
                            HalfSpace(N[-0.44, 0.06], N(0.78))])
        @assert isbounded(Pbdd)

        Mnotinv = N[1 0; 2 0]
        @assert !isinvertible(Mnotinv)

        Minv = N[1 2; -1 0.4]
        @assert isinvertible(Minv)

        # invertible matrix times a bounded polyhedron
        L = linear_map(Minv, Pbdd)
        @test L isa HPolyhedron{N}

        # invertible matrix times an unbounded polyhedron
        L = linear_map(Minv, Punbdd)
        @test L isa HPolyhedron{N}

        # not invertible matrix times a bounded polyhedron
        L = linear_map(Mnotinv, Pbdd, algorithm="vrep") # Requires Polyhedra because it works on vertices
        @test L isa VPolytope

        L = linear_map(Mnotinv, Pbdd, algorithm="elimination")
        @test L isa HPolyhedron

        # test default
        L = linear_map(Mnotinv, Pbdd)
        @test L isa HPolyhedron

        # not invertible matrix times an unbounded polyhedron
        @test linear_map(Mnotinv, Punbdd) isa HPolyhedron

        # check that we can use sparse matrices as well ; Requires SparseArrays
        L = linear_map(sparse(Minv), Pbdd, algorithm="inv_right")
        @test L isa HPolyhedron{N}
        L = linear_map(sparse(Minv), Punbdd, algorithm="inv_right")
        @test L isa HPolyhedron{N}
        L = linear_map(sparse(Mnotinv), Pbdd, algorithm="vrep") # Requires Polyhedra because it works on vertices
        @test L isa VPolytope
        # breaks because "inv_right" requires an invertible matrix
        @test_throws ArgumentError linear_map(sparse(Mnotinv), Punbdd, algorithm="inv_right")

        # remove a repeated constraint (#909)
        Q = HPolyhedron([HalfSpace([-1.29817, 1.04012], 6.07731),
                         HalfSpace([-1.29348, -0.0920708], 1.89515)])
        addconstraint!(Q, Q.constraints[2])
        remove_redundant_constraints!(Q)
        @test length(constraints_list(Q)) == 2

        # concrete projection of an unbounded set
        # P = {x, y, z : x >= 0, y >= 0, z >= 0} is the positive orthant
        P = HPolyhedron([HalfSpace(N[-1, 0, 0], N(0)),
                         HalfSpace(N[0, -1, 0], N(0)),
                         HalfSpace(N[0, 0, -1], N(0))])
        πP = project(P, [1, 2])
        @test πP isa HPolyhedron{N}
        @test ispermutation(constraints_list(πP), [HalfSpace(N[-1, 0], N(0)), HalfSpace(N[0, -1], N(0))])
    end

    # tests that require ModelingToolkit
    @static if VERSION >= v"1.3" && isdefined(@__MODULE__, :ModelingToolkit)
        vars = @variables x y
        p1 = HPolyhedron([x + y <= 1, x + y >= -1,  x - y <= 1, x - y >= -1], vars)
        b1 = Ball1(zeros(2), 1.0)
        @test isequivalent(p1, b1)

        p2 = HPolyhedron([x == 0, y <= 0], vars)
        h2 = HPolyhedron([HalfSpace([1.0, 0.0], 0.0), HalfSpace([-1.0, 0.0], 0.0), HalfSpace([0.0, 1.0], 0.0)])
        @test p2 ⊆ h2 && h2 ⊆ p2 # isequivalent(p2, h2) see #2370
    end
end


text:
cholcov = function(x::Array{Float64, 2})
    tol = 1e-6
    if(isposdef(x))
        out = Array(cholesky(x).U)
    else
        eigX = eigen(Hermitian(x))
        val = sqrt.(eigX.values[eigX.values .> tol])
        vec = eigX.vectors[:,eigX.values .> tol]
        out = Diagonal(val) * vec'
    end
    out
end

function makeSVDdict(A::Array{Float64, 2}, B::Array{Float64, 2})
    U,V,Q,C,S,R = svd(cholcov(A), cholcov(B))
    #svd_Design_Ainv.Q = inv(q)'
    #svd_Design_Ainv.
    #svd_Design_Ainv.s2 = diag(S2'*S2)
    #Qt_Design = svd_Design_Ainv.Q'*Design;
    H = R * Q'
    q = H'
    Q = ("Q", inv(q)')
    s1 = ("s1", Array(diag(C' * C)))
    s2 = ("s2", Array(diag(S' * S)))
    Dict([Q, s1, s2])
end

function cov2cor(x::Array{Float64, 2})
    sds = sqrt.(diag(x))
    x ./ (sds * sds')
end


text:
#=
    Contents: Utilities for authorization.
=#

function has_role(req::Request, role::AbstractString)
    result = false
    if is_logged_in(req)
	username = get_username(req)
	result   = has_role(username, role)
    end
    result
end

# EOF


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct CreateDataLakeStoreAccountProperties <: SwaggerModel
    defaultGroup::Any # spec type: Union{ Nothing, String } # spec name: defaultGroup
    encryptionConfig::Any # spec type: Union{ Nothing, EncryptionConfig } # spec name: encryptionConfig
    encryptionState::Any # spec type: Union{ Nothing, String } # spec name: encryptionState
    firewallRules::Any # spec type: Union{ Nothing, Vector{CreateFirewallRuleWithAccountParameters} } # spec name: firewallRules
    virtualNetworkRules::Any # spec type: Union{ Nothing, Vector{CreateVirtualNetworkRuleWithAccountParameters} } # spec name: virtualNetworkRules
    firewallState::Any # spec type: Union{ Nothing, String } # spec name: firewallState
    firewallAllowAzureIps::Any # spec type: Union{ Nothing, String } # spec name: firewallAllowAzureIps
    trustedIdProviders::Any # spec type: Union{ Nothing, Vector{CreateTrustedIdProviderWithAccountParameters} } # spec name: trustedIdProviders
    trustedIdProviderState::Any # spec type: Union{ Nothing, String } # spec name: trustedIdProviderState
    newTier::Any # spec type: Union{ Nothing, String } # spec name: newTier

    function CreateDataLakeStoreAccountProperties(;defaultGroup=nothing, encryptionConfig=nothing, encryptionState=nothing, firewallRules=nothing, virtualNetworkRules=nothing, firewallState=nothing, firewallAllowAzureIps=nothing, trustedIdProviders=nothing, trustedIdProviderState=nothing, newTier=nothing)
        o = new()
        validate_property(CreateDataLakeStoreAccountProperties, Symbol("defaultGroup"), defaultGroup)
        setfield!(o, Symbol("defaultGroup"), defaultGroup)
        validate_property(CreateDataLakeStoreAccountProperties, Symbol("encryptionConfig"), encryptionConfig)
        setfield!(o, Symbol("encryptionConfig"), encryptionConfig)
        validate_property(CreateDataLakeStoreAccountProperties, Symbol("encryptionState"), encryptionState)
        setfield!(o, Symbol("encryptionState"), encryptionState)
        validate_property(CreateDataLakeStoreAccountProperties, Symbol("firewallRules"), firewallRules)
        setfield!(o, Symbol("firewallRules"), firewallRules)
        validate_property(CreateDataLakeStoreAccountProperties, Symbol("virtualNetworkRules"), virtualNetworkRules)
        setfield!(o, Symbol("virtualNetworkRules"), virtualNetworkRules)
        validate_property(CreateDataLakeStoreAccountProperties, Symbol("firewallState"), firewallState)
        setfield!(o, Symbol("firewallState"), firewallState)
        validate_property(CreateDataLakeStoreAccountProperties, Symbol("firewallAllowAzureIps"), firewallAllowAzureIps)
        setfield!(o, Symbol("firewallAllowAzureIps"), firewallAllowAzureIps)
        validate_property(CreateDataLakeStoreAccountProperties, Symbol("trustedIdProviders"), trustedIdProviders)
        setfield!(o, Symbol("trustedIdProviders"), trustedIdProviders)
        validate_property(CreateDataLakeStoreAccountProperties, Symbol("trustedIdProviderState"), trustedIdProviderState)
        setfield!(o, Symbol("trustedIdProviderState"), trustedIdProviderState)
        validate_property(CreateDataLakeStoreAccountProperties, Symbol("newTier"), newTier)
        setfield!(o, Symbol("newTier"), newTier)
        o
    end
end # type CreateDataLakeStoreAccountProperties

const _property_map_CreateDataLakeStoreAccountProperties = Dict{Symbol,Symbol}(Symbol("defaultGroup")=>Symbol("defaultGroup"), Symbol("encryptionConfig")=>Symbol("encryptionConfig"), Symbol("encryptionState")=>Symbol("encryptionState"), Symbol("firewallRules")=>Symbol("firewallRules"), Symbol("virtualNetworkRules")=>Symbol("virtualNetworkRules"), Symbol("firewallState")=>Symbol("firewallState"), Symbol("firewallAllowAzureIps")=>Symbol("firewallAllowAzureIps"), Symbol("trustedIdProviders")=>Symbol("trustedIdProviders"), Symbol("trustedIdProviderState")=>Symbol("trustedIdProviderState"), Symbol("newTier")=>Symbol("newTier"))
const _property_types_CreateDataLakeStoreAccountProperties = Dict{Symbol,String}(Symbol("defaultGroup")=>"String", Symbol("encryptionConfig")=>"EncryptionConfig", Symbol("encryptionState")=>"String", Symbol("firewallRules")=>"Vector{CreateFirewallRuleWithAccountParameters}", Symbol("virtualNetworkRules")=>"Vector{CreateVirtualNetworkRuleWithAccountParameters}", Symbol("firewallState")=>"String", Symbol("firewallAllowAzureIps")=>"String", Symbol("trustedIdProviders")=>"Vector{CreateTrustedIdProviderWithAccountParameters}", Symbol("trustedIdProviderState")=>"String", Symbol("newTier")=>"String")
Base.propertynames(::Type{ CreateDataLakeStoreAccountProperties }) = collect(keys(_property_map_CreateDataLakeStoreAccountProperties))
Swagger.property_type(::Type{ CreateDataLakeStoreAccountProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_CreateDataLakeStoreAccountProperties[name]))}
Swagger.field_name(::Type{ CreateDataLakeStoreAccountProperties }, property_name::Symbol) =  _property_map_CreateDataLakeStoreAccountProperties[property_name]

const _allowed_CreateDataLakeStoreAccountProperties_encryptionState = ["Enabled", "Disabled"]

const _allowed_CreateDataLakeStoreAccountProperties_firewallState = ["Enabled", "Disabled"]

const _allowed_CreateDataLakeStoreAccountProperties_firewallAllowAzureIps = ["Enabled", "Disabled"]

const _allowed_CreateDataLakeStoreAccountProperties_trustedIdProviderState = ["Enabled", "Disabled"]

const _allowed_CreateDataLakeStoreAccountProperties_newTier = ["Consumption", "Commitment_1TB", "Commitment_10TB", "Commitment_100TB", "Commitment_500TB", "Commitment_1PB", "Commitment_5PB"]

function check_required(o::CreateDataLakeStoreAccountProperties)
    true
end

function validate_property(::Type{ CreateDataLakeStoreAccountProperties }, name::Symbol, val)
    if name === Symbol("encryptionState")
        Swagger.validate_param(name, "CreateDataLakeStoreAccountProperties", :enum, val, _allowed_CreateDataLakeStoreAccountProperties_encryptionState)
    end
    if name === Symbol("firewallState")
        Swagger.validate_param(name, "CreateDataLakeStoreAccountProperties", :enum, val, _allowed_CreateDataLakeStoreAccountProperties_firewallState)
    end
    if name === Symbol("firewallAllowAzureIps")
        Swagger.validate_param(name, "CreateDataLakeStoreAccountProperties", :enum, val, _allowed_CreateDataLakeStoreAccountProperties_firewallAllowAzureIps)
    end
    if name === Symbol("trustedIdProviderState")
        Swagger.validate_param(name, "CreateDataLakeStoreAccountProperties", :enum, val, _allowed_CreateDataLakeStoreAccountProperties_trustedIdProviderState)
    end
    if name === Symbol("newTier")
        Swagger.validate_param(name, "CreateDataLakeStoreAccountProperties", :enum, val, _allowed_CreateDataLakeStoreAccountProperties_newTier)
    end
end


text:
export decimation, masterTrigger, currentFrame, ramWriterMode, connectADC, startADC, stopADC, samplesPerPeriod, periodsPerFrame,
     currentWP, slowDACInterpolation, bufferSize, keepAliveReset, triggerMode, numSlowADCChan,
     ADCPerformanceData, RPPerformance, RPStatus, RPInfo, startPipelinedData, PerformanceData, numChan, dataRate

struct ADCPerformanceData
  deltaRead::UInt64
  deltaSend::UInt64
  chunkSize::UInt64
end

struct RPStatus
  overwritten::Bool
  corrupted::Bool
  stepsLost::Bool
  adcEnabled::Bool
  dacEnabled::Bool
end

struct PerformanceData
  wpRead::UInt64
  adc::ADCPerformanceData
  dac::DACPerformanceData
  status::RPStatus
end

struct RPPerformance
  data::Vector{PerformanceData}
end

struct RPInfo
  performances::Vector{RPPerformance}
end

function RPInfo()
  return RPInfo([RPPerformance([])])
end

function dataRate(chunkSize, deltaSend, decimation; unit="Mbits")
  bitsPerSample = (chunkSize * 4 * 8) / deltaSend # 4 bytes per samples
  freq = div(125e6, decimation)
  bitsPerSec = bitsPerSample * freq
  return bitsPerSec/1e6 #TODO unit conversion
end

function dataRate(adc::ADCPerformanceData, decimation; unit ="Mbits")
  return dataRate(adc.chunkSize, adc.deltaSend, decimation)
end

decimation(rp::RedPitaya) = query(rp,"RP:ADC:DECimation?", Int64)
function decimation(rp::RedPitaya, dec)
  rp.decimation = Int64(dec)
  send(rp, string("RP:ADC:DECimation ", rp.decimation))
end

numChan(rp::RedPitaya) = 2

function slowDACInterpolation(rp::RedPitaya, enable::Bool)
  enableI = Int32(enable)
  send(rp, string("RP:ADC:SlowDACInterpolation ", enableI))
end

numSlowADCChan(rp::RedPitaya) = query(rp,"RP:ADC:SlowADC?", Int64)
function numSlowADCChan(rp::RedPitaya, value)
  send(rp, string("RP:ADC:SlowADC ", Int64(value)))
end

function samplesPerPeriod(rp::RedPitaya)
  return rp.samplesPerPeriod
end
function samplesPerPeriod(rp::RedPitaya, value)
  rp.samplesPerPeriod = value
end

function periodsPerFrame(rp::RedPitaya)
  return rp.periodsPerFrame
end
function periodsPerFrame(rp::RedPitaya, value)
  rp.periodsPerFrame = value
end

function currentFrame(rp::RedPitaya)
  return Int64(floor(currentWP(rp) / (rp.samplesPerPeriod * rp.periodsPerFrame)))
end

function currentPeriod(rp::RedPitaya)
  return Int64(floor(currentWP(rp) / (rp.samplesPerPeriod)))
end

currentWP(rp::RedPitaya) = query(rp,"RP:ADC:WP:CURRENT?", Int64)
bufferSize(rp::RedPitaya) = query(rp,"RP:ADC:BUFFER:SIZE?", Int64)

function masterTrigger(rp::RedPitaya, val::Bool)
  valStr = val ? "ON" : "OFF"
  send(rp, string("RP:MasterTrigger ", valStr))
end
masterTrigger(rp::RedPitaya) = occursin("ON", query(rp,"RP:MasterTrigger?"))

function keepAliveReset(rp::RedPitaya, val::Bool)
  valStr = val ? "ON" : "OFF"
  send(rp, string("RP:KeepAliveReset ", valStr))
end
keepAliveReset(rp::RedPitaya) = occursin("ON", query(rp,"RP:KeepAliveReset?"))

# "TRIGGERED" or "CONTINUOUS"
function ramWriterMode(rp::RedPitaya, mode::String)
  send(rp, string("RP:RamWriterMode ", mode))
end

# "INTERNAL" or "EXTERNAL"
function triggerMode(rp::RedPitaya, mode::String)
  send(rp, string("RP:Trigger:Mode ", mode))
end

connectADC(rp::RedPitaya) = send(rp, "RP:ADC:ACQCONNect")
function startADC(rp::RedPitaya)
  send(rp, "RP:ADC:ACQSTATUS ON")
end
stopADC(rp::RedPitaya) = send(rp, "RP:ADC:ACQSTATUS OFF")

wasOverwritten(rp::RedPitaya) = query(rp, "RP:STATus:OVERwritten?", Bool)
wasCorrupted(rp::RedPitaya) = query(rp, "RP:STATus:CORRupted?", Bool)

function serverStatus(rp::RedPitaya)
  send(rp, "RP:STATus?")
  return readServerStatus(rp)
end

function readServerStatus(rp::RedPitaya)
  statusRaw = read!(rp.dataSocket, Array{Int8}(undef, 1))[1]
  status = RPStatus((statusRaw & 1) != 0, # overwritten
   (statusRaw & (1 << 1)) != 0, # corrupted
   (statusRaw & (1 << 2)) != 0, # stepsLost
   (statusRaw & (1 << 3)) != 0, # adcEnabled
   (statusRaw & (1 << 4)) != 0) # dacEnabled
  return status
end

function performanceData(rp::RedPitaya, numSamples = 0)
  send(rp, "RP:PERF?")
  return readPerformanceData(rp, numSamples)
end

function readPerformanceData(rp::RedPitaya, numSamples = 0)
  adc = readADCPerformanceData(rp, numSamples)
  dac = readDACPerformanceData(rp)
  return (adc, dac)
end

function readADCPerformanceData(rp::RedPitaya, numSamples = 0)
  perf = read!(rp.dataSocket, Array{UInt64}(undef, 2))
  return ADCPerformanceData(perf[1], perf[2], numSamples)
end

# Low level read. One has to take care that the numFrames are available
function readSamples_(rp::RedPitaya, reqWP, numSamples)
  command = string("RP:ADC:DATA? ",Int64(reqWP),",",Int64(numSamples))
  send(rp, command)

  @debug "read data ..."
  u = read!(rp.dataSocket, Array{Int16}(undef, 2 * Int64(numSamples)))
  @debug "read data!"
  return u
end

# Low level read, reads samples, error and perf. Values need to be already requested
function readSamplesChunk_(rp::RedPitaya, reqWP::Int64, numSamples::Int64, into=nothing)
  @debug "read samples chunk ..."
  if isnothing(into)
    into = Array{Int16}(undef, 2 * Int64(numSamples))
  end
  data = read!(rp.dataSocket, into)
  status = readServerStatus(rp)
  (adc, dac) = readPerformanceData(rp, numSamples)
  @debug "read samples chunk ..."
  perf = PerformanceData(UInt64(reqWP), adc, dac, status)
  return (data, perf)
end

# Low level read, that includes performance and error data
function readDetailedSamples_(rp::RedPitaya, reqWP::Int64, numSamples::Int64)
  command = string("RP:ADC:DATA:DETAILED? ",Int64(reqWP),",",Int64(numSamples))
  send(rp, command)
  return readSamplesChunk_(rp, reqWP, numSamples)
end

function readSamplesIntermediate_(rp::RedPitaya, reqWP::Int64, numSamples::Int64)
  data = readSamples_(rp, reqWP, numSamples)
  status = serverStatus(rp)
  (adc, dac) = performanceData(rp)
  perf = PerformanceData(UInt64(reqWP), adc, dac, status)
  return (data, perf)
end

function readSamplesOld_(rp::RedPitaya, reqWP::Int64, numSamples::Int64)
  data = readSamples_(rp, reqWP, numSamples)
  overwritten = wasOverwritten(rp)
  corrupted = wasCorrupted(rp)
  status = RPStatus(overwritten, corrupted, false, true, true)
  (adc, dac) = performanceData(rp)
  perf = PerformanceData(UInt64(reqWP), adc, dac, status)
  return (data, perf)
end

function startPipelinedData(rp::RedPitaya, reqWP::Int64, numSamples::Int64, chunkSize::Int64)
  command = string("RP:ADC:DATA:PIPELINED? ", reqWP, ",", numSamples, ",", chunkSize)
  send(rp, command)
end

# Low level read. One has to take care that the numFrames are available
function readDataSlow_(rp::RedPitaya, startFrame, numFrames)
  numPeriods = rp.periodsPerFrame
  numChan = numSlowADCChan(rp)

  command = string("RP:ADC:SLOW:FRAMES:DATA ",Int64(startFrame),",",Int64(numFrames))
  #send(rp, command)

  @debug "read data ..."
  #u = read!(rp.dataSocket, Array{Float32}(undef, numChan * numFrames * numPeriods))
  @debug "read data!"
  #return reshape(u, numChan, numPeriods, numFrames)
  return zeros(Float32, numChan, numPeriods, numFrames)
end

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: tknopp/RedPitayaDAQServer

max_forks_repo_name: tknopp/RedPitayaDAQServer