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

using Markdown
using InteractiveUtils

# ╔═╡ 14d058ad-8ef9-4929-8975-c63f5b4b45ac
begin
	using Pkg
	Pkg.add("ForwardDiff")
end

# ╔═╡ 631af3ab-cf2d-47c5-bdb9-014ec30bf2b7
md"""
# Intermediate Julia for scientific computing
"""

# ╔═╡ adadf28b-fc27-487a-8264-92747c0df902
md"""
This workshop is designed to introduce two fundamental concepts in Julia: **types** and **metaprogramming**.

In order to cover various key uses of types in Julia, we have chosen to frame the discussion around a concrete topic in scientific computing, namely **root-finding**.
The goal is *not* to learn algorithms for root finding *per se*, but rather to have a (pseudo-)real context in which to explore various concepts centered around types and how they arise naturally in real applications of Julia, in particular applications of **multiple dispatch**, which is one of the core choices in Julia that differentiate it from other common languages.

We will implement a couple of root-finding algorithms just to have something to work with. These will just be toy implementations that are far away from the best implementations.

Instead we should use one of the high-quality packages that are available in Julia for this purpose. The large number of them shows the importance of root finding. The ones that I am aware of are the following (in alphabetical order):
"""

# ╔═╡ 63574266-dcf4-4430-994e-001f20bad887
md"""
- Single root of a nonlinear function:
    - [`NLsolve.jl`](https://github.com/JuliaNLSolvers/NLsolve.jl)
    - [`Roots.jl`](https://github.com/JuliaMath/Roots.jl)

- All roots of polynomial:
    - [`HomotopyContinuation.jl`](https://www.juliahomotopycontinuation.org)
    - [`PolynomialRoots.jl`](https://github.com/giordano/PolynomialRoots.jl)
    - [`Polynomials.jl`](https://github.com/JuliaMath/Polynomials.jl)

- All roots of a nonlinear function:
    - [`ApproxFun.jl`](https://github.com/JuliaApproximation/ApproxFun.jl)
    - [`IntervalRootFinding.jl`](https://github.com/JuliaIntervals/IntervalRootFinding.jl)
    - [`MDBM.jl`](https://github.com/bachrathyd/MDBM.jl)
    - [`Roots.jl`](https://github.com/JuliaMath/Roots.jl)
"""

# ╔═╡ 1cf19f60-7e52-406b-b2c9-baad3544aea8
md"""
Each of these uses different techniques, with different advantages and disadvantages.
"""

# ╔═╡ 77899e06-1642-4336-b5ca-daed74fad1b4
md"""
The challenge exercise for the workshop is: develop a package which integrates all of these disparate packages into a coherent whole!
"""

# ╔═╡ c235a33d-8e07-4018-ad5e-8617e084f573
md"""
### Logistics of the workshop
"""

# ╔═╡ c9f7f7e3-b5c3-4296-a56a-d941fae5efb0
md"""
The workshop is based around a series of exercises to be done during the workshop. We will pause to work on the exercises and then I will discuss possible solutions during the workshop.
"""

# ╔═╡ f5ec31bb-d2d4-4372-aeb6-940a6d9ccaa9
md"""
These techniques are useful for both users and developers; indeed, in Julia the distinction between users and developers is not useful, since it's much easier than in other languages to join the two categories together.
"""

# ╔═╡ 31266442-3d7a-449e-8135-88301250f5b3
md"""
### Outline
"""

# ╔═╡ 77261e8f-a2f5-48bb-a2ac-67808cfb78b4
md"""
We will start by quickly reviewing roots of functions and quickly reviewing one of the standard algorithms, **Newton's algorithm**. We will restrict to finding roots of 1D functions for simplicity.

Newton's algorithm requires the calculation of derivatives, for which several choices of algorithm are available. We will see how to encode the choice of algorithm using dispatch.

Then we will define types which will contain all information about a root-finding problem.
"""

# ╔═╡ 572b5bcc-5443-4194-bd3a-d735aa7363a4
md"""
## Roots
"""

# ╔═╡ 641b531f-a2f1-486f-a2c9-b5d566ac323a
md"""
Given a function ``f: \mathbb{R} \to \mathbb{R}`` (i.e. that accepts a single real number as argument and returns another real number), recall that a **root** or **zero** of the function is a number $x^*$ such that

$$f(x^*) = 0,$$

i.e. it is a solution of the equation $f(x) = 0$.

In general it is impossible to solve this equation exactly for $x^*$, so we use iterative numerical algorithms instead.
"""

# ╔═╡ 204fafc2-4005-4c97-8835-8c51dd6ba43b
md"""
#### Example
"""

# ╔═╡ f269a3f7-491e-476a-8251-bf7d892526ec
md"""
Recall that the function $f$ given by $f(x) := x^2 - 2$ has exactly two roots, at $x^*_1 = +\sqrt{2}$ and $x^*_2 = -\sqrt{2}$. Note that it is impossible to represent these values exactly using floating-point arithmetic.
"""

# ╔═╡ 8a123283-a780-4ec4-8ca1-693c661087c3
md"""
## Newton algorithm
"""

# ╔═╡ 6e8e8929-d22e-45d6-afd0-fd56a2317efd
md"Based on the following first order approximation:"

# ╔═╡ e44493cf-f7fa-40fa-8b33-3ee76052d24b
md"""
$${f(x_{n+1})} = f(x_n) + (x_{n+1} - x_n) {f'(x_n)}$$
"""

# ╔═╡ 45e4574f-f352-437b-92f8-2bbf3dc650cd
md"""
Hypothesizing ``x_{n+1}`` is a root

``{f(x_{n+1})} = f(x_n) + (x_{n+1} - x_n) {f'(x_n)} = 0``
"""

# ╔═╡ c4aa101d-7f7c-4aff-ada3-f1253962859c
md"""
The Newton algorithm for (possibly) finding a root of a nonlinear function $f(x)$ in 1D is the following iteration:
"""

# ╔═╡ 76159eac-ae3f-4999-af01-19c3617aeaab
md"""
$$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)},$$

where $f'$ is the derivative of $f$. We start from an initial guess $x_0$ that can be almost anything (except points for which $f'(x_0) = 0$).
"""

# ╔═╡ 1f78d275-7b6f-48cf-8358-60d7ef2c51d1
md"""
#### Exercise 1
"""

# ╔═╡ 1a5f54d4-434e-4a3d-98df-b2b70f023d10
md"""
1. Implement the Newton algorithm for a fixed number $n$ of steps in a function `newton`, starting from a given starting point $x_0$.

    Hint: Which information does the function require?

2. Does your function work with other number types, such as `BigFloat`? What do you need in order to run it with those types? Use it to calculate $\sqrt{2}$. How many decimal places are correct with the standard precision of `BigFloat`?
"""

# ╔═╡ e67d1b09-aa18-4914-811e-b695d256d801
f(x) = x^2 - 2

# ╔═╡ 5e738a26-487a-4879-8806-785dab351683
df(x) = 2x

# ╔═╡ ef85ca3b-f7af-4205-8dbe-b928a6df3f9a
"""
	newton_steps(f::Function, df::Function, x0::Float64, nsteps::Int64=10)

Newton root-finding of function `f`
with explicit derivative fn `df`.
Starting from `x₀`,
Returns `nsteps` toward root x: f(x) == 0.0
"""
function newton_steps(f::Function, df::Function, x₀::Float64, nsteps::Int64=10)
	x = x₀
	for n in 1:nsteps
		x = x - f(x)/df(x)
	end
	return x
end

# ╔═╡ 15e80f75-f150-43a5-95ec-1eded7a2ee3d
"""
	newton(f::Function, df::Function, x0::Float64, n::Int64=10)

Newton root-finding of function `f`
with explicit derivative fn `df`.
Returns root x: f(x) == 0.0
"""
function newton(f::Function, df::Function, x0::Float64, n::Int64=10)
	x = x0
	x_new = -3.5

	num_iterations = 0

	for i in 1:n
		x_new = x - f(x) / df(x)

		x = x_new
	end

	return x_new
end

# ╔═╡ c2de19e1-cb00-4435-89a1-dedd7248fe66

# ╔═╡ 68e23b8d-e74e-495d-acec-5f608cea36fe
md"This algorithm oscillates between two near solutions, at numerical precision."

# ╔═╡ 3121853d-5bd1-4c16-9a0e-e5239bb450b7
newton_steps(f, df, 3.0, 7)

# ╔═╡ 3cb7e4a1-fecd-4a55-a0fa-759819c96fee
newton_steps(f, df, 3.0, 8)

# ╔═╡ de25cfec-8080-4cce-acf8-56e0028a3a3d
newton_steps(f, df, 3.0)

# ╔═╡ f0b69144-d578-4808-b7e2-13d16abf4abc
1e-128

# ╔═╡ 3417b888-c922-4e00-9053-afbda8038fcf
f(newton_steps(f, df, 3.0, 9)) # ≈ 0.0

# ╔═╡ 784b5ef2-e53d-431b-bedf-71c85045bfb0
sqrt(2)

# ╔═╡ a892e53e-9f96-408e-87b2-affc1340554b
md"""
## Calculating derivatives
"""

# ╔═╡ a5c74608-47c6-4034-b046-2a9885c0da9d
md"""
The Newton algorithm requires us to specify the derivative of a function. If $f$ is a complicated function, we certainly don't want to do that by hand.

One standard solution is to use a *finite-difference approximation* of the first derivative:

$$f'(a) \simeq \frac{f(a + h) - f(a - h)}{2h}.$$
Here we use the two-sided finite-difference approximation.
"""

# ╔═╡ 83f25a9b-24ca-45e9-bed6-c81fa283761f
md"""
#### Exercise 2
"""

# ╔═╡ 48bb1b5c-63d7-4fc7-a75b-13f71105095f
md"""
1. Implement a function `finite_difference` with a default value $h = 0.001$.

2. Implement a version of `newton` that does not take the derivative as argument and uses `finite_difference` to calculate the derivative. This version of `newton` should **re-use** the previous version by defining the function `fp` and calling that version.
"""

# ╔═╡ 3b84e81a-3e36-41c2-b16b-20a4e8c95f5c
"""
	finite_difference(f::Function, a::Float64, h::Float64=0.001)

two-side finite difference of `f` at `a`, step-size `h`.
"""
function finite_difference(f::Function, a::Float64, h::Float64=0.001)
    return ( f(a + h) - f(a - h) ) / (2h)
end

# ╔═╡ df21662e-3f94-4765-8925-9ab648e82752
"""
	newton(f::Function, x0::Float64, n::Int64=10)

Newton root-finding of function `f`
without explict derivative (use finite difference approx).
Returns root x: f(x) == 0.0
"""
function newton(f::Function, x0::Float64, n::Int64=10)
	x = x0
	x_new = -3.5

	num_iterations = 0

	for i in 1:n
		df = finite_difference(f, x)

		x_new = x - f(x) / df

		x = x_new
	end

	return x_new
end

# ╔═╡ aad415a8-bdd1-4698-a097-4f1a4ad4e3d9
begin
	g(x) = x^2 - 2
	dg(x) = 2x
end

# ╔═╡ 9e72e71f-cb58-4f88-a997-f69dc7b8ca6f
begin
	a = 3.0
	finite_difference(g, a)
end

# ╔═╡ 6d76b715-960f-4d27-bdb4-6da34af3959f
dg(a)

# ╔═╡ 8a03afb2-76ed-4a73-a909-22f6f0f0b56c
# finite diff with a larger step-size (less accurate)
finite_difference(g, a, 0.01)

# ╔═╡ 97151bb9-6905-48d7-a193-6b25985b2986
md"""
Mathematically: $x \mapsto f\prime(x)$
"""

# ╔═╡ a9368f4d-350e-4978-bec8-1bc873fa32a1
# this version defines the derivative fn as the finite diff,
# and then calls the newton method with the fn and its (approx) derivative.

function newton(f::Function, x0::Float64)
    # we need to calculate the derivative function dg
    df = x -> finite_difference(f, x)

    return newton(f, df, x0)   # using the previous method we defined
end

# ╔═╡ f3232e7a-2e69-4337-9352-4a4a9cc2e7ff
md"""
`newton` is a **generic function**
"""

# ╔═╡ e7d1b1be-ebc1-4f62-83a5-4e758047e1f1
md"""
### Algorithmic differentiation
"""

# ╔═╡ 85104806-f773-4785-885c-9d32cf51b93d
md"""
An alternative way to calculate derivatives is by using [**algorithmic differentiation**](https://en.wikipedia.org/wiki/Automatic_differentiation) (also called **automatic differentiation** or **computational differentiation**). This gives exact results (up to rounding error).

We will implement this algorithm in the next notebook, but for now let's just use the implementation in the excellent [`ForwardDiff.jl` package](https://github.com/JuliaDiff/ForwardDiff.jl).

"""

# ╔═╡ d394a4e8-c08c-4086-95a8-39901994b2ca
md"""
#### Exercise 3
"""

# ╔═╡ 6b5aa387-8f48-45fe-b9be-b0d427e833f1
md"""
1. Install `ForwardDiff.jl` if necessary.

2. Import it.

3. Define a function `forwarddiff` that uses the `ForwardDiff.derivative` function to calculate a derivative.
"""

# ╔═╡ 481a94d9-bc4f-4b0b-b808-65b4c80ec0ac
# ]add ForwardDiff

# ╔═╡ bdcce296-bffa-45fe-a658-1116331d9bb6
import ForwardDiff

# ╔═╡ 63ccc7f4-8cad-4cb3-b815-deff9f3f9a59
forwarddiff(f, x) = ForwardDiff.derivative(f, x)

# ╔═╡ 6cf3f90d-9c30-4cfe-a3eb-940f51b49dea
forwarddiff(g, a)

# ╔═╡ 064994f6-3823-4e43-a71d-00774804fd20
md"""
### Choosing between algorithms
"""

# ╔═╡ 718cec78-ed34-455f-9c24-c75b8b7d1f71
md"""
We now have two different algorithms available to calculate derivatives. This kind of situation is common in scientific computing; for example, the [`DifferentialEquations.jl`](http://docs.juliadiffeq.org/latest/) ecosystem has some 300 algorithms for solving differential equations. One of the techniques we will learn is how to easily be able to specify different algorithms.
"""

# ╔═╡ 92ecf5cd-d25d-45f9-95dc-504c354183cf
md"""
One possible solution is just by specifying the *function* to use as an argument to another function:
"""

# ╔═╡ 6189f7b2-9331-4f2b-bb58-760e46b03d82
md"""
#### Exercise 4
"""

# ╔═╡ a36de957-9ba5-45ad-8c82-3a21d3090fb7
md"""
1. Make a version of the Newton algorithm that takes an argument which is the algorithm to use to calculate the derivative.

    Hint: Here we will just pass in the function as an argument by giving its name as a parameter.
"""

# ╔═╡ c6fbb9e8-91a3-4c67-981b-40043198a1bf
md"""
Make this a keyword argument
"""

# ╔═╡ 200bd89d-a14a-4e4f-ac5c-a0b5e769e2ad
function newton(f, x, derivative::Function)
     df = x -> derivative(f, x)

    return newton(f, df, x)
end

# ╔═╡ e2d57494-2109-41c7-afea-4e4b7f014538
newton(f, df, 3.0)

# ╔═╡ 90ca5a0b-41ca-49f7-b673-31a0c5bd97cc
newton(f, df, 3.0) ≈ sqrt(2)

# ╔═╡ d4dca943-58fe-40d1-b23a-eeb3c21a273a
methods(newton)

# ╔═╡ 651df39c-60ca-4e0d-a6fd-0f13722a7572
newton(g, a)

# ╔═╡ 81a95d09-e0f7-418e-bcf7-dca908c6384c
@which newton(g, a)

# ╔═╡ a9eca302-7939-4d9b-9b3c-48680c734379
methods(newton)

# ╔═╡ c266a1bf-94e6-43fd-ae7f-a5ccdef7cf4b
methods(newton)

# ╔═╡ d5bf6a90-e938-43dd-a852-8f5b2f13bb49
newton(g, a; derivative=forwarddiff)

# ╔═╡ b2c9ee3b-6805-4bd1-9bee-06714190cbe8
newton(g, a; derivative=finite_difference)

# ╔═╡ ef59ee59-1364-47aa-af21-7742e90a108e
newton(g, a, forwarddiff)

# ╔═╡ 6034cc3d-6153-4a71-a200-9f15a643677c
@which newton(g, a; derivative=forwarddiff)

# ╔═╡ 9b639977-9310-4682-a5a6-32ea7c4f7ae3
@which newton(g, a, forwarddiff)

# ╔═╡ bc4461c8-e563-4cf9-a53b-2a70b343efd7
md"""
For some purposes this technique of passing functions may be sufficient. In a later notebook we will see a more powerful technique to specify different algorithms.
"""

# ╔═╡ Cell order:
# ╟─631af3ab-cf2d-47c5-bdb9-014ec30bf2b7
# ╟─adadf28b-fc27-487a-8264-92747c0df902
# ╟─63574266-dcf4-4430-994e-001f20bad887
# ╟─1cf19f60-7e52-406b-b2c9-baad3544aea8
# ╟─77899e06-1642-4336-b5ca-daed74fad1b4
# ╟─c235a33d-8e07-4018-ad5e-8617e084f573
# ╟─c9f7f7e3-b5c3-4296-a56a-d941fae5efb0
# ╟─f5ec31bb-d2d4-4372-aeb6-940a6d9ccaa9
# ╟─31266442-3d7a-449e-8135-88301250f5b3
# ╟─77261e8f-a2f5-48bb-a2ac-67808cfb78b4
# ╟─572b5bcc-5443-4194-bd3a-d735aa7363a4
# ╟─641b531f-a2f1-486f-a2c9-b5d566ac323a
# ╟─204fafc2-4005-4c97-8835-8c51dd6ba43b
# ╟─f269a3f7-491e-476a-8251-bf7d892526ec
# ╟─8a123283-a780-4ec4-8ca1-693c661087c3
# ╟─6e8e8929-d22e-45d6-afd0-fd56a2317efd
# ╟─e44493cf-f7fa-40fa-8b33-3ee76052d24b
# ╟─45e4574f-f352-437b-92f8-2bbf3dc650cd
# ╟─c4aa101d-7f7c-4aff-ada3-f1253962859c
# ╟─76159eac-ae3f-4999-af01-19c3617aeaab
# ╟─1f78d275-7b6f-48cf-8358-60d7ef2c51d1
# ╟─1a5f54d4-434e-4a3d-98df-b2b70f023d10
# ╠═e67d1b09-aa18-4914-811e-b695d256d801
# ╠═5e738a26-487a-4879-8806-785dab351683
# ╠═ef85ca3b-f7af-4205-8dbe-b928a6df3f9a
# ╠═15e80f75-f150-43a5-95ec-1eded7a2ee3d
# ╠═c2de19e1-cb00-4435-89a1-dedd7248fe66
# ╠═df21662e-3f94-4765-8925-9ab648e82752
# ╠═e2d57494-2109-41c7-afea-4e4b7f014538
# ╟─68e23b8d-e74e-495d-acec-5f608cea36fe
# ╠═3121853d-5bd1-4c16-9a0e-e5239bb450b7
# ╠═3cb7e4a1-fecd-4a55-a0fa-759819c96fee
# ╠═de25cfec-8080-4cce-acf8-56e0028a3a3d
# ╠═f0b69144-d578-4808-b7e2-13d16abf4abc
# ╠═3417b888-c922-4e00-9053-afbda8038fcf
# ╠═784b5ef2-e53d-431b-bedf-71c85045bfb0
# ╠═90ca5a0b-41ca-49f7-b673-31a0c5bd97cc
# ╟─a892e53e-9f96-408e-87b2-affc1340554b
# ╟─a5c74608-47c6-4034-b046-2a9885c0da9d
# ╟─83f25a9b-24ca-45e9-bed6-c81fa283761f
# ╟─48bb1b5c-63d7-4fc7-a75b-13f71105095f
# ╠═3b84e81a-3e36-41c2-b16b-20a4e8c95f5c
# ╠═aad415a8-bdd1-4698-a097-4f1a4ad4e3d9
# ╠═9e72e71f-cb58-4f88-a997-f69dc7b8ca6f
# ╠═6d76b715-960f-4d27-bdb4-6da34af3959f
# ╠═8a03afb2-76ed-4a73-a909-22f6f0f0b56c
# ╟─97151bb9-6905-48d7-a193-6b25985b2986
# ╠═a9368f4d-350e-4978-bec8-1bc873fa32a1
# ╠═d4dca943-58fe-40d1-b23a-eeb3c21a273a
# ╠═651df39c-60ca-4e0d-a6fd-0f13722a7572
# ╠═81a95d09-e0f7-418e-bcf7-dca908c6384c
# ╟─f3232e7a-2e69-4337-9352-4a4a9cc2e7ff
# ╠═a9eca302-7939-4d9b-9b3c-48680c734379
# ╟─e7d1b1be-ebc1-4f62-83a5-4e758047e1f1
# ╟─85104806-f773-4785-885c-9d32cf51b93d
# ╟─d394a4e8-c08c-4086-95a8-39901994b2ca
# ╟─6b5aa387-8f48-45fe-b9be-b0d427e833f1
# ╠═14d058ad-8ef9-4929-8975-c63f5b4b45ac
# ╠═481a94d9-bc4f-4b0b-b808-65b4c80ec0ac
# ╠═bdcce296-bffa-45fe-a658-1116331d9bb6
# ╠═63ccc7f4-8cad-4cb3-b815-deff9f3f9a59
# ╠═6cf3f90d-9c30-4cfe-a3eb-940f51b49dea
# ╟─064994f6-3823-4e43-a71d-00774804fd20
# ╟─718cec78-ed34-455f-9c24-c75b8b7d1f71
# ╟─92ecf5cd-d25d-45f9-95dc-504c354183cf
# ╟─6189f7b2-9331-4f2b-bb58-760e46b03d82
# ╟─a36de957-9ba5-45ad-8c82-3a21d3090fb7
# ╠═c266a1bf-94e6-43fd-ae7f-a5ccdef7cf4b
# ╟─c6fbb9e8-91a3-4c67-981b-40043198a1bf
# ╠═d5bf6a90-e938-43dd-a852-8f5b2f13bb49
# ╠═b2c9ee3b-6805-4bd1-9bee-06714190cbe8
# ╠═200bd89d-a14a-4e4f-ac5c-a0b5e769e2ad
# ╠═ef59ee59-1364-47aa-af21-7742e90a108e
# ╠═6034cc3d-6153-4a71-a200-9f15a643677c
# ╠═9b639977-9310-4682-a5a6-32ea7c4f7ae3
# ╟─bc4461c8-e563-4cf9-a53b-2a70b343efd7

max_stars_repo_path: 1. Root finding.jl

max_issues_repo_path: 1. Root finding.jl

max_forks_repo_path: 1. Root finding.jl


text:
using JuMP, AmplNLWriter, Bonmin_jll, Couenne_jll
using Plots

J = 50;
B = 10.0;
θ = 0.0024;
M = zeros(J).+0.1;
Gt = 1.0;
g = 9.81;
L = 1000.0;
n0 = 0.033;
nt = 0.0133;
Q = 20;
ηu = 0.0;
ηd = 0.0;
f = 1;

I = tan(θ);
Δx = L/(J-1);
x = LinRange(0,L,J);
E = (L.-x)*I;

model = Model(() -> AmplNLWriter.Optimizer(Bonmin_jll.amplexe));
# model = Model(() -> AmplNLWriter.Optimizer(Couenne_jll.amplexe));

# MOI.set(model, MOI.RawParameter("bonmin.algorithm"), "B-BB")
# MOI.set(model, MOI.RawParameter("bonmin.num_resolve_at_root"), 5)
MOI.set(model, MOI.RawParameter("max_iter"), 5000)

@variable(model, 0.0 <= H[1:J] <= Q, start = 2.0);
@variable(model, 0.0 <= A[1:J] <= Q*Q, start = 20.0);
@variable(model, 0.0 <= SM[1:J] <= Q, start = 0.01);
@variable(model, 0.0 <= V[1:J] <= Q, start = 1.0);
# @variable(model, 0.0 <= W[1:J] <= Q, start = 0.01);
@variable(model, z[1:J], Bin);
@NLexpression(model,W[j in 1:J], 5*(V[j]/2.8)^(3/2)/8);

@NLconstraint(model, [j in 1:J-1], V[j]/(2*g) + H[j] + E[j] == V[j+1]/(2*g) + H[j+1] + E[j+1] + SM[j]*Δx + W[j]*z[j]);
@constraint(model, sum(z) <= f);
@constraint(model, z[1] == 0);
@constraint(model, z[J] == 0);
@constraint(model,[j in 1:J], A[j] == B*H[j]);
@NLconstraint(model,[j in 1:J], V[j]*A[j]^2 == Q^2);
@NLexpression(model, R[j in 1:J], A[j]/(B+2*H[j]));
@constraint(model, [j in 1:J], V[j] <= g*H[j]);
@NLexpression(model, S[j in 1:J], V[j]*(n0^2+nt^2*z[j]+2*n0*nt*z[j])/(R[j]^(4/3)));
@NLconstraint(model,[j in 1:J-1], SM[j] == (S[j] + S[j+1])/2);
@constraint(model, [j in 1:J], H[j] >= Gt*z[j]);
@constraint(model, [j in 1:J-1], -M[j]*H[j] <= H[j+1] - H[j]);
@constraint(model, [j in 1:J-1], H[j+1] - H[j] <= M[j]*H[j]);

@NLobjective(model, Max, -H[1]*ηu + H[J]*ηd + sum(W[j]*z[j] for j in 1:J) )

@time status = optimize!(model)

W = value.(model[:W]);
H = value.(model[:H]);
V = value.(model[:V]);
A = value.(model[:A]);
z = value.(model[:z]);
S = value.(model[:S]);
SM = value.(model[:SM]);
U = sqrt.(V)

plot(x,[H,U] ,title=["Altura" "Velocidade"],xlabel="X", ylabel=["H" "U"],layout=(2,1))
# scatter!(x,H)
# scatter!(x,U)
# annotate!(x[4],H[4],"Turbine")
# annotate!(x[4],U[4],"Turbine")
# for i in 1:J
#     if z[i] == 1.0
#         @show i
#         annotate!(x[i],H[i],"Turbine")
#         annotate!(x[i],U[i],"Turbine")
#     end
# end

ext: jl

max_stars_repo_name: Thiago-NovaesB/alpheus

max_issues_repo_name: Thiago-NovaesB/alpheus

max_forks_repo_name: Thiago-NovaesB/alpheus


text:
using DocSeeker
import Markdown

handle("searchdocs") do data
  @destruct [mod || Main, nameOnly || false, exportedOnly || false, allPackages || false, query] = data
  items = @errs DocSeeker.searchdocs(query, mod = mod, exportedonly = exportedOnly,
                                     loaded = !allPackages, name_only = nameOnly)

  if items isa EvalError
    errstr = sprint(showerror, items.err)
    err = startswith(errstr, "Please regenerate the") ?
            """
            Please regenerate the documentation cache with the button above.
            Note that this might take a few minutes. You can track the progress with the progress bar
            in the lower left corner, and should be able to use Juno during normally during that time.
            """ : errstr
    Dict(
      :error => true,
      :errmsg => err
    )
  else
    Dict(
      :error => false,
      :items => [renderitem(i[2]) for i in items], :scores => [i[1] for i in items]
    )
  end
end

function renderitem(x)
  r = Dict(f => getfield(x, f) for f in fieldnames(DocSeeker.DocObj))
  r[:html] = view(renderMD(x.html))

  mod = getmodule(x.mod)
  name = Symbol(x.name)
  r[:typ], r[:icon], r[:nativetype] = if (name !== :ans || mod === Base) && iskeyword(name)
    "keyword", "k", x.typ
  else
    val = getfield′(mod, name)
    # @NOTE: DocSeeker can show docs for non-loaded packages via `createdocsdb()`
    nativetype = val isa Undefined ? "Undefined or not loaded yet" : x.typ
    wstype(mod, name, val), wsicon(mod, name, val), nativetype
  end
  r
end

handle("moduleinfo") do data
  @destruct [mod] = data
  d, items = getmoduleinfo(mod)
  items = [renderitem(i) for i in items]
  Dict(:doc => view(d), :items => items)
end

getmoduleinfo(mod) = ispackage(mod) ? packageinfo(mod) : moduleinfo(mod)
ispackage(mod) = Base.find_package(mod) ≠ nothing

function packageinfo(mod)
  path = DocSeeker.readmepath(mod)
  readme = ispath(path) ? String(read(path)) : ""
  description = Markdown.parse(mod ∈ values(Pkg.Types.stdlib()) ? "## Standard library package `$(mod)`" : readme)

  return  Hiccup.div(
            renderMD(description),
            renderMD("\n---\n## Defined symbols in `$(mod)`:")
          ), modulesymbols(mod)
end

function moduleinfo(mod)
  header = if mod ∈ ("Core", "Base", "Main") || first(split(mod, '.')) == "Base"
    "## Standard module `$(mod)`"
  else
    "## Module `$mod`"
  end * "\n---\n## Defined symbols:" |> renderMD

  header, modulesymbols(mod)
end

function modulesymbols(mod)
  syms = filter(x -> x.mod == mod, DocSeeker.alldocs())
  sort(syms, by = x -> x.name)[1:min(100, length(syms))]
end

using Logging: with_logger
using .Progress: JunoProgressLogger

handle("regeneratedocs") do
  with_logger(JunoProgressLogger()) do
    @errs DocSeeker.createdocsdb()
  end
end


text:
module test_test_logs

using Test
using Logging: Info

@test_logs (Info,)  @info "foo"
@test_logs (:info,) @info "foo"

end # module test_test_logs


text:
@testset "Function testing" begin

@testset ":Options" begin
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;traverse_strategy=:DBFS,obj_epsilon=0.5)...
        )
    )

    v = [10,20,12,23,42]
    w = [12,45,12,22,21]
    @variable(m, x[1:5], Bin)

    @objective(m, Max, dot(v,x))

    @NLconstraint(m, sum(w[i]*x[i]^2 for i=1:5) <= 45)

    MOIU.attach_optimizer(m)
    bm = JuMP.backend(m)
    options = bm.optimizer.model.optimizer.options

    nd_options = Juniper.get_non_default_options(options)
    @test nd_options[:obj_epsilon] == 0.5
    @test nd_options[:traverse_strategy] == :DBFS
    @test length(nd_options[:log_levels]) == 0
end

function option_not_available_t()
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;traverse_strategy=:DBS,obj_epsilon=0.5)...
        )
    )
end

function option_not_available_b()
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;branch_strategy=:Pseudo,obj_epsilon=0.5)...
        )
    )
end

function option_not_available()
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;branch=:Pseudo,obj_epsilon=0.5)...
        )
    )
end

function option_no_mip_solver()
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;branch=:Pseudo,obj_epsilon=0.5,feasibility_pump=true)...
        )
    )

    v = [10,20,12,23,42]
    w = [12,45,12,22,21]
    @variable(m, x[1:5], Bin)

    @objective(m, Max, dot(v,x))

    @NLconstraint(m, sum(w[i]*x[i]^2 for i=1:5) <= 45)
    # there it is detected that no mip solver is present
    optimize!(m)

    bm = JuMP.backend(m)
    return bm.optimizer.model.optimizer.options
end

@testset "Silent/TimeLimitSec" begin
    optimizer = Juniper.Optimizer(;nl_solver=optimizer_with_attributes(Ipopt.Optimizer))
    @test MOI.supports(optimizer, MOI.Silent()) === true
    @test MOI.supports(optimizer, MOI.TimeLimitSec()) === true
    MOI.set(optimizer, MOI.RawOptimizerAttribute("mip_gap"), 1.0)
    # parameter doesn't exist
    @test_logs (:error, r"doesn't exist") MOI.set(optimizer, MOI.RawOptimizerAttribute("mip_gap_1"), 1.0)
    # wrong parameter format
    @test_logs (:error, r"different type") MOI.set(optimizer, MOI.RawOptimizerAttribute("mip_gap"), "abc")
    MOI.set(optimizer, MOI.Silent(), true)
    MOI.set(optimizer, MOI.TimeLimitSec(), nothing)
    @test MOI.get(optimizer, MOI.Silent()) === true
    @test MOI.get(optimizer, MOI.RawOptimizerAttribute("mip_gap")) == 1.0
    @test_logs (:error, r"doesn't exist") MOI.get(optimizer, MOI.RawOptimizerAttribute("mip_gap_1"))
    @test isinf(MOI.get(optimizer, MOI.TimeLimitSec()))
    MOI.set(optimizer, MOI.TimeLimitSec(), 12.0)
    @test MOI.get(optimizer, MOI.TimeLimitSec()) == 12.0
end

@testset ":Option not available" begin
    @test_logs (:error, r"not supported") option_not_available_t()
    @test_logs (:error, r"not supported") option_not_available_b()
    opts = option_no_mip_solver()
    @test opts.feasibility_pump == false
    @test !isa(try option_not_available() catch ex ex end, Exception)
end

@testset "Info/Table" begin
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;branch_strategy=:StrongPseudoCost,processors=2, strong_restart=true)...
        )
    )

    v = [10,20,12,23,42]
    w = [12,45,12,22,21]
    @variable(m, x[1:5], Bin)

    @objective(m, Max, dot(v,x))

    @NLconstraint(m, sum(w[i]*x[i]^2 for i=1:5) <= 45)

    MOIU.attach_optimizer(m)
    bm = JuMP.backend(m)
    JuMP.optimize!(m)
    innermodel = internalmodel(m)

    options = innermodel.options
    jp = innermodel
    println("typeof(m): ", typeof(m))
    println("typeof(jp): ", typeof(jp))

    @test !isa(try Juniper.print_info(jp) catch ex ex end, Exception)
    @test !isa(try Juniper.print_options(jp;all=true) catch ex ex end, Exception)
    @test !isa(try Juniper.print_options(jp;all=false) catch ex ex end, Exception)

    fields, field_chars = Juniper.get_table_config(options)
    ln = Juniper.get_table_header_line(fields, field_chars)
    @test length(ln) == sum(field_chars)

    start_time = time()
    tree = Juniper.init(start_time,jp)
    tree.best_bound = 100
    node = Juniper.new_default_node(1,1,zeros(Int64,5),ones(Int64,5),zeros(Int64,5))
    step_obj = Juniper.new_default_step_obj(jp,node)
    step_obj.counter = 1
    tab_ln, tab_arr = Juniper.get_table_line(2,tree,node,step_obj,start_time,fields,field_chars;last_arr=[])
    @test length(fields) == length(field_chars)
    for i in 1:length(fields)
        fc = field_chars[i]
        f = string(fields[i])
        @test fc >= length(f)
    end
    # Test with incumbent
    tree.incumbent = Juniper.Incumbent(42,[0,0,0,0,1],false)
    tree.limit = :Time
    tab_ln, tab_arr = Juniper.get_table_line(2,tree,node,step_obj,start_time,fields,field_chars;last_arr=[])
    @test length(fields) == length(field_chars)
    for i in 1:length(fields)
        fc = field_chars[i]
        f = string(fields[i])
        @test fc >= length(f)
    end
    @test length(fields) == length(tab_arr)
    @test length(tab_ln) == sum(field_chars)

    # test for diff
    tab_arr_new = copy(tab_arr)
    tab_arr_new[1] = 100
    @test Juniper.is_table_diff(fields, tab_arr, tab_arr_new) == true
    @test Juniper.is_table_diff(fields, tab_arr, tab_arr) == false

    # test if :Time not exists
    i = 1
    for f in fields
        if f == :Time
            deleteat!(fields,i)
            deleteat!(field_chars,i)
        end
        i+= 1
    end
    tab_ln, tab_arr = Juniper.get_table_line(2,tree,node,step_obj,start_time,fields,field_chars;last_arr=[])
    tab_arr_new = copy(tab_arr)
    tab_arr_new[1] = 100
    @test Juniper.is_table_diff(field_chars, tab_arr, tab_arr_new) == true
    @test Juniper.is_table_diff(field_chars, tab_arr, tab_arr) == false

    # produce table line with high counter for nrestarts
    step_obj.counter = 100
    tab_ln, tab_arr = Juniper.get_table_line(2,tree,node,step_obj,start_time,fields,field_chars;last_arr=[])
    idx_restarts = findfirst(fields .== :Restarts)
    @test tab_arr[idx_restarts] == "-"

    # normal gain gap
    step_obj.gain_gap = 0.05
    tab_ln, tab_arr = Juniper.get_table_line(2,tree,node,step_obj,start_time,fields,field_chars;last_arr=[])
    idx_gain_gap = findfirst(fields .== :GainGap)
    @test tab_arr[idx_gain_gap] == "5.0%" || tab_arr[idx_gain_gap] == "5.00%"

    # Inf gain gap
    step_obj.gain_gap = Inf
    tab_ln, tab_arr = Juniper.get_table_line(2,tree,node,step_obj,start_time,fields,field_chars;last_arr=[])
    idx_gain_gap = findfirst(fields .== :GainGap)
    @test tab_arr[idx_gain_gap] == "∞"

    # Huge gain gap
    step_obj.gain_gap = 123456789000
    tab_ln, tab_arr = Juniper.get_table_line(2,tree,node,step_obj,start_time,fields,field_chars;last_arr=[])
    idx_gain_gap = findfirst(fields .== :GainGap)
    @test tab_arr[idx_gain_gap] == ">>"

    # very large gap
    tree.best_bound = 1000000
    tree.incumbent.objval = 1
    tab_ln, tab_arr = Juniper.get_table_line(2,tree,node,step_obj,start_time,fields,field_chars;last_arr=[])
    idx_gap = findfirst(fields .== :Gap)
    @test tab_arr[idx_gap] == ">>"

    # way too long
    fields, field_chars = Juniper.get_table_config(options)
    start_time = time()-123456789
    tab_ln, tab_arr = Juniper.get_table_line(2,tree,node,step_obj,start_time,fields,field_chars;last_arr=[])
    idx_time = findfirst(fields .== :Time)
    @test tab_arr[idx_time] == "t.l."
end

@testset "FP: Table config" begin
    mip_obj,nlp_obj,t, fields, field_chars = 5,2,2, [:MIPobj,:NLPobj,:Time], [20,20,10]
    ln, arr = Juniper.get_fp_table(mip_obj,nlp_obj,t, fields, field_chars, 1e-6)
    @test length(ln) == 50
    @test arr == ["5.0","2.0","2.0"]
end

@testset "Random restarts" begin
    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;branch_strategy=:StrongPseudoCost,processors=2, strong_restart=true)...
        )
    )

    v = [10,20,12,23,42]
    w = [12,45,12,22,21]
    @variable(m, x[1:5], Int)
    JuMP.set_lower_bound(x[1], 0)
    JuMP.set_lower_bound(x[2], 0.5)
    JuMP.set_upper_bound(x[3], 1)
    JuMP.set_upper_bound(x[2], 1)

    @objective(m, Max, dot(v,x))

    @NLconstraint(m, sum(w[i]*x[i]^2 for i=1:5) <= 45)

    MOIU.attach_optimizer(m)
    bm = JuMP.backend(m)
    JuMP.optimize!(m)
    model = internalmodel(m)

    cont_restart = Juniper.generate_random_restart(model)
    @test length(cont_restart) == 5
    @test 0 <= cont_restart[1] <= 20
    @test 0.5 <= cont_restart[2] <= 1
    @test -19 <= cont_restart[3] <= 1
    @test -10 <= cont_restart[4] <= 10
    @test -10 <= cont_restart[5] <= 10

    disc_restart = Juniper.generate_random_restart(model; cont=false)
    @test length(disc_restart) == 5
    for i=1:5
        @test isapprox(round(disc_restart[i])-disc_restart[i],0,atol=1e-6)
    end

    # test a different seed

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;branch_strategy=:StrongPseudoCost,processors=2, strong_restart=true, seed=92301)...
        )
    )

    v = [10,20,12,23,42]
    w = [12,45,12,22,21]
    @variable(m, x[1:5], Int)
    JuMP.set_lower_bound(x[1], 0)
    JuMP.set_lower_bound(x[2], 0.5)
    JuMP.set_upper_bound(x[3], 1)
    JuMP.set_upper_bound(x[2], 1)

    @objective(m, Max, dot(v,x))

    @NLconstraint(m, sum(w[i]*x[i]^2 for i=1:5) <= 45)

    MOIU.attach_optimizer(m)
    bm = JuMP.backend(m)
    JuMP.optimize!(m)
    model = internalmodel(m)

    cont_restart2 = Juniper.generate_random_restart(model)
    # all random numbers should be different then before
    for i=1:5
        @test cont_restart[i] != cont_restart2[i]
    end

    # test a different seed

    m = Model(optimizer_with_attributes(
        Juniper.Optimizer,
        DefaultTestSolver(;branch_strategy=:StrongPseudoCost,processors=2, strong_restart=true, seed=92301)...
        )
    )

    v = [10,20,12,23,42]
    w = [12,45,12,22,21]
    @variable(m, x[1:5], Int)
    JuMP.set_lower_bound(x[1], 0)
    JuMP.set_lower_bound(x[2], 0.5)
    JuMP.set_upper_bound(x[3], 1)
    JuMP.set_upper_bound(x[2], 1)

    @objective(m, Max, dot(v,x))

    @NLconstraint(m, sum(w[i]*x[i]^2 for i=1:5) <= 45)

    MOIU.attach_optimizer(m)
    bm = JuMP.backend(m)
    JuMP.optimize!(m)
    model = internalmodel(m)

    cont_restart3 = Juniper.generate_random_restart(model)
    # all random numbers should be the same as before
    for i=1:5
        @test cont_restart2[i] == cont_restart3[i]
    end
end

end


text:
function Base.show(io::IO, p::Polytope)
    diff = 0
    print(io, rank(p), "-polytope with")
    for (i, j) ∈ enumerate(vcat(start_indices(p), nv(p.graph)))
        n = j - diff
        k = i - 1
        noun = n == 1 && k == 0 ? "vertex" : n == 1 && k == 1 ? "edge" : n == 1 && k == 2 ? "face" : k == 0 ? "vertices" : k == 1 ? "edges" : k == 2 ? "faces" : n == 1 ? string(k, "-face") : string(k, "-faces")
        i ≠ 1 && print(io, ',')
        print(io, ' ', n, ' ', noun)
        diff = j
    end
end

Base.show(io::IO, gp::GeometricPolytope) = print(io, embeddim(gp), "-dimensional geometric ", string(gp.polytope))


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BlockCrossValidation(side)

Cross-validation with blocks of given `side`.

## References

* Roberts et al. 2017. Cross-validation strategies for data with
  temporal, spatial, hierarchical, or phylogenetic structure.
"""
struct BlockCrossValidation{T<:Real} <: AbstractErrorEstimator
  side::T
end

function Base.error(solver::AbstractLearningSolver,
                    problem::LearningProblem,
                    eestimator::BlockCrossValidation)
  sdata = sourcedata(problem)
  ovars = outputvars(task(problem))
  side  = eestimator.side

  blocks    = partition(sdata, BlockPartitioner(side))
  bsubsets  = subsets(blocks)
  neighbors = metadata(blocks)[:neighbors]
  allblocks = 1:length(blocks)

  solutions = pmap(allblocks) do b
    # source and target blocks
    sblocks = setdiff(allblocks, vcat(neighbors[b], b))
    tblocks = [b]

    # source and target indices
    sinds = reduce(vcat, bsubsets[sblocks])
    tinds = reduce(vcat, bsubsets[tblocks])

    # setup and solve learning sub-problem
    subproblem = LearningProblem(view(sdata, sinds),
                                 view(sdata, tinds),
                                 task(problem))
    solve(subproblem, solver)
  end

  result = pmap(ovars) do var
    𝔏 = defaultloss(sdata[1,var])
    losses = map(allblocks) do b
      dview = view(sdata, bsubsets[b])
      ŷ = solutions[b][var]
      y = dview[var]
      𝔏(ŷ, y)
    end
    var => mean(losses)
  end

  Dict(result)
end


text:
using Dates
# using GitCommand
using GitHub
using JSON
using Pkg
using Printf
using RegistryCI
using Test
using TimeZones

const AutoMerge = RegistryCI.AutoMerge

include("automerge-integration-utils.jl")

AUTOMERGE_INTEGRATION_TEST_REPO = ENV["AUTOMERGE_INTEGRATION_TEST_REPO"]::String
TEST_USER_GITHUB_TOKEN = ENV["BCBI_TEST_USER_GITHUB_TOKEN"]::String
INTEGRATION_TEST_READ_ONLY_TOKEN = ENV["INTEGRATION_TEST_READ_ONLY_TOKEN"]::String
GIT = "git"
auth = GitHub.authenticate(TEST_USER_GITHUB_TOKEN)
whoami = RegistryCI.AutoMerge.username(GitHub.DEFAULT_API, auth)
repo_url_without_auth = "https://github.com/$(AUTOMERGE_INTEGRATION_TEST_REPO)"
repo_url_with_auth = "https://$(whoami):$(TEST_USER_GITHUB_TOKEN)@github.com/$(AUTOMERGE_INTEGRATION_TEST_REPO)"
repo = GitHub.repo(AUTOMERGE_INTEGRATION_TEST_REPO; auth=auth)
@test success(`$(GIT) --version`)
@info("Authenticated to GitHub as \"$(whoami)\"")

_delete_branches_older_than = Dates.Hour(3)
delete_old_pull_request_branches(repo_url_with_auth, _delete_branches_older_than; GIT=GIT)

requires_commit = "1c843ed4d4d568345ca557fea48f43efdcd0271f"
hello_world_commit1 = "197e0e03f3f840f830cb5095ab6407d00fbee61c"
hello_world_commit2 = "57b0aec49622faa962c6752d4bc39a62b91fe37c"

@testset "Integration tests" begin
    for (
        test_number,
        (master_dir, feature_dir, public_dir, title, check_license, pass, commit),
    ) in enumerate([
        (
            "master_1",
            "feature_1",
            "",
            "New package: Requires v1.0.0",
            true,
            true,
            requires_commit,
        ), # OK: new package
        (
            "master_1",
            "feature_1",
            "",
            "New package: Requires v1.0.0",
            true,
            false,
            "659e09770ba9fda4a503f8bf281d446c9583ff3b",
        ), # FAIL: wrong commit!
        (
            "master_2",
            "feature_2",
            "",
            "New version: Requires v2.0.0",
            false,
            true,
            requires_commit,
        ),            # OK: new version
        (
            "master_1",
            "feature_3",
            "",
            "New package: Req v1.0.0",
            false,
            false,
            requires_commit,
        ),                # FAIL: name too short
        (
            "master_2",
            "feature_4",
            "",
            "New version: Requires v2.0.1",
            false,
            false,
            requires_commit,
        ),           # FAIL: skips v2.0.0
        (
            "master_3",
            "feature_5",
            "",
            "New version: Requires v2.0.0",
            false,
            false,
            requires_commit,
        ),           # FAIL: modifies extra file
        (
            "master_1",
            "feature_6",
            "",
            "New package: HelloWorldC_jll v1.0.6+0",
            false,
            true,
            hello_world_commit1,
        ),   # OK: new JLL package
        (
            "master_4",
            "feature_7",
            "",
            "New version: HelloWorldC_jll v1.0.8+0",
            false,
            true,
            hello_world_commit2,
        ),   # OK: new JLL version
        (
            "master_1",
            "feature_8",
            "",
            "New package: HelloWorldC_jll v1.0.6+0",
            false,
            false,
            hello_world_commit1,
        ),  # FAIL: unallowed dependency
        (
            "master_1",
            "feature_1",
            "public_1",
            "New package: Requires v1.0.0",
            false,
            true,
            requires_commit,
        ),    # OK: no UUID conflict
        (
            "master_1",
            "feature_1",
            "public_2",
            "New package: Requires v1.0.0",
            false,
            false,
            requires_commit,
        ),  # FAIL: UUID conflict, name differs
        (
            "master_1",
            "feature_1",
            "public_3",
            "New package: Requires v1.0.0",
            false,
            false,
            requires_commit,
        ),  # FAIL: UUID conflict, repo differs
        (
            "master_1",
            "feature_1",
            "public_4",
            "New package: Requires v1.0.0",
            false,
            true,
            requires_commit,
        ),   # OK: UUID conflict but name and repo match
    ])
        @info "Performing integration tests with settings" test_number master_dir feature_dir title pass
        with_master_branch(
            templates(master_dir), "master"; GIT=GIT, repo_url=repo_url_with_auth
        ) do master
            with_feature_branch(
                templates(feature_dir), master; GIT=GIT, repo_url=repo_url_with_auth
            ) do feature
                public_registries = String[]
                if public_dir != ""
                    public_git_repo = generate_public_registry(public_dir, GIT)
                    push!(public_registries, "file://$(public_git_repo)/.git")
                end
                head = feature
                base = master
                body = """
                - Foo: Bar
                - Commit: $commit
                - Hello: World
                """
                params = Dict(
                    "title" => title, "head" => head, "base" => base, "body" => body
                )

                sleep(1)
                pr = GitHub.create_pull_request(repo; auth=auth, params=params)
                pr = wait_pr_compute_mergeability(GitHub.DEFAULT_API, repo, pr; auth=auth)
                @test pr.mergeable
                sleep(1)
                with_pr_merge_commit(pr, repo_url_without_auth; GIT=GIT) do build_dir
                    withenv(
                        "AUTOMERGE_GITHUB_TOKEN" => TEST_USER_GITHUB_TOKEN,
                        "TRAVIS_BRANCH" => master,
                        "TRAVIS_BUILD_DIR" => build_dir,
                        "TRAVIS_EVENT_TYPE" => "pull_request",
                        "TRAVIS_PULL_REQUEST" => string(pr.number),
                        "TRAVIS_PULL_REQUEST_SHA" =>
                            string(AutoMerge.pull_request_head_sha(pr)),
                        "TRAVIS_REPO_SLUG" => AUTOMERGE_INTEGRATION_TEST_REPO,
                    ) do
                        sleep(1)
                        run_thunk =
                            () -> AutoMerge.run(;
                                merge_new_packages=true,
                                merge_new_versions=true,
                                new_package_waiting_period=Minute(typemax(Int32)),
                                new_jll_package_waiting_period=Minute(typemax(Int32)),
                                new_version_waiting_period=Minute(typemax(Int32)),
                                new_jll_version_waiting_period=Minute(typemax(Int32)),
                                registry=AUTOMERGE_INTEGRATION_TEST_REPO,
                                authorized_authors=String[whoami],
                                authorized_authors_special_jll_exceptions=String[whoami],
                                error_exit_if_automerge_not_applicable=true,
                                master_branch=master,
                                master_branch_is_default_branch=false,
                                check_license=check_license,
                                public_registries=public_registries,
                            )
                        @info "Running integration test for " test_number master_dir feature_dir public_dir title pass
                        if pass
                            run_thunk()
                        else
                            @test_throws(
                                RegistryCI.AutoMerge.AutoMergeGuidelinesNotMet, run_thunk()
                            )
                        end
                    end
                    withenv(
                        "AUTOMERGE_GITHUB_TOKEN" => TEST_USER_GITHUB_TOKEN,
                        "TRAVIS_BRANCH" => master,
                        "TRAVIS_BUILD_DIR" => build_dir,
                        "TRAVIS_EVENT_TYPE" => "cron",
                        "TRAVIS_PULL_REQUEST" => "false",
                        "TRAVIS_PULL_REQUEST_SHA" => "",
                        "TRAVIS_REPO_SLUG" => AUTOMERGE_INTEGRATION_TEST_REPO,
                    ) do
                        sleep(1)
                        AutoMerge.run(;
                            merge_new_packages=true,
                            merge_new_versions=true,
                            new_package_waiting_period=Minute(typemax(Int32)),
                            new_jll_package_waiting_period=Minute(typemax(Int32)),
                            new_version_waiting_period=Minute(typemax(Int32)),
                            new_jll_version_waiting_period=Minute(typemax(Int32)),
                            registry=AUTOMERGE_INTEGRATION_TEST_REPO,
                            authorized_authors=String[whoami],
                            authorized_authors_special_jll_exceptions=String[whoami],
                            error_exit_if_automerge_not_applicable=true,
                            master_branch=master,
                            master_branch_is_default_branch=false,
                        )
                        sleep(1)
                        AutoMerge.run(;
                            merge_new_packages=true,
                            merge_new_versions=true,
                            new_package_waiting_period=Minute(0),
                            new_jll_package_waiting_period=Minute(0),
                            new_version_waiting_period=Minute(0),
                            new_jll_version_waiting_period=Minute(0),
                            registry=AUTOMERGE_INTEGRATION_TEST_REPO,
                            authorized_authors=String[whoami],
                            authorized_authors_special_jll_exceptions=String[whoami],
                            error_exit_if_automerge_not_applicable=true,
                            master_branch=master,
                            master_branch_is_default_branch=false,
                        )
                    end
                end
            end
        end
    end
end


text:
# Test @assume and @predict macros on a model without conditioning.

using Turing
using Test

@model test_assume() = begin
  x ~ Bernoulli(1)
  y ~ Bernoulli(x / 2)
  x, y
end

smc = SMC(10000)
pg = PG(10,1000)

res = sample(test_assume(), smc)

@test reduce(&, res[:x]) == 1 # check that x is always 1
@test mean(res[:y]) ≈ 0.5 atol=0.1 # check that the mean of y is between 0.4 and 0.6

res = sample(test_assume(), pg)

@test reduce(&, res[:x]) == 1 # check that x is always 1
@test mean(res[:y]) ≈ 0.5 atol=0.1 # check that the mean of y is between 0.4 and 0.6


text:
using Solaris
using Flux, Statistics
using Flux.Data: DataLoader
using Flux: onehotbatch, onecold, logitcrossentropy, throttle, @epochs
using Base.Iterators: repeated
using MLDatasets

mutable struct Args
    η::Float64      # learning rate
    batchsize::Int   # batch size
    epochs::Int        # number of epochs
    device::Function
end

function prepare_data(args)
    MLDatasets.MNIST.download(i_accept_the_terms_of_use = true)

    # Loading Dataset
    xtrain, ytrain = MLDatasets.MNIST.traindata(Float32)
    xtest, ytest = MLDatasets.MNIST.testdata(Float32)

    # Reshape Data in order to flatten each image into a linear array
    xtrain = Flux.flatten(xtrain)
    xtest = Flux.flatten(xtest)

    # One-hot-encode the labels
    ytrain, ytest = onehotbatch(ytrain, 0:9), onehotbatch(ytest, 0:9)

    # Batching
    train_data = DataLoader(xtrain, ytrain, batchsize = args.batchsize, shuffle = true)
    test_data = DataLoader(xtest, ytest, batchsize = args.batchsize)

    return train_data, test_data
end

function build_model(; imgsize = (28, 28, 1), nclasses = 10)
    return Chain(Affine(prod(imgsize), 32, relu), Affine(32, nclasses))
end

function loss_all(dataloader, model)
    l = 0.0f0
    for (x, y) in dataloader
        l += logitcrossentropy(model(x), y)
    end
    l / length(dataloader)
end

function accuracy(data_loader, model)
    acc = 0
    for (x, y) in data_loader
        acc += sum(onecold(cpu(model(x))) .== onecold(cpu(y))) * 1 / size(x, 2)
    end
    acc / length(data_loader)
end

# ------------------------------------------------------------
# Training
# ------------------------------------------------------------
# initialize parameters
args = Args(1e-3, 1024, 10, cpu)

# load Data
train_data, test_data = prepare_data(args)
train_data = args.device.(train_data)
test_data = args.device.(test_data)

# construct model
m = build_model()
m = args.device(m)

# optimization
loss(x, y) = logitcrossentropy(m(x), y)
evalcb = () -> @show(loss_all(train_data, m))
opt = ADAM()

# training
@epochs args.epochs Flux.train!(loss, params(m), train_data, opt, cb = evalcb)

@show accuracy(train_data, m)
@show accuracy(test_data, m)


text:
module BifurcationInference

	using BifurcationKit: ContIterable, newton, ContinuationPar, NewtonPar, DeflationOperator
	using BifurcationKit: BorderedArray, AbstractLinearSolver, AbstractEigenSolver, BorderingBLS
	using BifurcationKit: ContState, detectBifucation

	using ForwardDiff: Dual,tagtype,derivative,gradient,jacobian
	using Flux: Momentum,update!

	using Setfield: @lens,@set,setproperties
	using Parameters: @unpack

	using InvertedIndices: Not
	using LinearAlgebra, StaticArrays

	include("Structures.jl")
	include("Utils.jl")

	include("Objectives.jl")
	include("Gradients.jl")
	include("Plots.jl")

	export plot,@unpack,BorderedArray,SizedVector
	export StateSpace,deflationContinuation,train!
	export getParameters,loss,∇loss,norm

	""" root finding with newton deflation method"""
	function findRoots!( f::Function, J::Function, roots::AbstractVector{<:AbstractVector},
		parameters::NamedTuple, hyperparameters::ContinuationPar;
		maxRoots::Int = 3, maxIter::Int=500, verbosity=0 )

		hyperparameters = @set hyperparameters.newtonOptions = setproperties(
			hyperparameters.newtonOptions; maxIter = maxIter, verbose = verbosity )

		# search for roots across parameter range
		pRange = range(hyperparameters.pMin,hyperparameters.pMax,length=length(roots))
		roots .= findRoots.( Ref(f), Ref(J), roots, pRange, Ref(parameters), Ref(hyperparameters); maxRoots=maxRoots )
	end

	function findRoots( f::Function, J::Function, roots::AbstractVector{V}, p::T,
		parameters::NamedTuple, hyperparameters::ContinuationPar{T, S, E}; maxRoots::Int = 3, converged = false
		) where { T<:Number, V<:AbstractVector{T}, S<:AbstractLinearSolver, E<:AbstractEigenSolver }

		Zero = zero(first(roots))
		inf = Zero .+ Inf

		# search for roots at specific parameter value
		deflation = DeflationOperator(one(T), dot, one(T), [inf] ) # dummy deflation at infinity
		parameters = @set parameters.p = p

        for u ∈ roots # update existing roots
    		u, residual, converged, niter = newton( f, J, u.+hyperparameters.ds, parameters,
				hyperparameters.newtonOptions, deflation)

				i = 0
				while any(isnan.(residual)) & (i<hyperparameters.newtonOptions.maxIter)
					u .= randn(length(u))

					u, residual, converged, niter = newton( f, J, u.+hyperparameters.ds, parameters,
						hyperparameters.newtonOptions, deflation)
					i += 1
				end

			@assert( !any(isnan.(residual)), "f(u,p) = $(residual[end]) at u = $u, p = $(parameters.p), θ = $(parameters.θ)")
    		if converged push!(deflation,u) else break end
        end

		u = Zero
		if converged || length(deflation)==1 # search for new roots
			while length(deflation)-1 < maxRoots

				u, residual, converged, niter = newton( f, J, u.+hyperparameters.ds, parameters,
					hyperparameters.newtonOptions, deflation)

				# make sure new roots are different from existing
				if any( isapprox.( Ref(u), deflation.roots, atol=2*hyperparameters.ds ) ) break end
				if converged push!(deflation,u) else break end
			end
		end

		filter!( root->root≠inf, deflation.roots ) # remove dummy deflation at infinity
		@assert( length(deflation.roots)>0, "No roots f(u,p)=0 found at p = $(parameters.p), θ = $(parameters.θ); try increasing maxIter")
		return deflation.roots
	end

	""" deflation continuation method """
	function deflationContinuation( f::Function, roots::AbstractVector{<:AbstractVector{V}},
		parameters::NamedTuple, hyperparameters::ContinuationPar{T, S, E};
		maxRoots::Int = 3, maxIter::Int=500, resolution=400, verbosity=0, kwargs...
		) where {T<:Number, V<:AbstractVector{T}, S<:AbstractLinearSolver, E<:AbstractEigenSolver}

		maxIterContinuation,ds = hyperparameters.newtonOptions.maxIter,hyperparameters.ds
		J(u,p) = jacobian(x->f(x,p),u)

		findRoots!( f, J, roots, parameters, hyperparameters; maxRoots=maxRoots, maxIter=maxIter, verbosity=verbosity)
		pRange = range(hyperparameters.pMin,hyperparameters.pMax,length=length(roots))
	    intervals = ([zero(T),step(pRange)],[-step(pRange),zero(T)])

		branches = Vector{Branch{V,T}}()
		hyperparameters = @set hyperparameters.newtonOptions.maxIter = maxIterContinuation
		linsolver = BorderingBLS(hyperparameters.newtonOptions.linsolver)

	    for (i,us) ∈ enumerate(roots)
	        for u ∈ us # perform continuation for each root

	            # forwards and backwards branches
	            for (pMin,pMax) ∈ intervals

					hyperparameters = setproperties(hyperparameters;
						pMin=pRange[i]+pMin, pMax=pRange[i]+pMax,
						ds=sign(hyperparameters.ds)*ds)

	                # main continuation method
					branch = Branch{V,T}()
					parameters = @set parameters.p = pRange[i]+hyperparameters.ds

					try
						iterator = ContIterable( f, J, u, parameters, (@lens _.p), hyperparameters, linsolver, verbosity=verbosity)
						for state ∈ iterator
							push!(branch,state)
						end

						midpoint = sum( s -> s.z.p, branch ) / length(branch)
						if minimum(pRange) < midpoint < maximum(pRange)
							push!(branches,pMin < 0 ? reverse(branch) : branch) end

					catch error
						printstyled(color=:red,"Continuation Error at f(u,p)=$(f(u,parameters))\nu=$u, p=$(parameters.p), θ=$(parameters.θ)\n")
						rethrow(error)
					end
	        		hyperparameters = @set hyperparameters.ds = -hyperparameters.ds
	            end
	    	end
	    end

		hyperparameters = @set hyperparameters.ds = ds
		updateParameters!(hyperparameters,branches;resolution=resolution)
		return unique(branches; atol=10*hyperparameters.ds)
	end
end


text:
import Base: parse

"""
 Parse complex numbers.

 This code doesn't deal with all possible forms of complex numbers, just those output by Mathematica

## Arguments
* `::Type{Complex{T}}`: the type to parse to, where T is a Real number type
* `s::AbstractString`: the string to parse

## Examples
```jldoctest; setup = :(using Polylogarithms)
julia> parse( Complex{Float64}, "1.2 - 3.1*I")
1.2 + 3.1im
```
"""
function parse(::Type{Complex{T}}, s::AbstractString) where {T<:Real}
    s = replace(s, "*^" => "e") # Mathematica seems to export in forms like "3.061616997868383*^-18"
    # real_r = r"([+-]?\d*(\.\d*)?(e[+-]?\d+)?)([\s+-]|$)"
    real_r = r"([+-]?\d*(\.\d*)?(e[+-]?\d+)?)([^*iI])"
    imag_r = r"([+-]?\d*(\.\d*)?(e[+-]?\d+)?)\s*\*\s*I"

    z = zero(Complex{T})

    m = match( imag_r, s )
    if m != nothing
        z += im*parse(T, m[1])
    end

    m = match( real_r, s )
    if m != nothing
        z += parse(T, m[1])
    end
    return z
end

"""
 Parse rational numbers.

    From https://github.com/JuliaLang/julia/issues/18328

## Arguments
* `::Type{Rational{T}}`: the type to parse to, where T is an Integer number type
* `s::AbstractString`: the string to parse

## Examples
```jldoctest; setup = :(using Polylogarithms)
julia> parse( Rational{Int64}, "1 / 2")
1//2
```
"""
function parse(::Type{Rational{T}}, x::AbstractString) where {T<:Integer}
    list = split(x, '/', keepempty=false)
    if length(list) == 1
        return parse(T, list[1]) // 1
    else
        @assert length(list) == 2
        return parse(T, list[1]) // parse(T, list[2])
    end
end


text:
module Corpora

using JSON

struct Corpus
    corpus_path::String
    meta_path::String
    num_documents::Integer
    num_words::Integer
    fields::Array
    index::Array

    function Corpus(path::String, meta_path=nothing)
        corpus_path = path
        base_path = path
        if (endswith(path, ".txt"))
            base_path = path[1:end-4]
        else
            corpus_path = "$path.txt"
        end

        corpus_path = abspath(corpus_path)
        if (!isfile(corpus_path))
            error("Corpus file: $corpus_path not found")
        end

        if (meta_path == nothing)
            meta_path = abspath("$base_path.meta.json")
        else
            meta_path = abspath(meta_path)
        end

        if (!isfile(meta_path))
            error("Metadata file: $meta_path not found")
        end

        meta = JSON.parsefile(meta_path)
        num_documents = meta["num_documents"]
        num_words = meta["num_words"]
        fields = [x for x in meta["fields"] if (x != "line" && x != "byte")]
        index = meta["index"]
        return new(corpus_path, meta_path, num_documents, num_words, fields,
                   index)
    end
end

function get_text(corpus::Corpus, idx::Integer)
    meta = corpus.index[idx]
    open(corpus.corpus_path) do f
        seek(f, meta["byte"])
        return readline(f)
    end
end

function get_text(io::IO, corpus::Corpus, idx::Integer)
    meta = corpus.index[idx]
    seek(io, meta["byte"])
    return readline(io)
end

function get_texts(corpus::Corpus, indices)
    open(corpus.corpus_path) do f
        return [get_text(f, corpus, idx) for idx in indices]
    end
end

end  # end module


text:
#####
### PSG_01 post production of player session gap data
### This produces the data for gap analysis between microsession. The goal is to
### produce more accurate player sessions by throttling the break points that
### player sessionization uses.
#####

using DataFrames, DataArrays

###
# this is the primary data frame - see the sql document in this repository
# for details on how it was made.
###

# df = readtable("C:/Users/Peter.Mancini/Documents/Datafiles/Gap_Analysis.csv")
df = readtable("C:/Users/Peter.Mancini/Documents/Datafiles/Gap_Analysis_Durant.csv")

#####
### This is critical - this is how we smooth out the results. All of the elements
### of the vectors need to be normalized to the same scale. If one is much bigger
### than the others then it will dominate the direction of the vector. Without
### normalized elements, you would just throw the small numbers away and compare
### the large numbers. So if a vector looked like [0.1, 15000, 98], you may as
### well just compare the middle values because the others won't make a difference.
### We normalize by multiplying the the smaller values until their median meets
### the median of the largest value.
#####

medians = [mean(df[:,:cashin]), mean(df[:,:ticketin]), mean(df[:,:bet_freq]), mean(df[:,:avg_bet])]

factors = 1 ./ (medians ./ maximum(medians))

#####
### This is it. Cosine similarity! This is the intelligence. Seems pretty basic!
#####

"""
    cossim(a,b)

compute the cosine similarity between vector a, and vector b.
Returns a Float 1.0 >= cossim >= 0.0
"""
function cossim(a,b)
  dot(a,b)/(vecnorm(a)*vecnorm(b))
end

"""
    session_similarity(x::Int)
Assumes df, the data frame, is global and checks for cosine similarity between
the record at index x and the one prior to it.
"""
function session_similarity(x::Int)
  if !(x>1 && x<=nrow(df))
    println(STDERR, "Index value must be > 1 and less than or equal to ", nrow(df))
    return -1.0
  end
  a = vec(Array(df[x-1,[:cashin, :ticketin, :bet_freq, :avg_bet]])) .* factors
  b = vec(Array(df[x,[:cashin, :ticketin, :bet_freq, :avg_bet]])) .* factors
  return cossim(a,b)
end

"""
     session_similarity(x::UnitRange{Int})
Takes a range and determines session similarity for each item in the range.

# Examples
```julia-repl
julia> session_similarity(510:515)
Float64[6]
0.953…
0.583…
0.666…
0.982…
0.995…
0.998…
```
"""
function session_similarity(x::UnitRange{Int})
  return_array = Array{Float64}(0)
  for y in x
    a = vec(Array(df[y-1,[:cashin, :ticketin, :bet_freq, :avg_bet]])) .* factors
    b = vec(Array(df[y,[:cashin, :ticketin, :bet_freq, :avg_bet]])) .* factors
    push!(return_array, cossim(a,b))
  end
  return return_array
end

"""
    day_hour(unixtime)
Just returns an array of the day of the week and the hour of the day as Ints.
Day 1 is Sunday.
"""
function day_hour(unixtime)
  dt = Dates.unix2datetime(unixtime)
  return ([Dates.dayofweek(dt) Dates.hour(dt)])
end

###
# Start calculationg the gaps, days, and hours for the data frame.
# We create a bunch of holder arrays to keep the values and initialize them.
###

gap = Array{Int32}(nrow(df))
gap[1] = 0
day = Array{Int32}(nrow(df))
day[1] = day_hour(df[1,:mss_begin])[1]
hour = Array{Int32}(nrow(df))
hour[1] = day_hour(df[1,:mss_begin])[2]

###
# this is where we store the similarity values.
###

similarity = Array{Float64}(nrow(df))
similarity[1] = 1.0

###
# this is our current status. These let us know when changes are made and
# what unix time we will use to measure the gap.
###

last_end = df[1,:mss_end]
current_machine = df[1,:machinekey]

###
# This is the main loop that does the work.
###

for i in 2:nrow(df)
  if df[i, :machinekey] != current_machine
    gap[i] = 0
    similarity[i] = 1.0
    day[i] = day_hour(df[i, :mss_begin])[1]
    hour[i] = day_hour(df[i, :mss_begin])[2]
    last_end = df[i,:mss_end]
    current_machine = df[i,:machinekey]
    continue
  end
  gap[i] = df[i,:mss_begin] - last_end
  similarity[i] = session_similarity(i)
  day[i] = day_hour(df[i, :mss_begin])[1]
  hour[i] = day_hour(df[i, :mss_begin])[2]
  if !(gap[i] >= 0)
    println(STDERR, "Error at row ", i, " gap is zero or negative.")
    exit(-1)
  end
  last_end = df[i,:mss_end]
end

###
# Update the data frame with the values we generated. They will then be columns
# we can access like any other column in the dataframe.
###

df[:gap] = gap
df[:similarity] = similarity
df[:day] = day
df[:hour] = hour

###
# this is where we save the results. If we run into a problem or want to do
# a new analysis, you can just reload this file into a df without having
# to re-run this entire program again.
###

writetable("C:/Users/Peter.Mancini/Documents/Datafiles/playersessiongap_post.csv", df)

###
# The next file you should run is DataFrameBinning.jl to get the code you need
# for the real analysis.
###


text:
module TerrainMap

export Terrain_Map, terrain_contour
type Terrain_Map
    north
    north_abs
    east
    east_abs
    down

    min_N
    max_N
    min_E
    max_E

    sig

function Terrain_Map()
    down = -10.0*ones(20,20);
    setO = Set();

    for e = 1:1:20
        for n = 1:1:20
            if(n == 15 && e == 8)
                down[n,e] = down[n,e] + 4;
            elseif(ceil(sqrt((n-15)^2+(e-8)^2)) == 1)
                down[n,e] = down[n,e] + 3;
            elseif(ceil(sqrt((n-15)^2+(e-8)^2)) == 2)
                down[n,e] = down[n,e] + 2;
            elseif(ceil(sqrt((n-15)^2+(e-8)^2)) == 3)
                down[n,e] = down[n,e] + 1;
            end

            if(floor(sqrt((n-8)^2+(e-15)^2))<= 1)
                down[n,e] = down[n,e] + 4;
            elseif(floor(sqrt((n-8)^2+(e-15)^2)) <= 2)
                down[n,e] = down[n,e] + 2;
            end

        end
    end
    this = new();
    this.down = down;
    this.north = [1.0:20.0];
    this.north_abs = this.north;
    this.east = [1.0:20.0];
    this.east_abs = this.east;
    this.min_N = 1.0;
    this.max_N = 20.0;
    this.min_E = 1.0;
    this.max_E = 20.0;
    this.sig = 0.1;
    return this
end

function Terrain_Map(filename)
    down = readdlm(filename,skipstart=2);
    down[down.==0]=NaN;
    (num_north,num_east) = size(down);
    temp = readdlm(filename);
    north = squeeze(float(temp[1,1:num_north]),1);
    east = squeeze(float(temp[2,1:num_east]),1);
    this = new();
    this.down = down;
    this.north = north-minimum(north)+1;
    this.north_abs = north;
    this.east = east-minimum(east)+1;
    this.east_abs = east;
    this.min_N = minimum(north);
    this.min_E = minimum(east);
    this.max_E = maximum(east)
    this.max_N = maximum(north)
    this.sig = 0.1;
    return this;
end
end

function terrain_countour(map)
    min_down = minimum(map.down)-1;
    max_down = maximum(map.down)+1;
    c = contour(map.east,map.north,map.down,25,cmap="ocean",vmin=min_down,vmax=max_down);
    axis("equal")
    return c
end

end # module


text:
# ? activate Splines
import Pkg
Pkg.activate("Splines")

using PlotlyJS
using Revise
Revise.revise()

using Splines

n, m = 8, 12
nodes1 = cat([0:n-1 zeros(n) sin.(range(0, stop=π, length=n))], dims=2)'
nodes2 = cat([zeros(m) 0:m-1 sin.(range(0, stop=π, length=m))], dims=2)'

# create two bezier curves
@info "N=1"
b₁ = Bezier(nodes1)
b₂ = Bezier(nodes2)

# create bezier surface
@info "N=2"
b2n = b₁ ⊗ b₂

@info "plots"
display(
    plot([
        plot_nodes(nodes1, color="red"),
        plot_curve(b₁.coordinates, color="red"),
        plot_nodes(nodes2, color="blue"),
        plot_curve(b₂.coordinates, color="blue"),
        plot_surface(b2n.coordinates, opacity=1)...,
]))

max_stars_repo_name: FedeClaudi/Splines

max_issues_repo_name: FedeClaudi/Splines

max_forks_repo_name: FedeClaudi/Splines


text:
# Test Knitro with MINLPTests
#
using MINLPTests, JuMP, Test

using KNITRO

const SOLVER = JuMP.with_optimizer(KNITRO.Optimizer, outlev=0)

const NLP_SOLVERS = [SOLVER]
const MINLP_SOLVERS = [SOLVER]
const POLY_SOLVERS = []
const MIPOLY_SOLVERS = []

@testset "JuMP Model Tests" begin
    @testset "$(solver.constructor): nlp" for solver in NLP_SOLVERS
        MINLPTests.test_nlp(solver, exclude = [
            "005_011",  # Uses the function `\`
            "008_011"   # Requires quadratic duals
        ])
        MINLPTests.test_nlp_cvx(solver)
    end
    @testset "$(solver.constructor): nlp_mi" for solver in MINLP_SOLVERS
        MINLPTests.test_nlp_mi(solver, exclude = [
            "005_011",  # Uses the function `\`
            "003_014",  # Bug in Knitro
            "003_015",  # Bug in Knitro
            "006_010"   # Bug in Juniper - handling of user-defined functions.
        ])
        MINLPTests.test_nlp_mi_cvx(solver)
    end
end


text:
# Definindo o conjunto de dados:
function get_protein_data(path) # Função que lê os arquivos e transforma em vetores com os caracteres.

    listfile = open(path*"/lista-de-proteinas.txt", "r")
    protein_list = readlines(listfile)
    close(listfile)

    #fastafile = open(path*"/proteinas.fa", "r")
    #fasta = readlines(fastafile)
    #close(file)

    Seq = Vector{Vector{Char}}(undef, 0) # Lista de Resíduos
    Str = Vector{Vector{Char}}(undef, 0) # Lista da Estrutura Secundária de Cada Resíduo

    println("Reading protein data files")
    @showprogress for i in 1:length(protein_list)

        # DSSP File
        dsspfile = open(path * "/DSSP/" * protein_list[i] * ".dssp")
        dssp_file = readlines(dsspfile)[29:end]
        close(dsspfile)

        # SEQUENCE
        x = Vector{Char}(undef, length(dssp_file))
        for j in 1:length(x)
            x[j] = dssp_file[j][14]
        end

        # STRUCTURE
        y = Vector{Char}(undef, length(dssp_file))
        for j in 1:length(y)
            y[j] = dssp_file[j][17]
        end

        bri = vcat(0, findall(f -> f == '!', x), length(x))
        for idx in 2:length(bri)
            x_vec = copy(x[bri[idx-1]+1:bri[idx]-1])
            y_vec = copy(y[bri[idx-1]+1:bri[idx]-1])
            Seq = vcat([x_vec], Seq)
            Str = vcat([y_vec], Str)
        end
    end
    println("End of the data reading")

    return Seq, Str
end


text:
# PassiveBranch
abstract type PassiveBranch{DF<:DeviceFormulation} <: Branch{DF} end

cost(d::PassiveBranch) = sum(d[:capital_cost] .* (d[:s_nom] .- getparam(d, :s_nom)))

# Split lines into AC and DC?
@adddevice(Line, PassiveBranch, LinearExpansionForm{LinearDispatchForm},
           :lines, (:L, :T=>:snapshots), joinpath(@__DIR__, "attrs", "lines.csv"))

@adddevice(Transformer, PassiveBranch, LinearExpansionForm{LinearDispatchForm},
           :transformers, (:L, :T=>:snapshots), joinpath(@__DIR__, "attrs", "transformers.csv"))

function addto!(jm::ModelView, m::EnergyModel, d::PassiveBranch{DF}) where
      {DDF, DF <: LinearExpansionForm{DDF}}

    addto!(jm, m, with_formulation(d, LinearExpansionInvestmentForm))
    addto!(jm, m, with_formulation(d, LinearExpansionDispatchForm{DDF}))
end

function addto!(jm::ModelView, m::EnergyModel, d::PassiveBranch{DF}) where
      {DF <: LinearExpansionInvestmentForm}

    L = axis(m, d)

    s_nom_min = get(d, :s_nom_min, L)
    s_nom_max = get(d, :s_nom_max, L)

    @variable(jm, s_nom_min[l] <= s_nom[l=L] <= s_nom_max[l])
end

function addto!(jm::ModelView, m::EnergyModel, d::PassiveBranch{DF}) where
      {DDF, DF <: LinearExpansionDispatchForm{DDF}}

    L = axis(m, d)
    T = axis(m, :snapshots)

    s_nom = get(d, :s_nom, L)
    s_max_pu = get(d, :s_max_pu, L, T)

    @variable(jm, p[L,T])

    @constraints jm begin
        p_upper[l=L,t=T], p[l,t] <= s_max_pu[l,t] * s_nom[l]
        p_lower[l=L,t=T], p[l,t] >= - s_max_pu[l,t] * s_nom[l]
    end
end

function addto!(jm::ModelView, m::EnergyModel, d::PassiveBranch{DF}) where {DF <: DispatchForm}

    L = axis(m, d)
    T = axis(m, :snapshots)

    s_nom = get(d, :s_nom, L)
    s_max_pu = get(d, :s_max_pu, L, T)

    @variable(jm, - s_max_pu[l,t] * s_nom[l] <= p[l=L,t=T] <= s_max_pu[l,t] * s_nom[l])
end

## Lines and Transformers
effectiveimpedance(d::PassiveBranch) = impedance(d)[d[:carrier] == :DC ? :r_pu_eff : :x_pu_eff]

phaseshift(d::PassiveBranch) = 0.
phaseshift(d::Transformer) = gettypeparams(d.model.data, d)[:phase_shift]

function impedance(d::Line)
    L = axis(d)
    bus_v_nom = d.model[Bus][:v_nom]
    v_nom = bus_v_nom[convert(Array{Int64}, indexin(get(d, :bus0, L), first(axisvalues(bus_v_nom))))]

    p = gettypeparams(d.model.data, d)
    if p !== nothing
        num_parallel = get(d, :num_parallel, L)
        length = get(d, :length, L)

        # line types
        x = AxisArray(p[:x_per_length] .* length ./ num_parallel, L)
        r = AxisArray(p[:r_per_length] .* length ./ num_parallel, L)
        # b = AxisArray(2pi*1e-9*p[:f_nom]*p[:d_per_length] .* length .* num_parallel, ax)
    else
        x = get(d, :x, L)
        r = get(d, :r, L)
        # b = d[:b]
    end

    x_pu = AxisArray(x ./ v_nom.^2, L)
    r_pu = AxisArray(r ./ v_nom.^2, L)
    # b_pu = AxisArray(r .* v_nom.^2, ax)

    Dict(:x => x, :r => r,
         # :b => b, :g => g,
         :x_pu => x_pu, :x_pu_eff => x_pu,
         :r_pu => r_pu, :r_pu_eff => r_pu,
         # :b_pu => b_pu, :g_pu => 0
         )
end

function impedance(d::Transformer)
    ax = axis(d)
    p = gettypeparams(d.model.data, d)

    if p !== nothing
        num_parallel = d[:num_parallel]
        s_nom = p[:s_nom]

        r = p[:vscr]/100 ./ num_parallel
        x = sqrt(p[:vsc]^2 - p[:vscr]^2)/100 ./ num_parallel
        g = p[:pfe]/(1000*p[:s_nom]) .* num_parallel
        b = (-sqrt(clamp((p[:i0])^2 - (p[:pfe]/(10*p[:s_nom]))^2, 0., Inf)) / 100
             .* num_parallel)
    else
        s_nom = d[:s_nom]

        x = d[:x]
        r = d[:r]
        g = d[:g]
        b = d[:b]
    end

    r_pu = AxisArray(r ./ s_nom, ax)
    x_pu = AxisArray(x ./ s_nom, ax)
    g_pu = AxisArray(g .* s_nom, ax)
    b_pu = AxisArray(b .* s_nom, ax)

    """Follows http://home.earthlink.net/~w6rmk/math/wyedelta.htm"""
    function wye_to_delta(z1,z2,z3)
        summand = z1*z2 + z2*z3 + z3*z1
        (summand/z2,summand/z1,summand/z3)
    end

    """Convert given T-model parameters to PI-model parameters using wye-delta transformation"""
    function apply_transformer_t_model!(r_pu, x_pu, g_pu, b_pu)
        for i in eachindex(r_pu, x_pu, g_pu, b_pu)
            z_series = r_pu[i] + im*x_pu[i]
            y_shunt  = g_pu[i] + im*b_pu[i]

            if y_shunt == 0 continue end

            za, zb, zc = wye_to_delta(z_series/2, z_series/2, 1/y_shunt)

            r_pu[i] = real(zc)
            x_pu[i] = imag(zc)
            g_pu[i] = real(2/za)
            b_pu[i] = imag(2/za)
        end
    end

    if d[:model] == "t"
        apply_transformer_t_model!(r_pu, x_pu, g_pu, b_pu)
    end

    tap_ratio = AxisArray(1 .+ (d[:tap_position] .- p[:tap_neutral]) .* (p[:tap_step]/100.), ax)
    Dict(:r_pu => r_pu, :r_pu_eff => AxisArray(r_pu .* tap_ratio, ax),
         :x_pu => x_pu, :x_pu_eff => AxisArray(x_pu .* tap_ratio, ax),
         :g_pu => g_pu, :b_pu => b_pu,
         :tap_ratio => tap_ratio,
         :phase_shift => p[:phase_shift],
         :s_nom => p[:s_nom])
end

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
module CUDNN

using ..CUDA
import ..CUDA: ndevices, getdevice
import Libdl

if Sys.iswindows()
    const libcudnn = Libdl.find_library(["cudnn64_7"])
else
    const libcudnn = Libdl.find_library("libcudnn")
end
isempty(libcudnn) && error("CUDNN cannot be found.")

const API_VERSION = Ref{Int}()

function __init__()
    API_VERSION[] = Int(ccall((:cudnnGetVersion,libcudnn),Cint,()))
    @info "CUDNN API $(API_VERSION[])"
end

macro cudnn(f, args...)
    quote
        status = ccall(($f,libcudnn), Cint, $(map(esc,args)...))
        if status != 0
            p = ccall((:cudnnGetErrorString,libcudnn), Ptr{UInt8}, (Cint,), status)
            throw(unsafe_string(p))
        end
    end
end

include("define.jl")
include("handle.jl")

function setstream(handle::Handle, stream)
    @cudnn :cudnnSetStream (Ptr{Cvoid},Ptr{Cvoid}) handle stream
end

const ALLOCATED = []

include("activation.jl")
include("add.jl")
include("convolution.jl")
include("filter.jl")
include("dropout.jl")
include("reduce.jl")
include("rnn.jl")
include("softmax.jl")
include("tensor.jl")

end


text:
using Pluto
notebook_path = joinpath(@__DIR__, "scenes.jl")
isfile(notebook_path)

options = Pluto.Configuration.from_flat_kwargs(; project=@__DIR__)
s = Pluto.ServerSession()
nb = Pluto.SessionActions.open(s, notebook_path; run_async=false)

open(joinpath(@__DIR__, "tutorials", "scenes.md"), "w") do io
    for (i, c) in enumerate(nb.cells)
        if startswith(c.code, "md\"\"\"")
            md = join(split(c.code, "\n")[2:end-1], "\n")
            println(io, md)
        else
            println(io, """
                       \\begin{examplefigure}{}
                        ```julia
                        """)
            println(io, c.code)
            println(io, """
                        ```
                        \\end{examplefigure}
                        """)
        end
    end
end


text:
clear
close all

test_case = 'test21'
dtoutput = 100;  # Seconds
modeloutputpathandfilename = [test_case '.mat']
CMM_driver(modeloutputpathandfilename,dtoutput)
CMM_graphics(modeloutputpathandfilename,10,test_case,dtoutput)

CMM_graphics(modeloutputpathandfilename,1,test_case,dtoutput)
CMM_graphics(modeloutputpathandfilename,3,test_case,dtoutput)
CMM_graphics(modeloutputpathandfilename,5,test_case,dtoutput)
CMM_graphics(modeloutputpathandfilename,6,test_case,dtoutput)
CMM_graphics(modeloutputpathandfilename,7,test_case,dtoutput)
CMM_graphics(modeloutputpathandfilename,8,test_case,dtoutput)

##
title_in = 'Perturbation Potential Temperature (K) at t = '
title_in = 'Pressure Perturbation (Pa) at t = '
title_in = 'W-wind (m/s) at t ='
title_in = 'U-wind (m/s) at t = '
title_in = 'Water Vapor Mixing Ratio (kg/kg) at t = '
title_in = 'Cloud Water Mixing Ratio (kg/kg) at t = '
title_in = 'Rain Water Mixing Ratio (kg/kg) at t = '

b_mat = [900 2400 3600];
for b=1:3
var_input = 'qrain'
test_case = 'test21'
modeloutputpathandfilename = [test_case '.mat']
n = b_mat(b)
CMM_graphics_v3(modeloutputpathandfilename,title_in,var_input,test_case,n)
end

clear
load test15.mat
x_interval = dx/1000;
x_end = (dx/1000)*(nx-1);
y_begin = dz/2000;
y_interval = dz/1000;
y_end = (dz*nz/1000)-(dz/2000);
x=0.0:x_interval:x_end;
y=y_begin:y_interval:y_end;
u_qv = fliplr(rot90(u_3600,3));
w_qv = fliplr(rot90(w_3600,3));
th_tot=zeros(nx,nz);
for k = 1 : nz
    for i = 1 : nx

th_tot(i,k) = tb(k) + th_3600(i,k);
    end
end

for i=1:3:nx-3
        u_qv(:,i) = 0;
        u_qv(:,i+2) = 0;
        w_qv(:,i) = 0;
        w_qv(:,i+2) = 0;
end

for i=1:2:nz-2
        u_qv(i,:) = 0;
        u_qv(i+2,:) = 0;
        w_qv(i,:) = 0;
        w_qv(i+2,:) = 0;
end

mask=zeros(nx,nz);
mask(260:340,1:3) = 5000;
mask(270:330,4:7) = 5000;
mask(280:320,8:11) = 5000;
mask(290:310,12:14) = 5000;

f = figure('Visible','on');
  [C, h]= contourf(x,y,fliplr(rot90(pprt_3600,3)),[-500:5:500]);
  caxis([-500 500])
  colorbar
  hold on
  quiver(x,y,u_qv,w_qv,3,'k')
  [C2,h2] = contourf(x,y,fliplr(rot90(mask,3)),[0 5000]);
  ch = get(h2,'child');alpha(ch,.9);
  grid on
  xlabel('x (km)')
  ylabel('z (km)')
  clabel(C,h)
  title_str = ['Wind Vectors and Perturbation Pressure (Pa) at t = 3600 s'];
  title(title_str)
  hold off
  str_print= ['print ' '-dpng ' '-f  forpaper_pressureandwinds_case1_3600s.png'];
  eval(str_print);
  close

##
# ## Quiver Plots
# clear
# load kessler.mat
# x_interval = dx/1000;
# x_end = (dx/1000)*(nx-1);
# y_begin = dz/2000;
# y_interval = dz/1000;
# y_end = (dz*nz/1000)-(dz/2000);
# x=0.0:x_interval:x_end;
# y=y_begin:y_interval:y_end;
#
# u_qv = fliplr(rot90(alltimes_u(:,:,670),3));
# w_qv = fliplr(rot90(alltimes_w(:,:,670),3));
#
# for i=1:3:nx-3
#         u_qv(:,i) = 0;
#         u_qv(:,i+2) = 0;
#         w_qv(:,i) = 0;
#         w_qv(:,i+2) = 0;
# end
#
# for i=1:2:nz-2
#         u_qv(i,:) = 0;
#         u_qv(i+2,:) = 0;
#         w_qv(i,:) = 0;
#         w_qv(i+2,:) = 0;
# end
#
#
# f = figure('Visible','on');
#   [C, h]= contour(x,y,fliplr(rot90(alltimes_th(:,:,670),3)),[-10:.25:10]);
#   colorbar
#   hold on
#   quiver(x,y,u_qv,w_qv,3)
#   grid on
#   xlabel('x (km)')
#   ylabel('z (km)')
#   clabel(C,h)
#   title_str = ['Perturbation Potential Temperature (K) at t = 670 s'];
#   title(title_str)
#   hold off
#   str_print= ['print ' '-dpng ' '-f' ' ' ...
#   nametosave '_theta_0secs.png'];
#   eval(str_print);
#   close
#
#   ## Hovmoller Plots
#
#  for n = 1:1200
#    reorg(:,n) = alltimes_th(:,4,n);
#  end
#
#   figure
#   contour(xs/1000,t(1:1200),flipud(rot90(reorg)))
#   colorbar
#
#
# ## Filled Contour Plots
# clear
# load kessler_test.mat
# x_interval = dx/1000;
# x_end = (dx/1000)*(nx-1);
# y_begin = dz/2000;
# y_interval = dz/1000;
# y_end = (dz*nz/1000)-(dz/2000);
# x=0.0:x_interval:x_end;
# y=y_begin:y_interval:y_end;
#
# f = figure('Visible','on');
#   [C, h]= contourf(x,y,fliplr(rot90(alltimes_th(:,:,1500),3)),[-18:1:18]);
#   colorbar
#   grid on
#   xlabel('x (km)')
#   ylabel('z (km)')
#
#
#
# f = figure('Visible','on');
#   [C, h]= contourf(x,y,fliplr(rot90(alltimes_qrain(:,:,1200),3)),[0.0:.001:015]);
#   colorbar
#   grid on
#   xlabel('x (km)')
#   ylabel('z (km)')
#
#   clabel(C,h)
#   title_str = ['Perturbation Potential Temperature (K) at t = 670 s'];
#   title(title_str)
#
#   str_print= ['print ' '-dpng ' '-f' ' ' ...
#   nametosave '_theta_0secs.png'];
#   eval(str_print);
#   close
#
#
#

max_stars_repo_name: jeffreycunn/toy-atmospheric-model

max_issues_repo_name: jeffreycunn/toy-atmospheric-model

max_forks_repo_name: jeffreycunn/toy-atmospheric-model


text:
function _build_disjunction()
    form = Cerberus.DMIPFormulation()
    Cerberus.add_variable(form)
    Cerberus.add_variable(form)
    Cerberus.set_bounds!(form, _CVI(1), _IN(-1.1, 1.1))
    Cerberus.set_bounds!(form, _CVI(2), _IN(-1.1, 1.1))
    x = [_SV(_VI(1)), _SV(_VI(2))]
    Cerberus.add_constraint(
        form,
        Cerberus.AffineConstraint(
            _CSAF([1.0, 1.0], [_CVI(1), _CVI(2)], 0.0),
            _LT(0.5),
        ),
    )
    Cerberus.add_constraint(
        form,
        Cerberus.AffineConstraint(
            _CSAF([1.0, -1.0], [_CVI(1), _CVI(2)], 0.0),
            _LT(0.6),
        ),
    )

    f_1 = 1.0 * x[1] + 1.0 * x[2]
    f_2 = 1.0 * x[1] - 1.0 * x[2]
    f_3 = 1.0 * x[1] + 0.5 * x[2]
    f_4 = 1.0 * x[1] - 0.5 * x[2]
    f = [convert(_CSAF, _f) for _f in (f_1, f_2, f_3, f_4)]

    lbs = [
        -Inf -Inf -Inf
        -Inf -Inf -Inf
        -Inf -Inf 0.5
        -Inf -Inf 0.5
    ]
    ubs = [
        0.0 Inf Inf
        Inf 0.0 Inf
        Inf Inf Inf
        Inf Inf Inf
    ]
    s = DisjunctiveConstraints.DisjunctiveSet(lbs, ubs)
    return form, Cerberus.Disjunction(f, s)
end

@testset "NaiveBigMFormulater" begin
    form, disjunction = _build_disjunction()
    formulater = Cerberus.DisjunctiveFormulater(
        disjunction,
        DisjunctiveConstraints.NaiveBigM(
            DisjunctiveConstraints.IntervalArithmetic(),
        ),
    )
    Cerberus.attach_formulater!(form, formulater)

    @testset "new_variables_to_attach" begin
        raw_indices = Cerberus.new_variables_to_attach(formulater)
        @test raw_indices == [_ZO(), _ZO(), _ZO()]
        for var_kind in raw_indices
            Cerberus.add_variable(form, var_kind)
        end
    end
    @testset "compute_disjunction_activity" begin
        let node = Cerberus.Node()
            activity = Cerberus.compute_disjunction_activity(
                form,
                formulater,
                node,
                CONFIG,
            )
            activity == [true, true, true]
        end
        let node = Cerberus.Node(
                [Cerberus.BoundUpdate(_CVI(5), _LT(0.0))],
                [Cerberus.BoundUpdate(_CVI(3), _GT(1.0))],
                2,
            )
            activity = Cerberus.compute_disjunction_activity(
                form,
                formulater,
                node,
                CONFIG,
            )
            @test activity == [true, false, false]
        end
        let node = Cerberus.Node(
                Cerberus.BoundUpdate{_LT}[],
                [
                    Cerberus.BoundUpdate(_CVI(3), _GT(1.0)),
                    Cerberus.BoundUpdate(_CVI(4), _GT(1.0)),
                ],
                2,
            )
            activity = Cerberus.compute_disjunction_activity(
                form,
                formulater,
                node,
                CONFIG,
            )
            @test activity == [false, false, false]
        end
    end
end


text:
# function MTiXsens_FillMissing(data)
# #fill MTiXsens missing datas with linear interpolation
#
#     # MTiXsens_temp=MTiXsens[filenameMTiXsens[1,1]][sensorXsensNames[1]][:,[1,6,7,8]]
#     MTiXsens_temp=data[:,:]
#     #search missing packet
#     MTiXsens_temp_missing=MTiXsens_temp.PacketCounter[2:end]-MTiXsens_temp.PacketCounter[1:end-1].-1
#
#     # as the Packcounter is coded on 16bits, the max value is 65535 ans return to 0 when reaching over
#     # we have to add 2^16=65536 to all PacketCounter after each 0 to code the PacketCounter with more bits
#     negative_row=findall(MTiXsens_temp_missing.<0)
#     if size(negative_row,1)!=0
#         for negative in negative_row
#             MTiXsens_temp.PacketCounter[negative+1:end]=MTiXsens_temp.PacketCounter[negative+1:end].+2^16
#         end
#         MTiXsens_temp_missing=MTiXsens_temp.PacketCounter[2:end]-MTiXsens_temp.PacketCounter[1:end-1].-1
#     end
#
#     #Loop to fill missing datas
#     while isnothing(findfirst(MTiXsens_temp_missing.!=0))==false
#         row_missing=findfirst(MTiXsens_temp_missing.!=0)
#         nb_missing=MTiXsens_temp_missing[row_missing]
#
#         #save temporary datas before and after missing datas
#         MTiXsens_temp_down=MTiXsens_temp[row_missing:row_missing,:]
#         MTiXsens_temp_up=MTiXsens_temp[row_missing+1:row_missing+1,:]
#
#         #linear coefficient to interpolate missing datas
#         MTiXsens_temp_frac=(MTiXsens_temp_up.-MTiXsens_temp_down)./(nb_missing+1)
#
#         #add missing datas row by row with interpolation
#         for k=1:nb_missing
#             MTiXsens_temp_final=MTiXsens_temp_down.+(MTiXsens_temp_frac).*k
#             MTiXsens_temp_final.PacketCounter=MTiXsens_temp_down.PacketCounter.+k
#
#             insert!.(eachcol(MTiXsens_temp),row_missing+k,
#                 [MTiXsens_temp_final[1,1],
#                 MTiXsens_temp_final[1,2],
#                 MTiXsens_temp_final[1,3],
#                 MTiXsens_temp_final[1,4]]
#                 );
#         end
#
#         # insert!.(eachcol(MTiXsens_temp),2,[5,6,7,8,9]);
#
#         MTiXsens_temp_missing=MTiXsens_temp.PacketCounter[2:end]-MTiXsens_temp.PacketCounter[1:end-1].-1
#     end
#
#     return MTiXsens_temp
#
# end

function MTiXsens_FillMissing(data)
#fill MTiXsens missing datas with linear interpolation

    # MTiXsens_temp=MTiXsens[filenameMTiXsens[1,1]][sensorXsensNames[1]][:,[1,6,7,8]]
    MTiXsens_temp=data[:,:]
    #search missing packet
    MTiXsens_temp_missing=MTiXsens_temp.PacketCounter[2:end]-MTiXsens_temp.PacketCounter[1:end-1].-1

    # as the Packcounter is coded on 16bits, the max value is 65535 ans return to 0 when reaching over
    # we have to add 2^16=65536 to all PacketCounter after each 0 to code the PacketCounter with more bits
    negative_row=findall(MTiXsens_temp_missing.<0)
    if size(negative_row,1)!=0
        for negative in negative_row
            MTiXsens_temp.PacketCounter[negative+1:end]=MTiXsens_temp.PacketCounter[negative+1:end].+2^16
        end
        MTiXsens_temp_missing=MTiXsens_temp.PacketCounter[2:end]-MTiXsens_temp.PacketCounter[1:end-1].-1
    end

    #Loop to fill missing datas
    while isnothing(findfirst(MTiXsens_temp_missing.!=0))==false
        row_missing=findfirst(MTiXsens_temp_missing.!=0)
        nb_missing=MTiXsens_temp_missing[row_missing]

        #save temporary datas before and after missing datas
        MTiXsens_temp_down=MTiXsens_temp[row_missing:row_missing,:]
        if row_missing+1==size(MTiXsens_temp,1)
            MTiXsens_temp_up=MTiXsens_temp[row_missing+1:row_missing+1,:]
        else
            MTiXsens_temp_up=MTiXsens_temp[row_missing+2:row_missing+2,:]
        end

        #linear coefficient to interpolate missing datas
        MTiXsens_temp_frac=(MTiXsens_temp_up.-MTiXsens_temp_down)./(nb_missing+2)

        #add missing datas row by row with interpolation
        for k=1:nb_missing
            MTiXsens_temp_final=MTiXsens_temp_down.+(MTiXsens_temp_frac).*k
            MTiXsens_temp_final.PacketCounter=MTiXsens_temp_down.PacketCounter.+k

            insert!.(eachcol(MTiXsens_temp),row_missing+k,
                convert(Array,MTiXsens_temp_final[1,:])
                # [MTiXsens_temp_final[1,1],
                # MTiXsens_temp_final[1,2]]
                );
        end
        MTiXsens_temp_final=MTiXsens_temp_down.+(MTiXsens_temp_frac).*(nb_missing+1)
        MTiXsens_temp_final.PacketCounter=MTiXsens_temp_down.PacketCounter.+(nb_missing+1)
        MTiXsens_temp[row_missing+nb_missing+1,:]=convert(Array,MTiXsens_temp_final[1,:])

        MTiXsens_temp_missing=MTiXsens_temp.PacketCounter[2:end]-MTiXsens_temp.PacketCounter[1:end-1].-1
    end

    return MTiXsens_temp

end


text:
module TestDeparse
    using Calculus2
    using Base.Test

    @test isequal(
        deparse(:(cos(x) + sin(x))),
        "cos(x) + sin(x)"
    )
    @test isequal(
        deparse(:(cos(x) + sin(x) + exp(-x))),
        "cos(x) + sin(x) + exp(-x)"
    )
    @test isequal(
        deparse(parse("x+y*z")),
        "x + y * z"
    )
    @test isequal(
        deparse(parse("(x+y)*z")),
        "(x + y) * z"
    )
    @test isequal(
        deparse(parse("1/(x/y)")),
        "1 / (x / y)"
    )
    @test isequal(
        deparse(parse("1/(x*y)")),
        "1 / (x * y)"
    )
    @test isequal(
        deparse(parse("z^(x+y)")),
        "z ^ (x + y)"
    )
    @test isequal(
        deparse(parse("z^(x*y)")),
        "z ^ (x * y)"
    )
end


text:
# Empty hcat defaults to row vector
hcat() = [].'
typed_hcat{T}(::Type{T}) = typed_vcat(T).'

hcat{T}(X::T...)         = T[ X[j] for j=1:length(X) ].'
hcat{T<:Number}(X::T...) = T[ X[j] for j=1:length(X) ].'

hcat(X::Number...) = hvcat_fill(Array{promote_typeof(X...),1}(length(X)), X).'
typed_hcat{T}(::Type{T}, X::Number...) = hvcat_fill(Array{T,1}(length(X)), X).'

hcat(X::TransposedVector...) = vcat(map(transpose, X)...).'
hcat(X::Union{TransposedVector,Number}...) = vcat(map(transpose, X)...).'


text:
using Knockoffs
using Test
using LinearAlgebra
using Random
using StatsBase
using Statistics

@testset "fixed equi knockoffs" begin
    Random.seed!(2021)

    # simulate matrix and normalize columns
    n = 3000
    p = 1000
    X = randn(n, p)
    zscore!(X, mean(X, dims=1), std(X, dims=1)) # center/scale Xj to mean 0 var 1
    normalize_col!(X) # normalize Xj to unit length

    # equi-correlated knockoff
    @time knockoff = fixed_knockoffs(X, :equi)
    X = knockoff.X
    X̃ = knockoff.X̃
    s = knockoff.s
    Σ = knockoff.Σ
    Σinv = knockoff.Σinv

    @test all(X' * X .≈ Σ)
    @test all(isapprox.(X̃' * X̃, Σ, atol=5e-2)) # numerical accuracy not good?
    @test all(s .≥ 0)
    @test all(1 .≥ s)
    # λ = eigvals(2Σ - Diagonal(s))
    # for λi in λ
    #     @test λi ≥ 0 || λi ≈ 0
    # end
    # @test all(isapprox.(Ũ' * X, 0, atol=1e-10))
    for i in 1:p, j in 1:p
        if i == j
            @test isapprox(dot(X[:, i], X̃[:, i]), Σ[i, i] - s[i])
            @test isapprox(dot(X[:, i], X̃[:, i]), 1 - s[i], atol=5e-2) # numerical accuracy not good?
        else
            @test dot(X[:, i], X̃[:, j]) ≈ dot(X[:, i], X[:, j])
        end
    end
end

@testset "fixed SDP knockoffs" begin
    Random.seed!(2021)

    # simulate matrix and normalize columns
    n = 1000
    p = 100
    X = randn(n, p)
    zscore!(X, mean(X, dims=1), std(X, dims=1)) # center/scale Xj to mean 0 var 1
    normalize_col!(X) # normalize Xj to unit length

    # SDP knockoff
    @time knockoff = fixed_knockoffs(X, :sdp)
    X = knockoff.X
    X̃ = knockoff.X̃
    s = knockoff.s
    Σ = knockoff.Σ
    Σinv = knockoff.Σinv

    @test all(X' * X .≈ Σ)
    @test all(isapprox.(X̃' * X̃, Σ, atol=5e-1)) # numerical accuracy not good?
    @test all(s .≥ 0)
    @test all(1 .≥ s)
    # λ = eigvals(2Σ - Diagonal(s))
    # for λi in λ
    #     @test λi ≥ 0 || λi ≈ 0
    # end
    # @test all(isapprox.(Ũ' * X, 0, atol=1e-10))
    for i in 1:p, j in 1:p
        if i == j
            @test isapprox(dot(X[:, i], X̃[:, i]), Σ[i, i] - s[i], atol=1e-8)
            @test isapprox(dot(X[:, i], X̃[:, i]), 1 - s[i], atol=1e-8)
        else
            @test dot(X[:, i], X̃[:, j]) ≈ dot(X[:, i], X[:, j])
        end
    end
end

@testset "Knockoff data structure" begin
    Random.seed!(2021)

    # simulate matrix and normalize columns
    n = 1000
    p = 100
    X = randn(n, p)
    zscore!(X, mean(X, dims=1), std(X, dims=1)) # center/scale Xj to mean 0 var 1
    normalize_col!(X) # normalize Xj to unit length

    # construct knockoff struct and the real [A Ã]
    @time A = fixed_knockoffs(X, :sdp)
    Atrue = [A.X A.X̃]

    # array operations
    @test size(Atrue) == size(A)
    @test eltype(Atrue) == eltype(A)
    @test getindex(Atrue, 127) == getindex(A, 127)
    @test getindex(Atrue, 2, 19) == getindex(A, 2, 19)
    @test getindex(Atrue, 900, 110) == getindex(A, 900, 110)
    @test all(@view(Atrue[:, 1]) .== @view(A[:, 1]))
    @test all(@view(Atrue[1:2:end, 1:5:end]) .== @view(A[1:2:end, 1:5:end]))

    # matrix-vector multiplication
    b = randn(2p)
    ctrue, c = zeros(n), zeros(n)
    mul!(ctrue, Atrue, b)
    mul!(c, A, b)
    @test all(ctrue .≈ c)

    # matrix-matrix multiplication
    B = randn(2p, n)
    Ctrue, C = zeros(n, n), zeros(n, n)
    mul!(Ctrue, Atrue, B)
    mul!(C, A, B)
    @test all(Ctrue .≈ C)
end

@testset "model X Guassian Knockoffs" begin
    Random.seed!(2021)

    # simulate matrix (but manually normalizing columns leads to rank deficiency for some reason)
    n = 300
    p = 600
    X = randn(n, p)
    # zscore!(X, mean(X, dims=1), std(X, dims=1))
    # normalize_col!(X)
    # @show rank(X)

    # generate knockoff
    @time knockoff = modelX_gaussian_knockoffs(X, :sdp, zeros(p));
    X = knockoff.X
    X̃ = knockoff.X̃
    s = knockoff.s
    Σ = knockoff.Σ
    Σinv = knockoff.Σinv

    # test properties
    @test all(X' * X .≈ Σ)
    @test all(isapprox.(X̃' * X̃, Σ, atol=0.5)) # numerical accuracy not good?
    @test all(s .≥ 0)
    @test all(1 .≥ s)
    for i in 1:p, j in 1:p
        if i == j
            @test isapprox(dot(X[:, i], X̃[:, i]), Σ[i, i] - s[i], atol=1.0) # numerical accuracy not good?
        else
            @test isapprox(dot(X[:, i], X̃[:, j]), dot(X[:, i], X[:, j]), atol=1.0) # numerical accuracy not good?
        end
    end
end

@testset "threshold functions" begin
    w = [0.1, 1.9, 1.3, 1.8, 0.8, -0.7, -0.1]
    @test threshold(w, 0.2) == 0.8
    @test threshold(w, 0.2, :knockoff_plus) == Inf

    w = [0.27, 0.76, 0.21, 0.1, -0.38, -0.01]
    @test threshold(w, 0.4) == 0.1
    @test threshold(w, 0.5, :knockoff_plus) == 0.1

    w = [0.74, -0.65, -0.83, -0.27, -0.19, 0.4]
    @test threshold(w, 0.25) == Inf
    @test threshold(w, 0.25, :knockoff_plus) == Inf
end

@testset "coefficient_diff" begin
    # knockoffs are at the end (e.g. [XX̃])
    β = [1.0, 0.2, -0.3, 0.8, -0.1, 0.5]
    w = coefficient_diff(β, :concatenated)
    @test length(w) == 3
    @test w[1] ≈ 0.2
    @test w[2] ≈ 0.1
    @test w[3] ≈ -0.2

    # knockoffs are interleaved (e.g. [x₁x̃₁x₂x̃₂...])
    β = [1.0, 0.2, -0.3, 0.8, -0.1, 0.5]
    w = coefficient_diff(β, :interleaved)
    @test length(w) == 3
    @test w[1] ≈ 0.8
    @test w[2] ≈ -0.5
    @test w[3] ≈ -0.4

    # multivariate: knockoffs are at the end (e.g. [XX̃])
    B = [[1.0, 2.0, -0.3, 0.8, -0.1, 0.5] [-0.4, 0.6, 1.8, -0.3, 1.4, 0.3]]
    w = coefficient_diff(B, :concatenated)
    @test length(w) == 3
    @test w[1] ≈ 0.3
    @test w[2] ≈ 1.1
    @test w[3] ≈ 1.3

    # multivariate: knockoffs are interleaved (e.g. [x₁x̃₁x₂x̃₂...])
    B = [[1.0, 2.0, -0.3, 0.8, -0.1, 0.5] [-0.4, 0.6, 1.8, -0.3, 1.4, 0.3]]
    w = coefficient_diff(B, :interleaved)
    @test length(w) == 3
    @test w[1] ≈ -1.2
    @test w[2] ≈ 1.0
    @test w[3] ≈ 0.7

    # knockoffs and original are randomly swapped
    β = [1.0, 0.2, -0.3, 0.8, -0.1, 0.5]
    original = [1, 4, 6]
    knockoff = [2, 3, 5]
    w = coefficient_diff(β, original, knockoff)
    @test length(w) == 3
    @test w[1] ≈ 0.8
    @test w[2] ≈ 0.5
    @test w[3] ≈ 0.4

    # group knockoffs
    β = [1.0, 0.2, -0.3, 0.8, -0.1, 0.5, 0.5, -0.1]
    groups = [1, 1, 1, 1, 2, 2, 2, 2]
    original = [1, 4, 6, 7]
    knockoff = [2, 3, 5, 8]
    w = coefficient_diff(β, groups, original, knockoff)
    @test length(w) == 2
    @test w[1] ≈ 1.8 - 0.5
    @test w[2] ≈ 1.0 - 0.2
end