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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 87e78056-fecc-4b15-a723-16f51fa19353
using DataFrames, PlutoUI,HypertextLiteral, Statistics

# ╔═╡ 0ebe3384-97aa-11ec-17b5-cbf23568df03
md" # 2 Mathematical characterization of systems

As a first approximatiomn we can focus on a steady state characterization.

## 2.1 Ecological Systems

In fact, diversity is related to robustness and adaptation, because if some part o the system perform different tasks, the sayste itself may respond better to sudden environmental changes and laugh changes, as occurs in ecosystems

How do we measure this diversity? is a system composed of parts that are almost identical a diverse system? Is a bag of apples diverse? Depends of what are we looking for, Therefore the measure of diversity is going to depend on the the of system we are dealing with.

Lets say we want to measure size diversity in a population of bacteria, We can just use standard deviation, Large std, very diverse, small std, no diversity. But lets say now that we focus on the genetic information of a system composed on foxes and  rabbits. ,, there is no way we can quantify the difference in the genes  based on standard deviation (standard deviation of what, the length of the genome,

"

# ╔═╡ 943c55a4-00d2-4cb7-bba0-a6c46b30e0fa
md"Suppose a biologist wants to measure the diversity of species in a local forest. She collects the following data: "

# ╔═╡ 8ad237fa-8a09-4bae-80f4-ba1298c09581
begin
	dog_slide = @bind 🐶 html"<input type=range min=0 max=100 step=1>"
	mouse_slide = @bind 🐭 html"<input type=range min=0 max=100 step=1>"
	cat_slide = @bind 😺 html"<input type=range min=0 max=100 step=1>"
	fox_slide = @bind 🦊 html"<input type=range min=0 max=100 step=1>"
	bear_slide = @bind 🐻 html"<input type=range min=0 max=100 step=1>"

	md"""
	**How many molecules do you have?**

	Amount of dogs: $(dog_slide)

	Amount of dogs: $(mouse_slide)

	Amount of cats: $(cat_slide)

	Amount of foxes: $(fox_slide)

	Amount of bears: $(bear_slide)
	"""
end

# ╔═╡ bee1dd58-c8e6-4fa6-8eae-7cd4cca0a466
values_animals=DataFrame("Species" => ["Dog", "Mouse", "Cat",  "Fox",  "Bear"],"Frequency" => [🐶, 🐭, 😺, 🦊, 🐻])

# ╔═╡ 68713a25-a962-47e2-8830-c5e5c0511b43
md"Next, she can calculate the total number of organisms"

# ╔═╡ 279b1402-f4ba-4a3b-9d7d-82ac81b515e4
 total_number_of_organisms=sum(values_animals.Frequency)

# ╔═╡ 4190734c-4c6e-4662-9be9-2c19c4956a38
md"A first measure is the Richness $R$, i.e., how many species do we have "

# ╔═╡ 9981448a-edee-435f-9d2b-ec3271c4016a
R=size(values_animals.Frequency,1)

# ╔═╡ 0dbc209e-91dd-4569-afe9-b6407bdab995
md" For ecosystems, or systems composed of a number of individuals from species, we often use the Dominance Index:

```math
D=\frac{\sum n_i(n_i-1)}{N(N-1)}
```

where $n_i$ is the number of organisms in a species i, and N is the total number of organisms

"

# ╔═╡ c3606f18-d301-44d8-8765-7746bb6493c0
N=sum(values_animals.Frequency)

# ╔═╡ 34a569b4-3f88-4b5f-928b-0edd40c4d6ef
md"the Dominace index is"

# ╔═╡ 3989ead3-2760-443f-80e8-21ac961e2ff5
D = (sum((values_animals.Frequency).* ((values_animals.Frequency) .- 1)) / (N*(N-1)))

# ╔═╡ 7085e737-82db-4879-9c2e-96c6f1f146f0
md"Another highly used idex is the Simpson’s Index of Diversity, which is simply

```math
S=1-D=1-\frac{\sum n_i(n_i-1)}{N(N-1)}
```

"

# ╔═╡ 7bfac697-d85b-4a60-b90a-2a69c1155f5c
S = 1 - D

# ╔═╡ c71435ea-4bca-4bf7-ba90-576563385b33
md" Species __evenness__ refers to how close in numbers each species in an environment is. Mathematically it is defined as a diversity index, a measure of biodiversity which quantifies how equal the community is numerically. So if there are 40 foxes and 1000 dogs, the community is not very even.

Simpson evenness is closely related to Simpson diversity. It is actually just the Inverse Simpson Diversity with the richness added to the denominator, to mitigate the effect of species count on the range of the metric.
```math
E = \frac{1}{D * R}
```
"

# ╔═╡ 8f31e31d-1080-4b22-b508-b71bf30c2ec0
E= 1/ D / R

# ╔═╡ aa33422d-a8cb-46ed-8e11-d70cfd7caf36
@htl("""

Exercise: Measure the Richness, the dominance, teh diversity and the evenness of the followin systems. Discuss and compare the values obtained between the three systems:

- Which one is more diverse?

- Which one is more even?

- which one is more rich?

""")

# ╔═╡ 721ad120-8999-4f9f-9bd3-b9d966246436
values_animals_2=DataFrame("Species" => ["A", "B", "C",  "D",  "E", "F", "G", "H", "I", "J"],"Frequency" => [10,9, 11, 10, 8, 12, 10,11,10,9])

# ╔═╡ 8de6461d-e9ba-46df-a347-a7648b6f6e6a
values_animals_3=DataFrame("Species" => ["A", "B", "C",  "D",  "E", "F", "G", "H", "I", "J"],"Frequency" => [72,9, 11, 10, 8, 12, 10,11,10,9])

# ╔═╡ 8a2c17b6-3f5c-4f9f-9714-09df2944dd50
values_animals_4=DataFrame("Species" => ["A", "B", "C",  "D",  "E", "F", "G", "H", "I", "J"],"Frequency" => [33,33, 34, 0, 0, 0, 0, 0,0,0])

# ╔═╡ b20cdd19-08c4-4837-947f-97594b8822d2
md" ## 2.2 Graph Theory
As we mentioned, the important feature of the system comes from the interactions between the parts. Now we will see some tools to study how the parts of a system are interacting and connected.

- Proteins can act as transcription factors for other genes, forming complex networks that ultimately set the identity of a cell.

- When some node of the network changes (for instance, due to mutations), the behavior of a cell can change (leading to genetic diseases, such as cancer).

- Therefore, to understand and treat genetic diseases, we need  to know the genes involved, how they interact with other genes, and how they ultimately affect cellular responses.

"

# ╔═╡ 73f4fd5b-4611-4d90-8f05-af7b279dd86a
md" ### - Type of Networks
The structure of a network determines many relevant aspects of its function.
Network architecture is established by topological and statistical analyses.
Networks can be studied based on its global properties, (size, connectivity, robustness…)

Depending of the nature of the links between nodes ewe can classify networks in:

- Undirected: links between nodes have no direction (Facebook)
- Directed: links between nodes have direction (Twitter)
- Directed with self-links: Feedback loops (Gene regulatory networks)

[![EC1wSj.md.png](https://iili.io/EC1wSj.md.png)](https://freeimage.host/i/EC1wSj)
"

# ╔═╡ bfa2aecf-b7a8-4d68-aa1b-82d29032bd70
tree_url = "https://er.yuvayana.org/wp-content/uploads/sites/11/2016/02/Binary-Tree-network-model.jpg"

# ╔═╡ 40f52fcb-4be1-40c1-8720-ee48433473f3
md" ### - Size of Networks
The size of a network can refer to the number of nodes N or the number of links E.
The number of links is minimal in a tree network

$(Resource(tree_url))

"

# ╔═╡ af0276a5-322a-4a7f-9457-7e58524c5cbd
md"The minimal number of links E in a network of N nodes $E_{min}$ = N-$1$
The maximum number of links depends of the type of network
For undirected

```math
E_{max}=\binom{N}{2}=\frac{N!}{2!(N-2)!)}=\frac{N(N-1)(N-2)...1}{2\cdot 1 \cdot (N-2)(N-3)...1}=\frac{N(N-1)}{2}
```

For directed

```math
\begin{align*}
E_{max}=2\binom{N}{2}=N+\frac{2 \cdot N!}{2!(N-2)!)}=\\
\\
E_{max}=\frac{2 \cdot N(N-1)(N-2)...1}{2\cdot 1 \cdot (N-2)(N-3)...1}=\\
\\
E_{max}=N(N-1)
\end{align*}
```

For directed with self-links

```math
\begin{align*}
E_{max}=N+2\binom{N}{2}=N+\frac{2 \cdot N!}{2!(N-2)!)}=\\
\\
E_{max}=N+\frac{2 \cdot N(N-1)(N-2)...1}{2\cdot 1 \cdot (N-2)(N-3)...1}=\\
\\
E_{max}=N+N(N-1)=N(1+N-1)=N^2
\end{align*}
```

"

# ╔═╡ d4c87e16-8dff-4818-9bfc-df89c372a498
md" ### - Density of Networks

The density D of a network of N nodes is defined as a ratio of the number of links E to the number of possible links $E_{max}$

```math
\begin{align*}
D=\frac{E-(N-1)}{E_{max}-(N-1)}
\end{align*}
```

"

# ╔═╡ 477ba0d8-1dd9-46c8-bff0-72e4bfa89238
md"Exercise, draw two networks (one directed and one undirected) and calculate the density, which one is mode dense? "

# ╔═╡ b5174a24-428e-4a41-9b85-fccc7bcf8d32
path_url = "https://cdncontribute.geeksforgeeks.org/wp-content/uploads/exampleFigure-1.png"

# ╔═╡ 6551021e-e6d4-4d34-af80-efd8a7b3010d
md" ### - Average shortest path

In a real network like the Internet, a short average path length facilitates the quick transfer of information and reduces costs

- $(Resource(path_url))

Input: source node = $0$ and destination node is = $7$.

Output: Shortest path length is:$2$ -> [0 3 7]

Input: source vertex is = $2$ and destination vertex is = $6$

Output: Shortest path length is:$5$ -> [2 1 0 3 4 6]
"

# ╔═╡ 7fd4ffe4-6143-43f6-b75c-c76eb831108e
@htl("""

Exercise: Measure the Richness, the dominance, teh diversity and the evenness of the followin systems. Discuss and compare the values obtained between the three systems:

""")

# ╔═╡ ac539912-861a-4779-bac6-79c001177faa
Distance_to_1=DataFrame("Species" => ["1->2", "1->0", "1->3",  "1->7",  "1->4", "1->6", "1->5"],"Distance" => [1,1,2,3,3,4,4])

# ╔═╡ f11db50b-6aa9-4b38-a1c5-7fc217c55e95
md"The average shortest distance to node 1 is $(mean(Distance_to_1.Distance))"

# ╔═╡ f5673e80-8f26-4e1b-adfc-e9ec84f18c5b
@htl("""

Exercise: Measure the Richness, the dominance, teh diversity and the evenness of the followin systems. Discuss and compare the values obtained between the three systems:

""")

# ╔═╡ 4ace57bb-d9be-4681-8f7e-61982bb34b1c
Distance_to_2=DataFrame("Species" => ["2->1", "2->0", "2->3",  "2->7",  "2->4", "2->5", "2->6"],"Distance" => [1,2,3,4,4,5,5])

# ╔═╡ d4d3ca6e-8a3b-42fd-afd3-0a6d1c40a58e
md"The average shortest distance to node 2 is $(mean(Distance_to_2.Distance))"

# ╔═╡ 45a3e3e1-53c0-4632-8eab-b67f691a0ac8
@htl("""

""")

# ╔═╡ 9179e620-b8ff-4559-b864-c9729a3a42ec
Distance_to_0=DataFrame("Species" => ["0->2", "0->1", "0->3",  "0->7",  "0->4", "0->5", "0->6"],"Distance" => [2,1,1,2,2,3,3])

# ╔═╡ 9c2a3452-8c7f-4f59-877f-a1c585107e0c
md"The average shortest distance to node 2 is $(mean(Distance_to_0.Distance)). Central nodes have a shortest distance to everyone else in teh network than nodes in the boundaries. It is likely that information travels through them."

# ╔═╡ 76ece02b-8170-487b-8e0b-8e32b579b131
md"## Adjacency Matrix: Matrix Representation of Networks

A = binary square matrix of $N$*$N$ elements.
$A_{ij}$ = 1 if there exists an interaction between $i$ and $j$
$A_{ij}$ = 0 if not.

- undirected graph
```math
A_{ij}  = A_{ji}, A_{ii} = 0
```
- directed graph
```math
A_{ij} \ne A_{ji}, A_{ii} = 0,
```

- directed graph with self-links
```math
A_{ij} \ne A_{ji}, A_{ii} \ne 0,
```

Example 1: Draw and characterize the network that correspodns to the following graph,:

```math
\begin{equation*}
A_{i,j} =
\begin{pmatrix}
0 & 1 & 0 & 1 \\
1 & 0 & 1 & 0 \\
0  & 1  & 0 & 1  \\
1 & 0 & 1 &0
\end{pmatrix}
\end{equation*}
```

Example 2: Draw and characterize the the network that correspodns to the following graph:

```math
\begin{equation*}
A_{i,j} =
\begin{pmatrix}
0 & 1 & 0 & 0 \\
0 & 0 & 0 & 1 \\
1  & 0  & 0 & 1  \\
0 & 1 & 0 & 0
\end{pmatrix}
\end{equation*}
```

Example 3: Draw and characterize the the network that correspodns to the following graph:

```math
\begin{equation*}
A_{i,j} =
\begin{pmatrix}
0 & 1 & 0 & 1 \\
1 & 0 & 1 & 1 \\
0  & 1  & 0 & 0  \\
1 & 1 & 0 & 1
\end{pmatrix}
\end{equation*}
```

"

# ╔═╡ ddd62a2f-b19d-41b7-8b56-3c279f9c61ab
md" ### - Degree of Networks
Another more direct measure of how connect are the different nodes is the number of links they received. The degree $k$ of a node is the number of connections the node has to other nodes. If a network is directed, such as GRN, nodes have two different degrees:
- in-degree $k_{in}$: the number of incoming edges, and the
- out-degree $k_{out}$:  the number of outgoing edges
- total-degree $k_{tot}$ = $k_{in}$ + $k_{out}$:

[![EC1wSj.md.png](https://iili.io/EC1wSj.md.png)](https://freeimage.host/i/EC1wSj)

_Undirected_:
- degree of node 1 is 2
- degree of node 2 is 2
- degree of node 3 is 3
- degree of node 4 is 1

In an undirected network the total number of links, E, can be expressed as the sum of the node degrees

```math
\begin{align*}
E=\frac{1}{2}\sum_{i=1}^{N} k_i
\end{align*}
```

the 1/2 factor corrects for the fact that in the sum (2.1) each link is counted twice. For example, the link connecting the nodes 2 and 3 will be counted once in the degree of node 2 and once in the degree of node 3

_Directed_:

- In-degree of node 1 is 0
- In-degree of node 2 is 2
- Out-degree of node 1 is 2
- Out degree of node 2 is 0

In a directed network the total number of links, E, can be expressed as the sum of all node in-degrees or out-degrees

```math
\begin{align*}
E=\sum_{i=1}^{N} k_{in} = \sum_{i=1}^{N} k_{out}
\end{align*}
```

The __in__ degrees are:

```math
\begin{align*}
k^{in}_1 = 0\\
k^{in}_2 = 2\\
k^{in}_3 = 2\\
k^{in}_4 = 1
\end{align*}
```

The _out_ degrees are:
```math
\begin{align*}
k^{out}_1 = 2\\
k^{out}_2 = 0\\
k^{out}_3 = 2\\
k^{out}_4 = 1
\end{align*}
```

So, for this particular network:
```math
\begin{align*}
E=\sum_{i=1}^{N} k_{in} = \sum_{i=1}^{N} k_{out} =5
\end{align*}
```
"

# ╔═╡ 2db33033-cb6c-4e5b-b952-c9f06c98def0
md"Using the adjacency matrix, the degree is the sum of the row

For undirected:
```math
\begin{align*}
k_1=\sum_{i=1}^{N} A_{ij}
\end{align*}
```

For directed
```math
\begin{align*}
k_{1}^{out}=\sum_{i=1}^{N} A_{ij}\\
k_{1}^{in}=\sum_{j=1}^{N} A_{ij}
\end{align*}
```

For weighted
```math
\begin{align*}
k_{1}^{out}=\sum_{i=1}^{N} A_{ij}\\
k_{1}^{in}=\sum_{j=1}^{N} A_{ij}
\end{align*}
```
"

# ╔═╡ 979d1f55-870f-424c-a488-71a369c1fd8b
degree_distribution_url = "https://www.researchgate.net/profile/Daniele-Condorelli/publication/242112667/figure/fig3/AS:669311584710678@1536587725639/An-important-graph-property-is-the-degree-distribution-function-P-k-that-describes_W640.jpg"

# ╔═╡ b41588c9-c52c-4e85-ad1d-2775d143bd06
md" ## Average degree

The average number of links per node. Gives you a basic measurement of how the network is connected in average.

In undirected networks:

```math
\begin{align*}
\left \langle k \right \rangle = \frac{k_1 + k_2+ k_3 + ... k_N}{N}=\frac{1}{N}\sum_{i=1}^{N}k_i=\frac{2E}{N}
\end{align*}
```

In directed networks:

```math
\begin{align*}
\left \langle k^{in}_i \right \rangle =\frac{1}{N}\sum_{i=1}^{N}k^{in}_i=\left \langle k^{out}_i \right \rangle =\frac{1}{N}\sum_{i=1}^{N}k^{out}_i=\frac{E}{N}
\end{align*}
```

using the notation of the adjacency matrix:
```math
\begin{align*}
\left \langle k \right \rangle =  \frac{ \sum_{i=j}^{N} \sum_{i=1}^{N} A_{ij}} {N}
\end{align*}
```

Exercise: find de average degree of the folliwng networks:

[![EC1wSj.md.png](https://iili.io/EC1wSj.md.png)](https://freeimage.host/i/EC1wSj)

For the undirected: $\left \langle k \right \rangle = (2+2+3+1)/ 4 = 2$

For the directed: $\left \langle k_{in} \right \rangle = (2+2+3+1)/ 4 = 2$

"

# ╔═╡ 010b6033-1546-478a-8a9f-c26cf3f194fe
@htl("""

""")

# ╔═╡ a44b48e7-5497-4d49-a4d7-64508ecb5c3a
adjacency_url="https://www.ebi.ac.uk/training/online/courses/network-analysis-of-protein-interaction-data-an-introduction/wp-content/uploads/sites/64/2020/08/new-fig-4.png"

# ╔═╡ f938eda7-cefc-4ea3-991a-094b918d0430
md" $(Resource(adjacency_url))"

# ╔═╡ f9f7a01d-fbcf-4e5d-9cab-152eff279e5a
types_of_networks_url="https://www.researchgate.net/publication/343839881/figure/fig2/AS:928266965774337@1598327499590/Four-types-of-networks-in-the-scale-free-network-the-white-and-striped-nodes-represent_W640.jpg"

# ╔═╡ d580711c-2aff-404c-9aa9-84d56759eac9
md" ## Degree Distribution

The degree distribution $p_k$ provides the probability that a randomly selected node in the network has degree $p_{k}$.
Since $p_{k}$ is a probability, it must be normalized, i.e.

```math
\begin{align*}
\sum_{k=1}^{N} p_k =1
\end{align*}
```

P(k) is the probability distribution of the degrees over the whole network.  pk of a network is defined as the fraction of nodes in the network with degree $k_i$. Thus if there are N nodes in total in a network and M of them have degree k, we have $p_{k}$ = M/N.

The degree distribution is probably the most important feature of a network, since it determines many properties of how the network acts. So, the number of nodes in a network with degree k can be obtained by:

```math
\begin{align*}
N_k = N \cdot p_k
\end{align*}
```

And the average degree can be also calculated as:

```math
\begin{align*}
\left \langle k \right \rangle =\frac{1}{N}\sum_{i=1}^{N}k_i=\sum_{k=0}^{\infty }k \cdot p_k
\end{align*}
```
We can characterize the networks based on the degree distribution:
$(Resource(types_of_networks_url))

- __Regular networks and trees__ (usually man-made, have the lowest heterogeneity (e.g. the number of connections each node has is more or less the same)

- __small-world networks__ higher clustering and almost the same average path than the random networks with the same number of nodes and edges. high modularity, high herogeneity

- __Scale-free networks__ a highly heterogeneous degree distribution (power-law degree distribution)

- __Random networks__, nodes paired with uniform probability (low heterogeneity , since most nodes have the same number of connections), the degree distribution will be a Gaussian bell-shaped curve.

$(Resource(degree_distribution_url))

The typical plot of the degree distribution of a random network as a maximum, which correspond to the characteristic scale of the network. For example, in this network, the characteristic scale is around 8, since there is a higher number of nodes with that number of connections.

Biological gene regulatory networks have a degree distribution very different from random. There is no characteristic scale (no node with a more common number of connections, no maximum in the degree distribution plot). Many nodes poorly connected, very few nodes (hubs) highly connected.

Degree distributions are often represented in log-log plot. These types of networks are called power-law or scale free networks

Some nonbiological networks also show a degree distribution that is close to scale free. Above we plot the degree distribution of some of them.

(A) Actor collaboration graph with N = 212,250,〈k〉 = 28.78.
(B) Small sample of the World Wide Web, N = 325,729, 〈k〉 = 5.46.
(C) Power grid data, N = 4941, 〈k〉 = 2.67

What does it mean that so many real networks have a power law (scale free) degree distribution?
As opposed to random networks, scale-free networks do not have a characteristic scale, meaning that there is no typical node in the network that represents the degree for the other nodes.

### - 2.3 Conclusion:
The study of Networks is a very big, active and interesting field. The characteristics of networks of interactions are quite useful to understand key properties of large systems where their parts are interacting (big data…). Biological networks are a perfect systems to study using this formalism of network theory. It is well known that most of biologically occurring networks  have a degree distribution close to scale free, characterized by the presence of large hubs:

These power law distributions have the same functional form at all scales.
Why so many different systems have the same scale free degree distribution?

This is quite difficult question to try to explain in such a basic course. But the current agreement is that Scale Free networks often apear in  in systems where you have two properties combined:

 They grow in size (new nodes and links added). They started from small networks and became large by adding new nodes (as oppose to networks were all the links are established at the same time, i.e., a group in WhatsApp)
There is some sort of “preferential attachment”, meaning that there is some popular nodes that every one wants to get information from (rich gets richer)

Both social networks and GRN have these two properties.

"

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[compat]
DataFrames = "~1.3.2"
HypertextLiteral = "~0.9.3"
PlutoUI = "~0.7.38"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.42.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "ae02104e835f219b8930c7664b8012c93475c340"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.3.2"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "670e559e5c8e191ded66fa9ea89c97f10376bb4c"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.38"

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "28ef6c7ce353f0b35d0df0d5930e0d072c1f5b9b"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.1"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "dfb54c4e414caa595a1f2ed759b160f5a3ddcba5"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.1"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
git-tree-sha1 = "5ce79ce186cc678bbb5c5681ca3379d1ddae11a1"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.7.0"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─0ebe3384-97aa-11ec-17b5-cbf23568df03
# ╟─943c55a4-00d2-4cb7-bba0-a6c46b30e0fa
# ╠═87e78056-fecc-4b15-a723-16f51fa19353
# ╟─8ad237fa-8a09-4bae-80f4-ba1298c09581
# ╟─bee1dd58-c8e6-4fa6-8eae-7cd4cca0a466
# ╟─68713a25-a962-47e2-8830-c5e5c0511b43
# ╟─279b1402-f4ba-4a3b-9d7d-82ac81b515e4
# ╟─4190734c-4c6e-4662-9be9-2c19c4956a38
# ╠═9981448a-edee-435f-9d2b-ec3271c4016a
# ╟─0dbc209e-91dd-4569-afe9-b6407bdab995
# ╠═c3606f18-d301-44d8-8765-7746bb6493c0
# ╟─34a569b4-3f88-4b5f-928b-0edd40c4d6ef
# ╠═3989ead3-2760-443f-80e8-21ac961e2ff5
# ╟─7085e737-82db-4879-9c2e-96c6f1f146f0
# ╠═7bfac697-d85b-4a60-b90a-2a69c1155f5c
# ╟─c71435ea-4bca-4bf7-ba90-576563385b33
# ╠═8f31e31d-1080-4b22-b508-b71bf30c2ec0
# ╠═aa33422d-a8cb-46ed-8e11-d70cfd7caf36
# ╟─721ad120-8999-4f9f-9bd3-b9d966246436
# ╟─8de6461d-e9ba-46df-a347-a7648b6f6e6a
# ╟─8a2c17b6-3f5c-4f9f-9714-09df2944dd50
# ╟─b20cdd19-08c4-4837-947f-97594b8822d2
# ╟─73f4fd5b-4611-4d90-8f05-af7b279dd86a
# ╟─40f52fcb-4be1-40c1-8720-ee48433473f3
# ╟─bfa2aecf-b7a8-4d68-aa1b-82d29032bd70
# ╟─af0276a5-322a-4a7f-9457-7e58524c5cbd
# ╟─d4c87e16-8dff-4818-9bfc-df89c372a498
# ╟─477ba0d8-1dd9-46c8-bff0-72e4bfa89238
# ╟─b5174a24-428e-4a41-9b85-fccc7bcf8d32
# ╟─6551021e-e6d4-4d34-af80-efd8a7b3010d
# ╟─7fd4ffe4-6143-43f6-b75c-c76eb831108e
# ╟─ac539912-861a-4779-bac6-79c001177faa
# ╟─f11db50b-6aa9-4b38-a1c5-7fc217c55e95
# ╟─f5673e80-8f26-4e1b-adfc-e9ec84f18c5b
# ╟─4ace57bb-d9be-4681-8f7e-61982bb34b1c
# ╟─d4d3ca6e-8a3b-42fd-afd3-0a6d1c40a58e
# ╟─45a3e3e1-53c0-4632-8eab-b67f691a0ac8
# ╟─9179e620-b8ff-4559-b864-c9729a3a42ec
# ╟─9c2a3452-8c7f-4f59-877f-a1c585107e0c
# ╟─76ece02b-8170-487b-8e0b-8e32b579b131
# ╟─ddd62a2f-b19d-41b7-8b56-3c279f9c61ab
# ╟─2db33033-cb6c-4e5b-b952-c9f06c98def0
# ╟─979d1f55-870f-424c-a488-71a369c1fd8b
# ╟─b41588c9-c52c-4e85-ad1d-2775d143bd06
# ╟─010b6033-1546-478a-8a9f-c26cf3f194fe
# ╟─a44b48e7-5497-4d49-a4d7-64508ecb5c3a
# ╟─f938eda7-cefc-4ea3-991a-094b918d0430
# ╟─f9f7a01d-fbcf-4e5d-9cab-152eff279e5a
# ╟─d580711c-2aff-404c-9aa9-84d56759eac9
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: 2 Networks as Complex Systems.jl

max_stars_repo_name: davidgmiguez/PlutoNotebooks

max_issues_repo_path: 2 Networks as Complex Systems.jl

max_issues_repo_name: davidgmiguez/PlutoNotebooks

max_forks_repo_path: 2 Networks as Complex Systems.jl

max_forks_repo_name: davidgmiguez/PlutoNotebooks


text:
const SEATCODES = Dict('F' => '0', 'B' => '1', 'L' => '0', 'R' => '1')

function seatid(line)
  parse(Int, map(c -> SEATCODES[c], line); base=2)
end

seatids = BitSet(seatid.(eachline("a5.txt")))

maxseatid = maximum(seatids)
#part1
println("$maxseatid")

#part2
println("$(maximum(setdiff(BitSet(1:maxseatid), seatids)))")

ext: jl


text:
module RocketTakeOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: take()" begin

    println("Testing: operator take()")

    run_proxyshowcheck("Take", take(1))

    run_testset([
        (
            source = from(1:5) |> take(3),
            values = @ts([ 1:3, c ])
        ),
        (
            source = from(1:5) |> take(0),
            values = @ts(c)
        ),
        (
            source = from(1:5) |> async(0) |> take(3),
            values = @ts([ 1 ] ~ [ 2 ] ~ [ 3, c ])
        ),
        (
            source = timer(100, 30) |> take(3),
            values = @ts(100 ~ [ 0 ] ~ 30 ~ [ 1 ] ~ 30 ~ [ 2, c ])
        ),
        (
            source = completed() |> take(10),
            values = @ts(c)
        ),
        (
            source = faulted("e") |> take(10),
            values = @ts(e("e"))
        ),
        (
            source = never() |> take(10),
            values = @ts()
        )
    ])

end

end


text:
function test_location_scale(
    rng::Union{AbstractRNG, Missing},
    μ::Real, σ::Real, ρ::UnivariateDistribution, dref::UnivariateDistribution,
)
    d = LocationScale(μ,σ,ρ)
    @test params(d) == (μ,σ,ρ)
    @test eltype(d) === eltype(dref)

    # Different ways to construct the LocationScale object
    if dref isa DiscreteDistribution
        # floating point division introduces numerical errors
        # Better: multiply with rational numbers
        d_dict = Dict(
            "original" => d,
            "sugar" => σ * ρ + μ,
        )
    else
        d_dict = Dict(
            "original" => d,
            "sugar" => σ * ρ + μ,
            "composed" =>  μ + ((2 * ρ) * σ - 2) / 2 + 1
        )
    end

    @test d == deepcopy(d)

    #### Support / Domain

    @testset "Support" begin
        function test_support(d, dref)
            @test minimum(d) == minimum(dref)
            @test maximum(d) == maximum(dref)
            @test extrema(d) == (minimum(d), maximum(d))
        end
        @testset "$k" for (k,dtest) in d_dict
            test_support(dtest, dref)
        end
    end

    #### Promotions and conversions

    @testset "Promotions and conversions" begin
        function test_promotions_and_conversions(d, dref)
            @test typeof(d.µ) === typeof(d.σ)
            @test location(d) ≈ μ atol=1e-15
            @test    scale(d) ≈ σ atol=1e-15
        end
        @testset "$k" for (k,dtest) in d_dict
            test_promotions_and_conversions(dtest, dref)
        end
    end

    #### Statistics

    @testset "Statistics" begin
        function test_statistics(d, dref)
            @test mean(d) ≈ mean(dref)
            @test median(d) ≈ median(dref)
            @test mode(d) ≈ mode(dref)
            @test modes(d) ≈ modes(dref)

            @test var(d) ≈ var(dref)
            @test std(d) ≈ std(dref)
            @test skewness(d) ≈ skewness(dref)
            @test kurtosis(d) ≈ kurtosis(dref)

            @test isplatykurtic(d) == isplatykurtic(dref)
            @test isleptokurtic(d) == isleptokurtic(dref)
            @test ismesokurtic(d) == ismesokurtic(dref)

            @test entropy(d) ≈ entropy(dref)
            @test mgf(d,-0.1) ≈ mgf(dref,-0.1)
        end
        @testset "$k" for (k,dtest) in d_dict
            test_statistics(dtest, dref)
        end
    end

    #### Evaluation & Sampling

    @testset "Evaluation & Sampling" begin
        function test_evaluation_and_sampling(rng, d, dref)
            xs = rand(dref, 5)
            x = first(xs)
            insupport(d, x) == insupport(dref, x)
            # might return `false` for discrete distributions
            insupport(d, -x) == insupport(dref, -x)

            @test pdf(d, x) ≈ pdf(dref, x)
            @test pdf.(d, xs) ≈ pdf.(dref, xs)
            @test logpdf(d, x) ≈ logpdf(dref, x)
            @test logpdf.(d, xs) ≈ logpdf.(dref, xs)
            @test loglikelihood(d, x) ≈ loglikelihood(dref, x)
            @test loglikelihood(d, xs) ≈ loglikelihood(dref, xs)

            @test cdf(d, x) ≈ cdf(dref, x)
            @test logcdf(d, x) ≈ logcdf(dref, x)
            @test ccdf(d, x) ≈ ccdf(dref, x) atol=1e-15
            @test logccdf(d, x) ≈ logccdf(dref, x) atol=1e-15

            @test quantile(d,0.1) ≈ quantile(dref,0.1)
            @test quantile(d,0.5) ≈ quantile(dref,0.5)
            @test quantile(d,0.9) ≈ quantile(dref,0.9)

            @test cquantile(d,0.1) ≈ cquantile(dref,0.1)
            @test cquantile(d,0.5) ≈ cquantile(dref,0.5)
            @test cquantile(d,0.9) ≈ cquantile(dref,0.9)

            @test invlogcdf(d,log(0.2)) ≈ invlogcdf(dref,log(0.2))
            @test invlogcdf(d,log(0.5)) ≈ invlogcdf(dref,log(0.5))
            @test invlogcdf(d,log(0.8)) ≈ invlogcdf(dref,log(0.8))

            @test invlogccdf(d,log(0.2)) ≈ invlogccdf(dref,log(0.2))
            @test invlogccdf(d,log(0.5)) ≈ invlogccdf(dref,log(0.5))
            @test invlogccdf(d,log(0.8)) ≈ invlogccdf(dref,log(0.8))

            r = Array{float(eltype(d))}(undef, 100000)
            if ismissing(rng)
                rand!(d,r)
            else
                rand!(rng,d,r)
            end
            @test mean(r) ≈ mean(dref) atol=0.02
            @test std(r) ≈ std(dref) atol=0.01
            @test cf(d, -0.1) ≈ cf(dref,-0.1)

            if dref isa ContinuousDistribution
                @test gradlogpdf(d, 0.1) ≈ gradlogpdf(dref, 0.1)
            end
        end
        @testset "$k" for (k,dtest) in d_dict
            test_evaluation_and_sampling(rng, dtest, dref)
        end
    end
end

function test_location_scale_normal(
    rng::Union{AbstractRNG, Missing}, μ::Real, σ::Real, μD::Real, σD::Real,
)
    ρ = Normal(μD, σD)
    dref = Normal(μ + σ * μD, σ * σD)
    return test_location_scale(rng, μ, σ, ρ, dref)
end

function test_location_scale_discretenonparametric(
    rng::Union{AbstractRNG, Missing}, μ::Real, σ::Real, support, probs,
)
    ρ = DiscreteNonParametric(support, probs)
    dref = DiscreteNonParametric(μ .+ σ .* support, probs)
    return test_location_scale(rng, μ, σ, ρ, dref)
end

@testset "LocationScale" begin
    rng = MersenneTwister(123)

    for _rng in (missing, rng)
        test_location_scale_normal(_rng, 0.3, 0.2, 0.1, 0.2)
        test_location_scale_normal(_rng, -0.3, 0.1, -0.1, 0.3)
        test_location_scale_normal(_rng, 1.3, 0.4, -0.1, 0.5)
    end
    test_location_scale_normal(rng, ForwardDiff.Dual(0.3), 0.2, 0.1, 0.2)

    probs = Distributions.pnormalize!(rand(10))
    for _rng in (missing, rng)
        test_location_scale_discretenonparametric(_rng, 1//3, 1//2, 1:10, probs)
        test_location_scale_discretenonparametric(_rng, -1//4, 1//3, (-10):(-1), probs)
        test_location_scale_discretenonparametric(_rng, 6//5, 3//2, 15:24, probs)
    end
end


text:
# divide-partition-and-search (DPS) heuristic
function divide_partition_search(
    Ct::Matrix{Float64},
    Cd::Matrix{Float64};
    n_groups::Int = 1,
    local_search_methods::Vector{Function} = Function[two_point_move, one_point_move, two_opt_move],
    flying_range::Float64 = MAX_DRONE_RANGE,
    time_limit::Float64 = MAX_TIME_LIMIT
)

    time0 = time()

    n1, n2 = size(Ct)
    n_nodes = n1 - 1

    @assert size(Ct) == size(Cd)

    tsp_tour = find_tsp_tour(Ct[1:end-1, 1:end-1])
    push!(tsp_tour, n_nodes+1) # adding the final depot

    total_tspd_len = 0.0
    total_t_route = Int[]
    total_d_route = Int[]
    # @show tsp_tour

    n = n_nodes + 1
    group_size = Int((floor(n / n_groups)))

    remaining_time_limit = time_limit - (time() - time0)
    time_limit_each_group = remaining_time_limit / n_groups

    for i in 1:n_groups
        start_idx = 1 + group_size * (i-1)
        end_idx = min(1 + group_size * i, n)

        group_idx = start_idx:end_idx
        group_nodes = tsp_tour[group_idx]

        Ct_ = Ct[group_nodes, group_nodes]
        Cd_ = Cd[group_nodes, group_nodes]
        init_tour = collect(1:length(group_nodes))

        tspd_len, t_route_idx, d_route_idx = tsp_ep_all(Ct_, Cd_, init_tour; local_search_methods=local_search_methods, flying_range=flying_range, time_limit=time_limit_each_group)

        total_tspd_len += tspd_len
        append!(total_t_route, group_nodes[t_route_idx])
        append!(total_d_route, group_nodes[d_route_idx])
    end

    unique!(total_t_route)
    unique!(total_d_route)

    obj_val = objective_value(total_t_route, total_d_route, Ct, Cd)
    @assert isapprox(obj_val, total_tspd_len)

    # @show total_tspd_len

    return total_tspd_len, total_t_route, total_d_route

end


text:
function objectivefunc(A1, A0, Y,Agrid1,EV1,intfV)

#-------------------------------------------------------------------------------%
# This function returns the following quantity:
# u(c) +  b V( A1)
# where c is calculated from today's assets and tomorrow's assets
#-------------------------------------------------------------------------------%

# ------------------------------------------------------------------------
# Declare global we need this file have access to
#global beta r interpMethod            % structural model parameters
#global Agrid1 EV1                     % tomorrow's asset grid and tomorrow's expected value function

# ------------------------------------------------------------------------
#Get tomorrow's consumption (cons), the value of left over assets (VA1) and
#total value (u(c) + b * VA1
cons = A0  + Y - (A1)/(1+r);

# In order to obatin the future value function, you need to interpolate on the grid
# of assets. We don't need that for Income as it was "integrated out" at the end
# of the solution loop

    #intfV  = InterpIrregular(Agrid1,EV1, BCnil, InterpLinear); #Already done outside this function
    VA1=interp1extra(Agrid1,EV1,A1,intfV);

value = utility(cons) + beta * VA1;

# ------------------------------------------------------------------------
#The optimisation routine that we will use searches for the minimum of the
#function. We want the maximum. So we multiply out function here by -1 so
#that the optimiser will fill the minimum of the negative of our function,
#i.e. the maximum of our functino

return value = - value;

end

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
using DifferentialEquations
using Plots

# Parameters
const ib = 2e-5; # g⋅m⁻²⋅day⁻¹ immigration rate of bream
const ip = 2e-5; # g⋅m⁻²⋅day⁻¹ immigration rate of pike
const r = 7.5e-3; # day⁻¹ maximum growth rate of bream
const H₁ = 0.5; # half saturation constant
const H₂ = 10; # % half saturation constant
const H₃ = 20; # g⋅m⁻² half saturation constant
const H₄ = 15; # g⋅m⁻² half saturation constant
const cb = 7.5e-5; # m⁻²⋅g⁻¹⋅day⁻¹ intraspecific competition constant for bream
const cp = 2.75e-4; # m⁻²⋅g⁻¹⋅day⁻¹ intraspecific competition constant for pike
const prmax = 5e-2; # day⁻¹ maximum predation rate of pike
const ce = 0.1; # pike food conversion efficiency to growth
const mp = 2.25e-3; # day⁻¹ mortality rate of pike
const K = 100; # % maximum vetetation coverage

function scheffer!(du,u,p,t)
    B, P = u # g⋅m⁻² bream/pike density
    ib, r, nutr, H₁, H₂, H₃, H₄, cb, prmax, ip, ce, mp, cp, K = p

    V = K*(H₃^2/(H₃^2 + B^2)); # % of lake covered vetetation
    FR = B^2/(B^2 + H₄^2); # fuctional response of pike

    du[1] = dB = ib + r*(nutr/(nutr + H₁))*B - cb*B^2 - prmax*FR*P
    du[2] = dP = ip + ce*prmax*FR*P*(V/(V + H₂)) - mp*P - cp*P^2
end

function run_ss_solver(;
             nutr = 0.5, #nutrient level
             u0 = [40.0, 2.6], #Initial populations
             p = [ib, r, nutr, H₁, H₂, H₃, H₄, cb, prmax, ip, ce, mp, cp, K])
    prob = SteadyStateProblem(scheffer!,u0,p);
    solve(prob)
end

equilibria = Set();
for i in range(0.0, stop=100.0, length=10)
    for j in range(0.0, stop=10.0, length=10)
        for k in [0.5, 1.5, 2.5];
            eq = run_ss_solver(u0=[i,j], nutr=k);
            push!(equilibria, (eq[1], eq[2]));
        end
    end
end

function dB(B, P)
    FR = B^2/(B^2 + H₄^2); # fuctional response of pike

    ib + r*(nutr/(nutr + H₁))*B - cb*B^2 - prmax*FR*P
end

function dP(B, P)
    V = K*(H₃^2/(H₃^2 + B^2)); # % of lake covered vetetation
    FR = B^2/(B^2 + H₄^2); # fuctional response of pike

    ip + ce*prmax*FR*P*(V/(V + H₂)) - mp*P - cp*P^2
end

bream = range(0.0, stop=100.0, length=100)
pike = range(0.0, stop=10.0, length=100)

nutr = 2.5;
contour(bream, pike, dB, levels=[0], color=:blue, cbar=false, legend=true, label="bream'=0, nutr=2.5");
nutr = 1.5;
contour!(bream, pike, dB, levels=[0], color=:red, cbar=false, label="bream'=0, nutr=1.5");
nutr = 0.5;
contour!(bream, pike, dB, levels=[0], color=:green, cbar=false, label="bream'=0, nutr=0.5");
contour!(bream, pike, dP, levels=[0], color=:black, cbar=false, label="pike'=0");

eq = collect(equilibria);
scatter!(first.(eq), last.(eq), label="equilibria");

gui()


text:
using RFCVGT
using Base.Test

# write your own tests here
@test 1 == 2

#using DataFrames,MLBase,DecisionTree

## 判別分類
### 学習用データと検証用データに分ける
using RDatasets
iris = dataset("datasets", "iris");
using SamplingDataFrame
dataindex=TestDataSetsNumbering(iris,PERCENTAGE=0.85)
#### 学習用
training=iris[dataindex["TrainingNumbering"] ,:];
#### 検証用
test=iris[dataindex["TestNumbering"] ,:];

### 目的変数の列番号を指定
yind=[5];
### 説明変数の列番号を指定
xind=convert(Array,1:4);

### パッケージがデータを扱えるように変換
Ytraining=convert(Array , training[:,yind] );
Xtraining=convert(Array , training[:,xind] );

### パラメータサーチ
r=searchparm(Ytraining,Xtraining,nfeature=[1,2,4],ntree=[250,500])

## 回帰
### テスト用データ
x1=[1,11,11,3,1,2,21,1,11,10];
x2=[29,31,56,71,31,54,47,40,66,68];
x3=[15,8,8,17,22,18,4,23,9,8];
x4=[52,47,20,6,44,22,26,34,12,12];
y=[74.3,87.6,104.3,102.7,72.5,93.1,115.9,83.8,113.3,109.4];

ds=DataFrame(y=y,x1=x1,x2=x2,x3=x3,x4=x4)

### 検証用データ
newx1=[7,7,11];
newx2=[26,52,55];
newx3=[6,6,9];
newx4=[60,33,22];
newy=[78.5,95.9,109.2];
newX = DataFrame(x1=newx1,x2=newx2,x3=newx3,x4=newx4)

### ランダムフォレスト用に変換
Ytraining=convert(Array,ds[:,:y])
Xtraining=convert(Array,ds[:,[:x1,:x2,:x3,:x4]])
newXrf=convert(Array,newX)

### パラメータサーチ
r=searchparm(Ytraining,Xtraining)

max_stars_repo_name: wakakusa/RFCVGT

max_issues_repo_name: wakakusa/RFCVGT

max_forks_repo_name: wakakusa/RFCVGT


text:
module FFNetwork

    include("ActFunction.jl")
    include("TestDatum.jl")
    include("Network.jl")
    include("Cost.jl")
    include("Tools.jl")
    include("BackPropagation.jl")
    include("Trainer.jl")
    include("MNIST.jl")

    export Network;
    export TestDatum;
    export input;
    export output;
    export inlength;
    export outlength;
    export feedforward!;
    export train!;

end # module


text:
for id_typ in (:BoardID, :ColumnID, :CardID, :LabelID, :CommentID)
    @eval $(
        quote
            struct $id_typ
                id::String
            end
        end
    )
end
## Board API
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "GET", "/boards", ["fields", "archived", "page", "per_page", "sort"]))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "POST", "/boards", []))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "GET", "/boards/{board_id}", ["fields"]))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "POST", "/boards/{board_id}", []))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "DELETE", "/boards/{board_id}", []))
## Column API
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "POST", "/boards/{board_id}/columns", []))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "POST", "/boards/{board_id}/columns/batch", []))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "POST", "/boards/{board_id}/columns/{column_id}", []))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "DELETE", "/boards/{board_id}/columns/{column_id}", []))
## Card API
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "GET", "/boards/{board_id}/cards", ["fields", "archived", "page", "per_page", "sort"]))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "POST", "/boards/{board_id}/cards", []))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "POST", "/boards/{board_id}/cards/batch", []))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "GET", "/boards/{board_id}/cards/{card_id}", ["fields"]))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "POST", "/boards/{board_id}/cards/{card_id}", []))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "DELETE", "/boards/{board_id}/cards/{card_id}", []))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "GET", "/boards/{board_id}/columns/{column_id}/cards", ["fields", "archived", "page", "per_page", "sort"]))
## Lables API
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "POST", "/boards/{board_id}/labels", []))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "POST", "/boards/{board_id}/labels/{label_id}", []))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "DELETE", "/boards/{board_id}/labels/{label_id}", []))
## Attachments API
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "GET", "/boards/{board_id}/cards/{card_id}/attachments", ["fields", "page", "per_page", "sort"]))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "POST", "/boards/{board_id}/cards/{card_id}/attachments", []))
## Comments API
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "GET", "/boards/{board_id}/cards/{card_id}/comments", ["fields", "page", "per_page", "sort"]))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "POST", "/boards/{board_id}/cards/{card_id}/comments", []))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "POST", "/boards/{board_id}/cards/{card_id}/comments/batch", []))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "POST", "/boards/{board_id}/cards/{card_id}/comments/{comment_id}", []))
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "DELETE", "/boards/{board_id}/cards/{card_id}/comments/{comment_id}", []))
## User API
@eval $(declare_api("https://gloapi.gitkraken.com/v1/glo", "GET", "/user", ["fields"]))


text:
struct MockIntegrator
    du::Vector
    u::Vector
    p::Any
end

function MockIntegrator(inputs)
    return MockIntegrator(
        zeros(PSID.get_variable_count(inputs)),
        zeros(PSID.get_variable_count(inputs)),
        inputs,
    )
end


text:
"""
  log_meixner_ortho_poly(a, b, x)

Compute the logarithm of the Meixner polynomial
``M_{b}(a; -(a + b - 1), exp(x))`` where `a` and `b` are positive integers, with
`a >= b`.

Define
``u_{n} = 1 + (1 + a (1 - exp(x)) - n) / ((c + 1 - n) exp(x))``
and
``v_{n} = (n - 1) / ((c + 1 - n) exp(x))``

If ``y_{n}`` represents the Meixner polynomial of degree `n`, then
``y_{n} = u_{n} y_{n - 1} + v_{n} y_{n - 2}``
where ``y_{0} = 1`` and ``y_{1} = 1 + a * (1 - exp(x)) / (c * exp(x))``
"""
function log_meixner_ortho_poly(
  a::F,
  b::F,
  x::F
)::F where {
  F <: AbstractFloat
}
  lom = log1mexp(-x)
  lem = logexpm1(x)

  # define H_{b} = y_{b} / y_{b - 1}, i. e. y_{b} = H_{b} * y_{b - 1}
  # y_{0} = 1, y_{1} = 1 + a * (1 - x) / ((a + b - 1) * x)
  tmp = log1mexp(log(a / (a + b - 1)) + lom)

  # v = [log(H), log(y), compensation_error]
  v = [tmp, tmp, F(0)]

  n = F(1)
  while n < b
    n += 1
    kbn_meixner_ortho_poly_iteration!(v, a, n, a + b - 1, lom, lem)
  end

  # finally add the compensation error
  v[2] + v[3]
end

"""
  log_meixner_ortho_poly_joint(a, b, x)

Compute the logarithm of the Meixner polynomials
``M_{b}(a; -(a + b - 1), exp(x))``,
``M_{b - 1}(a - 1; -(a + b - 2), exp(x))``, and
``M_{b - 2}(a - 2; -(a + b - 3), exp(x))`` at the same time. `a` and `b` are
positive integers, `a >= b`. We assume that `a` and `b` are both greater than 1.
"""
function log_meixner_ortho_poly_joint(
  a::F,
  b::F,
  x::F
)::Tuple{F, F, F} where {
  F <: AbstractFloat
}
  lom = log1mexp(-x)
  lem = logexpm1(x)

  # define H_{b} = y_{b} / y_{b - 1}, i. e. y_{b} = H_{b} * y_{b - 1}
  # y_{0} = 1, y_{1} = 1 + a * (1 - x) / (c * x)

  # M_{b}(a; -(a + b - 1), x)
  c1 = a + b - 1
  tmp = log1mexp(log(a / c1) + lom)
  v1 = [tmp, tmp, F(0)]

  # M_{b - 1}(a - 1; -(a + b - 2), x)
  c2 = a + b - 2
  tmp = log1mexp(log((a - 1) / c2) + lom)
  v2 = [tmp, tmp, F(0)]

  if b > 2
    # M_{b - 2}(a - 2; -(a + b - 3), x)
    c3 = a + b - 3
    tmp = log1mexp(log((a - 2) / c3) + lom)
    v3 = [tmp, tmp, F(0)]

    n = F(1)
    while n < b - 2
      n += 1
      kbn_meixner_ortho_poly_iteration!(v1, a, n, c1, lom, lem)
      kbn_meixner_ortho_poly_iteration!(v2, a - 1, n, c2, lom, lem)
      kbn_meixner_ortho_poly_iteration!(v3, a - 2, n, c3, lom, lem)
    end

    # We still need few iterations
    kbn_meixner_ortho_poly_iteration!(v1, a, b - 1, c1, lom, lem)
    kbn_meixner_ortho_poly_iteration!(v2, a - 1, b - 1, c2, lom, lem)

    kbn_meixner_ortho_poly_iteration!(v1, a, b, c1, lom, lem)

    # add compensation errors and return the final values
    v1[2] + v1[3], v2[2] + v2[3], v3[2] + v3[3]
  else
    # M_{2}(a; -(a + 1), x)
    # M_{1}(a - 1; -(a + 2), x)
    # M_{0}(a - 2; -(a + 3), x)
    kbn_meixner_ortho_poly_iteration!(v1, a, F(2), c1, lom, lem)
    v1[2], v2[2], 0
  end
end

"""
    kbn_meixner_ortho_poly_iteration!(v, a, n, c, x)

Basic Kahan iteration operations for computing Meixner orthogonal polynomial.
"""
@inline function kbn_meixner_ortho_poly_iteration!(
  v::Vector{F},
  a::F,
  n::F,
  c::F,
  # lom = log(1 - exp(-x))
  lom::F,
  # lem = log(exp(x) - 1)
  lem::F
) where {
  F <: AbstractFloat
}
  # new value for log(y_{n} / y_{n - 1})
  v[1] = log1mexp(log(a / (c + 1 - n)) + lom +
                  log1mexp(log((n - 1) / a) + logexpm1(-v[1]) - lem))

  # running sum including error
  s = v[2] + v[1]

  # Neumaier modification to Kahan summation
  v[3] += if abs(v[2]) >= abs(v[1])
    (v[1] - s) + v[2]
  else
    (v[2] - s) + v[1]
  end

  v[2] = s

  nothing
end


text:
using Pkg
Pkg.activate("RunDSO")

println("Loading required packages: ProgressMeter, CSV, DataFrames, Dates, JuMP, Gurobi.")
using ProgressMeter
using CSV
using DataFrames
using Dates
using JuMP
using Gurobi
println("Loading of packages complete.")

###################
# Terminal ########
###################

    usedayperiod = false
    startday = 1
    lastday = 2
    starthour = 1
    lasthour = 8760
	x_ag = 0.3
	x_re = 0.3
	x_ind = 0.4
	s_ag = 0.05
	s_re = 0.05
	s_ind = 0.1
    workhours = 9:16

    ### Create Dict for average RD Price
    datapath = "data/"
    inputpath = "Processed_output/BAU/"
    output = "output/"
    outputpath = joinpath(output, "output_LoadShift/")

    isworkhour = Dict()
    for i in 0:24
        if i in workhours isworkhour[i] = 1
        else isworkhour[i] = 0
        end
    end
    if usedayperiod == true
        timeperiod = string((startday-1)*24+1)*":"*string((lastday*24))
    else
        timeperiod = string(starthour)*":"*string(lasthour)
    end

    timeperiod_str = replace(string(timeperiod), ":" => "-")

    println()
    # println("Case >" * case * "< selected.")
    println("Calculating for new load profiles.")
    println()

###################
# Run CM ##########
###################

    include("model/src/Sesam_ptdf_LoadShift.jl")

    using .Sesam

###################
# Data ############
###################
# Read nodes and load from .csv
    println("Reading data from csv...")

	nodes_df = CSV.read(joinpath(datapath, "nodes.csv"), DataFrame)

	load_df = CSV.read(joinpath(datapath, "load.csv"), DataFrame)

	exchange_df = CSV.read(joinpath(datapath, "exchange.csv"), DataFrame)

    nodal_hourlycost_df = CSV.read(joinpath(inputpath, "hourly_nodal_RD_cost.csv"), DataFrame)

    nodal_hourlyvolume_df = CSV.read(joinpath(inputpath, "hourly_nodal_RD_volume.csv"), DataFrame)

    marketprice_df = CSV.read(joinpath(inputpath, "ED_price.csv"), DataFrame)

    println("Reading of data complete.")
    println()
    println("Creating dictionaries for costs, volumes and marketprice...")

    nod_vol_dict = Dict((nodal_hourlyvolume_df[i, 1] => Dict((j => nodal_hourlyvolume_df[i, j+1]) for j in 1:8760)) for i in 1:nrow(nodal_hourlyvolume_df))
    nod_cost_dict = Dict((nodal_hourlycost_df[i, 1] => Dict((j => nodal_hourlycost_df[i, j+1]) for j in 1:8760)) for i in 1:nrow(nodal_hourlycost_df))
    av_el_price_dict = Dict((i => (sum(marketprice_df[(i-1)*24+1:i*24, 2]))/24) for i in 1:365)
###################
# Model setup #####
###################

    println("Preparing nodes...")

    # Create sets
        # Create set of nodes
        nodes = Nodes(nodes_df, load_df, exchange_df, 1000, 1)

        N = nodes.id

    println("Preparation complete.")
    println()

    session_time = Dates.now()
    session = Dates.format(session_time, "yyyy-mm-dd__HH-MM-SS")*"_newload"

Newload_df = DataFrame(node = N)
LoadShift_df = DataFrame(hour = Int64[], node=Int64[], OriginalLoad=Float64[], NewLoad=Float64[], ShiftUp=Float64[], ShiftDn=Float64[])
CSV.write(joinpath(outputpath, "Loadshift.csv"), LoadShift_df)

prog = Progress(length(dayslicer(timeperiod)))

for i in 1:365
    hours = dayslicer(timeperiod)[i]
    println("Shifting load for hours " * string(hours) * "...")

    av_el_price = av_el_price_dict[i]
    av_rd_price_dict = Dict((n => Dict()) for n in N)
    LoadShift_df = DataFrame(hour = Int64[], node=Int64[], OriginalLoad=Float64[], NewLoad=Float64[], ShiftUp=Float64[], ShiftDn=Float64[])

    for n in N
        for t in hours
            if nod_vol_dict[n][t] == 0
                av_rd_price_dict[n][t] = 0
            else
                av_rd_price_dict[n][t] = nod_cost_dict[n][t]/nod_vol_dict[n][t]
            end
        end
    end

        global    Load_shifted_opt, ShiftUp_opt, ShiftDn_opt, ShiftAgUp_opt, ShiftAgDn_opt, ShiftReUp_opt, ShiftReDn_opt, ShiftIndUp_opt, ShiftIndDn_opt =
            DSO_min_LoadShift(hours, nodes, av_el_price, av_rd_price_dict, x_ag, x_re, x_ind, s_ag, s_re, s_ind, isworkhour)

    for t in hours
        Newload_df[!, string(t)] = ((Load_shifted_opt[t, n] for n in N) |> collect)
        for n in N
            newload = Load_shifted_opt[t, n]
            origload = nodes.load[n][t]
            if (newload - origload > 0)
                shiftup = newload - origload
                shiftdn = 0
            else
                shiftup = 0
                shiftdn = origload - newload
            end
            newload = Load_shifted_opt[t, n]
            origload = nodes.load[n][t]
            push!(LoadShift_df, [t n origload newload shiftup shiftdn])
        end
    end
    CSV.write(joinpath(outputpath, "Loadshift.csv"), LoadShift_df, append = true)
    println()
    next!(prog)
    println()
end

CSV.write(joinpath(outputpath, "Newload.csv"), Newload_df)

Newload_df
#=
sum(nodes.load[6][1:24])
sum((Load_shifted_opt[i, 6]) for i in 1:24)
nodes.load[6][1]
Load_shifted_opt[1, 6]
=#


text:
using FastClosures
using SpecialFunctions
using StaticArrays

#Definite integral of the envelope function f (devided by 2).
function cosine_pulse_integral(t, tf, θ)
	x = t/tf
	π2 = convert(eltype(t), pi)*2
	ct = (θ/2)*(x .- sin.(π2*x)/π2)
	return ct
end
#Envelope function.
function cosine_pulse(t, tf, θ; derivative::Bool = false)
	x = t/tf
	π2 = convert(eltype(t), pi)*2
	if !derivative
		return (θ/tf)*(1 .- cos.(π2*x))
	else
		return (θ*π2/(tf*tf))*sin.(π2*x)
	end
end

const normalized_integral = erf(3/sqrt(2))
function gaussian_pulse(t, tf, θ; derivative::Bool = false)
	σ = tf/6
	μ = tf/2
	n = normalized_integral*σ*sqrt(2*pi)
	x = (t .- μ)/σ
	ft = exp.(-x.*x/2)
	if !derivative
		ft = (θ/n)*ft
	else
		ft = -(θ/(n*σ))*x.*ft
	end
	return ft
end

function gaussian_pulse_integral(t, tf, θ)
	σ = tf/6
	μ = tf/2
	if t < μ
		x = (μ - t)/(sqrt(2)*σ)
		s = (normalized_integral - erf(x))/2
	else
		x = (t - μ)/(sqrt(2)*σ)
		s = (erf(x) + normalized_integral)/2
	end
	return s*θ/(2*normalized_integral)
end

function pre_pulse_basis!(pulse_vector, t::Real, tf, npulse; derivative::Bool=false)
	k0 = 2*pi/tf
	if derivative == false
		for n = 1: npulse
			φ = k0*n*t
			pulse_vector[2*n - 1] = 1 - cos(φ)
			pulse_vector[2*n] = sin(φ)
		end
	else
		for n = 1: npulse
			kn = k0*n
			pulse_vector[2*n - 1] = kn*sin(kn*t)
			pulse_vector[2*n] = kn*cos(kn*t)
		end
	end
	return nothing
end

function pre_Δbasis!(Δbasis, t, tf, nΔ; derivative::Bool=false)
	Δbasis[1] = 1.0
	k0 = 2*pi/tf
	for n = 1: nΔ
		kn = k0*n
		Δbasis[2*n] = cos(kn*t)
		Δbasis[2*n + 1] = sin(kn*t)
	end
	return nothing
end

function fourier_basis(tf, nΔ, npulse; basis_derivative::Bool=false)
	Δvector = zeros(1 + 2*nΔ)
	Δbasis! = @closure (v, t)-> pre_Δbasis!(v, t, tf, nΔ; derivative=false)
	Δbasis = (function! =  Δbasis!, vector = Δvector)

	pulse_vector = zeros(2*npulse)
	pulse_basis! = @closure (v, t)-> pre_pulse_basis!(v, t, tf, npulse, derivative=false)
	pulse_basis = (function! =  pulse_basis!, vector = pulse_vector)
	basis = (pulse = pulse_basis, detuning = Δbasis)

	if basis_derivative == true
		dΔvector = zeros(1 + 2*nΔ)
		dΔbasis! = @closure (v, t)-> pre_Δbasis!(v, t, tf, nΔ; derivative=true)
		dΔbasis = (function! =  dΔbasis!, vector = dΔvector)

		dpulse_vector = zeros(2*npulse)
		dpulse_basis! = @closure (v, t)-> pre_pulse_basis!(v, t, tf, npulse; derivative=true)
		dpulse_basis = (function! =  dpulse_basis!, vector = dpulse_vector)
		basis_derivative = (pulse= dpulse_basis, detuning= dΔbasis)
		return basis, basis_derivative
	else
		return basis
	end
end

max_stars_repo_name: thalesfr/MagnusCorrection

max_issues_repo_name: thalesfr/MagnusCorrection

max_forks_repo_name: thalesfr/MagnusCorrection


text:
# Support for the DVS128 camera.

const DVS128_SIZE = (128,128)

""" Event class for the DVS128 """
struct DVS128
    pol::Bool
    x::Int8
    y::Int8
    ext::Bool
end

function Base.convert(::Type{DVS128}, a::Union{UInt16,UInt32})
    # Unpack the DVS128 fields from a 16-bit or 32-bit address:
    pol = a & 0x1
    x = (a >> 1) & 0x7F
    y = (a >> 8) & 0x7F
    ext = (a >> 15) & 0x1
    DVS128(pol, x, y, ext)
end

event_coord(e::DVS128) = (e.x, e.y)
event_pol(e::DVS128) = e.pol
event_location(e::DVS128, imsize) = (DVS128_SIZE[1]-e.y, e.x+1)

image_size(e::Type{DVS128}) = DVS128_SIZE

# There's no positive way to identify a DVS128 event.
# These are normally not mixed with other event types
# anyways. They can be found in 1.0 and 2.0 files.
isevent(::Type{DVS128}, a::Event{UInt16}) = true
isevent(::Type{DVS128}, a::Event{UInt32}) = true


text:
function xticks!(simple_plot::SimplePlot, cur_ticks)
  if isa(cur_ticks, Tuple)
    cur_ticks = collect.(cur_ticks)
  else
    cur_ticks = tuple(collect(cur_ticks))
  end

  simple_plot.xticks = cur_ticks
  validate!(simple_plot)
end

function yticks!(simple_plot::SimplePlot, cur_ticks)
  if isa(cur_ticks, Tuple)
    cur_ticks = collect.(cur_ticks)
  else
    cur_ticks = tuple(collect(cur_ticks))
  end

  simple_plot.yticks = cur_ticks
  validate!(simple_plot)
end

xticks!(cur_ticks) = xticks!(_plot, cur_ticks)
yticks!(cur_ticks) = yticks!(_plot, cur_ticks)

export xticks!
export yticks!


text:
########################################################################
# Global simulation of ACis curves from Medlyn data
# Simon Treillou, 2022
########################################################################
using DataFrames
using Plots
using RCall
using CSV
using PlantBiophysics
using Statistics
using Revise
constants = Constants()
using CairoMakie
using Dates
using MonteCarloMeasurements
unsafe_comparisons(true)
using LeafGasExchange
using Cropbox
using Makie

# Add absolute path for PlantBiophysics-ComplCode
path = "/Users/simon/Code/PlantBiophysics-paper/data/"

# Did you already have fitted the parameters ?
fitted = false

# Do you want to save the simulations ?
saving_simulations = false

# Do you want to save the figure ?
saving_figure = false

function RMSE(obs,sim)
    return sqrt(sum((obs .- sim).^2)/length(obs))
end

function EF(obs,sim)
    SSres = sum((obs - sim).^2)
    SStot = sum((obs .- mean(obs)).^2)
    return 1-SSres/SStot
end

# Loading R packages
R"""
library(plantecophys)
library(dplyr)
library(readr)
library(microbenchmark)
"""

########################################################################
# Parameters : fit them or download already fitted parameters
########################################################################
if !fitted
    df = read_licor6400(path*"TumbarumbaGasex_ACis_Medlyn.csv")
    df.Date = Date.(df.Date, Dates.DateFormat("yyyy/mm/dd"))

    df.VcMaxRef .= df.JMaxRef .= df.RdRef .= df.TPURef .= df.g0 .= df.g1 .= 0.
    df.VcMaxRefPE .= df.JMaxRefPE .= df.RdRefPE .= df.TPURefPE .= df.g0PE .= df.g1PE .= 0.
    for i in unique(df.Curve)
        dfi = filter(x -> x.Curve == i, df)
        sort!(dfi,:Cᵢ)

        g0,g1 = PlantBiophysics.fit(Medlyn,dfi)
        df.g0[df.Curve.==i] .= g0
        df.g1[df.Curve.==i] .= g1

        filter!(x->x.PPFD > 1400., dfi)

        VcMaxRef,JMaxRef,RdRef,TPURef = PlantBiophysics.fit(Fvcb,dfi)
        df.VcMaxRef[df.Curve.==i] .= VcMaxRef
        df.JMaxRef[df.Curve.==i] .= JMaxRef
        df.RdRef[df.Curve.==i] .= RdRef
        df.TPURef[df.Curve.==i] .= TPURef

        dfiPE = dfi[:,3:end]
        P = mean(dfi.P)
        rename!(dfiPE, :Tₗ=>:Tleaf, :Cᵢ=>:Ci, :Cₐ=>:Ca, :Gₛ=>:gs)
        @rput P
        @rput dfiPE
        R"""
        fit = fitaci(
            dfiPE,
            varnames = list(ALEAF = "A", Tleaf = "Tleaf", Ci = "Ci", PPFD = "PPFD"),
            Tcorrect = TRUE,
            Patm = P,
            citransition = NULL,
            quiet = TRUE,
            startValgrid = TRUE,
            fitmethod = "bilinear",
            algorithm = "nls",
            fitTPU = TRUE,
            alphag = 0,
            useRd = FALSE,
            PPFD = dfiPE$PPFD,
            Tleaf = dfiPE$Tleaf,
            alpha = 0.24,
            theta = 0.7,
            gmeso = NULL,
            EaV = 58550.0,
            EdVC = 2e+05,
            delsC = 629.26,
            EaJ = 29680.0,
            EdVJ = 2e+05,
            delsJ = 631.88,
            #GammaStar = NULL,
            #Km = NULL,
            id = NULL)
        coefs = coef(fit)
        """
        @rget coefs
        dfi = filter(x -> x.Curve == i, df)
        dfiPE = dfi[:,3:end]
        rename!(dfiPE, :Tₗ=>:Tleaf, :Cᵢ=>:Ci, :Cₐ=>:Ca, :Gₛ=>:gs, :Dₗ=>:VpdL)
        filter!(x->x.A ./ (x.Ca)>0.,dfiPE)
        @rput dfiPE
        R"""
        dfiPE$Rh = dfiPE$Rh*100
        myfit = fitBB(
            dfiPE,
            varnames = list(ALEAF = "A", GS = "gs", VPD = "VpdL", Ca = "Ca", RH = "Rh"),
            gsmodel = c("BBOpti"),
            fitg0 = TRUE,
            D0 = NULL
            )
        gP = coef(myfit)
        """
        @rget gP
        df.VcMaxRefPE[df.Curve.==i] .= coefs[1]
        df.JMaxRefPE[df.Curve.==i] .= coefs[2]
        df.RdRefPE[df.Curve.==i] .= coefs[3]
        if !ismissing(coefs[4])
            df.TPURefPE[df.Curve.==i] .= coefs[4]
        else
            df.TPURefPE[df.Curve.==i] .= 1000000.
        end
        df.g0PE[df.Curve.==i] .= gP[1]
        df.g1PE[df.Curve.==i] .= gP[2]
    end
    CSV.write(path*"../tutorials/Medlyn_ACis_param.csv",df)

else
    df  = CSV.File(path*"../tutorials/Medlyn_ACis_param.csv", dateformat="yyyy/mm/dd") |> DataFrame
end
########################################################################
# Simulation
########################################################################
# Define characteristic length (here we have the Licor6400 area, we deduce a
# characteristic length in m
d = sqrt(df.Area[1])/100#sqrt(df.Area[1]/pi) / 100
Wind = 20.
@rput Wind
@rput d

########################################################################
# PlantBiophysics.jl

df.AsimPB .= df.EsimPB .= df.TlsimPB .= df.GssimPB .= 0.
for i in unique(df.Curve)
    dfi = filter(x->x.Curve == i,df)

    cols = fieldnames(Atmosphere)
    dfiMeteo = select(dfi, names(dfi, x -> Symbol(x) in cols))
    dfiMeteo.Wind .= Wind
    # Note that as we only use A-Ci curves, there is no NIR in the Licor6400
    dfiMeteo.Ri_SW_f .= dfi.PPFD ./ (4.57)
    meteo = Weather(dfiMeteo)

    leaf = LeafModels(
            energy = Monteith(aₛₕ = 2, aₛᵥ = 1,ε = 0.955 ,maxiter=100),
            photosynthesis = Fvcb(VcMaxRef = dfi.VcMaxRef[1], JMaxRef = dfi.JMaxRef[1], RdRef = dfi.RdRef[1], TPURef = dfi.TPURef[1]),
            stomatal_conductance = Medlyn(dfi.g0[1], dfi.g1[1]),
            Rₛ = meteo[:Ri_SW_f],
            sky_fraction = 1.0,
            PPFD = dfi.PPFD,
            d = d
    )

    energy_balance!(leaf,meteo)
    df.AsimPB[df.Curve .== i, :] = DataFrame(leaf).A
    df.EsimPB[df.Curve .== i, :] = DataFrame(leaf).λE ./(meteo[:λ] * constants.Mₕ₂ₒ)*1000
    df.TlsimPB[df.Curve .== i, :] = DataFrame(leaf).Tₗ
    df.GssimPB[df.Curve .== i, :] = DataFrame(leaf).Gₛ
end

########################################################################
# LeafGasExchange.jl

df.AsimLG .= df.EsimLG .= df.TlsimLG .= df.GssimLG .= 0.
for i in unique(df.Curve)
    dfi = filter(x->x.Curve == i,df)
    #dfi = filter(x->x.PPFD > 1400.,dfi)
    #dfi = filter(x->x.Cₐ > 150.,dfi)

    configs = []
    for i in 1:length(dfi.T)
        config = :Weather => (
                PFD = dfi.PPFD[i],
                CO2 = dfi.Cₐ[i],
                RH = dfi.Rh[i] * 100,
                T_air = dfi.T[i],
                wind = Wind,
                P_air = dfi.P[i],
                g0 = dfi.g0PE[i],
                g1 = 1.57*dfi.g1PE[i],
                Vc25 = dfi.VcMaxRefPE[i],
                Jm25 = dfi.JMaxRefPE[i],
                Rd25 = dfi.RdRefPE[i],
                Tp25 = dfi.TPURefPE[i],
                Ear = 46.39,
                Haj = 29.68,
                w = d,
                d = d,
                EaVc = 58.55,
                ϵ = 0.955,
                Dh = 21.5
        )
        push!(configs,config)
    end
    res = simulate(ModelC3MD;configs=configs,nounit=true)
    df.AsimLG[df.Curve .== i, :] = res.A_net
    df.EsimLG[df.Curve .== i, :] = res.E
    df.TlsimLG[df.Curve .== i, :] = res.T
    df.GssimLG[df.Curve .== i, :] = res.gsc
end

########################################################################
# plantecophys

R"""
A_sim = c()
E_sim = c()
Tl_sim = c()
Gs_sim = c()
failed =c()
"""

df.AsimPE .= df.EsimPE .= df.TlsimPE .= df.GssimPE .= df.PEfailed .= 0.
for i in unique(df.Curve)
    dfi = filter(x->x.Curve == i,df)
    dfi = dfi[:,3:end]
    @rput dfi
    Ca = dfi.Cₐ
    @rput Ca
    R"""
    VcMaxRef = dfi$VcMaxRefPE[1]
    JMaxRef = dfi$JMaxRefPE[1]
    TPURef = dfi$TPURefPE[1]
    RdRef = dfi$RdRefPE[1]
    g0 = dfi$g0PE[1]
    g1 = dfi$g1PE[1]
    res = PhotosynEB(Tair = dfi$T,Wind = Wind,VPD=dfi$VPD,
             Wleaf = d,Ca = Ca,  StomatalRatio = 1,
             LeafAbs = 1.,gsmodel = "BBOpti",g0 = g0, g1 = g1,
             EaV = 58550.0,EdVC = 2e+05, delsC = 629.26,
             EaJ = 29680.0,EdVJ = 2e+05,delsJ = 631.88,
             alpha = 0.24,theta = 0.7, Jmax = JMaxRef,
             Vcmax = VcMaxRef, TPU = TPURef,Rd = RdRef,
             RH = dfi$Rh*100,PPFD=dfi$PPFD,
             Patm = dfi$P,gk=0.,
             Tcorrect = FALSE)
    A_sim = append(A_sim,res$ALEAF)
    Gs_sim = append(Gs_sim,res$GS)
    failed = append(failed,res$failed)
    Tl_sim = append(Tl_sim,res$Tleaf)
    E_sim = append(E_sim,res$ELEAF)
    """
end
@rget A_sim

@rget failed
@rget Gs_sim
@rget E_sim
@rget Tl_sim

df.AsimPE .= A_sim
df.EsimPE .= E_sim
df.TlsimPE .= Tl_sim
df.GssimPE .= Gs_sim
df.PEfailed .= failed

if saving_simulations
    CSV.write(path*"../tutorials/Medlyn_ACis_simulations.csv",df)
end

########################################################################
# Plotting with Makie.jl
########################################################################
CairoMakie.activate!(type = "svg")
noto_sans = assetpath("fonts", "NotoSans-Regular.ttf")
fig = Makie.Figure(backgroundcolor = RGBf(1, 1, 1),resolution = (1000, 1000),font=noto_sans,dpi=600,size=(1000,1000),xminorgridstyle=true)
sideinfo1 = Label(fig[1:2, 0], "Simulations", rotation = pi/2,labelsize=25)
sideinfo2 = Label(fig[3,2:3], "Observations")

########################################################################
# Assimilation
axa = Axis(fig[1,2],title="(a) Net CO₂ assimilation Aₙ")
index = 1:length(df.PPFD)
index = ifelse.(df.PPFD .> 1400, index, 9999.)
filter!(x->x!=9999.,index)
indexPE = 1:length(df.PPFD)
indexPE = ifelse.((df.PPFD .> 1400).&(df.PEfailed .==0.), indexPE, 9999.)
filter!(x->x!=9999.,indexPE)
RMSPB=round(RMSE(df.A[index],df.AsimPB[index]),digits=3)
RMSLG=round(RMSE(df.A[index],df.AsimLG[index]),digits=3)
RMSPE=round(RMSE(df.A[indexPE],df.AsimPE[indexPE]),digits=3)
EF(df.A[index],df.AsimPB[index])
EF(df.A[index],df.AsimLG[index])
EF(df.A[indexPE],df.AsimPE[indexPE])
Makie.lines!(axa,df.AsimPB[index],df.AsimPB[index],color=(:grey,0.4),linewidth=4)
LG = Makie.scatter!(axa,df.A[index],df.AsimLG[index],color=(:black,0.5),markersize=8,marker='□')
PE = Makie.scatter!(axa,df.A[indexPE],df.AsimPE[indexPE],color=(:blue,0.5),markersize=8,marker='∆')
PB = Makie.scatter!(axa,df.A[index],df.AsimPB[index],color=(:red,0.8),markersize=8,marker='o')
axislegend(axa, [PB,PE,LG], ["PlantBiophysics.jl (RMSE="*string(RMSPB)*")","plantecophys (RMSE="*string(RMSPE)*")","LeafGasExchange.jl (RMSE="*string(RMSLG)*")"], "", position = :rb,
    orientation = :vertical,labelsize=13)

########################################################################
# Transpiration
index = 1:length(df.PPFD)
index = ifelse.(df.Cₐ .> 140., index, 9999.)
filter!(x->x!=9999.,index)
indexPE = 1:length(df.PPFD)
indexPE = ifelse.((df.Cₐ .> 140.).&(df.PEfailed .==0.), indexPE, 9999.)
filter!(x->x!=9999.,indexPE)

axb = Axis(fig[1,3],title="(b) Transpiration rate E")

RMSPB=round(RMSE(df.Trmmol[index],df.EsimPB[index]),digits=3)
RMSLG=round(RMSE(df.Trmmol[index],df.EsimLG[index]),digits=3)
RMSPE=round(RMSE(df.Trmmol[indexPE],df.EsimPE[indexPE]),digits=3)
EF(df.Trmmol[index],df.EsimPB[index])
EF(df.Trmmol[index],df.EsimLG[index])
EF(df.Trmmol[indexPE],df.EsimPE[indexPE])
Makie.lines!(axb,df.EsimPB[index],df.EsimPB[index],linewidth=4,color=(:grey,0.4))
LG = Makie.scatter!(axb,df.Trmmol[index],df.EsimLG[index],color=(:black,0.5),markersize=8,marker='□')
PE = Makie.scatter!(axb,df.Trmmol[indexPE],df.EsimPE[indexPE],color=(:blue,0.5),markersize=8,marker='∆')
PB = Makie.scatter!(axb,df.Trmmol[index],df.EsimPB[index],color=(:red,0.8),markersize=8,marker='o')
axislegend(axb, [PB,PE,LG], ["PlantBiophysics.jl (RMSE="*string(RMSPB)*")","plantecophys (RMSE="*string(RMSPE)*")","LeafGasExchange.jl (RMSE="*string(RMSLG)*")"], "", position = :rb,
    orientation = :vertical,labelsize=13)

########################################################################
# Leaf temperature
indexPE = 1:length(df.PPFD)
indexPE = ifelse.((df.PEfailed .==0.), indexPE, 9999.)
filter!(x->x!=9999.,indexPE)

axd = Axis(fig[2,3],title="(d) Leaf temperature Tₗ")

RMSPB=round(RMSE(df.Tₗ,df.TlsimPB),digits=3)
RMSLG=round(RMSE(df.Tₗ,df.TlsimLG),digits=3)
RMSPE=round(RMSE(df.Tₗ[indexPE],df.TlsimPE[indexPE]),digits=3)
EF(df.Tₗ[index],df.TlsimPB[index])
EF(df.Tₗ[index],df.TlsimLG[index])
EF(df.Tₗ[indexPE],df.TlsimPE[indexPE])
Makie.lines!(axd,df.TlsimPB,df.TlsimPB,linewidth=4,color=(:grey,0.4))
LG = Makie.scatter!(axd,df.Tₗ,df.TlsimLG,color=(:black,0.5),markersize=8,marker='□')
PE = Makie.scatter!(axd,df.Tₗ[indexPE],df.TlsimPE[indexPE],color=(:blue,0.5),markersize=8,marker='∆')
PB = Makie.scatter!(axd,df.Tₗ,df.TlsimPB,color=(:red,0.8),markersize=8,marker='o')
axislegend(axd, [PB,PE,LG], ["PlantBiophysics.jl (RMSE="*string(RMSPB)*")","plantecophys (RMSE="*string(RMSPE)*")","LeafGasExchange.jl (RMSE="*string(RMSLG)*")"], "", position = :rb,
    orientation = :vertical,labelsize=13)

########################################################################
# Stomatal conductance
index = 1:length(df.PPFD)
index = ifelse.(df.Cₐ .> 140., index, 9999.)
filter!(x->x!=9999.,index)
indexPE = 1:length(df.PPFD)
indexPE = ifelse.((df.Cₐ .> 140.).&(df.PEfailed .==0.), indexPE, 9999.)
filter!(x->x!=9999.,indexPE)
axc = Axis(fig[2,2],title="(c) CO₂ stomatal conductance Gₛ")

RMSLG=round(RMSE(df.Gₛ[index],df.GssimLG[index]),digits=3)
RMSPB=round(RMSE(df.Gₛ[index],df.GssimPB[index]),digits=3)
RMSPE=round(RMSE(df.Gₛ[indexPE],df.GssimPE[indexPE]),digits=3)
EF(df.Gₛ[index],df.GssimPB[index])
EF(df.Gₛ[index],df.GssimLG[index])
EF(df.Gₛ[indexPE],df.GssimPE[indexPE])
Makie.lines!(axc,df.GssimPB,df.GssimPB,linewidth=4,color=(:grey,0.4))
LG = Makie.scatter!(axc,df.Gₛ[index],df.GssimLG[index],color=(:black,0.5),markersize=8,marker='□')
PE = Makie.scatter!(axc,df.Gₛ[indexPE],df.GssimPE[indexPE],color=(:blue,0.5),markersize=8,marker='∆')
PB = Makie.scatter!(axc,df.Gₛ[index],df.GssimPB[index],color=(:red,0.8),markersize=8,marker='o')
Makie.ylims!(axc,(0.,0.4))
Makie.xlims!(axc,(0.,0.4))
axislegend(axc, [PB,PE,LG], ["PlantBiophysics.jl (RMSE="*string(RMSPB)*")","plantecophys (RMSE="*string(RMSPE)*")","LeafGasExchange.jl (RMSE="*string(RMSLG)*")"], "", position = :rb,
    orientation = :vertical,labelsize=13)
fig

if saving_figure
    save(path*"../tutorials/figure_global_simulation.png", fig, px_per_unit = 2) # output size = 1600 x 1200 pixels
end

max_stars_repo_name: VEZY/PlantBiophysics-paper

max_issues_repo_name: VEZY/PlantBiophysics-paper

max_forks_repo_name: VEZY/PlantBiophysics-paper


text:
module Examples

using ..Pages

export examples

function examples(;port=8000,launch=false)
    Endpoint("/examples") do request::HTTP.Request
        read(joinpath(@__DIR__,"examples.html"),String)
    end

    Endpoint("/examples/requests") do request::HTTP.Request
        read(joinpath(@__DIR__,"requests","requests.html"),String)
    end

    Endpoint("/examples/requests/echo") do request::HTTP.Request
        params = HTTP.queryparams(HTTP.URI(request.target).query)
        println("Body: $(params)")
        response = JSON3.write(params)
    end

    Endpoint("/examples/requests/echo",POST) do request::HTTP.Request
        data = String(request.body)
        println("Parameters: $(data)")
        response = JSON3.write(Dict(:data => data))
    end

    Endpoint("/examples/blank") do request::HTTP.Request
        read(joinpath(@__DIR__,"blank","blank.html"),String)
    end

    include(joinpath(@__DIR__,"blank","blank.jl"))
    include(joinpath(@__DIR__,"randomping","randomping.jl"))

    Pages.start(port=port)
    launch && Pages.launch("http://localhost:$(port)/examples")
    return
end

end


text:
include("RangeInfiniteVector.jl")
include("DiscreteInfiniteVector.jl")


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# Run the function on the given items, using multi-threading.
# Kwarg `numThreads` can be set to reduce cpu usage, though this is only approximate.
# Requires environment variable JULIA_NUM_THREADS to be set before starting julia.
function runParallel!(f::Function, items...; numThreads::Int = maxNumThreads)
	@assert(numThreads >= 1)
	@threads for i = 1:numThreads
		for j = i:numThreads:length(items)
			f(items[j])
		end
	end
end


text:
struct WasmtimeInstance
    # wasmtime instances are "just" a 64 bits identifier associated to a store, as such
    # they don't have an associated destructor of type wasmtime_delete_instance()
    # see the definition of wasmtime_instance_t in LibWasmtime for more info.
    identifier::wasmtime_instance
    store::WasmtimeStore
end

Base.unsafe_convert(::Type{Ptr{wasmtime_instance_t}}, instance::WasmtimeInstance) =
    Base.convert(Ptr{wasmtime_instance}, Base.pointer_from_objref(instance.identifier))

function WasmtimeInstance(store::WasmtimeStore, mod::WasmtimeModule)
    module_imports = imports(mod)
    n_expected_imports = length(module_imports.wasm_imports)
    @assert n_expected_imports == 0 "No imports provided, expected $n_expected_imports"

    instance = Ref(wasmtime_instance_t(0, 0))
    wasm_trap_ptr = Ref(Ptr{wasm_trap_t}())
    @wt_check GC.@preserve instance wasmtime_instance_new(
        store,
        mod,
        C_NULL,
        0,
        instance,
        wasm_trap_ptr,
    )
    if wasm_trap_ptr[] != C_NULL
        trap_msg = WasmByteVec()
        wasm_trap_message(wasm_trap_ptr[], trap_msg)
        wasm_trap_delete(wasm_trap_ptr[])

        error_message = unsafe_string(trap_msg.data, trap_msg.size)
        error(error_message)
    end

    WasmtimeInstance(instance[], store)
end

function WasmtimeInstance(linker::WasmtimeLinker, store::WasmtimeStore, mod::WasmtimeModule)
    instance = Ref(wasmtime_instance_t(0, 0))
    wasm_trap_ptr = Ref(Ptr{wasm_trap_t}())
    @wt_check GC.@preserve instance wasmtime_linker_instantiate(
        linker,
        store,
        mod,
        instance,
        wasm_trap_ptr,
    )
    if wasm_trap_ptr[] != C_NULL
        trap_msg = WasmByteVec()
        wasm_trap_message(wasm_trap_ptr[], trap_msg)
        wasm_trap_delete(wasm_trap_ptr[])

        error_message = unsafe_string(trap_msg.data, trap_msg.size)
        error(error_message)
    end

    WasmtimeInstance(instance[], store)
end

mutable struct WasmtimeInstanceExport <: AbstractWasmExport
    wasm_export_ptr::Ptr{wasm_exporttype_t}
    extern::Ref{wasmtime_extern}

    wasmtime_instance::WasmtimeInstance
    name::String
end

function wasmtime_valkind_to_julia(valkind::wasmtime_valkind_t)::Type
    if valkind == WASMTIME_I32
        Int32
    elseif valkind == WASMTIME_I64
        Int64
    elseif valkind == WASMTIME_F32
        Float32
    elseif valkind == WASMTIME_F64
        Float64
    elseif valkind == WASMTIME_FUNCREF
        wasmtime_func_t
    elseif valkind == WASMTIME_EXTERNREF
        Ptr{wasmtime_externref_t}
    elseif valkind == WASMTIME_V128
        wasmtime_v128
    else
        error("Invalid value kind $type")
    end
end

function (wasmtime_export::WasmtimeInstanceExport)(args...)
    extern = wasmtime_export.extern[]
    @assert extern.kind == WASM_EXTERN_FUNC "Expected an exported function but got type $(extern.kind)"

    func = Ref(extern.of.func)
    functype = wasmtime_func_type(wasmtime_export.wasmtime_instance.store, func)

    wasm_params = Base.unsafe_convert(Ptr{WasmVec{wasm_valtype_vec_t,Ptr{wasm_valtype_t}}}, wasm_functype_params(functype)) |> Base.unsafe_load
    wasm_results = Base.unsafe_convert(Ptr{WasmVec{wasm_valtype_vec_t,Ptr{wasm_valtype_t}}}, wasm_functype_results(functype)) |> Base.unsafe_load

    @assert length(args) == length(wasm_params) "Expected $(wasm_params.size) arguments but got $(length(args))"
    params_kind = wasm_valtype_kind.(wasm_params)

    wasmtime_params = wasmtime_val[]
    for (i, (param, kind)) in enumerate(zip(args, params_kind))
        jtype = typeof(param)
        etype = wasmtime_valkind_to_julia(kind)
        @assert jtype == etype "Parameter #$i is of type $jtype, expected $etype"

        valunion = Ref{wasmtime_valunion}(wasmtime_valunion(Tuple(zero(UInt8) for _ in 1:16)))
        ptr = Base.unsafe_convert(Ptr{wasmtime_valunion}, valunion)
        GC.@preserve valunion ptr.i32 = etype(param)
        push!(wasmtime_params, wasmtime_val(WASM_I32, valunion[]))
    end

    wasmtime_results = Vector{wasmtime_val}(undef, length(wasm_results))

    trap = Ref(Ptr{wasm_trap_t}())
    @wt_check GC.@preserve wasmtime_params wasmtime_results wasmtime_func_call(
        wasmtime_export.wasmtime_instance.store,
        func,
        pointer(wasmtime_params),
        length(wasmtime_params),
        pointer(wasmtime_results),
        length(wasmtime_results),
        trap
    )

    results = map(wasmtime_results) do result
        if result.kind == WASMTIME_I32
            result.of.i32
        elseif result.kind == WASMTIME_I64
            result.of.i64
        elseif result.kind == WASMTIME_F32
            result.of.f32
        elseif result.kind == WASMTIME_F64
            result.of.f64
        elseif result.kind == WASMTIME_V128
            result.of.v128
        elseif result.kind == WASMTIME_FUNCREF
            result.of.funcref
        elseif result.kind == WASMTIME_EXTERNREF
            result.of.externref
        else
            error("Unknown value kind $(result.kind)")
        end
    end

    length(results) == 1 ?
        first(results) : results
end

function exports(instance::WasmtimeInstance)
    instance_type = wasmtime_instance_type(instance.store, instance)
    @assert instance_type != C_NULL "Failed to get module type from WasmtimeModule"

    wasm_exports = WasmPtrVec(wasm_exporttype_t)
    wasmtime_instancetype_exports(instance_type, wasm_exports)

    wasmtime_exports = map(enumerate(wasm_exports)) do (i, wasm_export_ptr)
        owned_wasm_export_ptr = wasm_exporttype_copy(wasm_export_ptr)
        @assert owned_wasm_export_ptr != C_NULL "Failed to copy WASM export"

        extern = wasmtime_extern(
            UInt8(0),
            wasmtime_extern_union(Tuple(UInt8(0xab) for _ = 1:16)),
        )
        name_vec = Ref(wasm_name(0, Ptr{wasm_byte_t}()))
        char_ptr =
            Base.unsafe_convert(
                Ptr{Nothing},
                Base.unsafe_convert(Ptr{wasm_byte_vec_t}, name_vec),
            ) + Base.sizeof(Csize_t)
        len_ptr = Base.unsafe_convert(
            Ptr{Csize_t},
            Base.unsafe_convert(Ptr{wasm_byte_vec_t}, name_vec),
        )
        exists = wasmtime_instance_export_nth(
            instance.store,
            instance,
            i - 1,
            char_ptr,
            len_ptr,
            Base.pointer_from_objref(extern),
        )
        exists || error("Export #$i does not exists")
        name = unsafe_string(name_vec[].data, name_vec[].size)

        WasmtimeInstanceExport(wasm_export_ptr, extern, instance, name)
    end

    WasmExports(instance, wasmtime_exports)
end


text:
@testset "CircularArrayBuffer" begin
    A = ones(2, 2)
    C = ones(Float32, 2, 2)
    @testset "1D Int" begin
        b = CircularArrayBuffer{Int}(3)

        @test eltype(b) == Int
        @test capacity(b) == 3
        @test isfull(b) == false
        @test isempty(b) == true
        @test length(b) == 0
        @test nframes(b) == 0
        @test size(b) == (0,)
        # element must has the exact same length with the element of buffer
        @test_throws Exception push!(b, [1, 2])

        for x in 1:3
            push!(b, x)
        end

        @test capacity(b) == 3
        @test isfull(b) == true
        @test length(b) == 3
        @test nframes(b) == 3
        @test size(b) == (3,)
        @test b[1] == 1
        @test b[end] == 3
        @test b[1:end] == [1, 2, 3]

        for x in 4:5
            push!(b, x)
        end

        @test capacity(b) == 3
        @test length(b) == 3
        @test nframes(b) == 3
        @test size(b) == (3,)
        @test b[1] == 3
        @test b[end] == 5
        @test b[1:end] == [3, 4, 5]

        empty!(b)
        @test isfull(b) == false
        @test isempty(b) == true
        @test length(b) == 0
        @test nframes(b) == 0
        @test size(b) == (0,)

        push!(b, 6)
        @test isfull(b) == false
        @test isempty(b) == false
        @test length(b) == 1
        @test nframes(b) == 1
        @test size(b) == (1,)
        @test b[1] == 6

        push!(b, 7)
        push!(b, 8)
        @test isfull(b) == true
        @test isempty(b) == false
        @test length(b) == 3
        @test nframes(b) == 3
        @test size(b) == (3,)
        @test b[[1, 2, 3]] == [6, 7, 8]

        push!(b, 9)
        @test isfull(b) == true
        @test isempty(b) == false
        @test length(b) == 3
        @test nframes(b) == 3
        @test size(b) == (3,)
        @test b[[1, 2, 3]] == [7, 8, 9]

        update!(b, 0)
        @test length(b) == 3
        @test b[[1, 2, 3]] == [7, 8, 0]

        update!(b, 1)
        @test length(b) == 3
        @test b[[1, 2, 3]] == [7, 8, 1]

        x = pop!(b)
        @test x == 1
        @test length(b) == 2
        @test b[[1, 2]] == [7, 8]

        x = pop!(b)
        @test x == 8
        @test length(b) == 1
        @test b[1] == 7

        x = pop!(b)
        @test x == 7
        @test length(b) == 0

        @test_throws ArgumentError pop!(b)
    end

    @testset "2D Float64" begin
        b = CircularArrayBuffer{Float64}(2, 2, 3)

        @test eltype(b) == Float64
        @test capacity(b) == 3
        @test isfull(b) == false
        @test length(b) == 0
        @test nframes(b) == 0
        @test size(b) == (2, 2, 0)

        for x in 1:3
            push!(b, x * A)
        end

        @test capacity(b) == 3
        @test isfull(b) == true
        @test nframes(b) == 3
        @test length(b) == 2 * 2 * 3
        @test size(b) == (2, 2, 3)
        for i in 1:3
            @test b[:, :, i] == i * A
        end
        @test b[:, :, end] == 3 * A

        for x in 4:5
            push!(b, x * ones(2, 2))  # collection is also OK
        end

        @test capacity(b) == 3
        @test length(b) == 2 * 2 * 3
        @test nframes(b) == 3
        @test size(b) == (2, 2, 3)
        @test b[:, :, 1] == 3 * A
        @test b[:, :, end] == 5 * A

        @test b == reshape([c for x in 3:5 for c in x * A], 2, 2, 3)

        push!(b, 6 * ones(Float32, 2, 2))
        push!(b, 7 * ones(Int, 2, 2))
        @test b == reshape([c for x in 5:7 for c in x * A], 2, 2, 3)
    end

    @testset "2D Float32" begin
        b = CircularArrayBuffer{Float32}(2, 2, 3)

        @test eltype(b) == Float32
        @test capacity(b) == 3
        @test isfull(b) == false
        @test length(b) == 0
        @test nframes(b) == 0
        @test size(b) == (2, 2, 0)

        for x in 1:3
            push!(b, x * C)
        end

        @test capacity(b) == 3
        @test isfull(b) == true
        @test nframes(b) == 3
        @test length(b) == 2 * 2 * 3
        @test size(b) == (2, 2, 3)
        for i in 1:3
            @test b[:, :, i] == i * C
        end
        @test b[:, :, end] == 3 * C

        for x in 4:5
            push!(b, x * ones(Float32, 2, 2))  # collection is also OK
        end

        @test capacity(b) == 3
        @test length(b) == 2 * 2 * 3
        @test nframes(b) == 3
        @test size(b) == (2, 2, 3)
        @test b[:, :, 1] == 3 * C
        @test b[:, :, end] == 5 * C

        @test b == reshape([c for x in 3:5 for c in x * C], 2, 2, 3)

        push!(b, 6 * ones(Float64, 2, 2))
        push!(b, 7 * ones(Int, 2, 2))
        @test b == reshape([c for x in 5:7 for c in x * C], 2, 2, 3)
    end
end


text:
using CalibrationErrors
using Literate: Literate
using Pkg: Pkg

if haskey(ENV, "GITHUB_ACTIONS")
    # Print `@debug` statements (https://github.com/JuliaDocs/Documenter.jl/issues/955)
    ENV["JULIA_DEBUG"] = "Documenter"
    # Bypass the accept download prompt
    ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"
end

const EXAMPLES = ("classification", "distribution")
const INPUT = joinpath(@__DIR__, "..", "examples")
const OUTPUT = joinpath(@__DIR__, "src", "examples")

ispath(OUTPUT) && rm(OUTPUT; recursive=true)
mkpath(joinpath(OUTPUT, "figures"))

# Add link to nbviewer below the first heading of level 1
function preprocess(content)
    sub = SubstitutionString(
        """
\\0
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/@__NAME__.ipynb)
#md #
# You are seeing the
#md # HTML output generated by [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl) and
#nb # notebook output generated by
# [Literate.jl](https://github.com/fredrikekre/Literate.jl) from the
# [Julia source file](@__REPO_ROOT_URL__/examples/@__NAME__/script.jl).
# The corresponding
#md # notebook can be viewed in [nbviewer](@__NBVIEWER_ROOT_URL__/examples/@__NAME__.ipynb).
#nb # HTML output can be viewed [here](https://devmotion.github.io/CalibrationErrors.jl/dev/examples/@__NAME__/).
#md #
#md # ```@setup @__NAME__
#md # using Pkg: Pkg
#md # Pkg.activate("$(INPUT)/@__NAME__")
#md # Pkg.instantiate()
#md # ```
#
        """,
    )
    return replace(content, r"^# # [^\n]*"m => sub; count=1)
end

for name in EXAMPLES
    file = joinpath(INPUT, name, "script.jl")

    # Activate project environment
    Pkg.activate(dirname(file)) do
        Pkg.instantiate()

        Literate.markdown(file, OUTPUT; name=name, documenter=true, preprocess=preprocess)
        Literate.notebook(file, OUTPUT; name=name, documenter=true, preprocess=preprocess)
    end
end

using Documenter

makedocs(;
    modules=[CalibrationErrors],
    authors="David Widmann <david.widmann@it.uu.se>",
    repo="https://github.com/devmotion/CalibrationErrors.jl/blob/{commit}{path}#L{line}",
    sitename="CalibrationErrors.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://devmotion.github.io/CalibrationErrors.jl",
        assets=String[],
    ),
    pages=[
        "index.md",
        "introduction.md",
        "ece.md",
        "kce.md",
        "others.md",
        "Examples" => [joinpath("examples", "$(name).md") for name in EXAMPLES],
    ],
    strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/devmotion/CalibrationErrors.jl.git",
    push_preview=true,
    devbranch="main",
)


text:
export net2dot

# visualize a net with GraphViz, by outputing dot description file
function net2dot(io::IO, net::Net)
  println(io, """digraph "$(net.name)" {""")
  println(io, """rankdir = BT;""")
  println(io, """node [shape=record,fontname="Fira Mono"];""")

  #------------------------------------------------------------
  # define layers
  for i = 1:length(net.layers)
    layer = net.layers[i]

    layer_type = replace("$(typeof(layer))", r"Layer$", "")
    layer_name = layer.name
    neuron_name = ""
    if in(:neuron, fieldnames(layer)) && !isa(layer.neuron, Neurons.Identity)
      neuron_name = "\\n($(typeof(layer.neuron)))"
    end

    # layer name
    print(io, """layer$i [label="$layer_type\\n$layer_name$neuron_name|{top:|bot:}|""")
    print(io, """{""")

    # top blobs
    if in(:tops, fieldnames(layer))
      print(io, """{""")
      for j = 1:length(layer.tops)
        if j != 1
          print(io, "|")
        end
        print(io, "<t$j> $(layer.tops[j])")
        print(io, " ($(join(map(x -> "$x", size(net.states[i].blobs[j])), "x")))")
      end
      print(io, """}""")
    end

    print(io, """|""")

    # bottom blobs
    if in(:bottoms, fieldnames(layer))
      print(io, """{""")
      for j = 1:length(layer.bottoms)
        if j != 1
          print(io, "|")
        end
        print(io, "<b$j> $(layer.bottoms[j])")
        print(io, " ($(join(map(x -> "$x", size(net.blobs_forward[i][j])), "x")))")
      end
      print(io, """}""")
    end

    println(io, """}"];""")
  end

  #------------------------------------------------------------
  # Connect blobs
  top_maps = Dict{Symbol,AbstractString}()
  for i = 1:length(net.layers)
    if in(:tops, fieldnames(net.layers[i]))
      for j = 1:length(net.layers[i].tops)
        top_maps[net.layers[i].tops[j]] = "layer$i:t$j"
      end
    end
  end
  for i = 1:length(net.layers)
    layer = net.layers[i]
    if in(:bottoms, fieldnames(layer))
      for j = 1:length(layer.bottoms)
        println(io, "$(top_maps[layer.bottoms[j]]) ->layer$i:b$j;")
      end
    end
  end

  println(io, "}") # digraph
end


text:
# copied from https://github.com/JuliaLang/julia/blob/v1.5.4/stdlib/Logging/src/ConsoleLogger.jl
using Logging
import Logging: min_enabled_level, shouldlog, catch_exceptions, handle_message

function shorten_str(message, max_len)
    suffix = "…"
    if length(message) > max_len
        message[1:min(end, max_len-length(suffix))] * suffix
    else
        message
    end
end

function make_log_truncated(max_var_len=5_000)
    function log_truncated(io, args)
        levelstr = args.level == Logging.Warn ? "Warning" : string(args.level)
        msglines = split(chomp(string(shorten_str(args.message, max_var_len))), '\n')
        println(io, "┌ ", levelstr, ": ", msglines[1])
        for i in 2:length(msglines)
            str_line = sprint(print, "│ ", msglines[i])
            println(io, shorten_str(str_line, max_var_len))
        end
        for (key, val) in args.kwargs
            str_line = sprint(print, "│   ", key, " = ", val)
            println(io, shorten_str(str_line, max_var_len))
        end
        println(io, "└ @ ", something(args._module, "nothing"), " ",
                something(args.file, "nothing"), ":", something(args.line, "nothing"))
        nothing
    end
    log_truncated
end


text:
## test VB GMM for the standard example data
x = readdlm("faithful.txt")
## only do k-means
g = GMM(8, x, kind=:full, nIter=0)
p = GMMprior(g.d, 1.0, 1.0)
v = VGMM(g, p)
em!(v, x)
## show the reduction in Gaussians
println(history(v))
## we have no pretty printing yet...
for f in [:α, :β, :m, :ν, :W]
    println(f, " = ", getfield(v, f))
end


text:
# module ModifiedNodesLagRefFEs
#
# using Gridap

function ScaledNodesLagrangianRefFE(::Type{T}, p::Polytope{D}, orders::Vector{Int}, scaling) where {D,T}
  @assert length(orders) == D
  @assert D > 0
  if !(all(extrusion(p).array.== HEX_AXIS) || all(orders.==orders[1]))
    error("One can consider anisotropic orders on n-cubes only")
  end
  if (all(orders.==1))
    nodes, nfacenodes = Gridap.RefFEs._linear_lagrangian_nodes_polytope(p)
  else
    nodes, nfacenodes = Gridap.RefFEs._high_order_lagrangian_nodes_polytope(p,orders)
  end
  nodes = nodes*scaling
  dofsb = Gridap.RefFEs.LagrangianDOFBasis(T,nodes)
  # @santiagobadia : This part is missing for anisotropic orders !
  prebasis = Gridap.RefFEs._monomial_basis(p,T,orders[1])
  aux = zeros(Float64,numlocaldofs(dofsb),numlocaldofs(dofsb))
  @assert numlocaldofs(dofsb) == length(prebasis)
  evaluate!(dofsb,prebasis,aux)
  changeofbasis=inv(aux)
  basis = change_basis(prebasis, changeofbasis)
  nfacedofs = Gridap.RefFEs._genereate_nface_to_dofs(nfacenodes,dofsb.node_and_comp_to_dof)
  return LagrangianRefFE{D,T}(p, dofsb, basis, nfacedofs, nfacenodes)
end

function ScaledNodesLagrangianRefFE(::Type{T}, p::Polytope{D}, order::Int, scaling) where {D,T}
  _order = order*ones(Int,D)
  return ScaledNodesLagrangianRefFE(T,p,_order,scaling)
end

# F = Gridap.RefFEs
#
# T = Float64
# D = 1
# model = CartesianDiscreteModel(partition=(1,))
# diri = [1]
# p = Polytope(1,)
# order = 2
# reffe1 = LagrangianRefFE(Float64,p,order)
# orders = order*ones(Int,D)
# reffe2 = ScaledLagrangianRefFE(Float64,p,orders,0.5)
#
#
# end # module