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

using Markdown
using InteractiveUtils

# ╔═╡ 166bbbb2-8f7a-11eb-19c6-19dfe8ee5123
using ChainRules

# ╔═╡ a8439cf0-8eab-11eb-3d6b-a3d65e46a635
md"""
# Automatic Differentiation (AD)

Automatic differentiation is a method to evaluate the derivatives of a given program, you may heard this in context of deep learning as back-propagation. However, in more general scientific machine learning, the situation can be more complicated, thus it'd be necessary to understand better about the mechanism of automatic differentiation.

In this section, we will introduce the automatic differentiation technique in a more systematic aspect
"""

# ╔═╡ 12cb96f8-8f48-11eb-212d-956d24e7a892
md"""
# Dual Numbers

Dual number is the simplest way of implementing AD. And thanks to Julia's powerful generic programming features, we can easily implement one ourselves!

**definition** define ``ϵ`` to satisfy ``ϵ^2 = 0``, then similar to a complex number we define dual number as ``x = a + b ⋅ ϵ``

then we can intermediately get the rule of multiplication

```math
(a + b⋅ϵ)(c + d⋅ϵ) = ac + (ad + bc)ϵ
```

thus we can implement our own dual number as following
"""

# ╔═╡ af43976a-8f48-11eb-1485-0b76c7ee1e79
begin
	struct Dual{T<:Real} <: Real
		a::T
		b::T
	end

	Dual(a) = Dual(a, one(a))
end

# ╔═╡ d1522092-8f48-11eb-30f2-f9c3183ae364
Base.show(io::IO, x::Dual) = print(io, x.a, "+", x.b, "ϵ")

# ╔═╡ e0a6fe78-8f48-11eb-152f-5b448f08fc49
Dual(2)

# ╔═╡ e4fbb1f8-8f48-11eb-220e-f33784f6e617
Base.:(*)(x::Dual, y::Dual) = Dual(x.a * y.a, x.a * y.b + x.b * y.a)

# ╔═╡ 25669170-8f4b-11eb-1cae-a78fc79be709
Base.:(+)(x::Dual, y::Dual) = Dual(x.a + y.a, x.b + y.b)

# ╔═╡ 0c7cdab8-8f49-11eb-3bac-2de9b80bb7d0
Dual(3, 1) * Dual(3, 1)

# ╔═╡ 82b578de-8f49-11eb-13dd-aba04a2c196d
md"now let's define a simply function see what happens if we instead of feeding a normal real number, we feed in a dual number"

# ╔═╡ 28b6ea34-8f49-11eb-2625-d58fbfba667f
f(x) = x^3

# ╔═╡ 9f7849bc-8f49-11eb-0f59-1bdf987ecdfe
f(Dual(2))

# ╔═╡ a7d73bb6-8f49-11eb-2c45-99fbd3a2c248
md"""
the prefactor of ϵ gives the derivative! let's see why

given a real polynomail

```math
P(x) = p_0 + p_1 x + p_2 x^2 + ⋯ + p_n x^n
```

if we put a dual number ``x = a + b ϵ``, we have

```math
\begin{aligned}
P(a+bϵ)
&= p_0 + p_1 a + p_2 a^2 + ⋯ + p_n a^n + p_1 bϵ + 2p_2 abϵ + ⋯ + n p_n a^{n-1}b ϵ\\
&= P(a) + b P'(a)ϵ
\end{aligned}
```

more generally, using Taylor series, this is also true for any real analytical functions.

However, if we plugin our dual number into functions like `sin`, it doesn't work
"""

# ╔═╡ 8a4b5ec8-8f4a-11eb-1155-f358a173daa2
sin(Dual(2))

# ╔═╡ 94bc3b16-8f4a-11eb-0871-59f5469f4c0e
md"""
this is because in practical numerical program, for efficiency, a lot mathematical functions are implemented specifically for floating point type, thus we don't have a generic function for these functions to make things "just work"

thus we will need to define the derivative of `sin` manually here
"""

# ╔═╡ cfdb25b6-8f4a-11eb-2248-3f7783956f35
Base.sin(x::Dual) = Dual(sin(x.a), cos(x.a) * x.b)

# ╔═╡ de6df808-8f4a-11eb-0f70-43400519fd6f
sin(Dual(2))

# ╔═╡ f2b1a960-8f4a-11eb-1276-b333b086605c
md"however if we happen to have a function definition on generic types, e.g the matrix multiplication, this will just work"

# ╔═╡ 0906c78e-8f4b-11eb-1d3b-dd8d6d891175
A = [Dual(1.0) Dual(2.0); Dual(1.0) Dual(2.0)]

# ╔═╡ 149734d0-8f4b-11eb-0af4-6368d8bb34c0
B = [Dual(1.0) Dual(2.0); Dual(1.0) Dual(3.0)]

# ╔═╡ 1a483212-8f4b-11eb-0fa5-5994548ea6b9
A * B

# ╔═╡ f4228cb2-8eab-11eb-1b98-fbb30f8eeac3
md"""
# Forward Mode AD

The forward mode AD is the simplest way of implementing an AD system. The above dual number is one way to implement forward mode AD with a generic type system.

More formally, we say forward mode AD is to apply the chain rule on primal derivatives
when we evaluate the value of an expression, e.g if we are going to evaluate the following function
"""

# ╔═╡ c6d79168-8f4f-11eb-1c4f-4bf4abbb2c72
function foo(x)
	x1 = sin(x)
	x2 = sin(x1)
	x3 = cos(x2) + x2
end

# ╔═╡ d2f7c68e-8f4f-11eb-0d90-4f2f22bd87b6
md"""
when we put in a dual number we actually do the following things:

1. evaluate `sin(x)` and its derivative, store them in `x1`
2. evaluate `sin(x1)` and its derivative, multiply the derivative of `sin(x1)` to previous derivative since we have ``P'(x)b ϵ`` term in the dual number
3. evaluate the derivative of `cos(x2)` and do the same thing as above, but let's call the intermediate value of `cos(x2)` to be `y`
4. evaluate the value of `y+x2` and derivative then multiply and store them in `x3`

thus we see every time we evaluate the derivative, we can throw away the variables we calculated before, but for every single scalar number we need to calculate the entire function, as a result, this gives us the time complexity $O(mn)$ where $n$ is the number of parameters and $m$ is the complexity of the original code.

the memory complexity is $O(m+n)$ where $m$ is the memory complexity of the original code
"""

# ╔═╡ f75c7e22-8f52-11eb-0b32-f9ed6267fb18
md"""
to formalize the above process in terms of chain rule, the forward mode means we evaluate the derivatives in the following accumulation:

```math
\begin{aligned}
\frac{\partial y_n}{\partial x}
&= \frac{\partial y_n}{\partial y_{n-1}} \frac{\partial y_{n-1}}{\partial x}\\
&= \frac{\partial y_n}{\partial y_{n-1}} (\frac{\partial y_{n-1}}{\partial y_{n-2}} \frac{\partial y_{n-2}}{\partial x})\\
&= \cdots\\
&= \frac{\partial y_n}{\partial y_{n-1}} (\frac{\partial y_{n-1}}{\partial y_{n-2}} \cdots (\frac{\partial y_2}{\partial y_1} \frac{\partial y_1}{\partial x}))
\end{aligned}
```

"""

# ╔═╡ 9ef95826-8f53-11eb-11ed-a35e16d20083
md"""
we can further generalize this to multi-variable case using Jacobians

```math
\begin{aligned}
\frac{\partial \vec{y}_n}{\partial x}
&= \frac{\partial \vec{y}_n}{\partial \vec{y}_{n-1}} \frac{\partial \vec{y}_{n-1}}{\partial x}\\
&= \frac{\partial \vec{y}_n}{\partial \vec{y}_{n-1}} (\frac{\partial \vec{y}_{n-1}}{\partial \vec{y}_{n-2}} \frac{\partial \vec{y}_{n-2}}{\partial x})\\
&= \cdots\\
&= \frac{\partial \vec{y}_n}{\partial \vec{y}_{n-1}} (\frac{\partial \vec{y}_{n-1}}{\partial \vec{y}_{n-2}} \cdots (\frac{\partial \vec{y}_2}{\partial \vec{y}_1} \frac{\partial \vec{y}_1}{\partial x}))
\end{aligned}
```
"""

# ╔═╡ 0e07a38e-8f55-11eb-0c46-ef5d11df07de
md"""
now we see to calculate the derivative of ``\frac{\partial \vec{y}_n}{\partial x}`` it is actually about calculatnig a chain of **Jacobian-vector product**
"""

# ╔═╡ 13c539e8-8eac-11eb-35b0-61108d5f3c5f
md"""
# Reverse Mode AD

For large number of parameters, reverse mode is more efficient. This method has been re-discovered many times in history in different fields, you may hear it called back-propagation in context of deep learning.
"""

# ╔═╡ 45826916-8f55-11eb-3f3c-95d6b601191b
md"""
similar to forward mode AD, the reverse mode AD is actually about calculating the chain rule in the following accumulation:

```math
\begin{aligned}
\frac{\partial \vec{y}_n}{\partial x}
&= (\frac{\partial \vec{y}_n}{\partial \vec{y}_{n-1}}) \frac{\partial \vec{y}_{n-1}}{\partial x}\\
&= ((\frac{\partial \vec{y}_n}{\partial \vec{y}_{n-1}}) \frac{\partial \vec{y}_{n-1}}{\partial \vec{y}_{n-2}}) \frac{\partial \vec{y}_{n-2}}{\partial x}\\
&= ⋯\\
&= (\cdots ((\frac{\partial \vec{y}_n}{\partial \vec{y}_{n-1}}) \frac{\partial \vec{y}_{n-1}}{\partial \vec{y}_{n-2}}\cdots \frac{\partial \vec{y}_2}{\partial \vec{y}_1}) \frac{\partial \vec{y}_1}{\partial x}\\
(\frac{\partial \vec{y}_n}{\partial x})^T &= (\frac{\partial \vec{y}_1}{\partial x})^T (\frac{\partial \vec{y}_2}{\partial \vec{y}_1})^{T} \cdots (\frac{\partial \vec{y}_{n-1}}{\partial \vec{y}_{n-2}})^T (\frac{\partial \vec{y}_n}{\partial \vec{y}_{n-1}})^T
\end{aligned}
```

thus the reverse mode automatic differentiation is acutally about **Jacobian-transpose-vector product**.

We can see an intermediate advantage of this accumulation is that if ``x`` is a vector, we can directly calculate the derivative of this vector by a chain of matrix multiplication.
"""

# ╔═╡ e026ffa6-8f5f-11eb-0651-dfb9e28fcc4f
md"""
this process will be easier to understand if we use a graphical language to describe it, this the **computational graph**, to demonstrate this better, I'll use a more complicated function

```math
y = \mathbf{x}^T\mathbf{A}\mathbf{x} + \mathbf{b}\cdot\mathbf{x} + c
```
"""

# ╔═╡ 10ba0d9c-8f5f-11eb-0744-6be46cfacebe
md"""
We will need to call several functions in Julia to get the result $y$, which is

1. ``\mathbf{z_1} = \mathbf{x}^T``: `transpose` function.
2. ``\mathbf{z_2} = \mathbf{z_1} A`` matrix-vector multiplication, which can be `gemv` in `LinearAlgebra.BLAS`, or just `*`.
3. ``y_1 = \mathbf{z_2} \mathbf{x}`` vector dot operation, which is `LinearAlgebra.dot` or the UTF-8 operator `x ⋅ y`
4. ``y_2 = \mathbf{b} \cdot \mathbf{x}`` another vector dot
5. ``y_1 + y_2 + c`` a scalar add function, one can calculate it by simply calling `+` operator in Julia.

In fact, we can draw a graph of this expression, which illustrates the relationship between each variable in this expression.
Each node in the graph with an output arrow represents a variable and each node with an input arrow represents a function/operator.

![](https://blog.rogerluo.me/images/comput-graph-forward.gif)

The evaluation of the math equation above can then be expressed as a process called **forward evaluation**, it starts from the leaf nodes, which represents the inputs of the whole expression, e.g they are $\mathbf{x}, \mathbf{A}, \mathbf{b}, c$ in our expression. Each time, we receive the value of a node in the graph, we mark the node with **green**.
"""

# ╔═╡ 4921f4f8-8f60-11eb-0e45-3f3034dbea83
md"""
the derivative calculation can be then visualized as following

![](https://blog.rogerluo.me/images/comput-graph-backward.gif)
"""

# ╔═╡ 29fa38b6-8f61-11eb-118a-616fc647715a
md"""
as we demonstrate above, although the reverse mode can calculate the derivatives for all the parameters in single evaluation. it will need to store the intermediate values.

thus we need a specific data structure to store these values， this data structure is often called the computational graph, tape or Wengert list.

Now the problem of implementing reverse mode AD becomes how to create the tape.
"""

# ╔═╡ 0580566c-8ff4-11eb-35a3-bf0505fe5a1c
md"""
# Relation with Differential Geometry

Since we have represented our AD algorithm as Jacobian or Jacobian transpose product, it is natural to have a geometric interpretation of them.

## Definitions

Let ``p`` be a point in the space ``\mathbf{M}`` defined by ``(x_1, x_2, ⋯, x_m)``, and we have a function ``f: \mathbf{M} → \mathbf{K}`` that is a smooth map (differentiable everywhere), then we can define the **tangent space** ``\mathbf{T}_p \mathbf{M}`` at point ``p`` as the space span by the derivatives ``\frac{d}{dx}``, the tangent space acts linearly on the space of functions.

Now we define **the collection of the tangent spaces** ``\mathbf{T}_p\mathbf{M}`` for ``p ∈ \mathbf{M}`` as the **tangent bundle** of ``\mathbf{M}``

![tangent space](https://upload.wikimedia.org/wikipedia/commons/thumb/e/e7/Tangentialvektor.svg/200px-Tangentialvektor.svg.png)

let ``df`` denotes the first order information of ``f`` at each point, which is called the **differential** of ``f``. If the derivative of ``f`` and ``g`` agree at ``p``, we say that ``df`` and ``dg`` represents the same cotangent at ``p``. The covectors ``dx_1, dx_2, ⋯, dx_m`` form the basis of the cotangent space ``\mathbf{T}^*_p \mathbf{M}`` for ``p∈\mathbf{M}``

**push-forward**

![](https://upload.wikimedia.org/wikipedia/commons/thumb/3/37/Pushforward.svg/330px-Pushforward.svg.png)

**pullback**

![](https://www.mathphysicsbook.com/wp-content/uploads/2013/10/39.pullback-v2.png)
"""

# ╔═╡ eb4df5e4-8f79-11eb-2429-a31bd0009c4b
md"""
## the manual AD

so we can do reverse mode "automatic" differentiation manually first, to have a feeling about it.

We will use the package `ChainRules` here.
"""

# ╔═╡ 72b98a42-8ff8-11eb-17f3-6d14da3964be
x = 3.0

# ╔═╡ 1a78b2be-8f7a-11eb-0470-e7febaac6272
y, sin_pullback = rrule(sin, x)

# ╔═╡ 627f170a-8ff8-11eb-21eb-47796905e10b
sin_pullback(1)

# ╔═╡ 67c7f3d0-8ff8-11eb-04dc-4f796cfb940e
cos(x)

# ╔═╡ 7cffb62a-8ff8-11eb-0288-f55b0d333571
md"""
now we can manually generate the pullback function for a given function defined on some primal functions, in AD, we usually call these functions **adjoint**s.

let's use the previous example `foo` here
"""

# ╔═╡ 98d147c4-8ff8-11eb-3ab8-21e26217e02a
function adjoint_foo(x)
	x1, x1_pullback = rrule(sin, x)
	x2, x2_pullback = rrule(sin, x1)
	y, y_pullback = rrule(cos, x2)
	x3, x3_pullback = rrule(+, y, x2)

	return x3, function pullback(Δ)
		_, ∂y, ∂x2_1 = x3_pullback(Δ)
		_, ∂x2_2 = y_pullback(∂y)
		_, ∂x1 = x2_pullback(∂x2_1 + ∂x2_2)
		_, ∂x = x1_pullback(∂x1)
		return Zero(), ∂x
	end
end

# ╔═╡ 8e593406-8ff9-11eb-101a-695b2b026a20
_, foo_pullback = adjoint_foo(x)

# ╔═╡ 9e31a492-8ff9-11eb-108d-1b7bbd152170
foo_pullback(1)

# ╔═╡ a21aa0c2-8ff9-11eb-1cc8-7b63a20d0791
(foo(x+1e-10)-foo(x-1e-10))/2e-10

# ╔═╡ c5f14d48-8ff9-11eb-1350-790162aa8b34
md"""
thus the simplest reverse mode AD is about how to create this tiny `pullback` function automatically, and to know how to create this `pullback` function we need to know what primal functions are called by the given function, so that we can simply reverse the order of calls and replace them with the `pullback` functions.
"""

# ╔═╡ 3fc036f8-8eac-11eb-15dc-017c0d5237bc
md"""
## Operator Overloading AD

The easiest way of creating a tape is using operator overloading. By operator overloading, it means the exactly same way we implemented our toy symbolic program in the previous section.

for example, we can define our own type to dispatch the functions to a `track` function so that we can store the function call into a tape.

```julia
mutable struct Variable{T} <: ADExpr
    value::T
    grad::T

    Variable(val::T) where T = new{T}(val)
    Variable(val::T, grad::T) where T = new{T}(val)
end

struct Node{FT <: Function, ArgsT <: Tuple, KwargsT <: NamedTuple} <: ADExpr
    f::FT
    args::ArgsT
    kwargs::KwargsT
end
```

then overload some primal functions to track the call into tape

```
Base.sin(x::ADExpr) = register(Base.sin, x)
```
"""

# ╔═╡ 35cfef70-8ffa-11eb-1a48-e906596eebe6
md"""
I will not demonstrate the details of implementing such AD engine, but once you understand the mechanism, you should be able to write one yourself in mins!

For the curious, checkout my blog post **Implementing Your Own Automatic Differentiation Engine in ONE day**: [https://blog.rogerluo.me/2018/10/23/write-an-ad-in-one-day/](https://blog.rogerluo.me/2018/10/23/write-an-ad-in-one-day/)
"""

# ╔═╡ 73ffb3fe-8ffd-11eb-2b71-038543ea5093
md"""

## What packages/software uses this method?

#### Python

- PyTorch
- autograd
- ...

#### Julia

- Tracker
- YAAD
- Nabla
- ...

"""

# ╔═╡ bd69872c-8f78-11eb-2197-2f760c33b197
md"""

## Limitations

however, the above method **won't work for general programs that contains control flows**. This is because in the representation of Wengert list, or the tape, we don't have the semantic of control flows, and we don't have a way to record control flows using operator overloading, such as `for` loop or `if else` statements.

This is because control flows are not functions in most of the programming languages, we cannot change the control flow behaviour by overloading them on different types, and it may not make sense to do so.

On the other hand, it **won't work for generic programming**, since we have used a specific type to dispatch functions, which kinda "abused" the type system.
"""

# ╔═╡ 46ad0fa2-8eac-11eb-28df-4f2c23f8d3ff
md"""
## Source Code Transformation AD

As you may noticed, in both forward and reverse mode AD, in order to represent intermediate derivatives, we were forced to use a new temperory variable for every function call.

This leads to a special intermediate representation for general programs called the **Static Single Assignment** intermediate representation (SSA IR), by **Static Single Assignment** we mean:

> every variable in the program is only assigned once statically

the control flows in a SSA IR are usually canonicalized to `goto` or `goto ifnot` if it has conditions, since every variable is only assigned once we use `%<name>` or `%<int>` to represent a variable.
"""

# ╔═╡ cbe1949e-8ffc-11eb-0747-e576545133b4
md"""
the SSA IR is a widely used representation in your daily used compilers, for example, your C/C++ compiler clang, or some of your new fancy langauges like Julia, rust and more. This is because SSA IR is a very convenient representation to optimize and analysis.

![](https://slideplayer.com/slide/13900843/85/images/7/C-like+IR+In+SSA+form+x1+%3D+f%28x0%29%3B+x+%3D+f%28x%29%3B+if+%28x1+%3E+y0%29.jpg)
"""

# ╔═╡ 6616ee24-8ffd-11eb-0d13-7f7ea01d3f73
md"""
this is a natural extension of the Wengert list (the tape) with control flows. Thus instead of doing operator overloading, we can directly write a compiler extension to transform the source encoded in SSA IR to the corresponding differential function in SSA IR.

And in Julia, this is extremely easy (in terms of implementing it), because Julia is a dynamic compiled language

![](https://blog.rogerluo.me/images/julia-compile-diagram.png)
"""

# ╔═╡ df1d9e58-90de-11eb-341a-bb942e20e713
md"""
For the curious, checkout my blog post **Build your own source to source automatic differentiation in ONE day!** :[https://blog.rogerluo.me/2019/07/27/yassad/](https://blog.rogerluo.me/2019/07/27/yassad/)
"""

# ╔═╡ c5f18fe2-8ffe-11eb-0a96-1f3ce397e943
md"""
# Software/Package uses this method

- **Julia**: Zygote, Enzyme, Yota, Diffractor (work-in-progress), ...
- **Python**: JAX, ...
"""

# ╔═╡ 420d72f4-8ffe-11eb-1b0f-43bdc9edf955
md"""
## Limitations

No limitations in theory.

We have the information about control flows, which means it **will work for control flows**.

We don't need to reply on type system to track the function calls anymore which means it works for **generic programming**

But covering the entire general purpose programming language with most of the corner cases may take years to accomplish.
"""

# ╔═╡ 84e5fe6e-8eac-11eb-04f8-cf7499b2c218
md"""
## Other type of AD engines

There are other type of AD engines that works under some context, such as polynomials, reversible programs, etc.

You can find a full list of AD engines in Julia at [https://juliadiff.org/](https://juliadiff.org/)
"""

# ╔═╡ 2e3fc494-90de-11eb-1c90-f9cc1357aab3
md"""

### Is there an algebra gives us reverse mode AD?

- reversible programming

![](https://giggleliu.github.io/NiLang.jl/dev/asset/revcomp.jpg)

"""

# ╔═╡ 8d13dbe0-90de-11eb-1bcf-dfbd8118e59d
md"""
## NiLang: reversible domain-specific programming language in Julia

![](https://github.com/GiggleLiu/NiLang.jl/raw/master/docs/src/asset/logo3.png)
"""

# ╔═╡ 4c23d3bc-8eac-11eb-2870-a37dfc123bc5
md"""
## Differentiable Programming in Condensed Matter Physics: Tensor Network

One of the most poplular method in numerical condensed matter physics is the tensor network methods.

Tensor network is also way to approximate the solution, in previous lectures, you have learned how to use neural network states to solve ground state problem using variational Monte Carlo.

More traditionally physicists have been using tensor network as the physical ansatz. However calculating the gradient of tensor networks can be complicated. Thus, by making use of automatic differentiation, one will be able to do gradient based optimization on tensor network easily.

![](https://github.com/Roger-luo/PIML/blob/master/notebooks/assets/tensor-network.png?raw=true)

[https://arxiv.org/pdf/1903.09650.pdf](https://arxiv.org/pdf/1903.09650.pdf)
"""

# ╔═╡ 619d6b8e-90d4-11eb-17fa-dd41790aa420
md"""
## Further Readings

- the ChainRules package documentation: [https://juliadiff.org/ChainRulesCore.jl/stable/](https://juliadiff.org/ChainRulesCore.jl/stable/)
- NiLang tutorial for reversible programming: [https://giggleliu.github.io/NiLang.jl/dev/why/](https://giggleliu.github.io/NiLang.jl/dev/why/)
- the differential and pullback: [https://www.mathphysicsbook.com/mathematics/manifolds/mapping-manifolds/the-differential-and-pullback/](https://www.mathphysicsbook.com/mathematics/manifolds/mapping-manifolds/the-differential-and-pullback/)
"""

# ╔═╡ Cell order:
# ╟─a8439cf0-8eab-11eb-3d6b-a3d65e46a635
# ╟─12cb96f8-8f48-11eb-212d-956d24e7a892
# ╠═af43976a-8f48-11eb-1485-0b76c7ee1e79
# ╠═d1522092-8f48-11eb-30f2-f9c3183ae364
# ╠═e0a6fe78-8f48-11eb-152f-5b448f08fc49
# ╠═e4fbb1f8-8f48-11eb-220e-f33784f6e617
# ╠═25669170-8f4b-11eb-1cae-a78fc79be709
# ╠═0c7cdab8-8f49-11eb-3bac-2de9b80bb7d0
# ╟─82b578de-8f49-11eb-13dd-aba04a2c196d
# ╠═28b6ea34-8f49-11eb-2625-d58fbfba667f
# ╠═9f7849bc-8f49-11eb-0f59-1bdf987ecdfe
# ╟─a7d73bb6-8f49-11eb-2c45-99fbd3a2c248
# ╠═8a4b5ec8-8f4a-11eb-1155-f358a173daa2
# ╟─94bc3b16-8f4a-11eb-0871-59f5469f4c0e
# ╠═cfdb25b6-8f4a-11eb-2248-3f7783956f35
# ╠═de6df808-8f4a-11eb-0f70-43400519fd6f
# ╟─f2b1a960-8f4a-11eb-1276-b333b086605c
# ╠═0906c78e-8f4b-11eb-1d3b-dd8d6d891175
# ╠═149734d0-8f4b-11eb-0af4-6368d8bb34c0
# ╠═1a483212-8f4b-11eb-0fa5-5994548ea6b9
# ╟─f4228cb2-8eab-11eb-1b98-fbb30f8eeac3
# ╠═c6d79168-8f4f-11eb-1c4f-4bf4abbb2c72
# ╟─d2f7c68e-8f4f-11eb-0d90-4f2f22bd87b6
# ╟─f75c7e22-8f52-11eb-0b32-f9ed6267fb18
# ╟─9ef95826-8f53-11eb-11ed-a35e16d20083
# ╟─0e07a38e-8f55-11eb-0c46-ef5d11df07de
# ╟─13c539e8-8eac-11eb-35b0-61108d5f3c5f
# ╟─45826916-8f55-11eb-3f3c-95d6b601191b
# ╟─e026ffa6-8f5f-11eb-0651-dfb9e28fcc4f
# ╟─10ba0d9c-8f5f-11eb-0744-6be46cfacebe
# ╟─4921f4f8-8f60-11eb-0e45-3f3034dbea83
# ╟─29fa38b6-8f61-11eb-118a-616fc647715a
# ╟─0580566c-8ff4-11eb-35a3-bf0505fe5a1c
# ╟─eb4df5e4-8f79-11eb-2429-a31bd0009c4b
# ╠═166bbbb2-8f7a-11eb-19c6-19dfe8ee5123
# ╠═72b98a42-8ff8-11eb-17f3-6d14da3964be
# ╠═1a78b2be-8f7a-11eb-0470-e7febaac6272
# ╠═627f170a-8ff8-11eb-21eb-47796905e10b
# ╠═67c7f3d0-8ff8-11eb-04dc-4f796cfb940e
# ╟─7cffb62a-8ff8-11eb-0288-f55b0d333571
# ╠═98d147c4-8ff8-11eb-3ab8-21e26217e02a
# ╠═8e593406-8ff9-11eb-101a-695b2b026a20
# ╠═9e31a492-8ff9-11eb-108d-1b7bbd152170
# ╠═a21aa0c2-8ff9-11eb-1cc8-7b63a20d0791
# ╟─c5f14d48-8ff9-11eb-1350-790162aa8b34
# ╟─3fc036f8-8eac-11eb-15dc-017c0d5237bc
# ╟─35cfef70-8ffa-11eb-1a48-e906596eebe6
# ╟─73ffb3fe-8ffd-11eb-2b71-038543ea5093
# ╟─bd69872c-8f78-11eb-2197-2f760c33b197
# ╟─46ad0fa2-8eac-11eb-28df-4f2c23f8d3ff
# ╟─cbe1949e-8ffc-11eb-0747-e576545133b4
# ╟─6616ee24-8ffd-11eb-0d13-7f7ea01d3f73
# ╟─df1d9e58-90de-11eb-341a-bb942e20e713
# ╟─c5f18fe2-8ffe-11eb-0a96-1f3ce397e943
# ╟─420d72f4-8ffe-11eb-1b0f-43bdc9edf955
# ╟─84e5fe6e-8eac-11eb-04f8-cf7499b2c218
# ╟─2e3fc494-90de-11eb-1c90-f9cc1357aab3
# ╟─8d13dbe0-90de-11eb-1bcf-dfbd8118e59d
# ╟─4c23d3bc-8eac-11eb-2870-a37dfc123bc5
# ╟─619d6b8e-90d4-11eb-17fa-dd41790aa420

max_stars_repo_name: Roger-luo/PIML

max_issues_repo_name: Roger-luo/PIML

max_forks_repo_name: Roger-luo/PIML


text:
importall MathProgBase.SolverInterface
importall SemidefiniteModels

export CSDPMathProgModel, CSDPSolver

immutable CSDPSolver <: AbstractMathProgSolver
    options::Dict{Symbol,Any}
end
CSDPSolver(;kwargs...) = CSDPSolver(Dict{Symbol,Any}(kwargs))

type CSDPMathProgModel <: AbstractSDModel
    C
    b
    As
    X
    y
    Z
    status::Cint
    pobj::Cdouble
    dobj::Cdouble
    options::Dict{Symbol,Any}
    function CSDPMathProgModel(; kwargs...)
        new(nothing, nothing, nothing, nothing, nothing, nothing,
            -1, 0.0, 0.0, Dict{Symbol, Any}(kwargs))
    end
end
SDModel(s::CSDPSolver) = CSDPMathProgModel(; s.options...)
ConicModel(s::CSDPSolver) = SDtoConicBridge(SDModel(s))

supportedcones(s::CSDPSolver) = [:Free,:Zero,:NonNeg,:NonPos,:SDP]
function setvartype!(m::CSDPMathProgModel, vtype, blk, i, j)
    if vtype != :Cont
        error("Unsupported variable type $vtype by CSDP")
    end
end

function loadproblem!(m::CSDPMathProgModel, filename::AbstractString)
    if endswith(filename,".dat-s")
       m.C, m.b, As = read_prob(filename)
       m.As = [ConstraintMatrix(As[i], i) for i in 1:length(As)]
    else
       error("unrecognized input format extension in $filename")
    end
end
#writeproblem(m, filename::String)
function loadproblem!(m::CSDPMathProgModel, blkdims::Vector{Int}, constr::Int)
    m.C = blockmatzeros(blkdims)
    m.b = zeros(Cdouble, constr)
    m.As = [constrmatzeros(i, blkdims) for i in 1:constr]
end

function setconstrB!(m::CSDPMathProgModel, val, constr::Integer)
    m.b[constr] = val
end
function setconstrentry!(m::CSDPMathProgModel, coef, constr::Integer, blk::Integer, i::Integer, j::Integer)
    m.As[constr][blk][i,j] = coef
end
function setobjentry!(m::CSDPMathProgModel, coef, blk::Integer, i::Integer, j::Integer)
    m.C[blk][i,j] = coef
end

function optimize!(m::CSDPMathProgModel)
    As = map(A->A.csdp, m.As)

    let wrt = string(get(m.options, :write_prob, ""))
        if length(wrt) > 0
            k = 1
            wrtf = "$wrt.$k"
            while isfile(wrtf)
                wrtf = "$wrt.$k"
                k += 1
            end
            info("Writing problem to $(pwd())/$(wrtf)")
            write_prob(wrtf, m.C, m.b, As)
        end
    end

    m.X, m.y, m.Z = initsoln(m.C, m.b, As)
    m.status, m.pobj, m.dobj = easy_sdp(m.C, m.b, As, m.X, m.y, m.Z)
end

function status(m::CSDPMathProgModel)
    status = m.status
    if status == 0
        return :Optimal
    elseif status == 1
        return :Infeasible
    elseif status == 2
        return :Unbounded
    elseif 3 <= status <= 7
        return :Unknown
    elseif 8 <= status <= 9
        return :Error
    elseif status == -1
        return :Uninitialized
    else
        error("Internal library error: status=$status")
    end
end

function getobjval(m::CSDPMathProgModel)
    m.pobj
end
function getsolution(m::CSDPMathProgModel)
    m.X
end
function getdual(m::CSDPMathProgModel)
    m.y
end
function getvardual(m::CSDPMathProgModel)
    m.Z
end

ext: jl


text:
using WhereTheWaterFlows
const WWF = WhereTheWaterFlows
using Test

"An artificial DEM"
function dem1(x, y; withpit=false, randfac=0.0)
    out = - (x^2 - 1)^2 - (x^2*y - x - 1)^2 + 6 + 0.1*x + 3*y
    if withpit
        out -= 2*exp(-(x^2 + y^2)*50)
    end
    out += randfac*randn()
    return out<0 ? 0.0 : out
end

const _offset = [0.2, 0.5]

"DEM with a few maxs and mins"
function peaks(n=100)
    coords = range(-pi+_offset[1], pi-_offset[2], length=n)
    return coords, sin.(coords) .* cos.(coords')
end

"DEM with a few more maxs and mins"
function peaks2(n=100, randfac=0.0)
    coords = range(-pi+_offset[1], pi-_offset[2], length=n)
    return coords, sin.(coords) .* cos.(coords') .-
        0.7*(sin.(coords.+1) .* cos.(coords')).^8 .+
        0.01 * (coords.+coords') .+
        randfac*randn(n,n) # 0.02
end

"DEM with a few more maxs and mins.  NaN masked middle."
function peaks2_nan(n=100, randfac=0.0)
    coords = range(-pi+_offset[1], pi-_offset[2], length=n)
    dem = sin.(coords) .* cos.(coords') .-
        0.7*(sin.(coords.+1) .* cos.(coords')).^8 .+
        0.01 * coords.*coords' .+
        randfac*randn(n,n) # 0.02
    dem[(n÷2:n÷2+2).-8, n÷2+1:n÷2+4] .= NaN
    return coords, dem
end

"DEM with a few more maxs and mins.  NaN masked edge."
function peaks2_nan_edge(n=100, randfac=0.0)
    coords = range(-pi+_offset[1], pi-_offset[2], length=n)
    dem = sin.(coords) .* cos.(coords') .-
        0.7*(sin.(coords.+1) .* cos.(coords')).^8 .+
        0.01 * (coords.+coords') .+
        randfac*randn(n,n) # 0.02
    outflow = dem[n÷2, 1]
    dem[1, :] .= NaN
    dem[end, :] .= NaN
    dem[:, 1] .= NaN
    dem[:, end] .= NaN
    dem[n÷2, 1] = outflow
    return coords, dem
end

# Test low-level fns

# flowsinto
# """
# Tests whether a cell `J` with flowdir `dirJ` flows into cell `I`.
# """
# flowsinto(J::CartesianIndex, dirJ, I::CartesianIndex) = ind2dir(I-J) == dirJ
@testset "flowsinto" begin
    J = CartesianIndex(5,6)
    dirs = 0:9
    for i=[3,7], j=4:8, dir in dirs
        @test_throws BoundsError WWF.flowsinto(J, dir, CartesianIndex(i,j))
    end
    for i=3:7, j=[4,8], dir in dirs
        @test_throws BoundsError WWF.flowsinto(J, dir, CartesianIndex(i,j))
    end
    trueflow = Dict(1 => J+CartesianIndex(-1,-1),
                    2 => J+CartesianIndex(0,-1),
                    3 => J+CartesianIndex(1,-1),
                    4 => J+CartesianIndex(-1,0),
                    5 => J+CartesianIndex(0,0),
                    6 => J+CartesianIndex(1,0),
                    7 => J+CartesianIndex(-1,1),
                    8 => J+CartesianIndex(0,1),
                    9 => J+CartesianIndex(1,1)
                    )
    for (td, I) in trueflow
        for d in dirs
            if d==td
                @test WWF.flowsinto(J, d, I)
            else
                @test !WWF.flowsinto(J, d, I)
            end
        end
    end
end

@test WWF.dirnums==reverse([ 7 8 9
                             4 5 6
                             1 2 3]', dims=2)
tmp = convert(Matrix, WWF.dirnums)

for i=1:9
    dem = ones(3,3)
    dem[i] = 0
    dir, nout, nin, pits = WWF.d8dir_feature(dem, false)

    @test dir[2,2] == WWF.dirnums[i]
    @test dem[WWF.dir2ind(dir[2,2])+CartesianIndex(2,2)]==0

    @test WWF.d8dir_feature(dem, false)[1][2,2] == i
    @test maximum(nout)<=1
    @test minimum(nout)>=0
    @test maximum(nin)<=8
    @test minimum(nin)>=0
    @test sum(nout)==sum(nin)

    @test sum(dem[WWF.iterate_D9(CartesianIndex(2,2), dem)]) == 8

end

@testset "dir2ind" begin
    @test WWF.dir2ind(1) == CartesianIndex(-1,-1)
    @test WWF.dir2ind(2) == CartesianIndex(0,-1)
    @test WWF.dir2ind(3) == CartesianIndex(1,-1)
    @test WWF.dir2ind(4) == CartesianIndex(-1,0)
    @test WWF.dir2ind(5) == CartesianIndex(0,0)
    @test WWF.dir2ind(6) == CartesianIndex(1,0)
    @test WWF.dir2ind(7) == CartesianIndex(-1,1)
    @test WWF.dir2ind(8) == CartesianIndex(0,1)
    @test WWF.dir2ind(9) == CartesianIndex(1,1)
    for i=1:9
        @test i==WWF.ind2dir(WWF.dir2ind(i))
    end
end

# Test high-level functions
@testset "DEM: dem1" begin
    dx = 0.9
    xs = -1.5:dx:1
    ys = -0.5:dx:3.0
    dem = dem1.(xs, ys')
    @test size(dem)==(length(xs), length(ys))

    area, slen, dir, nout, nin, pits = WWF.waterflows(dem, bnd_as_pits=false)
    @test area == [1.0 1.0 5.0 1.0; 1.0 1.0 1.0 1.0; 5.0 3.0 2.0 1.0]
    @test slen == [1 1 2 1; 1 1 1 1; 4 3 2 1]
    @test dir == Int8[5 8 5 5; 6 7 4 1; 5 2 2 2]
    @test nout == Bool[false true false false; true true true true; false true true true]
    @test nin == Int8[0 0 4 0; 0 0 0 0; 2 1 1 0]
    @test pits ==CartesianIndex{2}[CartesianIndex(1, 1), CartesianIndex(3, 1), CartesianIndex(1, 3), CartesianIndex(1, 4)]
end

@testset "DEM: peaks2" begin
    xs, dem = peaks2()
    ys = xs
    @test size(dem)==(length(xs), length(ys))
    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=false);
    #plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 6
    @test maximum(slen)==67
    @test maximum(area)==4831
    @test length(unique(c))==6
    @test sum(c) == 33876
    @test sum(diff(c[:])) == 5
    @test sort(unique(c))[[1,end]] ==[1,6]

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 4
    @test maximum(slen)==118
    @test maximum(area)==7714
    @test length(unique(c))==4
    @test sum(c) == 21990
    @test sum(diff(c[:])) == 2
    @test sort(unique(c))[[1,end]] ==[1,4]
    @test bnds isa Array{Array{CartesianIndex,1},1}
end

@testset "DEM: peaks2_nan" begin
    xs, dem = peaks2_nan()
    ys = xs
    @test size(dem)==(length(xs), length(ys))

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=false, bnd_as_pits=false);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 7
    @test maximum(slen)==67
    @test maximum(area[.!isnan.(area)])==4625
    @test length(unique(c))==8
    @test sum(c) == 36607
    @test sum(diff(c[:])) == 6
    @test sort(unique(c))[[1,end]] ==[0,7]

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=false, bnd_as_pits=true);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 408
    @test maximum(slen)==66
    @test maximum(area[.!isnan.(area)])==4114
    @test length(unique(c))==409
    @test sum(c) == 2083959
    @test sum(diff(c[:])) == 407
    @test sort(unique(c))[[1,end]] ==[0,408]

    # fill pits
    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=false);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 4
    @test maximum(slen)==116
    @test maximum(area[.!isnan.(area)])==7707
    @test length(unique(c))==5
    @test sum(c) == 24117
    @test sum(diff(c[:])) == 3
    @test sort(unique(c))[[1,end]] ==[1,5]
    @test bnds isa Array{Array{CartesianIndex,1},1}

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=true);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 406
    @test maximum(slen)==75
    @test maximum(area[.!isnan.(area)])==4115
    @test length(unique(c))==407
    @test sum(c) == 1444498
    @test sum(diff(c[:])) == 406
    @test sort(unique(c))[[1,end]] ==[1,407]
    @test bnds isa Array{Array{CartesianIndex,1},1}
end

@testset "DEM: peaks2_nan_edge" begin
    xs, dem = peaks2_nan_edge()
    ys = xs
    @test size(dem)==(length(xs), length(ys))

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=false, bnd_as_pits=false);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 6
    @test maximum(slen)==66
    @test maximum(area[.!isnan.(area)])==4727
    @test length(unique(c))==7
    @test sum(c) == 32498
    @test sum(diff(c[:])) == 0
    @test sort(unique(c))[[1,end]] ==[0,6]

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=false, bnd_as_pits=true);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 390
    @test maximum(slen)==65
    @test maximum(area[.!isnan.(area)])==4372
    @test length(unique(c))==391
    @test sum(c) == 2003347
    @test sum(diff(c[:])) == -1
    @test sort(unique(c))[[1,end]] ==[0,390]

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=false);
    #plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 1
    @test maximum(slen)==165
    @test maximum(area[.!isnan.(area)])==9605
    @test length(unique(c))==2
    @test sum(c) == 19605
    @test sum(diff(c[:])) == 0
    @test sort(unique(c))[[1,end]] ==[1,2]
    @test bnds isa Array{Array{CartesianIndex,1},1}

    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, drain_pits=true, bnd_as_pits=true);
    # plotarea_dem(xs, ys, dem, area, pits)
    @test length(pits) == 388
    @test maximum(slen)==81
    @test maximum(area[.!isnan.(area)])==4831
    @test length(unique(c))==389
    @test sum(c) == 1372917
    @test sum(diff(c[:])) == 98
    @test sort(unique(c))[[1,end]] ==[1,389]
    @test bnds isa Array{Array{CartesianIndex,1},1}
end

@testset "fill_dem" begin
    # TODO: need more fill tests
    dx = 0.1
    xs = -1.5:dx:1
    ys = -0.5:dx:3.0
    dem = dem1.(xs, ys', withpit=true)
    area, slen, dir, nout, nin, pits = WWF.waterflows(dem)
    demf = WWF.fill_dem(dem, pits, dir)
    @test sum(demf.-dem) ≈ 2.1499674517313414
    @test sum(demf.-dem .> 0) == 5
end


text:
using StableRNGs
using LinearAlgebra
using SparseArrays
using GraphRecipes
using Plots
using LightGraphs
using AbstractTrees

cd(@__DIR__)
include("parse_readme.jl")
include("functions.jl")
cd("../assets")

julia_logo_pun()
savefig("readme_julia_logo_pun.png")

random_labelled_graph()[1]
savefig("random_labelled_graph.png")

random_3d_graph()
savefig("random_3d_graph.png")

light_graphs()
savefig("light_graphs.png")

directed()
savefig("directed.png")

edgelabel()
savefig("edgelabel.png")

selfedges()
savefig("selfedges.png")

multigraphs()
savefig("multigraphs.png")

arc_chord_diagrams()
savefig("arc_chord_diagrams.png")

marker_properties()
savefig("marker_properties.png")

ast_example()
savefig("ast_example.png")

julia_type_tree()
savefig("julia_type_tree.png")

julia_dict_tree()
savefig("julia_dict_tree.png")

funky_edge_and_marker_args()
savefig("funky_edge_and_marker_args.png")

custom_nodeshapes_single()
savefig("custom_nodeshapes_single.png")

custom_nodeshapes_various()
savefig("custom_nodeshapes_various.png")


text:
"""
    PyTable

Abstract type for Python wrappers which may be interpretable as Tables.jl-compatible tables.

If `pyconvert(PyTable, x::Py)` is a table, then `x` is also a table.
"""
abstract type PyTable end
export PyTable

PyTable(x) = pyconvert(PyTable, x)

Tables.istable(x::Py) = Tables.istable(@pyconvert(PyTable, x, return false))
Tables.rowaccess(x::Py) = Tables.rowaccess(@pyconvert(PyTable, x, return false))
Tables.rows(x::Py) = Tables.rows(pyconvert(PyTable, x))
Tables.columnaccess(x::Py) = Tables.columnaccess(@pyconvert(PyTable, x, return false))
Tables.columns(x::Py) = Tables.columns(pyconvert(PyTable, x))
Tables.materializer(x::Py) = Tables.materializer(pyconvert(PyTable, x))


text:
module scVIDE

using StatsBase
using Conda
using PyCall
using Random
using Distributed

# Load python packages
const random = PyNULL()
const os = PyNULL()
const np = PyNULL()
const pd = PyNULL()
const scvi = PyNULL()
#const h5py = PyNULL()
const scvi_dataset = PyNULL()
const scvi_models = PyNULL()
const scvi_inference = PyNULL()
const torch = PyNULL()
const scvi_inference_autotune = PyNULL()
const hyperopt = PyNULL()

function __init__()
    copy!(random, pyimport("random"))
    copy!(os, pyimport("os"))
    copy!(np, pyimport("numpy"))
    copy!(pd, pyimport("pandas"))
    copy!(scvi, pyimport("scvi"))
    copy!(scvi_dataset, pyimport("scvi.dataset"))
    copy!(scvi_models, pyimport("scvi.models"))
    copy!(scvi_inference, pyimport("scvi.inference"))
    copy!(torch, pyimport("torch"))
    copy!(scvi_inference_autotune, pyimport("scvi.inference.autotune"))
    copy!(hyperopt, pyimport("hyperopt"))
end

#function __init__()
#    copy!(random, pyimport_conda("random","random"))
#    copy!(os, pyimport_conda("os","os"))
#    copy!(np, pyimport_conda("numpy","numpy"))
#    copy!(pd, pyimport_conda("pandas","pandas"))
#    copy!(scvi, pyimport_conda("scvi=0.6.0","scvi","bioconda"))
#    copy!(h5py, pyimport_conda("h5py=3.1","h5py"))
#    copy!(scvi_dataset, pyimport_conda("scvi.dataset"))
#    copy!(scvi_models, pyimport_conda("scvi.models"))
#    copy!(scvi_inference, pyimport_conda("scvi.inference"))
#    copy!(torch, pyimport_conda("torch"))
#    copy!(scvi_inference_autotune, pyimport_conda("scvi.inference.autotune"))
#    copy!(hyperopt, pyimport_conda("hyperopt"))
#end

export
    alternative_model_scVI,
    null_model_scVI,
    subsample_scVI_cells!,
    subsample_scVI_cells_idxs,
    jackstraw_scVI

include("scVIDE_algorithm.jl")

end # of module scVIDE


text:
export get_exp_dir, get_exp_path, get_bootstrap_paths
export get_last_dualcgp_paths, get_last_ind_path
export init_backup, fetch_backup
export find_yaml, cfg_from_exp_dir, log_from_exp_dir

export LOG_HEADER, NSGA2_LOG_HEADER

using CSV
using YAML
# using DataFrames
using Dates
using Plots

LOG_HEADER = ["date", "lib", "type", "gen_number", "best", "mean", "std"]
NSGA2_LOG_HEADER = ["date", "lib", "type", "gen_number", "rank", "fitness",
                    "normalized_fitness", "reached_frames", "dna_id"]

function init_backup(logid::String, resdir::String, cfg_path::String)
    #new_resu_dir = joinpath(resdir, "results", logid)
    #new_logs_dir = joinpath(new_resu_dir, "logs")
    #new_gens_dir = joinpath(new_resu_dir, "gens")
    #mkdir(new_resu_dir)
    #mkdir(new_logs_dir)
    #mkdir(new_gens_dir)
    #new_cfg_path = joinpath(new_resu_dir, cfg_path[length("cfg/")+1:end])
    expdir = joinpath(resdir, logid)
    mkpath(joinpath(expdir, "logs"))
    mkpath(joinpath(expdir, "gens"))
    new_cfg_path = joinpath(expdir, "logs", basename(cfg_path))
    cp(cfg_path, new_cfg_path, force=true)
end

logsdir_from_logid(logid::String, resdir::String) = joinpath(resdir, "results", logid, "logs")
gensdir_from_logid(logid::String, resdir::String) = joinpath(resdir, "results", logid, "gens")

"""
    fetch_backup(rootdir::String)

Fetch all the data in `logs/` and `gens/` and copy them to `results/`.
Assumes the function `init_backup` to be run previously to the evolution.
Corresponding files are found using unique logid.
"""
function fetch_backup(rootdir::String; clean::Bool=false)
    each_log_name = Array{String,1}()
    each_log_id = Array{String,1}()
    each_log_dir = Array{String,1}()
    each_gen_dir = Array{String,1}()
    each_log_new_name = Array{String,1}()
    logsdir = joinpath(rootdir, "logs")
    gensdir = joinpath(rootdir, "gens")
    for l in setdiff(readdir(logsdir), ["placeholder.txt"])
        is_csv = endswith(l, ".csv")
        id = is_csv ? l[1:end-length(".csv")] : l
        if is_csv
            push!(each_log_name, l)
            push!(each_log_id, id)
            push!(each_log_new_name, "controller.csv")
        else
            log_dir = joinpath(logsdir, l)
            push!(each_log_dir, log_dir)
            for p in readdir(log_dir)
                push!(each_log_name, joinpath(l, p))
                push!(each_log_id, id)
                push!(each_log_new_name, p)
            end
        end
    end
    # Move all logs
    for i in eachindex(each_log_name)
        old_logfile = joinpath(logsdir, each_log_name[i])
        new_logdir = logsdir_from_logid(each_log_id[i], rootdir)
        new_logfile = joinpath(new_logdir, each_log_new_name[i])
        mkpath(new_logdir)  # Make path if not created
        cp(old_logfile, new_logfile, force=true)
    end
    # Delete all log files
    if clean
        for d in each_log_dir
            rm(d, recursive=true, force=true)
        end
    end

    for id in setdiff(readdir(gensdir), ["placeholder.txt"])
        gen_dir = joinpath(gensdir, id)
        push!(each_gen_dir, gen_dir)
        for g in readdir(gen_dir)
            old_genpath = joinpath(gensdir, id, g)
            if startswith(g, "encoder_") || startswith(g, "controller_")
                new_filename = g
            else
                new_filename = string("controller_", g)
            end
            new_gendir = gensdir_from_logid(id, rootdir)
            mkpath(new_gendir)  # Make path if not created
            new_genpath = joinpath(new_gendir, new_filename)
            cp(old_genpath, new_genpath, force=true)
        end
    end
    # Delete all gen files
    if clean
        for d in each_gen_dir
            rm(d, recursive=true, force=true)
        end
    end
end

"""
    get_exp_dir(game_name::String)

Get a new experiment directory.
"""
function get_exp_path(game_name::String, resdir::String)
    exp_dir = string(Dates.now(), "_", game_name)
    joinpath(resdir, exp_dir)
end

"""
    function get_exp_dir(
        resdir::String;
        min_date::DateTime=DateTime(0),
        max_date::DateTime=DateTime(0),
        games::Array{String,1}=Array{String,1}(),
        reducers::Array{String,1}=Array{String,1}()
    )

Fetch the experiments directories corresponding to a time lapse, specific
games and specific reducers. If no game list / no reducers list / no dates are
provided, return everything.

Example:

    using Dates
    exp_dirs, games = exp_dir(
        resdir,
        min_date=DateTime(2021, 07, 12),
        max_date=DateTime(2021, 07, 15),
        games=["freeway"],
        reducers=["pooling"]
    )
"""
function get_exp_dir(
    resdir::String;
    min_date::DateTime=DateTime(0),
    max_date::DateTime=DateTime(0),
    games::Array{String,1}=Array{String,1}(),
    reducers::Array{String,1}=Array{String,1}()
)
    no_time = min_date == max_date == DateTime(0)
    no_specified_games = length(games) == 0
    no_specified_reducers = length(reducers) == 0
    existing_res = readdir(resdir)
    filtered_res = Array{String,1}()
    filtered_games = Array{String,1}()
    for exp_dir in existing_res
        splt = findnext('_', exp_dir, 1)
        exp_date = DateTime(exp_dir[1:splt-1])
        exp_game = exp_dir[splt+1:end]
        if no_time || (min_date < exp_date < max_date)
            if no_specified_games || (exp_game in games)
                exp_full_path = string(resdir, exp_dir)
                cfg = cfg_from_exp_dir(exp_full_path)
                if no_specified_reducers || (cfg["reducer"]["type"] in reducers)
                    push!(filtered_res, exp_full_path)
                    push!(filtered_games, exp_game)
                end
            end
        end
    end
    filtered_res, filtered_games
end

function are_same_cfg(cfg_a::NamedTuple, cfg_b::NamedTuple)
    if keys(cfg_a) != keys(cfg_b)
        false
    else
        for k in keys(cfg_a)
            if (k != :id) && (cfg_a[k] != cfg_b[k])
                false
            end
        end
    end
    true
end

"""
    get_last_ind_path(exp_dir::String, indname::String)

Given the experiment directory, return the path of the last saved individual.
The argument `indname` specifies the name of the individual, which should have
the following form: `string(indname, "_", max_gen_str)``.
"""
function get_last_ind_path(exp_dir::String, indname::String)
    gens_dir = joinpath(exp_dir, "gens")
    max_gen = 0
    max_gen_str = ""
    for g in readdir(gens_dir)
        if startswith(g, indname)
            gen_str = g[length(indname)+2:end]
            gen = parse(Int64, gen_str)
            if gen > max_gen
                max_gen = gen
                max_gen_str = gen_str
            end
        end
    end
    dna_dir = joinpath(gens_dir, string(indname, "_", max_gen_str))
    n_ind = length(readdir(dna_dir))
    dna_filename_length = maximum([length(p) for p in readdir(dna_dir)]) - length(".dna")
    dna_file = string(lpad(n_ind, dna_filename_length, "0"), ".dna")
    dna_path = joinpath(dna_dir, dna_file)
    dna_path
end

"""
    get_best_individuals_paths(exp_dir::String)

Given the experiment directory, return the paths of the last saved encoder and
controller.
"""
function get_last_dualcgp_paths(exp_dir::String)
    enco_dna_path = get_last_ind_path(exp_dir, "encoder")
    cont_dna_path = get_last_ind_path(exp_dir, "controller")
    enco_dna_path, cont_dna_path
    #=
    gens_dir = joinpath(exp_dir, "gens")
    n_gen = readdir(gens_dir)[end]  # Last subdir in alphabetical order
    n_gen = n_gen[length("encoder_")+1:end]  # Gen number as a String
    enco_dna_path = joinpath(exp_dir, "gens", string("encoder_", n_gen))
    cont_dna_path = joinpath(exp_dir, "gens", string("controller_", n_gen))
    n_indiv = length(readdir(enco_dna_path))
    dna_file_length = maximum([length(p) for p in readdir(enco_dna_path)]) - length(".dna")
    dna_file = string(lpad(n_indiv, dna_file_length, "0"), ".dna")
    enco_dna_path = joinpath(enco_dna_path, dna_file)
    cont_dna_path = joinpath(cont_dna_path, dna_file)
    enco_dna_path, cont_dna_path
    =#
end

"""
    get_bootstrap_paths(
        resdir::String,
        enco_cfg::NamedTuple,
        cont_cfg::NamedTuple,
        game::String
    )

Deprecated.
Get path to bootstrap individuals for re-starting evolution from a previously
reached stage.
"""
function get_bootstrap_paths(
    resdir::String,
    enco_cfg::NamedTuple,
    cont_cfg::NamedTuple,
    game::String
)
    enco_bs_path = resdir
    cont_bs_path = resdir

    dirs = readdir(resdir)
    filter!(e -> e[end-length("nobootstrap")+1:end] != "nobootstrap", dirs)

    for dir in dirs
        exp_dir = joinpath(resdir, dir)
        cfg_dir = cfg_from_exp_dir(exp_dir)
        _, ecfg, ccfg, _, _ = dualcgp_config(cfg_dir, game)
        println()
        println("________________________________________________")
        println(are_same_cfg(enco_cfg, ecfg))
        println(are_same_cfg(cont_cfg, ccfg))
        println("________________________________________________")
        println(exp_dir)
        println("________________________________________________")
    end
    enco_bs_path, cont_bs_path
end

"""
    find_yaml(path::String)

Find first `.yaml` file at given path.
"""
function find_yaml(path::String)
    for d in readdir(path)
        if endswith(d, ".yaml")
            return d
        end
    end
end

"""
    cfg_from_exp_dir(exp_dir::String)

Get configuration file at specified experiment directory.
"""
function cfg_from_exp_dir(exp_dir::String)
    cfg_dir = joinpath(exp_dir, "logs")
    yaml = find_yaml(cfg_dir)
    if yaml == nothing
        cfg_dir = exp_dir
        yaml = find_yaml(cfg_dir)
    end
    cfg_path = joinpath(cfg_dir, yaml)
    YAML.load_file(cfg_path)
end

"""
    log_from_exp_dir(
        exp_dir::String;
        log_file::String="logs/controller.csv",
        log_header::Vector{String}=LOG_HEADER
    )

Get log file at specified experiment directory.
Default: return the controller's log file.
"""
function log_from_exp_dir(
    exp_dir::String;
    log_file::String="logs/controller.csv",
    header::Union{Int64,Vector{String}}=LOG_HEADER,
    sep::String=","
)
    log_file = joinpath(exp_dir, log_file)
    CSV.File(log_file, header=header, delim=sep)
end


text:
using Pkg
using SafeTestsets

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")

function activate_downstream_env()
    Pkg.activate("GPU")
    Pkg.develop(PackageSpec(path=dirname(@__DIR__)))
    Pkg.instantiate()
end

if GROUP == "All" || GROUP == "Core"
    @safetestset "Dispatch" begin include("core_dispatch.jl") end
    @safetestset "ODE tests" begin include("core_odes.jl") end
    @safetestset "Resizing" begin include("core_resizing.jl") end
    @safetestset "Nested Duals" begin include("core_nesteddual.jl") end
end

if !is_APPVEYOR && GROUP == "GPU"
    activate_downstream_env()
    @safetestset "GPU tests" begin include("gpu_all.jl") end
end


text:
module BenchAvgFilter1DConstProp

import Random
using BenchmarkTools
using TapirBenchmarks

Random.seed!(1234)

const SUITE = BenchmarkGroup()

for n in [2^25, 2^27]
    SUITE[:n=>n] = s0 = BenchmarkGroup()
    ys0, xs = avgfilter1d_constprop_setup(n)
    s0[:impl=>:seq] = @benchmarkable avgfilter1d_constprop_seq!($(copy(ys0)), $xs)
    s0[:impl=>:threads] = @benchmarkable avgfilter1d_constprop_threads!($(copy(ys0)), $xs)
    s0[:impl=>:tapir_dac] = @benchmarkable avgfilter1d_constprop_tapir_dac!($(copy(ys0)), $xs)
    # s0[:impl=>:tapir_seq] = @benchmarkable avgfilter1d_constprop_tapir_seq!($(copy(ys0)), $xs)
end

end
BenchAvgFilter1DConstProp.SUITE


text:
#Compute reduced quantities for error checking without resorting to the "ncdiff" tool
function reductions(model, grid)
    mass = 0.0
    te = 0.0
    nx, nz = grid.nx, grid.nz
    dx, dz = grid.dx, grid.dz

    for k = 1:nz, i = 1:nx
        r = model.state[hs+i, hs+k, ID_DENS] + model.hy_dens_cell[hs+k]       # Density
        u = model.state[hs+i, hs+k, ID_UMOM] / r                           # U-wind
        w = model.state[hs+i, hs+k, ID_WMOM] / r                           # W-wind
        th = (model.state[hs+i, hs+k, ID_RHOT] + model.hy_dens_theta_cell[hs+k]) / r # Potential Temperature (theta)
        p = C0 * (r * th)^gamma      # Pressure
        t = th / (p0 / p)^(rd / cp)  # Temperature
        ke = r * (u * u + w * w)           # Kinetic Energy
        ie = r * cv * t                # Internal Energy
        mass = mass + r * dx * dz # Accumulate domain mass
        te = te + (ke + ie) * dx * dz # Accumulate domain total energy
        #println(r, u, w, th)
      end

    #    double glob[2], loc[2];
    #    loc[0] = mass;
    #    loc[1] = te;
    #    int ierr = MPI_Allreduce(loc,glob,2,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    #    mass = glob[0];
    #    te   = glob[1];
    return mass, te
end

function createSnapshot(model, grid)

    snapshot = Array{Float64}(undef, 4, grid.nx, grid.nz)
    #Store perturbed values in the temp arrays for output
    for k = 1:grid.nz
        for i = 1:grid.nx
            snapshot[1, i, k] =  model.state[i, hs+k, ID_DENS]
            snapshot[2, i, k] =
                model.state[i, hs+k, ID_UMOM] /
                (model.hy_dens_cell[hs+k] + model.state[i, hs+k, ID_DENS])
            snapshot[3, i, k] =
                model.state[i, hs+k, ID_WMOM] /
                (model.hy_dens_cell[hs+k] + model.state[i, hs+k, ID_DENS])
            snapshot[4, i, k] =
                (model.state[i, hs+k, ID_RHOT] + model.hy_dens_theta_cell[hs+k]) /
                (model.hy_dens_cell[hs+k] + model.state[i, hs+k, ID_DENS]) -
                model.hy_dens_theta_cell[hs+k] / model.hy_dens_cell[hs+k]
        end
    end

    return snapshot
end

#
#Output the fluid state (state) to a NetCDF file at a given elapsed model time (etime)
#The file I/O uses parallel-netcdf, the only external library required for this mini-app.
#If it's too cumbersome, you can comment the I/O out, but you'll miss out on some potentially cool graphics
function output(snapshots, grid, etimes, ncfile="output.nc")
  #Create new file
  isfile(ncfile) && rm(ncfile)
  _, nx, nz, _ = size(snapshots)
  nt = length(etimes)
  nccreate(ncfile, "dens", "x", collect(range(0,stop=grid.nx*grid.dx, length=grid.nx)), Dict("units"=>"m"),
                           "z", collect(range(0,stop=grid.nz*grid.dz, length=grid.nz)), Dict("units"=>"m"),
                           "time", etimes, Dict("units"=>"s"))
  nccreate(ncfile, "uwnd", "x", collect(range(0,stop=grid.nx*grid.dx, length=grid.nx)), Dict("units"=>"m"),
                           "z", collect(range(0,stop=grid.nz*grid.dz, length=grid.nz)), Dict("units"=>"m"),
                           "time", etimes, Dict("units"=>"s"))
  nccreate(ncfile, "wwnd", "x", collect(range(0,stop=grid.nx*grid.dx, length=grid.nx)), Dict("units"=>"m"),
                           "z", collect(range(0,stop=grid.nz*grid.dz, length=grid.nz)), Dict("units"=>"m"),
                           "time", etimes, Dict("units"=>"s"))
  nccreate(ncfile, "theta", "x", collect(range(0,stop=grid.nx*grid.dx, length=grid.nx)), Dict("units"=>"m"),
                           "z", collect(range(0,stop=grid.nz*grid.dz, length=grid.nz)), Dict("units"=>"m"),
                           "time", etimes, Dict("units"=>"s"))

                           #nccreate(ncfile, "uwnd", "x", nx, Dict("units"=>"m"), "z", nz, Dict("units"=>"m"), "time", nt, Dict("units"=>"s"))
  #nccreate(ncfile, "wwnd", "x", nx, Dict("units"=>"m"), "z", nz, Dict("units"=>"m"), "time", nt, Dict("units"=>"s"))
  #nccreate(ncfile, "theta", "x", nx, Dict("units"=>"m"), "z", nz, Dict("units"=>"m"), "time", nt, Dict("units"=>"s"))
  #exit()
  #Write the grid data to file with all the processes writing collectively
#  ncwrite(etimes, ncfile, "time")
  ncwrite(snapshots[1, :, :, 1:nt], ncfile, "dens")
  ncwrite(snapshots[2, :, :, 1:nt], ncfile, "uwnd")
  ncwrite(snapshots[3, :, :, 1:nt], ncfile, "wwnd")
  ncwrite(snapshots[4, :, :, 1:nt], ncfile, "theta")

end


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

mutable struct NetworkInterfacePropertiesFormat2 <: SwaggerModel
    virtualMachine::Any # spec type: Union{ Nothing, SubResource } # spec name: virtualMachine
    networkSecurityGroup::Any # spec type: Union{ Nothing, NetworkSecurityGroup } # spec name: networkSecurityGroup
    privateEndpoint::Any # spec type: Union{ Nothing, PrivateEndpoint } # spec name: privateEndpoint
    ipConfigurations::Any # spec type: Union{ Nothing, Vector{NetworkInterfaceIPConfiguration2} } # spec name: ipConfigurations
    tapConfigurations::Any # spec type: Union{ Nothing, Vector{NetworkInterfaceTapConfiguration2} } # spec name: tapConfigurations
    dnsSettings::Any # spec type: Union{ Nothing, NetworkInterfaceDnsSettings2 } # spec name: dnsSettings
    macAddress::Any # spec type: Union{ Nothing, String } # spec name: macAddress
    primary::Any # spec type: Union{ Nothing, Bool } # spec name: primary
    enableAcceleratedNetworking::Any # spec type: Union{ Nothing, Bool } # spec name: enableAcceleratedNetworking
    enableIPForwarding::Any # spec type: Union{ Nothing, Bool } # spec name: enableIPForwarding
    hostedWorkloads::Any # spec type: Union{ Nothing, Vector{String} } # spec name: hostedWorkloads
    resourceGuid::Any # spec type: Union{ Nothing, String } # spec name: resourceGuid
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function NetworkInterfacePropertiesFormat2(;virtualMachine=nothing, networkSecurityGroup=nothing, privateEndpoint=nothing, ipConfigurations=nothing, tapConfigurations=nothing, dnsSettings=nothing, macAddress=nothing, primary=nothing, enableAcceleratedNetworking=nothing, enableIPForwarding=nothing, hostedWorkloads=nothing, resourceGuid=nothing, provisioningState=nothing)
        o = new()
        validate_property(NetworkInterfacePropertiesFormat2, Symbol("virtualMachine"), virtualMachine)
        setfield!(o, Symbol("virtualMachine"), virtualMachine)
        validate_property(NetworkInterfacePropertiesFormat2, Symbol("networkSecurityGroup"), networkSecurityGroup)
        setfield!(o, Symbol("networkSecurityGroup"), networkSecurityGroup)
        validate_property(NetworkInterfacePropertiesFormat2, Symbol("privateEndpoint"), privateEndpoint)
        setfield!(o, Symbol("privateEndpoint"), privateEndpoint)
        validate_property(NetworkInterfacePropertiesFormat2, Symbol("ipConfigurations"), ipConfigurations)
        setfield!(o, Symbol("ipConfigurations"), ipConfigurations)
        validate_property(NetworkInterfacePropertiesFormat2, Symbol("tapConfigurations"), tapConfigurations)
        setfield!(o, Symbol("tapConfigurations"), tapConfigurations)
        validate_property(NetworkInterfacePropertiesFormat2, Symbol("dnsSettings"), dnsSettings)
        setfield!(o, Symbol("dnsSettings"), dnsSettings)
        validate_property(NetworkInterfacePropertiesFormat2, Symbol("macAddress"), macAddress)
        setfield!(o, Symbol("macAddress"), macAddress)
        validate_property(NetworkInterfacePropertiesFormat2, Symbol("primary"), primary)
        setfield!(o, Symbol("primary"), primary)
        validate_property(NetworkInterfacePropertiesFormat2, Symbol("enableAcceleratedNetworking"), enableAcceleratedNetworking)
        setfield!(o, Symbol("enableAcceleratedNetworking"), enableAcceleratedNetworking)
        validate_property(NetworkInterfacePropertiesFormat2, Symbol("enableIPForwarding"), enableIPForwarding)
        setfield!(o, Symbol("enableIPForwarding"), enableIPForwarding)
        validate_property(NetworkInterfacePropertiesFormat2, Symbol("hostedWorkloads"), hostedWorkloads)
        setfield!(o, Symbol("hostedWorkloads"), hostedWorkloads)
        validate_property(NetworkInterfacePropertiesFormat2, Symbol("resourceGuid"), resourceGuid)
        setfield!(o, Symbol("resourceGuid"), resourceGuid)
        validate_property(NetworkInterfacePropertiesFormat2, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type NetworkInterfacePropertiesFormat2

const _property_map_NetworkInterfacePropertiesFormat2 = Dict{Symbol,Symbol}(Symbol("virtualMachine")=>Symbol("virtualMachine"), Symbol("networkSecurityGroup")=>Symbol("networkSecurityGroup"), Symbol("privateEndpoint")=>Symbol("privateEndpoint"), Symbol("ipConfigurations")=>Symbol("ipConfigurations"), Symbol("tapConfigurations")=>Symbol("tapConfigurations"), Symbol("dnsSettings")=>Symbol("dnsSettings"), Symbol("macAddress")=>Symbol("macAddress"), Symbol("primary")=>Symbol("primary"), Symbol("enableAcceleratedNetworking")=>Symbol("enableAcceleratedNetworking"), Symbol("enableIPForwarding")=>Symbol("enableIPForwarding"), Symbol("hostedWorkloads")=>Symbol("hostedWorkloads"), Symbol("resourceGuid")=>Symbol("resourceGuid"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_NetworkInterfacePropertiesFormat2 = Dict{Symbol,String}(Symbol("virtualMachine")=>"SubResource", Symbol("networkSecurityGroup")=>"NetworkSecurityGroup", Symbol("privateEndpoint")=>"PrivateEndpoint", Symbol("ipConfigurations")=>"Vector{NetworkInterfaceIPConfiguration2}", Symbol("tapConfigurations")=>"Vector{NetworkInterfaceTapConfiguration2}", Symbol("dnsSettings")=>"NetworkInterfaceDnsSettings2", Symbol("macAddress")=>"String", Symbol("primary")=>"Bool", Symbol("enableAcceleratedNetworking")=>"Bool", Symbol("enableIPForwarding")=>"Bool", Symbol("hostedWorkloads")=>"Vector{String}", Symbol("resourceGuid")=>"String", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ NetworkInterfacePropertiesFormat2 }) = collect(keys(_property_map_NetworkInterfacePropertiesFormat2))
Swagger.property_type(::Type{ NetworkInterfacePropertiesFormat2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkInterfacePropertiesFormat2[name]))}
Swagger.field_name(::Type{ NetworkInterfacePropertiesFormat2 }, property_name::Symbol) =  _property_map_NetworkInterfacePropertiesFormat2[property_name]

function check_required(o::NetworkInterfacePropertiesFormat2)
    true
end

function validate_property(::Type{ NetworkInterfacePropertiesFormat2 }, name::Symbol, val)
end


text:
################################################################################

import ArgParse: ArgParseSettings, @add_arg_table!, parse_args

################################################################################

function shArgParser(args)
  # minimal argument parsing
  s = ArgParseSettings(description = "MindReader command line utility.")
  @add_arg_table! s begin
      "--file", "-f"
        arg_type             = String
        required             = true
        help                 = "`edf` file to read"
      "--indir", "-i"
        arg_type             = String
        required             = true
        help                 = "`edf` file directory"
      "--outdir", "-o"
        arg_type             = String
        required             = false
        default              = "."
        help                 = "output directory"
      "--outsvg", "-s"
        arg_type             = String
        required             = false
        help                 = "output directory svg. Defined as `file`.svg at `outdir` if not specified"
      "--outcsv", "-c"
        arg_type             = String
        required             = false
        help                 = "output directory csv. Defined as `file`.csv at `outdir` if not specified"
      "--window-size", "-w"
        nargs                = '?'
        arg_type             = Int
        default              = 256
        help                 = "window size along raw signal"
      "--bin-overlap", "-b"
        nargs                = '?'
        arg_type             = Int
        default              = 4
        help                 = "bin overlap along raw signal"
  end
  parsed_args = parse_args(s)
  return parsed_args
end

################################################################################

#  parse shell arguments
shArgs = shArgParser(ARGS)

################################################################################

# create directory if does not exist
if !isdir(shArgs["outdir"]) && shArgs["outdir"] != "."
  @info "Creating directory `$(shArgs["outdir"])`"
  mkdir(shArgs["outdir"])
end

################################################################################

max_stars_repo_name: DanielRivasMD/MindReader

max_issues_repo_name: DanielRivasMD/MindReader

max_forks_repo_name: DanielRivasMD/MindReader


text:
module RNG

using StableRNGs
using Colors
using Random

const STABLE_RNG = StableRNG(123)

rand(args...) = Base.rand(STABLE_RNG, args...)
randn(args...) = Base.randn(STABLE_RNG, args...)

seed_rng!() = Random.seed!(STABLE_RNG, 123)

function Base.rand(r::StableRNGs.LehmerRNG, ::Random.SamplerType{T}) where T <: ColorAlpha
    return T(Base.rand(r), Base.rand(r), Base.rand(r), Base.rand(r))
end

function Base.rand(r::StableRNGs.LehmerRNG, ::Random.SamplerType{T}) where T <: AbstractRGB
    return T(Base.rand(r), Base.rand(r), Base.rand(r))
end

end

using .RNG


text:
export MPIIODriver

import JSON3

# Version of internal MPIIO format.
# If the version is updated, it should match the upcoming PencilArrays version.
const MPIIO_VERSION = 0.3

"""
    MPIIODriver(; sequential = false, uniqueopen = false, deleteonclose = false)

MPI-IO driver using the MPI.jl package.

Keyword arguments are passed to
[`MPI.File.open`](https://juliaparallel.github.io/MPI.jl/latest/io/#MPI.File.open).

This driver writes binary data along with a JSON file containing metadata.
When reading data, this JSON file is expected to be present along with the raw
data file.
"""
Base.@kwdef struct MPIIODriver <: ParallelIODriver
    sequential    :: Bool = false
    uniqueopen    :: Bool = false
    deleteonclose :: Bool = false
end

const MetadataDict = Dict{Symbol,Any}
const DatasetKey = Symbol
const DatasetDict = Dict{DatasetKey,Any}

"""
    MPIFile

Wraps a `MPI.FileHandle`, also including file position information and metadata.

File position is updated when reading and writing data, and is independent of
the individual and shared file pointers defined by MPI.
"""
mutable struct MPIFile
    file       :: MPI.FileHandle
    comm       :: MPI.Comm
    filename   :: String
    meta       :: MetadataDict
    position   :: Int  # file position in bytes
    write_mode :: Bool
    MPIFile(file, comm, filename, meta; write) =
        new(file, comm, filename, meta, 0, write)
end

function MPIFile(comm::MPI.Comm, filename; kws...)
    flags, other_kws = keywords_to_open(; kws...)
    meta = if flags.write && !flags.append
        mpiio_init_metadata()
    else
        mpiio_load_metadata(filename_meta(filename))
    end
    file = MPIFile(
        MPI.File.open(comm, filename; kws...),
        comm, filename, meta, write=flags.write,
    )
    if flags.append
        # Synchronise position in file.
        pos = MPI.File.get_position_shared(parent(file))
        seek(file, pos)
    end
    file
end

mpiio_init_metadata() = MetadataDict(
    :driver => (type = "MPIIODriver", version = MPIIO_VERSION),
    :datasets => DatasetDict(),
)

function mpiio_load_metadata(filename)
    isfile(filename) || error("metadata file not found: $filename")
    meta = open(JSON3.read, filename, "r")
    # Convert from specific JSON3 type to Dict, so that datasets can be appended.
    MetadataDict(
        :driver => meta.driver,
        :datasets => DatasetDict(meta.datasets),
    )
end

function Base.close(ff::MPIFile)
    if should_write_metadata(ff)
        write_metadata(ff, metadata(ff))
    end
    close(parent(ff))
end

function write_metadata(ff, meta)
    MPI.Comm_rank(ff.comm) == 0 || return
    open(filename_meta(ff), "w") do io
        JSON3.write(io, meta)
        write(io, '\n')
    end
    nothing
end

filename_meta(fname) = string(fname, ".json")
filename_meta(ff::MPIFile) = filename_meta(get_filename(ff))
metadata(ff::MPIFile) = ff.meta
should_write_metadata(ff::MPIFile) = ff.write_mode
get_filename(ff::MPIFile) = ff.filename
Base.parent(ff::MPIFile) = ff.file
Base.position(ff::MPIFile) = ff.position
Base.skip(ff::MPIFile, offset) = ff.position += offset
Base.seek(ff::MPIFile, pos) = ff.position = pos

"""
    open([f::Function], driver::MPIIODriver, filename, comm::MPI.Comm; keywords...)

Open parallel file using the MPI-IO driver.

See [`open(::ParallelIODriver)`](@ref) for common options for all drivers.

Driver-specific options may be passed via the `driver` argument. See
[`MPIIODriver`](@ref) for details.

## Driver notes

- the `truncate` keyword is ignored.
"""
function Base.open(::MPIIODriver) end

Base.open(D::MPIIODriver, filename::AbstractString, comm::MPI.Comm; keywords...) =
    MPIFile(
        comm, filename;
        sequential=D.sequential, uniqueopen=D.uniqueopen,
        deleteonclose=D.deleteonclose, keywords...,
    )

"""
    setindex!(file::MPIFile, x::MaybePencilArrayCollection,
              name::AbstractString; chunks = false, collective = true, infokws...)

Write [`PencilArray`](@ref) to binary file using MPI-IO.

# Optional arguments

- if `chunks = true`, data is written in contiguous blocks, with one block per
  process.
  Otherwise, each process writes to discontiguous sections of disk, using
  `MPI.File.set_view!` and custom datatypes.
  Note that discontiguous I/O (the default) is more convenient, as it allows to
  read back the data using a different number or distribution of MPI processes.

- if `collective = true`, the dataset is written collectivelly. This is
  usually recommended for performance.

- when writing discontiguous blocks, additional keyword arguments are passed via
  an `MPI.Info` object to `MPI.File.set_view!`. This is ignored if `chunks = true`.

"""
function Base.setindex!(
        ff::MPIFile, x::MaybePencilArrayCollection, name::AbstractString;
        collective=true, chunks=false, kw...,
    )
    file = parent(ff)
    offset = position(ff)
    for u in collection(x)
        # TODO write all collection components at once (should be easier in the
        # discontiguous case)
        if chunks
            write_contiguous(file, u; offset=offset, collective=collective, kw...)
        else
            write_discontiguous(file, u; offset=offset, collective=collective, kw...)
        end
        offset += sizeof_global(u)
    end
    add_metadata(ff, x, name, chunks)
    skip(ff, sizeof_global(x))
    x
end

eltype_collection(x::PencilArray) = eltype(x)
eltype_collection(x::PencilArrayCollection) = eltype(first(x))

function add_metadata(file::MPIFile, x, name, chunks::Bool)
    meta = metadata(file)
    size_col = collection_size(x)
    size_log = size_global(x, LogicalOrder())
    size_mem = size_global(x, MemoryOrder())
    meta[:datasets][DatasetKey(name)] = (
        metadata(x)...,
        element_type = eltype_collection(x),
        dims_logical = (size_log..., size_col...),
        dims_memory = (size_mem..., size_col...),
        chunks = chunks,
        offset_bytes = position(file),
        size_bytes = sizeof_global(x),
    )
    nothing
end

"""
    read!(file::MPIFile, x::PencilArray, name::AbstractString;
          collective = true, infokws...)

Read binary data from an MPI-IO stream, filling in [`PencilArray`](@ref).

See [`setindex!`](@ref setindex!(::MPIFile)) for details on keyword arguments.
"""
function Base.read!(ff::MPIFile, x::MaybePencilArrayCollection, name::AbstractString;
                    collective=true, kw...)
    meta = get(metadata(ff)[:datasets], DatasetKey(name), nothing)
    meta === nothing && error("dataset '$name' not found")
    file = parent(ff)
    offset = meta.offset_bytes :: Int
    chunks = meta.chunks :: Bool
    check_metadata(x, meta.element_type, Tuple(meta.dims_memory), meta.size_bytes)
    chunks && check_read_chunks(x, meta.process_dims, name)
    for u in collection(x)
        if chunks
            read_contiguous!(file, u; offset=offset, collective=collective, kw...)
        else
            read_discontiguous!(file, u; offset=offset, collective=collective, kw...)
        end
        offset += sizeof_global(u)
    end
    x
end

function check_metadata(x, file_eltype, file_dims, file_sizeof)
    T = eltype_collection(x)
    if string(T) != file_eltype
        error("incompatible type of file and array: $file_eltype ≠ $T")
    end
    sz = (size_global(x, MemoryOrder())..., collection_size(x)...)
    if sz !== file_dims
        error("incompatible dimensions of dataset in file and array: $file_dims ≠ $sz")
    end
    @assert sizeof_global(x) == file_sizeof
    nothing
end

function check_read_chunks(x, pdims_file, name)
    pdims = size(topology(x))
    if length(pdims) != length(pdims_file) || any(pdims .!= pdims_file)
        error("dataset '$name' was written in chunks with a different MPI topology" *
              " ($pdims ≠ $pdims_file)")
    end
    nothing
end

function write_discontiguous(ff::MPI.FileHandle, x::PencilArray;
                             offset, collective, infokws...)
    to = timer(pencil(x))
    @timeit_debug to "Write MPI-IO discontiguous" begin
        set_view!(ff, x, offset; infokws...)
        A = parent(x)
        if collective
            MPI.File.write_all(ff, A)
        else
            MPI.File.write(ff, A)
        end
    end
    nothing
end

function read_discontiguous!(ff::MPI.FileHandle, x::PencilArray;
                             offset, collective, infokws...)
    to = timer(pencil(x))
    @timeit_debug to "Read MPI-IO discontiguous" begin
        set_view!(ff, x, offset; infokws...)
        A = parent(x)
        if collective
            MPI.File.read_all!(ff, A)
        else
            MPI.File.read!(ff, A)
        end
    end
end

function set_view!(ff, x::PencilArray, offset; infokws...)
    etype = MPI.Datatype(eltype(x))
    filetype = create_discontiguous_datatype(x, MemoryOrder())  # TODO cache datatype?
    MPI.File.set_view!(ff, offset, etype, filetype; infokws...)
    nothing
end

function create_discontiguous_datatype(x::PencilArray, order)
    sizes = size_global(x, order)
    subsizes = size_local(x, order)
    offsets = map(r -> first(r) - 1, range_local(x, order))
    oldtype = MPI.Datatype(eltype(x), commit=false)
    dtype = MPI.Types.create_subarray(sizes, subsizes, offsets, oldtype)
    MPI.Types.commit!(dtype)
    dtype
end

function write_contiguous(ff::MPI.FileHandle, x::PencilArray;
                          offset, collective, infokws...)
    to = timer(pencil(x))
    @timeit_debug to "Write MPI-IO contiguous" begin
        offset += mpi_io_offset(x)
        A = parent(x)
        if collective
            MPI.File.write_at_all(ff, offset, A)
        else
            MPI.File.write_at(ff, offset, A)
        end
    end
    nothing
end

function read_contiguous!(ff::MPI.FileHandle, x::PencilArray;
                          offset, collective, infokws...)
    to = timer(pencil(x))
    @timeit_debug to "Read MPI-IO contiguous" begin
        offset += mpi_io_offset(x)
        A = parent(x)
        if collective
            MPI.File.read_at_all!(ff, offset, A)
        else
            MPI.File.read_at!(ff, offset, A)
        end
    end
    nothing
end

function mpi_io_offset(x::PencilArray)
    topo = topology(pencil(x))
    # Linear index of this process in the topology.
    # (TODO This should be stored in MPITopology...)
    n = LinearIndices(topo)[coords_local(topo)...]
    off = 0
    for m = 1:(n - 1)
        r = range_remote(x, m)
        off += length(CartesianIndices(r))  # length of data held by remote process
    end
    T = eltype(x)
    off * sizeof(T)
end


text:
# Poisson

"""
poipmf!(lambda, prob; left = 0, right = length(prob)-1+left)
poipmf(lambda, right; left = 0)

Compute the p.m.f. of Poisson distribution with mean lambda.
`left` and `right` are the domain of Poisson distribution.
They should be choosen so that the total probability in the domain becomes 1.
In `poipmf!``, the p.m.f. is saved to the vector `prob`.
The right can be obtained from `rightbound`.
The retuen value is the normalizing constant so that the total sum of `prob` is 1, called the weight.
`poipmf` returns a tuple of (weight, prob).
"""

@origin (prob => left) function poipmf!(lambda::Tv, prob::Vector{Tv};
    left::Ti = 0, right::Ti = length(prob)-1+left) where {Tv, Ti}
    @inbounds begin
        log2piOver2::Tv = log(2*pi) / 2
        mode::Ti = floor(Ti, lambda)
        if mode >= 1
            prob[mode] = exp(-lambda + mode * log(lambda)
                - log2piOver2 - (mode + 1/2) * log(mode) + mode)
        else
            prob[mode] = exp(-lambda)
        end
        # down
        for j = mode:-1:left+1
            prob[j-1] = j / lambda * prob[j]
        end
        # up
        for j = mode:right-1
            prob[j+1] = lambda / (j+1) * prob[j]
        end
        # compute W
        weight::Tv = 0
        s::Ti = left
        t::Ti = right
        while s < t
            if prob[s] <= prob[t]
                weight += prob[s]
                s += 1
            else
                weight += prob[t]
                t -= 1
            end
        end
        weight += prob[s]
    end
end

function poipmf(lambda::Tv, right::Ti; left::Ti = 0) where {Tv, Ti}
    prob = Vector{Tv}(undef, right-left+1)
    weight = poipmf!(lambda, prob, left=left, right=right)
    (weight, prob)
end

"""
cpoipmf!(lambda, prob, cprob; left = 0, right = length(prob)-1+left)
cpoipmf(lambda, right; left = 0)

Compute the p.m.f. and complementary c.d.f. of Poisson distribution with mean lambda.
`left` and `right` are the domain of Poisson distribution.
They should be choosen so that the total probability in the domain becomes 1.
In `cpoipmf!``, the p.m.f. and c.d.f. are stored to `prob` and `cprob`, respectively.
The right can be obtained from `rightbound`.
The retuen value is the normalizing constant so that the total sum of `prob` is 1, called the weight.
`cpoipmf` returns a tuple of (weight, prob, cprob).
"""

@origin (prob => left, cprob => left) function cpoipmf!(lambda::Tv, prob::Vector{Tv}, cprob::Vector{Tv}; left::Ti = 0, right::Ti = length(prob)-1+left) where {Tv, Ti}
    weight::Tv = poipmf!(lambda, prob, left=left, right=right)
    @inbounds begin
        cprob[right] = 0
        for k = right-1:-1:left
            cprob[k] = cprob[k+1] + prob[k+1]
        end
        weight
    end
end

function cpoipmf(lambda::Tv, right::Ti; left::Ti = 0) where {Tv, Ti}
    prob = Vector{Tv}(undef, right-left+1)
    cprob = Vector{Tv}(undef, right-left+1)
    weight = cpoipmf!(lambda, prob, cprob, left=left, right=right)
    (weight, prob, cprob)
end

"""
rightbound(::Type{Ti} = Int, lambda::Tv, q::Tv = Tv(1.0e-8))

Compute the rightbound of Poisson distribution with mean lambda.
The rightbound finds a quantile so that the complementary c.d.f. becomes `q`.
"""

function rightbound(lambda::Tv, q::Tv = Tv(1.0e-8))::Int where {Tv}
    rightbound(Int, lambda, q)
end

function rightbound(::Type{Ti}, lambda::Tv, q::Tv = Tv(1.0e-8))::Ti where {Tv, Ti}
    z = cquantile(Normal(), q)
    if lambda < 3.0
        ll = exp(-lambda)
        total = ll
        right::Ti = 0
        while true
            right += 1
            ll *= lambda / right
            total += ll
            if total + q >= 1.0
                break
            end
        end
        right
    else
        right = floor(Ti, (z + sqrt(4.0 * lambda - 1.0))^2 / 4.0 + 1.0)
    end
end


text:
include("misc.jl")

function findMin(funObj,w;maxIter=100,epsilon=1e-2,derivativeCheck=false,verbose=true)
	# funObj: function that returns (objective,gradient)
	# w: initial guess
	# maxIter: maximum number of iterations
	# epsilon: stop if the gradient gets below this
	# derivativeCheck: whether to check against numerical gradient

	# Evalluate the intial objective and gradient
	(f,g) = funObj(w)

	# Optionally check if gradient matches finite-differencing
	if derivativeCheck
		g2 = numGrad(funObj,w)

		if maximum(abs.(g-g2)) > 1e-4
			@show([g g2])
			@printf("User and numerical derivatives differ\n")
			sleep(1)
		else
			@printf("User and numerical derivatives agree\n")
		end
	end

	# Initial step size and sufficient decrease parameter
	gamma = 1e-4
	alpha = 1
	for i in 1:maxIter

		# Try out the current step-size
		wNew = w - alpha*g
		(fNew,gNew) = funObj(wNew)

		# Decrease the step-size if we increased the function
		gg = dot(g,g)
		while fNew > f - gamma*alpha*gg

			if verbose
				@printf("Backtracking\n")
			end

			# Fit a degree-2 polynomial to set step-size
			alpha = alpha^2*gg/(2(fNew - f + alpha*gg))

			# Try out the smaller step-size
			wNew = w - alpha*g
			(fNew,gNew) = funObj(wNew)
		end

		# Guess the step-size for the next iteration
		y = gNew - g
		alpha *= -dot(y,g)/dot(y,y)

		# Sanity check on the step-size
		if (!isfinitereal(alpha)) | (alpha < 1e-10) | (alpha > 1e10)
			alpha = 1
		end

		# Accept the new parameters/function/gradient
		w = wNew
		f = fNew
		g = gNew

		# Print out some diagnostics
		gradNorm = norm(g,Inf)
		if verbose
			@printf("%6d %15.5e %15.5e %15.5e\n",i,alpha,f,gradNorm)
		end

		# We want to stop if the gradient is really small
		if gradNorm < epsilon
			if verbose
				@printf("Problem solved up to optimality tolerance\n")
			end
			return w
		end
	end
	if verbose
		@printf("Reached maximum number of iterations\n")
	end
	return w
end

function findMinL1(funObj,w,lambda;maxIter=100,epsilon=1e-2)
	# funObj: function that returns (objective,gradient)
	# w: initial guess
	# lambda: value of L1-regularization parmaeter
	# maxIter: maximum number of iterations
	# epsilon: stop if the gradient gets below this

	# Evalluate the intial objective and gradient
	(f,g) = funObj(w)

	# Initial step size and sufficient decrease parameter
	gamma = 1e-4
	alpha = 1
	for i in 1:maxIter

		# Gradient step on smoooth part
		wNew = w - alpha*g
		# Proximal step on non-smooth part
		wNew = sign.(wNew).*max.(abs.(wNew) - lambda*alpha,0)
		(fNew,gNew) = funObj(wNew)

		# Decrease the step-size if we increased the function
		gtd = dot(g,wNew-w)
		while fNew + lambda*norm(wNew,1) > f + lambda*norm(w,1) - gamma*alpha*gtd
			@printf("Backtracking\n")
			alpha /= 2

			# Try out the smaller step-size
			wNew = w - alpha*g
			wNew = sign.(wNew).*max.(abs.(wNew) - lambda*alpha,0)
			(fNew,gNew) = funObj(wNew)
		end

		# Guess the step-size for the next iteration
		y = gNew - g
		alpha *= -dot(y,g)/dot(y,y)

		# Sanity check on the step-size
		if (!isfinitereal(alpha)) | (alpha < 1e-10) | (alpha > 1e10)
			alpha = 1
		end

		# Accept the new parameters/function/gradient
		w = wNew
		f = fNew
		g = gNew

		# Print out some diagnostics
		optCond = norm(w-sign.(w-g).*max.(abs.(w-g)-lambda,0),Inf)
		@printf("%6d %15.5e %15.5e %15.5e\n",i,alpha,f+lambda*norm(w,1),optCond)

		# We want to stop if the gradient is really small
		if optCond < epsilon
			@printf("Problem solved up to optimality tolerance\n")
			return w
		end
	end
	@printf("Reached maximum number of iterations\n")
	return w
end


text:
using Test, EngThermBase

# EngTherm top-level tests
include("abstract.test.jl")


text:
module PlayingRoster

include("shared.jl")

function get_roster(game_id::String)::String
    A = game_id[1:4] # Year start = 2014
    B = parse(Int16, game_id[1:4]) + 1 # year start + 1
    C = game_id[5:end] #non-date game id
    request = HTTP.request("GET", "http://www.nhl.com/scores/htmlreports/$(A)$(B)/RO$C.HTM")
    response = String(request.body)
    return response
end

function get_coaches(response::String)::Dict{Symbol, String}
    html = parsehtml(response)
    coaches = eachmatch(Selector("tr#HeadCoaches"), html.root)
    if coaches == nothing
        return coaches
    end

    tds = eachmatch(Selector("td"), coaches[1])

    return Dict{Symbol, String}(
        :away => uppercase(nodeText(tds[2])),
        :home => uppercase(nodeText(tds[4]))
    )
end

"""
    fix_name(player::Array)
    Get rid of (A) or (C) when a player has it attached to their name

    :param player: list of player info -> [number, position, name]

    :return: fixed list
"""
function fix_name(player::Array{String})::Array{String}
    fixed = player
    fixed[3] = replace(player[3], "(A)" => "")
    fixed[3] = replace(fixed[3], "(C)" => "")
    fixed[3] = strip(fixed[3])
    return fixed
end

function get_players(response::String)::Dict{Symbol, Array}
    html = parsehtml(response)

    tables = eachmatch(Selector("td.border > table"), html.root)
    player_info = [eachmatch(Selector("td"), table) for table in tables]
    player_info = [[nodeText(x) for x in group] for group in player_info]
    # Make list of list of 3 each. The three are: number, position, name (in that order)
    player_info = [[group[i:i+2] for i in 1:3:length(group) - 3] for group in player_info]
    # Get rid of header column
    player_info = [[player for player in group if player[1] != "#"] for group in player_info]

    for i in eachindex(player_info)
        for j in eachindex(player_info[i])
            # Need to use strings here as the array is typed to Strings
            was_scratched = i == 3 || i == 4 ? "true" : "false"
            push!(player_info[i][j], was_scratched)
        end
    end

    players = Dict{Symbol, Array}(:away => player_info[1], :home => player_info[2])

    # Merge scratches back into player dataframe
    # Sometimes scratches aren't included
    if length(player_info) > 2
        append!(player_info[3], players[:away])
        append!(player_info[4], players[:home])
    end

    # Get rid when just whitespace
    # For those with (A) or (C) in name field get rid of it
    # First condition is to control when we get whitespace as one of the indices
    players[:away] = [p[1] != "\xa0" ? fix_name(p) : p for p in players[:away]]
    players[:home] = [p[1] != "\xa0" ? fix_name(p) : p for p in players[:home]]
    players[:away] = [p for p in players[:away] if p[1] != "\xa0"]
    players[:home] = [p for p in players[:home] if p[1] != "\xa0"]

    return players
end

function scrape_roster_and_coaches(game_id::AbstractString)::Tuple{Dict{Symbol, Array}, Dict{Symbol, String}}
    response = get_roster(game_id)
    players = get_players(response)
    coaches = get_coaches(response)
    return players, coaches
end
#
# response = get_roster("2016020426")
# print(get_coaches(response))
# get_players(response)
#
# players, coaches = get_content(response)

end


text:
# This file is a part of FaceCraker. License is MIT

using Documenter, FaceCraker

makedocs(
	modules = [FaceCraker],
	sitename = "FaceCraker.jl",
	pages = Any[
	  "index.md"
	],
	versions = ["v#.#", "dev" => "dev"],
	assets = [""],
)

deploydocs(
	repo = "github.com/fetaxyu/FaceCraker.jl",
)


text:
"""
    Convenience function to build JuMP model.
"""
function build_model(fileNm::String; startT::Int=1, T::Int=-1, verb::Bool=true,
    cold_storage::Float64=90.e3, hot_storage::Float64=600.,
    carbonPrice::Float64=100., capCostTh::Bool=false, capCostMa::Bool=false,
    elecStorage::Bool=false, elec_storage::Float64=0.,nHRC::Float64=3.,
    nChiller::Float64=4., nHeater::Float64=3., dc_ratio::Float64=1.,
    startDate::Date=Date(), endDate::Date=Date(),
    umLoadsL2::Bool=true)

    if T != -1
        data = CSV.read(fileNm, datarow=1+startT, rows=T, allowmissing=:none)
    else
        data = CSV.read(fileNm, datarow=1+startT, allowmissing=:none)
    end

    data[:Date] = DateTime.(data[:Date])

    p = Sesi.design_base!()
    Sesi.design_water!(p)
    p[:coldTankMax] = cold_storage
    p[:cwTankInit] = 0.1 * cold_storage
    p[:hotTankMax] = hot_storage
    p[:hwTankInit] = 0.1 * hot_storage
    p[:HRCCAP] *= nHRC / 3.
    p[:chillerCAP] *= nChiller / 4.
    p[:heaterCAP] *= nHeater / 3.
    p[:umLoadsL2] = umLoadsL2

    Sesi.design_prices!(p)
    p[:demandCharge] = Dict(k => v*dc_ratio for (k,v) in p[:demandCharge])
    Sesi.design_carbon!(p)
    if elecStorage
        Sesi.design_storage!(p)
        p[:elecStorageMax] = elec_storage
        p[:elecStorageInit] = elec_storage
    end

    f = Sesi.forecast_base!(data, p; startT=startT, T=T, startDate=startDate,
        endDate=endDate)
    Sesi.forecast_water!(f, data, p)
    Sesi.forecast_prices!(f, data, p)
    Sesi.forecast_carbon!(f, data, p)
    f[:carbonPrice] = carbonPrice
    if elecStorage
        Sesi.forecast_storage!(f, data, p)
    end

    sesi = Sesi.model_base!(f, p)
    if capCostTh
        Sesi.model_cap_thstorage!(sesi, f, p)
    end
    if capCostMa
        Sesi.model_cap_machines!(sesi, f, p)
    end
    Sesi.model_water!(sesi, f, p)
    Sesi.model_water_switching!(sesi, f, p)
    if elecStorage
        Sesi.model_storage!(sesi, f, p)
    end
    Sesi.model_demand!(sesi, f, p)
    Sesi.model_energy_price!(sesi, f, p)
    Sesi.model_carbon!(sesi, f, p)
    Sesi.model_finalize!(sesi, f, p)

    return sesi, p, f
end

function carbon_cost(res, sc, base, curr)
    ((res[sc][curr]["agg_dollar"]-res[sc][base]["agg_dollar"])
        /(res[sc][base]["agg_carbon"]-res[sc][curr]["agg_carbon"]))
end

function carbon_reduction(res, sc, base, curr)
    100*(1-res[sc][curr]["CEP_carbon"]/res[sc][base]["CEP_carbon"])
end

function sanity(res, sc, base, curr)
    ((res[sc][curr]["CEP_carbon"]-res[sc][base]["CEP_carbon"])/
    (res[sc][curr]["agg_carbon"]-res[sc][base]["agg_carbon"]))
end

max_stars_repo_name: jdechalendar/sesi

max_issues_repo_name: jdechalendar/sesi

max_forks_repo_name: jdechalendar/sesi


text:
using Breadth
using Test
using SafeTestsets

@testset "Breadth.jl" begin
    @time @safetestset "Basics" begin include("Types.jl") end
end


text:
using Documenter, BioSequences

makedocs(
    format = Documenter.HTML(),
    sitename = "BioSequences.jl",
    pages = [
        "Home" => "index.md",
        "Biological Symbols" => "symbols.md",
        "Sequence types" => [
            "Overview" => "sequences/sequences.md",
            "BioSequence" => "sequences/bioseq.md",
            "Reference Sequences" => "sequences/refseq.md",
            "Nucleic acid k-mers" => "sequences/kmer.md"
        ],
        "IO" => [
            "FASTA formatted files" => "io/fasta.md",
            "FASTQ formatted files" => "io/fastq.md",
            "2bit formatted files" => "io/twobit.md"
        ],
        "Searching" => "searching.md",
        "Sequence Composition" => "composition.md",
        "Demultiplexing" => "demultiplexer.md",
        "Contributing" => "contributing.md"
    ],
    authors = "D. C. Jones, Kenta Sato, Ben J. Ward, The BioJulia Organisation and other contributors."
)

deploydocs(
    repo = "github.com/BioJulia/BioSequences.jl.git",
    julia = "1.0",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
__precompile__()

module Ranges

include("linspace.jl")

end # module


text:
#ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля
#РЕЗУЛЬТАТ: Робот - в исходном положении, и клетки поля промакированы так: нижний ряд - полностью,
#следующий - весь, за исключением одной последней клетки на Востоке, следующий - за исключением двух
#последних клеток на Востоке, и т.д.

function kray(r,side) #функция перемещает робота к краю
    while isborder(r,side) == false
        move!(r,side)
    end
end

function moverobot(r,side)
    while isborder(r,side) == false
        move!(r,side)
        putmarker!(r)
    end
end

function moveback(r,side)
    while ismarker(r) == true
        move!(r,side)
    end
    putmarker!(r)
end

function lesenka(r)
    kray(r,Sud)
    kray(r,West)
    putmarker!(r)
    moverobot(r,Ost)
    while isborder(r,Nord) == false
        move!(r,Nord)
        move!(r,West)
        moverobot(r,West)
        moveback(r,Ost)
    end
end

max_stars_repo_name: gardend/mirea-progs

max_issues_repo_name: gardend/mirea-progs

max_forks_repo_name: gardend/mirea-progs


text:
using ITensorGPU, Test, CUDA

CUDA.allowscalar(false)
@testset "ITensorGPU.jl" begin
  #@testset "$filename" for filename in ("test_cucontract.jl",)
  #  println("Running $filename with autotune")
  #  cmd = `$(Base.julia_cmd()) -e 'using Pkg; Pkg.activate(".."); Pkg.instantiate(); include("test_cucontract.jl")'`
  #  run(pipeline(setenv(cmd, "CUTENSOR_AUTOTUNE" => 1); stdout=stdout, stderr=stderr))
  #end
  @testset "$filename" for filename in (
    "dmrg.jl",
    "test_cuitensor.jl",
    "test_cudiag.jl",
    "test_cudense.jl",
    "test_cucontract.jl",
    "test_cumpo.jl",
    "test_cumps.jl",
    "test_cuiterativesolvers.jl",
    "test_cutruncate.jl",
  )
    println("Running $filename")
    include(filename)
  end
end


text:
using LazyArrays, FillArrays, Test
import LazyArrays: CachedArray, CachedMatrix, CachedVector

@testset "Cache" begin
    @testset "basics" begin
        A = 1:10
        C = cache(A)
        @test size(C) == (10,)
        @test axes(C) == (Base.OneTo(10),)
        @test all(Vector(C) .=== Vector(A))
        @test cache(C) isa CachedArray{Int,1,Vector{Int},UnitRange{Int}}
        C2 = cache(C)
        @test C2.data !== C.data
        @test C[:] == C[Base.Slice(Base.OneTo(10))] == C

        @test cache(A)[2,1] == 2
        @test_throws BoundsError cache(A)[2,2]

        A = reshape(1:10^2, 10,10)
        C = cache(A)
        @test size(C) == (10,10)
        @test axes(C) == (Base.OneTo(10),Base.OneTo(10))
        @test all(Array(C) .=== Array(A))
        @test cache(C) isa CachedArray{Int,2,Matrix{Int},typeof(A)}
        C2 = cache(C)
        @test C2.data !== C.data
        @test C[1:2,1:2] == [1 11; 2 12]

        A = reshape(1:10^3, 10,10,10)
        C = cache(A)
        @test size(C) == (10,10,10)
        @test axes(C) == (Base.OneTo(10),Base.OneTo(10),Base.OneTo(10))
        @test all(Array(C) .=== Array(A))

        A = reshape(1:10^3, 10,10,10)
        C = cache(A)
        LazyArrays.resizedata!(C,5,5,5)
        LazyArrays.resizedata!(C,8,8,8)
        @test all(C.data .=== Array(A)[1:8,1:8,1:8])

        @test C[1:3,1,:] == A[1:3,1,:]
    end

    @testset "Matrix cache" begin
        A = collect(1:5)
        C = cache(A)
        @test C isa Vector{Int}
        C[1] = 2
        @test A[1] ≠ 2

        A = cache(Matrix(reshape(1:6,2,3)))
        C = cache(A)
        @test C isa Matrix{Int}
        C[1,1] = 2
        @test A[1,1] ≠ 2
        C[1] = 3
        @test A[1,1] ≠ 3
    end

    @testset "setindex!" begin
        A = (1:5)
        C = cache(A)
        C[1] = 2
        @test C[1] == 2

        A = cache(reshape(1:6,2,3))
        C = cache(A)
        C[1,1] = 2
        @test A[1,1] ≠ 2
        @test C[1,1] == 2
        C[1] = 3
        @test C[1,1] == 3

        @test_throws BoundsError C[3,1]
        @test_throws BoundsError C[7]
    end

    @testset "colsupport past size" begin
        C = cache(Zeros(5,5)); C[5,1];
        @test colsupport(C,1) == Base.OneTo(5)
        @test colsupport(C,3) == 1:0
        @test rowsupport(C,1) == Base.OneTo(1)
    end

    @testset "broadcast" begin
        x = CachedArray([1,2,3],1:8);
        y = 1:8;
        f = (x,y) -> cos(x*y)
        @test f.(x,y) isa CachedArray
        @test @inferred(broadcast(f,x,y)) == f.(Vector(x), Vector(y))

        @test (x + y) isa CachedArray
        @test (x + y).array isa AbstractRange
        @test (x + y) == Vector(x) + Vector(y)

        z = CachedArray([1,4],Zeros{Int}(8));
        @test (x .+ z) isa CachedArray
        @test (x + z) isa CachedArray
        @test Vector( x .+ z) == Vector( x + z) == Vector(x) + Vector(z)

        # Lazy mixed with Static treats as Lazy
        s = SVector(1,2,3,4,5,6,7,8)
        @test f.(x , s) isa CachedArray
        @test f.(x , s) == f.(Vector(x), Vector(s))
    end

    @testset "padded CachedVector getindex" begin
        v = CachedArray([1,2,3],Zeros{Int}(1000))
        @test v[3:100] == [3; zeros(Int,97)]
        @test v[3:100] isa CachedArray{Int,1,Vector{Int},Zeros{Int,1,Tuple{Base.OneTo{Int}}}}
        @test v[4:end] isa CachedArray{Int,1,Vector{Int},Zeros{Int,1,Tuple{Base.OneTo{Int}}}}
        @test all(iszero,v[4:end])
        @test isempty(v[4:0])
        v = CachedArray([1,2,3],Fill{Int}(1,1000))
        @test v[3:100] == [3; ones(97)]
        @test norm(v) == norm(Array(v)) == norm(v,2)
        @test norm(v,1) == norm(Array(v),1)
        @test norm(v,Inf) == norm(Array(v),Inf)
        @test norm(v,3) ≈ norm(Array(v),3)
    end

    @testset "ambiguity broadcast" begin
        c = cache(1:100)
        v = Vcat([1,2,3],0:96)
        z = Zeros(100)
        @test v .+ c == c .+ v == Array(c) + Array(v)
        @test z .+ c == c .+ z == Array(c)
    end

    @testset "Fill" begin
        c = cache(Fill(1,10))
        @test c[2:3] isa LazyArrays.CachedVector{Int,Vector{Int},Fill{Int,1,Tuple{Base.OneTo{Int}}}}
        @test c[[2,4,6]] isa LazyArrays.CachedVector{Int,Vector{Int},Fill{Int,1,Tuple{Base.OneTo{Int}}}}
    end

    @testset "linalg" begin
        c = cache(Fill(3,3,3))
        @test fill(2,1,3) * c == fill(18,1,3)
        @test ApplyMatrix(exp,fill(3,3,3)) * c == exp(fill(3,3,3)) * fill(3,3,3)
        @test BroadcastMatrix(exp,fill(3,3,3)) * c == exp.(fill(3,3,3)) * fill(3,3,3)
        @test fill(2,3)' * c == fill(18,1,3)
        @test fill(2,3,1)' * c == fill(18,1,3)
    end

    @testset "broadcast" begin
        c = cache(Fill(3.0,10)); c[1] = 2;
        @test exp.(c) isa typeof(c)
        @test exp.(c) == exp.(Vector(c))
        @test c .+ 1 isa typeof(c)
        @test c .+ 1 == Vector(c) .+ 1
        @test 1 .+ c isa typeof(c)
        @test 1 .+ c == 1 .+ Vector(c)
        @test c .+ Ref(1) isa typeof(c)
        @test Ref(1) .+ c isa typeof(c)
    end

    @testset "*" begin
        A = randn(2,2)
        B = ApplyMatrix(exp,A)
        C = BroadcastMatrix(exp,A)
        x = cache(Fill(3,2))
        @test A*x ≈ A*Vector(x)
        @test B*x ≈ Matrix(B)*Vector(x)
        @test C*x ≈ Matrix(C)*Vector(x)
        @test A'x ≈ Matrix(A)'Vector(x)
        @test B'x ≈ Matrix(B)'Vector(x)
        @test C'x ≈ Matrix(C)'Vector(x)
    end
end


text:
module NoveltyLazy

using Scruff
using Scruff.SFuncs
using Scruff.Algorithms
import Scruff.Operators: expectation

struct NoveltySetup
    known_sfs::Vector{Dist{Float64}}
    known_probs::Vector{Float64}
    novelty_prob::Float64
    novelty_prior_mean::Float64
    novelty_prior_sd::Float64
    novelty_sd::Float64
end

function novelty_network(setup::NoveltySetup, numobs::Int)::InstantNetwork
    known = Cat(setup.known_sfs, setup.known_probs)()(:known)
    is_novel = Flip(setup.novelty_prob)()(:is_novel)
    novelty_mean = Normal(setup.novelty_prior_mean, setup.novelty_prior_sd)()(:novelty_mean)
    novelty = Det(Tuple{Float64}, Dist{Float64}, m -> Normal(m[1], setup.novelty_sd))()(:novelty)
    behavior = If{Dist{Float64}}()()(:behavior)

    variables = [known, is_novel, novelty_mean, novelty, behavior]
    graph = VariableGraph(novelty => [novelty_mean], behavior => [is_novel, novelty, known])
    for i in 1:numobs
        obs = Generate{Float64}()()(obsname(i))
        push!(variables, obs)
        graph[obs] = [behavior]
    end

    return InstantNetwork(variables, graph)
end

obsname(i) = Symbol("obs", i)

function do_experiment(setup, obs)
    net = novelty_network(setup, length(obs))
    is_novel = get_node(net, :is_novel)
    novelty_mean = get_node(net, :novelty_mean)
    evidence = Dict{Symbol, Score}()
    for (i,x) in enumerate(obs)
        evidence[obsname(i)] = HardScore(x)
    end

    alg = LSFI([is_novel, novelty_mean]; start_size = 5, increment = 5)
    runtime = Runtime(net)
    prepare(alg, runtime, evidence)

    for i = 1:10
        println("Range size: ", alg.state.next_size)
        refine(alg, runtime)
        is_novel_lb = probability_bounds(alg, runtime, is_novel, [false, true])[1]
        println("Probability of novel = ", is_novel_lb[2])
        println("Posterior mean of novel behavior = ", mean(alg, runtime, novelty_mean))
    end
end

function setup(generation_sd::Float64, prob_novel::Float64)::NoveltySetup
    known = [Normal(0.0, generation_sd), Normal(generation_sd, generation_sd)]
    return NoveltySetup(known, [0.75, 0.25], prob_novel, 0.0, 10.0, generation_sd)
end
setup1 = setup(1.0, 0.1)
setup2 = setup(4.0, 0.1)
obs = [5.0, 6.0, 7.0, 8.0, 9.0]

println("Lazy Inference")
println("Narrow generation standard deviation")
do_experiment(setup1, obs)
println("\nBroad generation standard deviation")
do_experiment(setup2, obs)

end