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

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 el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 89fcfe1c-913a-11eb-1b04-0d87fd8dbcbe
#Set up the working environment
begin
	using Pkg
	Pkg.activate("..")
	using Plots, PlutoUI,LaTeXStrings
end

# ╔═╡ 9f7c1628-9139-11eb-0d69-f910a995eca7
md"""## Differential Equations

When modelling growth, we rely on **differential equations** that describe the change in some variable $y$ with respect to some other variable $x$.

$$\frac{dy}{dx} = f(y)$$

Here, $\frac{dy}{dx}$ represents the **rate of change in $y$ with respect to $x$**, and $f(y)$ is some function that depends on $y$.

Solving differential equations to obtain a function that expresses $y$ in terms of $x$ can be relatively straight-forward, but this is generally not the case. So why bother? It turns out that in the natural world for many processes *it tends to be much easier to describe what the change in some quantity depends on*, rather than coming up with an *explicit* formula how $y$ changes as $x$ varies proves difficult.

To make the mathematical language a bit less abstract, let us rewrite the equation above, so that the the variables more directly correspond to what we are talking about when modelling growth:

$$(1)\frac{dP}{dt} = f(P)$$

Here, $\frac{dP}{dt}$ denotes the growth rate of some population $P$ over time, which depends on the population size $P$ in a yet-to-be specified manner $f(P)$."""

# ╔═╡ babc9748-9139-11eb-0da1-0d942d591481
md"""## The simplest possible case

To familiarise ourselves with growth models, let's look at a very simple model, while ignoring the fact that the model does not make sense when modelling biological populations.

$$(2)\frac{dP}{dt} = r$$

Here, our growth rate $\frac{dP}{dt}$  is constant and equal to a growth rate constant $r$.

### Rate Plots
A useful and simple trick to understand how systems described by differential equations behave involves looking at the **rate plot**, where we plot how $\frac{dP}{dt}$  changes as the population size $P$ varies.

> **Exercise:** Sketch a rate plot for the system described by equation (2)."""

# ╔═╡ ee2bb580-9145-11eb-2bc2-dbd9aae68be9
md"""
Since the rate of growth is constant,  changes by the same amount in the same timespan, regardless of how much time has already passed. Thus the graph of  versus  is a straight line that intersects the y-Axis at  and grows with a slope of .
"""

# ╔═╡ 001fc24a-9146-11eb-39b9-955dac6bcf58
md"""
### Solving a simple differential equation
Nice - We've analysed a simple differential equation using a rate plot and deduced an approximate plot of $P$ over time from it. Now on to the somewhat harder part - solving the differential equation mathematically...

In general, differential equations can be solved by separating the variables ($P$ and $t$ in our case) and then integrating both sides with respect to the variables. If you're unfamiliar with integration, but are familiar with differentiation, as taught in highschool, it should suffice to say that integration is essentially differentiation in reverse.

Upon separating our variables in equation (2), we get

$$(3) \int dP = r \int dt $$

Doing differentiation in reverse for each side, we find

$$(4) P + C_1 = rt + C_2 $$

We find the derivative of a variable by lowering its exponent by 1.

$$(5) f(x) = x^n, f'(x) = x^{n-1}$$

Hence, integrating the expression above (which can be viewed as containing $P^0=1$ and $t^0=1$, respectively) with respect to some variable amounts to raising the exponents of our "invisible" variables by 1. You may wonder where the two constants come from. They need to be included, because constants disappear in differentiation and this needs to be accounted for by creating them in the course of integration. However, we don't need to worry too much about the constants for now. In fact, they behave a bit like plasticine and can be liberally lumped together for now, as we only need to determine the constant we're left with once we're done with our arithmetic.

By absorbing $C_1$ into $C_2$, we get a new constant $C = C_1 - C_2$, and equation (4) rearranges to

$$(6) P(t) = rt + C $$

Voilá, here's an equation for a straight line!

To figure out the value of our constant $C$, we need one more trick: We set $P(0)=P_0$ and $t=0$, and we get

$$(7) P_0 = C$$

Hence our final equation is

$$(8) P(t) = rt + P_0$$
"""

# ╔═╡ ad08de50-913d-11eb-3f87-fb27a282bdb3
md"""### Interactive Simulation
*Change the sliders to adjust the corresponding values*:
- growth rate constant ($r$) $(@bind r_0 Slider(-0.5:0.01:0.5, default=0.01, show_value=true))
- initial population size($P_0$) $(@bind P0_0 Slider(0:5:100, default=20, show_value=true))
"""

# ╔═╡ f58d273e-9139-11eb-0fe0-b11e90166682
begin
	# Our differential equation
	dPdt_0(P, r) = r
	# A vector holding different population size values
	population_0 = 0:0.5:100
	# A vector holding our growth rate values for different population sizes
	growthrate_0 = [dPdt_0(p,r_0) for p in population_0]

	# The solution to our differential equation
	P_0(t, r, P0) = P0 + r*t
	# A vector holding values for different timepoints
	time_0 = 0:0.5:50
	# A vector holding our simulated population sizes
	sim_pop_0 = [P_0(t,r_0,P0_0) for t in time_0]

	rateplot_0 = plot(population_0, growthrate_0,
				xlims = (0,100),
				ylims = (-1.2, 1.2),
				title = "Rateplot",
				label = "",
				xlabel = "P", ylabel = L"$\frac{dP}{dt}$")
	simplot_0 = plot(time_0, sim_pop_0,
				color = "red",
				ylims = (0,100),
				title = "Simulation",
				label = "",
				xlabel = "t", ylabel = "P")
	plot(rateplot_0, simplot_0, layout = (1, 2),fontfamily="Computer Modern")
end

# ╔═╡ aa53e7e0-913d-11eb-32ab-4fc1ce2b72ab
md"""
## 2. A more realistic case

The growth of most biological populations will at the very least depend on the number of individuals already present. After all, having twice as many bacteria that could potentially divide, we'd expect the growth rate to increase accordingly. To capture this relationship, we can write

$$(9) \frac{dP}{dt} = rP$$

> **Exercise:** Sketch a rate plot for the system described by equation (9).

From our rate plot, we can see that in this case, the growth rate increases linearly as our population grows.

> **Exercise:** Try to sketch a plot that shows how $P$ changes over time $t$, starting from an arbitrarily chosen initial population size $P_0$.

Maybe you were able to figure out that since the growth rate increases as $P$ grows, the growth of the population should accelerate in exponential manner. To get an explicit formula for $P(t)$, and double-check our intuition, let us use the methods we've learned above.

Separation of variables in equation (9) yields

$$(10) \int \frac{1}{P} dP = r \int dt$$

Integrating equation 10 (and maybe consulting a trusty resource for common antiderivatives to find $\int \frac{1}{x} dx = ln(x) + C$), we get

$$(11) ln(P) + C_1 = rt + C_2$$

Again, we can absorn $C_1$ and $C_2$ into a new constant $C = C_2 - C_1$ to obtain

$$(12) ln(P) = rt + C $$

Now, to extract our population size $P$ from the somewhat unyieldy $ln(P)$, we make use of the following identities $e^{ln(x)} = x$, $e^{a+b} = e^a \cdot e^b$ and obtain:

$$(13) P(t) = e^{rt + C} = e^C \cdot e^{rt}$$

Exploiting the plasticine-like qualities of our integration constant and settin $e^C = A$, we clean up equation (13) and get

$$(14) P(t) = A \cdot e^{rt}$$

To figure out the value of $A$, we employ the same trick as above and set $t = 0$ and $P(0) = P_0$.

$$(15) P_0 = A \cdot e^{r \cdot 0} = A$$

Hence our solution to equation (9) is

$$(16) P(t) = P_0 \cdot e^{rt}$$

This equation ought to look familiar. It is the equation for exponential growth (and decay, if $r<0$) and abounds in science (e.g. radioactive decay). Yet, to model populations in biology, it is of very limited use, as it suggests a world of unlimited resources. Eventually, any population is going to encounter limiting factors - it may run out of space or food, or may find itself fall prey to predators or parasites.
"""

# ╔═╡ 008851e2-9147-11eb-3b1c-0ff07fd4e84f
md"""### Interactive Simulation
*Change the sliders to adjust the corresponding values*:
- growth rate constant ($r$) $(@bind r_1 Slider(-0.25:0.01:0.25, default=0.01, show_value=true))
- initial population size($P_0$) $(@bind P0_1 Slider(0:5:100, default=20, show_value=true))
"""

# ╔═╡ 20245898-9147-11eb-16cb-7320a2ffbe0d
begin
	# Our differential equation
	dPdt_1(P, r) = r*P
	# A vector holding different population size values
	population_1 = 0:0.5:100
	# A vector holding our growth rate values for different population sizes
	growthrate_1 = [dPdt_1(p,r_1) for p in population_1]

	# The solution to our differential equation
	P_1(t, r, P0) = P0*exp(r*t)
	# A vector holding values for different timepoints
	time_1 = 0:0.5:50
	# A vector holding our simulated population sizes
	sim_pop_1 = [P_1(t,r_1,P0_1) for t in time_1]

	rateplot_1 = plot(population_1, growthrate_1,
				xlims = (0,100),
				ylims = (-1.2, 1.2),
				title = "Rateplot",
				label = "",
				xlabel = "P", ylabel = L"$\frac{dP}{dt}$")
	simplot_1 = plot(time_1, sim_pop_1,
				color = "red",
				ylims = (0,100),
				title = "Simulation",
				label = "",
				xlabel = "t", ylabel = "P")
	plot(rateplot_1, simplot_1, layout = (1, 2),fontfamily="Computer Modern")
end

# ╔═╡ 24f0552e-9148-11eb-1226-c326de038a70
md"""
## 3. The Logistic Model

Without being to explicit about the precise nature of our limiting factor, let us assume the population size be limited to a maximum number that can be supported in a given habitat. We call this the **carrying capacity, $K$**.

$$(17) \frac{dP}{dt} = rP(1-\frac{P}{K})$$

### Visual Analysis
If our population is close to $0$, it will won't "feel" any limiting effect and grow in a similar manner to the exponential model discussed above. As the population gets closer to the carrying capacity, the growth will begin to slow down.
"""

# ╔═╡ 3d7b3688-9148-11eb-3224-9b832c4692d9
md"""
*Change the sliders to adjust the corresponding values*:
- growth rate constant ($r$) $(@bind r_2 Slider(0:0.001:0.25, default=0.01, show_value=true))
- carrying capacity ($K$) $(@bind K_2 Slider(1:5:80, default=30, show_value=true))
- initial population size($P_0$) $(@bind P0_2 Slider(1:5:100, default=20, show_value=true))
"""

# ╔═╡ 5fbb5af8-9148-11eb-0ade-b3fd39af29c1
begin
	# Our differential equation
	dPdt_2(P, r, K) = r*P*(1- P/K)
	# A vector holding different population size values
	population_2 = 0:0.5:100
	# A vector holding our growth rate values for different population sizes
	growthrate_2 = [dPdt_2(p,r_2, K_2) for p in population_2]

	# The solution to our differential equation
	P_2(t, r, K,  P0) = K/((K-P0)/P0*exp(-r*K*t) + 1)
	# A vector holding values for different timepoints
	time_2 = 0:0.5:50
	# A vector holding our simulated population sizes
	sim_pop_2 = [P_2(t,r_2/K_2,K_2, P0_2) for t in time_2]

	rateplot_2 = plot(population_2, growthrate_2,
				xlims = (0,100),
				ylims = (-1.2, 2.5),
				title = "Rateplot",
				label = "",
				xlabel = "P", ylabel = L"$\frac{dP}{dt}$")
		plot!(rateplot_2, [K_2 for p in -2:3] , -2:3,
				color = "darkgrey",
				linestyle = :dash,
				label = "";
				annotatios = (K_2, 1, "P=K"))
		plot!(rateplot_2, population_2, [r_2 * p for p in population_2],
				color = "grey",
				linestyle = :dot,
				label = L"$\frac{dP}{dt} = rP$")

	simplot_2 = plot(time_2, sim_pop_2,
				color = "red",
				ylims = (0,100),
				title = "Simulation",
				label = "",
				xlabel = "t", ylabel = "P")
		plot!(simplot_2, time_2 , [K_2 for t in time_2],
				color = "darkgrey",
				linestyle = :dash,
				label = "",
				annotatios = ([40], [K_2], "P=K"))
	plot(rateplot_2, simplot_2, layout = (1, 2),fontfamily="Computer Modern")
end

# ╔═╡ 64b31f38-9149-11eb-11a0-718a61330e76
md"""
The corresponding rate plot looks like a prabola that crosses the x-axis at $P=0$ and $P=K$ and has a maximum at $P=\frac{K}{2}$. The x-crossings are important. Here, the growth rate is $0$, so they represent **equilibria**.  However, the two equilibria are not the same: Just to the right of $P=0$, the growth rate increases, so the equilibrium is **unstable**. The case for the equilibrium at $P=K$ is different: if $P<K$, the growth rate is positive, whereas the growth rate is negative if $P>K$. Consequently, the population will be attracted to $P=K$ from both sides. This is called a **stable equilibrium** - small fluctuations around $P=K$ will result in the population falling back to $P=K$.

Using the insight gained from the rate plot, we can infer the shape of the curve for $P(t)$. It will be *s-shaped* or *sigmoidal*, reflecting the change in growth rate around $P=\frac{K}{2}$, and level off around $P=K$.
"""

# ╔═╡ 6f8ed9ce-9149-11eb-2c08-1b1e1f3e3f5c
md"""
### Analytic solution
To obtain an explicit formula for $P(t)$, we need to marshall what we've learned so far, integrate, and brace ourselves for quite a bit of algebra...

To make our lives a bit easier, we rearrange the equation as follows

$$(18) \frac{dP}{dt} = r_{max}P(1-\frac{P}{K}) = r_{max}P(\frac{K}{K}-\frac{P}{K})= \frac{r_{max}}{K}P(K-P) = rP(K-P)$$

Next, we separate our variables:

$$(19) \int \frac{1}{P(K-P)} dP = r \int dt$$

To solve the integral on the left hand side, it helps to simplify the fraction using a *partial fraction decomposition*. Here we want to decompose our fraction $\frac{1}{P(K-P)}$ into two partial fractions of the form $\frac{A}{P} + \frac{B}{K-P}$, where $A$ and $B$ are constants that can be easily determined as follows:

$$\frac{A}{P} + \frac{B}{K-P} = \frac{1}{P(K-P)}$$

Multiply the two fractions on the left by $P$ and $(K-P)$ respectively to obtain a common denominator with the fraction on the right

$$\frac{A(K-P)}{P(K-P)} + \frac{BP}{P(K-P)} = \frac{1}{P(K-P)}$$
Multiply the whole expression by $P(K-P)$ and rearrange to find

$$AK + P(B-A) = 1 \rightarrow A,B = \frac{1}{K}$$

Hence,

$$\frac{1}{K} (\frac{1}{P} + \frac{1}{K-P}) = \frac{1}{P(K-P)}$$

Substituting the partial fraction decomposition into equation (19), we obtain
$$(20) \frac{1}{K} \int \frac{1}{P} + \frac{1}{K-P} dP = r \int dt$$

This way of rearranging our expression allows us to treat each fraction separately

$$(21) \frac{1}{K} (\int \frac{1}{P} dp + \int \frac{1}{K-P} dP) = r \int dt$$

###### Finding integrals by substitution
We already know how to integrate the first fraction - $\int \frac{1}{P} dP = ln(P) + C_1$, as well as the right hand side - $r \int dt = rt + C_2$. However, to integrate the second fraction, we need to use substitution.
Let

$$(22) u = K-P$$

It follows that

$$(23)\frac{du}{dP} = -1 \rightarrow  du = -dP \rightarrow dP = - du$$

We substitute equations (21) and (23) into $\int \frac{1}{K-P} dP$ and obtain$-\int \frac{1}{u} du$, which can be easily integrated and gives $-ln(u) +C =-ln(K-P) + C_3$. Plugging our solved integrals into equation (21) and absorbing our various constants into a new constant $B$, we get

$$(24) \frac{1}{K}[ln(P) - ln(K-P)] = rt + B$$

Using the identity $ln(a) - ln(b) = ln(\frac{a}{b})$, this rearranges to

$$(25) ln(\frac{P}{K-P}) = rKt + B$$

Applying the exponential function to both sides yields

$$(26) \frac{P}{K-P} = e^{rKt + B} = e^{rKt} \cdot e^{B}$$

Let $A = e^B$:

$$(27) \frac{P}{K-P} = Ae^{rKt}$$

###### Finding the value of the constant A

Let $P_0 = P(0)$ and $t=0$:

$$(28) \frac{P_0}{K-P_0} = A$$

Hence,

$$(29) \frac{P}{K-P} = \frac{P_0}{K-P_0}e^{rKt}$$

###### Solve expression for P

$$(30) P(K-P_0) = P_0(K-P)e^{rKt} = P_0Ke^{rKt} - P_0Pe^{rKt}$$

Next, we collect the terms that depend on $P$:

$$(31) P(K-P_0) + P_0Pe^{rKt} = P_0Ke^{rKt}$$

$$(32) P[(K-P_0) + P_0e^{rKt}] = P_0Ke^{rKt}$$

$$(33) P = \frac{P_0Ke^{rKt}}{(K-P_0) + P_0e^{rKt}}$$

This can be rewritten as

$$(34) P(t) = \frac{K}{\frac{K-P_0}{P_0}e^{-rKt} + 1}$$

## Conclusion
And that is that. Solving differential equations this way is possible for some, but not all. The algebra gets hairy quite quickly, particularly, when looking at systems of differential equations, where differential equations for two or more variables depend on each other. In such instances, we obtain timecourses by simulating these systems numerically. More on this soon...

"""

# ╔═╡ Cell order:
# ╟─89fcfe1c-913a-11eb-1b04-0d87fd8dbcbe
# ╟─9f7c1628-9139-11eb-0d69-f910a995eca7
# ╟─babc9748-9139-11eb-0da1-0d942d591481
# ╟─ee2bb580-9145-11eb-2bc2-dbd9aae68be9
# ╟─001fc24a-9146-11eb-39b9-955dac6bcf58
# ╟─ad08de50-913d-11eb-3f87-fb27a282bdb3
# ╟─f58d273e-9139-11eb-0fe0-b11e90166682
# ╟─aa53e7e0-913d-11eb-32ab-4fc1ce2b72ab
# ╟─008851e2-9147-11eb-3b1c-0ff07fd4e84f
# ╟─20245898-9147-11eb-16cb-7320a2ffbe0d
# ╟─24f0552e-9148-11eb-1226-c326de038a70
# ╟─3d7b3688-9148-11eb-3224-9b832c4692d9
# ╟─5fbb5af8-9148-11eb-0ade-b3fd39af29c1
# ╟─64b31f38-9149-11eb-11a0-718a61330e76
# ╟─6f8ed9ce-9149-11eb-2c08-1b1e1f3e3f5c

max_stars_repo_name: el-uhu/intro-to-modelling

max_issues_repo_name: el-uhu/intro-to-modelling

max_forks_repo_name: el-uhu/intro-to-modelling


text:
export increase_difficulty, decrease_difficulty, shuffle!
function shuffle!(specset::DataFrame)

    # Set all priorities to be random

    for i = 1:size(specset)[1]
        specset[i,:priority] = rand()
    end

    # Sort
    sort!(specset, :priority, rev=true)
end

function shuffle!(superset)
    for i = 1:size(superset)[1]
        superset[i][2] = rand()
        shuffle!(superset[i][1])
    end

    SwoLeios.sort_sets!(superset)

    output_set(superset)
end

function increase_difficulty(set::NamedTuple)
    new_set = (exercise = set[:exercise], reps_0 = set[:reps_0],
               reps_10 = set[:reps_10], difficulty = set[:difficulty] + 1,
               priority = set[:priority])
    output_set(new_set, select_reps(new_set))
    return new_set
end

function decrease_difficulty(set::NamedTuple)
    new_set = (exercise = set[:exercise], reps_0 = set[:reps_0],
               reps_10 = set[:reps_10], difficulty = set[:difficulty] - 1,
               priority = set[:priority])
    output_set(new_set, select_reps(new_set))
    return new_set
end

function increase_difficulty!(set::DataFrameRow)
    set[:difficulty] += 1
end

function decrease_difficulty!(set::DataFrameRow)
    set[:difficulty] -= 1
end

ext: jl

max_stars_repo_name: leios/SwoLeios

max_issues_repo_name: leios/SwoLeios

max_forks_repo_name: leios/SwoLeios


text:
@testset "Parallel.Bellman Ford" begin
    g4 = PathDiGraph(5)

    d1 = float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])
    d2 = sparse(float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]))
    for g in testdigraphs(g4)
        y = @inferred(Parallel.bellman_ford_shortest_paths(g, [2], d1))
        z = @inferred(Parallel.bellman_ford_shortest_paths(g, [2], d2))
        @test y.dists == z.dists == [Inf, 0, 6, 17, 33]
        @test @inferred(enumerate_paths(z))[2] == []
        @test @inferred(enumerate_paths(z))[4] == enumerate_paths(z, 4) == [2, 3, 4]
        @test @inferred(!Parallel.has_negative_edge_cycle(g, d1))

        y = @inferred(Parallel.bellman_ford_shortest_paths(g, [2], d1))
        z = @inferred(Parallel.bellman_ford_shortest_paths(g, [2], d2))
        @test y.dists == z.dists == [Inf, 0, 6, 17, 33]
        @test @inferred(enumerate_paths(z))[2] == []
        @test @inferred(enumerate_paths(z))[4] == enumerate_paths(z, 4) == [2, 3, 4]
        z = @inferred(Parallel.bellman_ford_shortest_paths(g, [2]))
        @test z.dists == [typemax(Int), 0, 1, 2, 3]
    end

    # Negative Cycle
    gx = CompleteGraph(3)
    for g in testgraphs(gx)
        d = [1 -3 1; -3 1 1; 1 1 1]
        @test_throws LightGraphs.NegativeCycleError Parallel.bellman_ford_shortest_paths(g, [1], d)
        @test Parallel.has_negative_edge_cycle(g, d)

        d = [1 -1 1; -1 1 1; 1 1 1]
        @test_throws LightGraphs.NegativeCycleError Parallel.bellman_ford_shortest_paths(g, [1], d)
        @test Parallel.has_negative_edge_cycle(g, d)
    end

    # Negative cycle of length 3 in graph of diameter 4
    gx = CompleteGraph(4)
    d = [1 -1 1 1; 1 1 1 -1; 1 1 1 1; 1 1 1 1]
    for g in testgraphs(gx)
        @test_throws LightGraphs.NegativeCycleError Parallel.bellman_ford_shortest_paths(g, [1], d)
        @test Parallel.has_negative_edge_cycle(g, d)
    end
end


text:
@testset "GravityEffect/layerResponse write" begin
	sensor=Dict(:x => 100., :y => 150., :z => 100., :sensHeight => 1.0)
	layers=Dict(:start => [0.0, 1.0],:stop  => [1.0, 2.0])
	dem_in = joinpath(dirname(@__DIR__),"test","input","dem_data.asc");
	zones =Dict(:dem   => [dem_in,dem_in], :radius=> [50.,200.],
				:resolution => [0.1,0.2],:interpAltitude => [true,false]);
	exclude = Dict(:cylinder=>Dict(:radius=>1.,:start=>0.5,:stop=>2.),
			   :prism=>Dict(:x=>[100.,106.],:y=>[155.,156.],
							:dx=>[2.,3.],:dy=>[1.5,1.8],
							:start=>[0.,1.],:stop=>[1.,2.]),
			   :polygon=>Dict(:file=>joinpath(dirname(@__DIR__),"test","input","exclusion_polygon.txt"),
							  :start=>0.,:stop=>2.))
	outputfile = joinpath(dirname(@__DIR__),"test","output","write_layerResponse.txt");
	#outdata = GravityEffect.layerResponse(sensor,layers,zones,
	#					exclude=exclude,nanheight=true,
	#					outfile=outputfile,def_density=10.);
	outdata = DataFrame(layer=[1,2],start=layers[:start],stop=layers[:stop],
	 		   	   total=[3.82615,3.76013].*1e-9,zone1=[3.76682,3.66139].*1e-9,
				   zone2=[0.05933,0.09874].*1e-9);
    FileTools.write_layerResponse(sensor,layers,zones,exclude,false,
			   		    outputfile,10.,outdata);
	t = readdlm(outputfile,comments=true,comment_char='%');
	@test size(t) == (2,6)
	@test all(isapprox.(t[:,4],outdata[:total].*1e+9))
	@test all(isapprox.(t[:,5],outdata[:zone1].*1e+9))
	@test all(isapprox.(t[:,6],outdata[:zone2].*1e+9))
	@test t[:,1]==[1.,2.]
	@test t[:,2]==layers[:start]
	@test t[:,3]==layers[:stop]
end

@testset "GravityEffect/layerResponse read" begin
	filein = joinpath(dirname(@__DIR__),"test","input","read_layerResponse.txt");
	t = read_layerResponse(filein,"results")
	@test t[:layer]==[1.,2.]
	@test t[:start]==[0.,1.]
	@test t[:stop]==[1.,2.]
	@test t[:total]==[4.,3.]
	@test t[:zone1]==[3.1,2.2]
	@test t[:zone2]==[0.9,0.8]

	t = read_layerResponse(filein,"layers")
	@test t == Dict(:start => [0.0, 1.0],:stop  => [1.0, 2.0])
	t = read_layerResponse(filein,"sensor")
	@test t == Dict(:x => 100., :y => 150., :z => 100., :sensHeight => 1.0)
	t = read_layerResponse(filein,"def_density")
	@test t == 10.
	t = read_layerResponse(filein,"nanheight")
	@test !t

	dem_in = joinpath(dirname(@__DIR__),"test","input","dem_data.asc");
	zones =Dict(:dem   => [dem_in,dem_in], :radius=> [50.,200.],
				:resolution => [0.1,0.2],:interpAltitude => [true,false]);
	t = read_layerResponse(filein,"zones")
	#@test t == zones

	exclude = Dict(:cylinder=>Dict(:radius=>1.,:start=>0.5,:stop=>2.),
			   :prism=>Dict(:x=>[100.,106.],:y=>[155.,156.],
							:dx=>[2.,3.],:dy=>[1.5,1.8],
							:start=>[0.,1.],:stop=>[1.,2.]),
			   :polygon=>Dict(:file=>joinpath(dirname(@__DIR__),"test","input","exclusion_polygon.txt"),
							  :start=>0.,:stop=>2.))
	t = read_layerResponse(filein,"exclude")
	#@test t == exclude
end


text:
module sync

import io
import struct

self.INFINITE = 2147483647
self.LOCK_ACQUIRED = 0

class AtomicInt {
  field mutex
  new(field i) {
    mutex = Mutex()
  }
  method set(v) {
    mutex.acquire()
    i = v
    mutex.release()
  }
  method get() {
    mutex.acquire()
    cpy = i
    mutex.release()
    return cpy
  }
  method inc(a) {
    mutex.acquire()
    i = i + a
    mutex.release()
  }
  method to_s() {
    mutex.acquire()
    tos = str(i)
    mutex.release()
    return tos
  }
}

class ThreadPool {
  field mutex, task_lock, threads, tasks
  new(field num_threads) {
    mutex = Mutex()
    task_lock = Semaphore(0, num_threads)
    threads = []
    tasks = []
    for i=0, i<num_threads, i=i+1 {
      threads.append(Thread(_execute_task, None))
    }
    for (_,t) in threads {
      t.start()
    }
  }
  method _execute_task() {
    task = None
    while True {
      task_lock.lock(INFINITE)
      mutex.acquire()
      if tasks.len > 0 {
        task = tasks.pop(1)[0]
      } else {
        task_lock.unlock()
      }
      mutex.release()

      if task{
        task.exec()
        task = None
      }
    }
  }
  method _execute_future(f) {
    mutex.acquire()
    if f.mark {
      mutex.release()
      return f
    }
    tasks.append(f)
    f.mark = True
    mutex.release()

    task_lock.unlock()
    return f
  }
  method execute(fn, args=None) {
    f = None
    if fn is TaskGraph {
      f = TaskGraphFuture(fn, self)
    } else {
      f = Future(fn, args, self)
    }
    return _execute_future(f)
  }
}

class Future {
  field listeners, listeners_mutex, result, has_result, read_mutex, mark
  new(field fn, field args, field ex) {
    listeners = []
    listeners_mutex = Mutex()
    result = None
    has_result = False
    read_mutex = Semaphore(0, 1)
    mark = False
  }
  method exec() {
    result = fn(args)
    listeners_mutex.acquire()
    has_result = True
    for (_,f) in listeners {
      f.args = result
      ex._execute_future(f)
    }
    listeners_mutex.release()
    read_mutex.unlock()
  }
  method get() {
    while True {
      read_mutex.lock(INFINITE)
      if ~has_result {
        read_mutex.unlock()
        continue
      }
      read_mutex.unlock()
      return result
    }
  }
  method then_do(fn) {
    return _then_do_future(Future(fn, None, ex))
  }
  method _then_do_future(f) {
    listeners_mutex.acquire()
    if has_result {
      f.args = result
      ex._execute_future(f)
    } else {
      listeners.append(f)
    }
    listeners_mutex.release()
    return f
  }
  method is_complete() {
    read_mutex.lock(0)
    res = has_result
    read_mutex.unlock()
    return res
  }
  method wait() {
    get()
  }
}

class TaskGraph {
  field next_id, nodes, roots
  new() {
    next_id = 0
    nodes = []
    roots = []
  }
  method add_node(work, deps=None) {
    id = _next_id()
    if ~deps {
      roots.append(id)
    }
    deps = _wrap_deps(deps)
    node = _WorkNode(id, deps, work)
    nodes.append(node)
    deps.each(node_id -> nodes[node_id]._add_dep_for(id))
    return id
  }
  method _next_id() {
    retval = next_id
    next_id = next_id + 1
    return retval
  }
  method _wrap_deps(deps) {
    if ~deps {
      return []
    } else if deps is Array {
      return deps
    } else {
      return [deps]
    }
  }
}

class _WorkNode {
  field deps_for
  new(field id, field depends_on, field work) {
    deps_for = []
  }
  method _add_dep_for(dep_id) {
    deps_for.append(dep_id)
  }
  method to_s() {
    return concat('_WorkNode(id=', id, ',depends_on=', depends_on, ',deps_for=', deps_for, ')')
  }
}

class TaskGraphFuture : Future {
  new(field task_graph, field ex) {
    self.Future.new(_execute_graph, None, ex)
  }
  ; Why do I have to do this?
  method get() {
    return self.Future.get().map(f -> f.get())
  }
  method exec() {
    return self.Future.exec()
  }
  method _execute_graph() {
    visited = []
    visited[task_graph.nodes.len - 1] = None
    _create_futures(visited)
    task_graph.roots.each(id -> ex._execute_future(visited[id]))
    return visited
  }
  method _create_futures(visited) {
    next_to_process = task_graph.roots.copy()
    while next_to_process.len > 0 {
      id = next_to_process.pop()
      ; Skip if already visited
      if visited[id] {
        continue
      }
      node = task_graph.nodes[id]
      future = _WorkNodeFuture(node, visited, task_graph, ex)
      visited[id] = future
      node.deps_for.each(next_to_process.append)
    }
    return visited
  }
}

class _WorkNodeFuture : Future {
  new(field node, field visited, field task_graph, field ex) {
    self.Future.new(node.work, None, ex)
  }
  method get() {
    return self.Future.get()
  }
  method exec() {
    result_map = {}
    for (_, dep) in node.depends_on {
      result_map[dep] = visited[dep].get()
    }
    self.Future.args = result_map
    self.Future.exec()
    node.deps_for.each((id) {
      child = task_graph.nodes[id]
      child_future = visited[id]
      all_deps_completed = child.depends_on.stream()
                                .map(id -> visited[id].is_complete())
                                .collect((a,b) -> a and b)
      if all_deps_completed {
        ex._execute_future(child_future)
      }
    })
  }
  method is_complete() {
    return self.Future.is_complete()
  }
}


text:
const continues = ["c", "Continue", "C", "continue"]

function c_to_continue()
  while true
    println("Press C/c and then Enter to continue:")
    ans = readline()
    ans ∈ continues && return
  end
end

# Visualize a contraction and then perform the contraction
function contract_visualize(As::ITensor...; contract = true, pause = false, kwargs...)
  scene = visualize_tensornetwork(As...; kwargs...)
  display(scene)
  if pause
    c_to_continue()
  end
  if contract
    return *(As...)
  end
  return nothing
end

expr_to_string(s::Symbol) = String(s)
expr_to_string(ex::Expr) = String(repr(ex))[3:end-1]

"""
    @visualize

Visualize a contraction of ITensors, returning the result of the contraction.

The contraction should be written in terms of a series of ITensors contracted with `*`.

# Examples
```julia
using ITensors
using ITensorsVisualization

i = Index(2, "index_i")
j = Index(10, "index_j")
k = Index(40, "index_k")
l = Index(40, "index_l")
m = Index(40, "index_m")
A = randomITensor(i, j, k)
B = randomITensor(i, j, l, m)
C = randomITensor(k, l)

# Contract the tensors over the common indices
# and visualize the results
ABC = @visualize A * B * C

# Pause to display intermediate results
AB = @visualize A * B pause = true
ABC = @visualize AB * C labels = ["A*B", "C"]
```

# Keyword arguments:
- `pause::Bool = false`: pause after displaying the tensor network visualization.
- `showtags::Bool = true`: show the Index tags on the edge labels of the network.
- `showplevs::Bool = true`: show the Index prime levels on the edge labels of the network.
- `showids::Bool = true`: show a shortened version of the Index id numbers in the edge labels of the network.
- `showdims::Bool = true`: show the Index dimensions on the edge labels of the network.
- `showqns::Bool = false`: show the quantum number sectors. Only available for ITensors with QNs.
- `showarrows::Bool = all(hasqns, tensors)`: show the arrow directions on the edges of the network which correspond to the Index directions for QN indices (corresponding to contravariant and covariant spaces). Only well defined for ITensors with QNs.
- `labels::Vector{String}`: custom tensor labels to display on the nodes of the digram. If not specified, they are determined automatically from the input to the macro.
"""
macro visualize(ex, kwargs...)
  # Must be a tensor contraction
  @assert ex.args[1] == :*
  expr_kwargs = [esc(a) for a in kwargs]
  res = if any(arg -> arg.args[1].args[1] == :labels, expr_kwargs)
    # If labels was passed, don't automatically generate the labels
    :(contract_visualize($(esc.(ex.args[2:end])...); $(expr_kwargs...)))
  else
    # labels was not passed, automatically generate them from the macro input
    :(contract_visualize($(esc.(ex.args[2:end])...); $(expr_kwargs...), labels = expr_to_string.($(ex.args[2:end]))))
  end
  return res
end

max_stars_repo_name: mtfishman/ITensorsVisualization

max_issues_repo_name: mtfishman/ITensorsVisualization

max_forks_repo_name: mtfishman/ITensorsVisualization


text:
"""
`function find_sound{T}(song::AbstractArray{T, 1}, cutoff::Real, maxsilence::Real; mindur::Real=0, maxdur::Real=-1)`

Automatically detect bird calls in a recording and return vectors of their start and end times

Paramaters:
 - `song`: An array containing the samples of the recording
 - `cutoff`: The cutoff (as a proportion of maximum volume) that is used to seperate sound from silence
 - `maxsilence:` The maximum length (in seconds) of silence that can occur before the current song is declared to have ended.
 - `mindur:` The minimum duration (in seconds) of sound with amplitude above `cutoff` that is deemed meaningful.
 - `maxdur`: The maximum duration (in seconds) of sound with amplitude above `cutoff` that is desired.
Sounds that are found but have duration less than `mindur` or greater than `maxdur` are not reported.

Algorithmic details:
1) Compute the STFT (spectrogram) of `song`.
2) For each piece of `song` in the spectrogram, find the centroid frequency
3) We now have a path through the spectrogram that will follow a bird song *if* it is sufficiently loud relative to background noise
4) Compute the amplitude along this path (actually along a small band centered on this path to compensate for windowing effects) We now have an estimate of amplitude vs. time that is insensitive to noise
5) The domain of amplitude(t) may be thought of as consisting of two disjoint subsets: "silence", all t such that amplitude(t) < `cutoff` and "sounds", all t such that amplitude(t) >= `cutoff`
6) However, we decide that silences of duration shorter than `maxsilence` don't count; when one of these silences occurs between two sounds, we just declare the two sounds to be one.
7) Repeat this process until all silences of duration < `maxsilence` have been eliminated.
8) Return two arrays containing start and end times of all the sounds

**Note**: The full spectrogram is never actually computed, because to do so  would be an egregious waste of memory. Instead, the centroid frequency and
amplitudes are computed for each transformed windowed piece of `song` in turn.
"""
function find_sound{T}(song::AbstractArray{T, 1}, cutoff::Real, maxsilence::Real; mindur::Real=0, maxdur::Real=-1)::Tuple{Array{Float64, 1}, Array{Float64, 1}, Array{T, 1}}
    if maxdur == -1
        maxdur = length(song)/samplerate
    end

    song_env = envelope_along_centroid(song, STFT.gaussian(windowlength), windowstep, lowharmonic, highharmonic, STFT.powSpec, envelope_bandwidth)
    song_env /= maximum(song_env)

    maxsilence *= samplerate/windowstep #Convert from seconds to indecies in song_env

    inasong = false
    current_silence_length = 0
    call_starts = Array{Float64, 1}(0)
    call_ends = Array{Float64, 1}(0)
    for (i, amp) in enumerate(song_env)
        if !inasong
            if amp >= cutoff
                current_silence_length = 0
                inasong = true
                push!(call_starts, i)
            end
        else
            if amp < cutoff
                if current_silence_length <= maxsilence
                    current_silence_length += 1
                    continue
                else
                    inasong = false
                    push!(call_ends, i-maxsilence)
                    current_silence_length = 0
                end
            else
                current_silence_length = 0
            end
        end
    end
    if length(call_starts) > length(call_ends)
        push!(call_ends, length(song))
    end

    call_starts *= windowstep
    call_starts += Int(floor(windowlength/2))
    call_starts /= samplerate

    call_ends *= windowstep
    call_ends += Int(floor(windowlength/2))
    call_ends /= samplerate

    indshift = 0
    for i in eachindex(call_starts)
        if call_ends[i-indshift] - call_starts[i-indshift] < mindur || call_ends[i-indshift]-call_starts[i-indshift] > maxdur
            deleteat!(call_starts, i-indshift)
            deleteat!(call_ends, i-indshift)
            indshift += 1
        end
    end

    return call_starts, call_ends, song_env
end

"""
`function find_sound(songname::AbstractString, cutoff::Real, maxsilence::Real; mindur::Real=0, maxdur::Real=-1, make_sonograms::Bool=false)`

find_sound method that takes a filename and writes out sonograms of each sound found
Parameters:
 - `make_sonograms`: Determines whether to create and save sonograms of what was found.
"""
function find_sound(songname::AbstractString, cutoff::Real, maxsilence::Real; mindur::Real=0, maxdur::Real=-1, make_sonograms::Bool=false)
    song = wavread(songname)[1][:, 1]
    call_starts, call_ends, song_env = find_sound(song, cutoff, maxsilence, mindur=mindur, maxdur=maxdur)

    if make_sonograms
        for i in eachindex(call_starts)
	    sonog = sonogram(view(song, Int(round(call_starts[i]*samplerate)):Int(round(call_ends[i]*samplerate)))).^.5
            gr = mglGraph()
            SetRanges(gr, 0, size(sonog, 2), 0, size(sonog, 1))
            SetRange(gr, 'c', minimum(sonog), maximum(sonog))
            SetDefScheme(gr, "wrcbk")
            Dens(gr, sonog)
            #Line(gr, MathGL.mglPoint(call_starts[i], 0, 0), MathGL.mglPoint(call_starts[i], maximum(song_env), 0), "r")
            #Line(gr, MathGL.mglPoint(call_ends[i], 0, 0), MathGL.mglPoint(call_ends[i], maximum(song_env), 0), "r")
            SetRange(gr, 'x', 0, size(sonog, 2)*samplerate/windowstep)
            Axis(gr)
            MathGL.Box(gr)
            WritePNG(gr, "$songname-$i.png")
	end
    end

    return call_starts, call_ends, song_env
end

export find_sound


text:
using Unitful

export ploss_coeff, ploss_coeff_nice, isunitless

function nounit(value)
    try
        uconvert(Unitful.NoUnits, value)
        true
    catch
        false
    end
end

# Source for notation and formulas:
# "Evaluating Gas Network Capacities", T. Koch et al., SIAM MO21, 2015
#
# Some constants are taken from that book, or the related source code of
# Lamatto++, see: http://www.mso.math.fau.de/edom/projects/lamatto.html

# pipe roughness, from REKO data (4390 * 0.012mm, 926 * 0.006mm, 166 * 0.002mm)
const k = 0.012u"mm"

# pipe (mean) diameter for friction factor
const D_m = 1.0u"m"

"friction factor by Nikoradse, for turbulent flow"
function nikoradse(diameter, roughness)
    (log10(diameter/roughness) + 1.138)^(-2)
end

# friction factor
const λ = nikoradse(D_m, k)

# universal gas constant (https://en.wikipedia.org/wiki/Gas_constant)
const R = 8.3144598u"J/K/mol"

# molar mass of methane
const molmass = 16.04u"g/mol"

# specific gas constant
const R_s = R/molmass

# critical pressure (default in Lamatto++)
const p_c = 46.4512u"bar"

# critical temperature (default in Lamatto++)
const T_c = 192.033u"K"

# mean pressure
const p_m = 60u"bar"

# mean temperature (default in Lamatto++)
const T_m = 283.15u"K"

"compressibility factor by AGA8"
function aga8(pressure, temperature)
    # "reduced" pressure and temperature
    p_r = pressure / p_c
    T_r = temperature / T_c
    z = 1.0 + 0.257 * p_r - 0.533 * p_r/T_r
end

# z_m: mean compressibility factor
const z_m = aga8(p_m, T_m)

# Pressure Loss equation ("Weymouth"), when we care about length and diameter
#   p_o^2 - p_i^2 = L/D^5 C q|q|
const weymouth = (π/4)^(-2) * λ * R_s * z_m * T_m
@assert !nounit(weymouth)

# alternative coefficient with same order of magnitude for nice, round numbers
const weymouth_nice = 3000u"m^2/s^2"
@assert nounit(weymouth/weymouth_nice) # same unit

# Unit-less Weymouth coefficient: Add types for length, diameters and flow, then
# divide by the expected type to get just a float
const ploss_coeff_unit = 1u"(km*m^(-5)*(kg/s)^2) / (bar^2)"
const ploss_coeff = weymouth * ploss_coeff_unit
const ploss_coeff_nice = weymouth_nice * ploss_coeff_unit
@assert nounit(ploss_coeff)
@assert nounit(ploss_coeff_nice)


text:
struct Term{T1 <: Real, T2 <: ProximableFunction, T3 <: AbstractExpression}
	lambda::T1
	f::T2
	A::T3
	Term(lambda::T1, f::T2, ex::T3) where {T1,T2,T3} = new{T1,T2,T3}(lambda,f,ex)
end

function Term(f::T, ex::AbstractExpression) where {T<:ProximableFunction}
	A = convert(Expression,ex)
	Term(1,f, A)
end

# Operations

# Define sum of terms simply as their vcat

import Base: +

(+)(a::Term,b::Term) = (a,b)
(+)(a::NTuple{N,Term},b::Term)    where {N} = (a...,b)
(+)(a::Term,b::NTuple{N,Term})    where {N} = (a,b...)
(+)(a::NTuple{N,Term},b::Tuple{}) where {N} = a
(+)(a::Tuple{},b::NTuple{N,Term}) where {N} = b
(+)(a::NTuple{N,Term},b::NTuple{M,Term}) where {N,M} = (a...,b...)

# Define multiplication by constant

import Base: *

function (*)(a::T1, t::Term{T,T2,T3}) where {T1<:Real, T, T2, T3}
	coeff = *(promote(a,t.lambda)...)
	Term(coeff, t.f, t.A)
end

function (*)(a::T1, t::T2) where {T1<:Real, N, T2 <: Tuple{Vararg{<:Term,N}} }
	return a.*t
end

# Properties

variables(t::Term) = variables(t.A)
operator(t::Term) = operator(t.A)
affine(t::Term) = affine(t.A)
displacement(t::Term) = displacement(t.A)

#importing properties from ProximalOperators
import ProximalOperators:
			  is_affine,
			  is_cone,
			  is_convex,
			  is_generalized_quadratic,
			  is_prox_accurate,
			  is_quadratic,
			  is_separable,
			  is_set,
			  is_singleton,
			  is_smooth,
			  is_strongly_convex

#importing properties from AbstractOperators
is_f = [:is_linear,
        :is_eye,
        :is_null,
        :is_diagonal,
        :is_AcA_diagonal,
        :is_AAc_diagonal,
        :is_orthogonal,
        :is_invertible,
        :is_full_row_rank,
        :is_full_column_rank,
        :is_sliced
       ]

for f in is_f
	@eval begin
		import AbstractOperators: $f
		$f(t::Term) = $f(operator(t))
		$f(t::NTuple{N,Term}) where {N} = all($f.(t))
	end
end

is_smooth(t::Term) = is_smooth(t.f)
is_convex(t::Term)    = is_convex(t.f) && is_linear(t)
is_quadratic(t::Term) = is_quadratic(t.f) && is_linear(t)
is_strongly_convex(t::Term) = is_strongly_convex(t.f) && is_full_column_rank(operator(t.A))

include("proximalOperators_bind.jl")

# other stuff, to make Term work with iterators
import Base: iterate, isempty
iterate(t::Term, state = true) = state ? (t, false) : nothing
isempty(t::Term) =  false


text:
a=(1,2,3.0, "abc") # 用圆括号构造元组
typeof(a)          # Tuple{Int64, Int64, Float64, String}
                   # 不同类型可以共处
a[2]               # 第二个元素是 2
a[2]=20            # 报错，元组内容不可变

b=[1,2,3.0,"abc"]  # 用方括号构造数组
typeof(b)          # Vector{Any} (alias for Array{Any, 1})
                   # 各类型被提升为共同的父类型
                   # 在 Julia 中，所有类型共同的父类型是 Any
b[2]="bcd"         # 可以改变数组的内容
b                  # [1, "bcd", 3.0, "abc"]

max_stars_repo_name: zhaiyusci/Julia-Notes

max_issues_repo_name: zhaiyusci/Julia-Notes

max_forks_repo_name: zhaiyusci/Julia-Notes


text:
# ---
# title: 493. Reverse Pairs
# id: problem493
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Binary Search, Divide and Conquer, Sort, Binary Indexed Tree, Segment Tree
# link: <https://leetcode.com/problems/reverse-pairs/description/>
# hidden: true
# ---
#
# Given an array `nums`, we call `(i, j)` an **_important reverse pair_** if `i
# < j` and `nums[i] > 2*nums[j]`.
#
# You need to return the number of important reverse pairs in the given array.
#
# **Example1:**
#
#
#
#     **Input** : [1,3,2,3,1]
#     **Output** : 2
#
#
# **Example2:**
#
#
#
#     **Input** : [2,4,3,5,1]
#     **Output** : 3
#
#
# **Note:**
#
#   1. The length of the given array will not exceed `50,000`.
#   2. All the numbers in the input array are in the range of 32-bit integer.
#
#
## @lc code=start
using LeetCode

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


text:
s = open("input.txt") do file
    readlines(file);
end

function FollowPath(instructions)
    pos = [0,0];
    dir = 1;
    directions = [
        [1,0], [0,-1], [-1, 0], [0, 1]
    ];
    for command in instructions
        action = command[1];
        num = parse(Int, command[2:end]);

        if action == 'F'
            pos .+= directions[dir].*num;
        elseif action == 'N'
            pos .+= directions[4]*num;
        elseif action == 'S'
            pos .+= directions[2]*num;
        elseif action == 'E'
            pos .+= directions[1]*num;
        elseif action == 'W'
            pos .+= directions[3]*num;
        elseif action == 'R'
            dir += Int(num/90);
            while dir > 4
                dir -= 4;
            end
        elseif action == 'L'
            dir -= Int(num/90);
            while dir < 1
                dir += 4;
            end
        end
    end

    return sum(abs.(pos));
end

FollowPath(s)


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using IterativeLQR
using LinearAlgebra
using Plots

# ## horizon
T = 51

# ## car
num_state = 3
num_action = 2
num_parameter = 0

function car(x, u, w)
    [u[1] * cos(x[3]); u[1] * sin(x[3]); u[2]]
end

function midpoint_explicit(x, u, w)
    h = 0.1 # timestep
    x + h * car(x + 0.5 * h * car(x, u, w), u, w)
end

# ## model
dynamics = Dynamics(midpoint_explicit, num_state, num_action, num_parameter)
model = [dynamics for t = 1:T-1]

# ## initialization
x1 = [0.0; 0.0; 0.0]
xT = [1.0; 1.0; 0.0]

# ## rollout
ū = [1.0e-2 * [1.0; 0.1] for t = 1:T-1]
w = [zeros(num_parameter) for t = 1:T]
x̄ = rollout(model, x1, ū, w)

# ## objective
ot = (x, u, w) -> 1.0 * dot(x - xT, x - xT) + 1.0e-2 * dot(u, u)
oT = (x, u, w) -> 1000.0 * dot(x - xT, x - xT)
ct = Cost(ot, num_state, num_action, num_parameter)
cT = Cost(oT, num_state, 0, num_parameter)
objective = [[ct for t = 1:T-1]..., cT]

# ## constraints
ul = -5.0 * ones(num_action)
uu = 5.0 * ones(num_action)

p_obs = [0.5; 0.5]
r_obs = 0.1

function stage_con(x, u, w)
    e = x[1:2] - p_obs
    [
    ul - u; # control limit (lower)
    u - uu; # control limit (upper)
    r_obs^2.0 - dot(e, e); # obstacle
    ]
end

function terminal_con(x, u, w)
    e = x[1:2] - p_obs
    [
    x - xT; # goal
    r_obs^2.0 - dot(e, e); # obstacle
    ]
end

cont = Constraint(stage_con, num_state, num_action, indices_inequality=collect(1:5))
conT = Constraint(terminal_con, num_state, num_action, indices_inequality=collect(3 .+ (1:1)))
constraints = [[cont for t = 1:T-1]..., conT]

# ## problem
prob = Solver(model, objective, constraints)
initialize_controls!(prob, ū)
initialize_states!(prob, x̄)

# ## solve
solve!(prob)

# ## solution
x_sol, u_sol = get_trajectory(prob)

# ## visualize
plot(hcat(x_sol...)')
plot(hcat(u_sol[1:end-1]...)', linetype=:steppost)

# ## benchmark allocations + timing
using BenchmarkTools
info = @benchmark solve!($prob, x̄, ū) setup=(x̄=deepcopy(x̄), ū=deepcopy(ū))


text:
# XXX this is not complete
function approx_path_leifer(op::OrienteeringProblem)
    m = Model(solver=Gurobi.GurobiSolver(OutputFlag=output))
    N = length(op)

    V = [1:N]
    without_start = [1:op.start-1; op.start+1:N]
    without_stop = [1:op.stop-1; op.stop+1:N]
    Vprime = intersect(without_start, without_stop)

    @defVar(m, x[1:N,1:N], Bin) # (7)
    @defVar(m, y[1:N], Bin) # (8)

    @setObjective(m, Max, sum{ op.r[j]*y[j] j=1:N }) # (1)

    @addConstraint(m, sum{x[op.start,j], j=without_start} == 1) # (2)
    @addConstraint(m, sum{x[i,op.stop], i=without_stop} == 1) # (3)

    @addConstraint(m, conservation[j=Vprime], sum{x[i,j],i=1:N} + sum{x[j,k],k=1:N} == 2*y[j]) # (4)

    @addConstraint(m, sum{ sum{op.distances[i,j]*x[i,j], j=1:N}, i=1:N } <= op.distance_limit) # (5)

    status = solve(m)

    # if status != :Optimal
    #     @bp
    # end
    if status != :Optimal
        warn("Not solved to optimality:\n$op")
    end

    # @show getObjectiveValue(m)
    # @show getValue(x)
    # @show getValue(u)

    path = [op.start]
    current = -1
    dist_sum = 0.0
    xopt = round(Int,getValue(x))
    while current != op.stop
        if current==-1
            current=op.start
        end
        current = findfirst(xopt[current,:])
        push!(path, current)
    end
    # @assert(distance(op, path) <= op.distance_limit)
    if distance(op,path) > op.distance_limit
        warn("Path Length: $(distance(op,path)), Limit: $(op.distance_limit)")
    end
    return path
end


text:
# --------------------------------------------------------------------------
# ACE.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

import JuLIP: evaluate!, evaluate_d!, SVec, write_dict, read_dict
import ACE: alloc_B, alloc_dB

# ------------------------------------------------------------
#   Fourier Basis evaluation

struct FourierBasis{T}
   deg::Int
   _fltt::Type{T}
   FourierBasis(deg::Integer, _fltt::Type{T}) where {T <: AbstractFloat} =
      new{T}(Int(deg), _fltt)
end

cyl_l2i(l, b::FourierBasis) = cyl_l2i(l, b.deg)
cyl_i2l(i, b::FourierBasis) = cyl_l2i(i, b.deg)

FourierBasis(deg::Integer) = FourierBasis(deg, Float64)

fltype(::FourierBasis{T}) where {T} = T
Base.length(fB::FourierBasis) = 2 * fB.deg + 1

write_dict(b::FourierBasis) =
      Dict("__id__" => "ACE_FourierBasis",
              "deg" => b.deg)
read_dict(::Val{:ACE_FourierBasis}, D::Dict) =
      FourierBasis(D["deg"], Float64)

# Dict(fB::FourierBasis) = Dict(
#       "__id__" => "ACE_FourierBasis",
#       "deg" => deg,
#       "fltt" => "$(fG._fltt)"
#    )
#
# FourierBasis(D::Dict) = FourierBasis(D["deg"], D["fltt"])
# FourierBasis(deg::Integer, fltt::AbstractString) =
#       FourierBasis(deg, eval(Meta.parse(fltt)))
#
# convert(::Val{:ACE_FourierBasis}, D::Dict) = FourierBasis(D)

alloc_B(fB::FourierBasis{T}) where {T}  = zeros(Complex{T}, length(fB))
alloc_dB(fB::FourierBasis{T}, args...) where {T} = zeros(SVec{2,Complex{T}}, length(fB))

# specify ordering
# 0, 1, -1, 2, -2, 3, -3
# cyl_l2i(l, maxL) = maxL + 1 + l  # = i
# cyl_i2l(i, maxL) = i - maxL - 1

cyl_l2i(l, args...) = 2*abs(l) + (sign(l)<=0)

function cyl_i2l(i, args...)
   absk = i ÷ 2
   r = mod(i, 2)
   sigk = 1-2*(r==1)
   return sigk * absk
end

function evaluate!(P, _::Nothing, fB::FourierBasis,
                   c::CylindricalCoordinates{T}) where {T}
   @assert length(P) >= length(fB)
   maxL = fB.deg
   z = c.cosθ + im * c.sinθ
   zl = one(T) + im * zero(T)
   P[cyl_l2i(0, maxL)] = zl
   for l = 1:maxL
      zl *= z
      P[cyl_l2i( l, maxL)] = zl
      P[cyl_l2i(-l, maxL)] = conj(zl)
   end
   return P
end

# we only return ∂Pl/∂x̂ since
#     ∂Pl/∂ŷ = ∂Pl/∂z * ∂z / ∂ŷ
#            = im * ∂Pl/∂z
#            = im * ∂Pl/∂x̂
function evaluate_d!(P, dP, _::Nothing, fB::FourierBasis,
                     c::CylindricalCoordinates{T}) where {T}
   @assert length(P) >= length(fB)
   maxL = fB.deg
   z = c.cosθ + im * c.sinθ    # z = x̂ + i ŷ
   zl = one(T) + im * zero(T)  # zl = z^l -> initialise to z^0 = 1
   P[cyl_l2i(0, maxL)] = zl
   dP[cyl_l2i(0, maxL)] =  SVec(zero(T), zero(T))

   for l = 1:maxL
      # ∂P_{ l}/∂x̂ = l z^(l-1)
      # ∂P_{-l}/∂x̂ = l z̄^(l-1)
      dP[cyl_l2i( l, maxL)] = SVec(l * zl,  im * l * zl)
      dP[cyl_l2i(-l, maxL)] = SVec(l * conj(zl), -im * l * conj(zl))

      zl *= z  # zl = z^l
      P[cyl_l2i( l, maxL)] = zl
      P[cyl_l2i(-l, maxL)] = conj(zl)

      # dP[cyl_l2i( l, maxL)] =  im * l * zl
      # dP[cyl_l2i(-l, maxL)] = -im * l * conj(zl)
   end
   return P
end


text:
export MagPhase, zerophase
export SignalElement

import Base: Complex, promote_rule, abs, angle

struct MagPhase{T<:Real} <: Number
    mag::T
    phi::T
    function MagPhase{T}(mag::T, phi::T) where {T <: Real}
        mag < 0 && throw(ArgumentError("negative magnitude!"))
        new(mag, phi)
    end
end

const SignalElement{T} = Union{Complex{T}, MagPhase{T}} where T <: Real

MagPhase(x::T, y::T) where {T <: Real} = MagPhase{T}(x, y)
MagPhase(x::Real, y::Real) = MagPhase(promote(x,y)...)
MagPhase(x::Real) = MagPhase(x, zero(x))
MagPhase(z::Complex) = MagPhase(abs(z), angle(z))
MagPhase{T}(z::C) where {T, C <: Complex{<:T}} = MagPhase{T}(abs(z), angle(z))
MagPhase{T}(x) where T = MagPhase(convert(T, x))

Base.Complex{T}(z::MagPhase{T}) where T = Complex{T}(z.mag * cos(z.phi), z.mag * sin(z.phi))
Base.Complex(z::MagPhase{T}) where T = Complex{T}(z)

Base.abs(z::MagPhase) = z.mag
Base.angle(z::MagPhase) = z.phi

zerophase(z::MagPhase) = MagPhase(abs(z))
zerophase(z::Complex) = Complex(abs(z))

Base.promote_rule(::Type{Complex{T}}, ::Type{MagPhase{T}}) where T = Complex{T}

function Base.show(io::IO, z::MagPhase)
    print(io, "(")
    show(io, z.mag)
    print(io, ", ")
    show(io, z.phi)
    print(io, ")")
end


text:
import ApproxFun
import ApproxFun: Fun, ⊗, (..)

## 2nd Poincaré Invariant for general two-form

struct PoincareInvariant2ndApproxFun{DT,ND,NC,NV,ET,ΩT,ϑT} <: AbstractPoincareInvariant2ndApproxFun{DT,ND,NC,NV}
    equ::ET
    ω::ΩT
    D²ϑ::ϑT
    Δt::DT
    nx::Int
    ny::Int
    ntime::Int
    nsave::Int
    nt::Int
    I::OffsetArray{Double64,1,Vector{Double64}}
    J::OffsetArray{Double64,1,Vector{Double64}}
    K::OffsetArray{Double64,1,Vector{Double64}}
    L::OffsetArray{Double64,1,Vector{Double64}}
    ΔI::OffsetArray{Double64,1,Vector{Double64}}
    ΔJ::OffsetArray{Double64,1,Vector{Double64}}
end

function PoincareInvariant2ndApproxFun(f_equ::Function, f_surface::Function, ω::ΩT, D²ϑ::ϑT, Δt::TT, nd::Int, nx::Int, ny::Int, ntime::Int, nsave::Int=1, DT=Float64) where {TT,ΩT,ϑT}

    # compute Chebyshev points
    c = ApproxFun.points(ApproxFun.Chebyshev(0..1)^2, nx*ny)

    # compute initial conditions
    q₀ = [f_surface(c[i][1], c[i][2]) for i in eachindex(c)]

    # initialise euation
    equ = f_equ(q₀)

    # create arrays for results
    nt = div(ntime, nsave)

    I  = OffsetArray(zeros(Double64, nt+1), 0:nt)
    J  = OffsetArray(zeros(Double64, nt+1), 0:nt)
    K  = OffsetArray(zeros(Double64, nt+1), 0:nt)
    L  = OffsetArray(zeros(Double64, nt+1), 0:nt)
    ΔI = OffsetArray(zeros(Double64, nt+1), 0:nt)
    ΔJ = OffsetArray(zeros(Double64, nt+1), 0:nt)

    # get size of coefficient and value vectors
    SC = ApproxFun.Chebyshev(0..1)^2
    SU = ApproxFun.Ultraspherical(1, 0..1)^2
    Dx = ApproxFun.Derivative(SC, [1,0])
    Dy = ApproxFun.Derivative(SC, [0,1])

    fqc = Fun(SC, ApproxFun.transform(SC, hcat(q₀...)[1,:]))
    fqu = Fun(Dx * fqc, SU)
    nc = ApproxFun.ncoefficients(fqu)
    nv = length(ApproxFun.values(fqu))

    if get_config(:verbosity) > 1
        println()
        println("Second Euler-Poincaré Integral Invariant (ApproxFun)")
        println("====================================================")
        println()
        println(" nx    = ", nx)
        println(" ny    = ", ny)
        println(" np    = ", length(c))
        println(" nc(CC)= ", ApproxFun.ncoefficients(fqc))
        println(" nv(CC)= ", length(ApproxFun.values(fqc)))
        println(" nc(US)= ", nc)
        println(" nv(US)= ", nv)
        println(" ntime = ", ntime)
        println(" nsave = ", nsave)
        println(" Δt    = ", Δt)
        println()
    end

    # initialise Poincare invariant
    PoincareInvariant2ndApproxFun{DT,nd,nc,nv,typeof(equ),ΩT,ϑT}(equ, ω, D²ϑ, DT(Δt), nx, ny, ntime, nsave, nt, I, J, K, L, ΔI, ΔJ)
end

function PoincareInvariant2ndApproxFun(f_equ::Function, f_surface::Function, ω::ΩT, Δt::TT, nd::Int, nx::Int, ny::Int, ntime::Int, nsave::Int=1, DT=Float64) where {TT,ΩT}
    PoincareInvariant2ndApproxFun(f_equ, f_surface, ω, (), Δt, nd, nx, ny, ntime, nsave, DT)
end

function evaluate_poincare_invariant(pinv::PoincareInvariant2ndApproxFun{DT}, sol::Solution) where {DT}

    local verbosity = get_config(:verbosity)

    # local q = permutedims(sol.q, [3,1,2])
    # if isdefined(sol, :p)
    #     local p = permutedims(sol.p, [3,1,2])
    # end
    # if isdefined(sol, :λ)
    #     local λ = permutedims(sol.λ, [3,1,2])
    # end

    verbosity ≤ 1 ? prog = Progress(size(sol.q,1), 5) : nothing

    for i in axes(sol.q,1)
        verbosity > 1 ? println("      it = ", i) : nothing
        pinv.I[i] = compute_noncanonical_invariant(pinv, sol.t[i], hcat(sol.q[i,:]...))
        pinv.ΔI[i] = abs(pinv.I[0]) < sqrt(eps()) ? pinv.I[i] : (pinv.I[i] .- pinv.I[0]) ./ pinv.I[0]
        verbosity > 1 ? println("           I_q = ", pinv.I[i], ",   ε_q = ", pinv.ΔI[i]) : nothing
    end

    if isdefined(sol, :p)
        for i in axes(sol.q,1)
            pinv.J[i] = compute_canonical_invariant(pinv, hcat(sol.q[i,:]...), hcat(sol.p[i,:]...))
            pinv.ΔJ[i] = abs(pinv.J[0]) < sqrt(eps()) ? pinv.J[i] : (pinv.J[i] .- pinv.J[0]) ./ pinv.J[0]
            verbosity > 1 ? println("           I_p = ", pinv.J[i], ",   ε_p = ", pinv.ΔJ[i]) : nothing
        end
    end

    return (DT.(pinv.I), DT.(pinv.J), DT.(pinv.ΔI), DT.(pinv.ΔJ))
end

function evaluate_poincare_invariant_correction(pinv::PoincareInvariant2ndApproxFun{DT}, sol::Solution) where {DT}

    local verbosity = get_config(:verbosity)

    if isdefined(sol, :λ)
        for i in axes(sol.q,1)
            pinv.K[i] = compute_noncanonical_correction(pinv, sol.t[i], hcat(sol.q[i,:]...), hcat(sol.λ[i,:]...))
            pinv.L[i] = pinv.I[i] - pinv.Δt^2 * pinv.K[i]
            verbosity > 1 ? println("           I_λ = ", pinv.L[i], ",   ε_λ = ", (pinv.L[i]-pinv.L[0])/pinv.L[0]) : nothing
            verbosity > 1 ? println("           K_λ = ", pinv.K[i]) : nothing
        end

        verbosity ≤ 1 ? next!(prog) : nothing
    end

    return (DT.(pinv.K), DT.(pinv.L))
end

function write_to_hdf5(pinv::PoincareInvariant2ndApproxFun{DT}, sol::Solution, output_file::String) where {DT}
    # h5open(output_file, isfile(output_file) ? "r+" : "w") do h5
    h5open(output_file, "w") do h5

        write(h5, "t", sol.t)
        write(h5, "I", pinv.I)

        isdefined(sol, :p) ? write(h5, "J", pinv.J) : nothing
        isdefined(sol, :λ) ? write(h5, "K", pinv.L) : nothing
        isdefined(sol, :λ) ? write(h5, "L", pinv.L) : nothing
    end
end

## 2nd Poincaré Invariant for canonical two-form

struct PoincareInvariant2ndApproxFunCanonical{DT,ND,NC,NV,ET} <: AbstractPoincareInvariant2ndApproxFun{DT,ND,NC,NV}
    equ::ET
    Δt::DT
    nx::Int
    ny::Int
    ntime::Int
    nsave::Int
    nt::Int
    I::OffsetArray{Double64,1,Vector{Double64}}
    ΔI::OffsetArray{Double64,1,Vector{Double64}}
end

function PoincareInvariant2ndApproxFunCanonical(f_equ::Function, f_surface_q::Function, f_surface_p::Function, Δt::TT, nd::Int, nx::Int, ny::Int, ntime::Int, nsave::Int=1, DT=Float64) where {TT}

    if get_config(:verbosity) > 1
        println()
        println("Second Canonical Euler-Poincaré Integral Invariant (ApproxFun)")
        println("==============================================================")
        println()
        println(" nx    = ", nx)
        println(" ny    = ", ny)
        println(" ntime = ", ntime)
        println(" nsave = ", nsave)
        println(" Δt    = ", Δt)
        println()
    end

    # compute Chebyshev points
    c = ApproxFun.points(ApproxFun.Chebyshev(0..1)^2, nx*ny)

    # compute initial conditions
    q₀ = [f_surface_q(c[i][1], c[i][2]) for i in eachindex(c)]
    p₀ = [f_surface_p(c[i][1], c[i][2]) for i in eachindex(c)]

    # initialise euation
    equ = f_equ(q₀, p₀)

    # create arrays for results
    nt = div(ntime, nsave)

    I  = OffsetArray(zeros(Double64, nt+1), 0:nt)
    ΔI = OffsetArray(zeros(Double64, nt+1), 0:nt)

    # get size of coefficient and value vectors
    SC = ApproxFun.Chebyshev(0..1)^2
    SU = ApproxFun.Ultraspherical(1, 0..1)^2
    Dx = ApproxFun.Derivative(SC, [1,0])
    Dy = ApproxFun.Derivative(SC, [0,1])

    fq = Fun(Dx * Fun(SC, ApproxFun.transform(SC, hcat(q₀...)[1,:])), SU)
    nc = ApproxFun.ncoefficients(fq)
    nv = length(ApproxFun.values(fq))

    # initialise Poincare invariant
    PoincareInvariant2ndApproxFunCanonical{DT,nd,nc,nv,typeof(equ)}(equ, DT(Δt), nx, ny, ntime, nsave, nt, I, ΔI)
end

function evaluate_poincare_invariant(pinv::PoincareInvariant2ndApproxFunCanonical{DT}, sol::Solution) where {DT}
    local verbosity = get_config(:verbosity)

    verbosity ≤ 1 ? prog = Progress(size(sol.q,1), 5) : nothing

    for i in axes(sol.q,1)
        verbosity > 1 ? println("      it = ", i) : nothing
        pinv.I[i] = compute_canonical_invariant(pinv, hcat(sol.q[i,:]...), hcat(sol.p[i,:]...))
        verbosity > 1 ? println("           I_p = ", pinv.I[i], ",   ε_p = ", (pinv.I[i]-pinv.I[0])/pinv.I[0]) : nothing

        verbosity ≤ 1 ? next!(prog) : nothing
    end

    return pinv.I
end

function write_to_hdf5(pinv::PoincareInvariant2ndApproxFunCanonical{DT}, sol::Solution, output_file::String) where {DT}
    # h5open(output_file, isfile(output_file) ? "r+" : "w") do h5
    h5open(output_file, "w") do h5
        write(h5, "t", sol.t)
        write(h5, "I", pinv.I)
    end
end

## Common functionality

@generated function compute_integral(Is)
    SU = ApproxFun.Ultraspherical(1, 0..1)^2
    Q  = ApproxFun.DefiniteIntegral(SU)

    quote
        return Number( $Q * Fun($SU, ApproxFun.transform($SU, Is)) )
    end
end

@generated function compute_derivative1(q)
    Dx = ApproxFun.Derivative(ApproxFun.Chebyshev(0..1)^2, [1,0])
    CU = ApproxFun.Conversion(ApproxFun.Ultraspherical(1, 0..1) ⊗ ApproxFun.Chebyshev(0..1), ApproxFun.Ultraspherical(1, 0..1)^2)

    quote
        return values($CU * ($Dx * q))
    end
end

@generated function compute_derivative2(q)
    Dy = ApproxFun.Derivative(ApproxFun.Chebyshev(0..1)^2, [0,1])
    CU = ApproxFun.Conversion(ApproxFun.Chebyshev(0..1) ⊗ ApproxFun.Ultraspherical(1, 0..1), ApproxFun.Ultraspherical(1, 0..1)^2)

    quote
        return values($CU * ($Dy * q))
    end
end

@generated function compute_canonical_invariant(pinv::AbstractPoincareInvariant2ndApproxFun{DT,ND,NC,NV}, q::AbstractArray{DT,2}, p::AbstractArray{DT,2}) where {DT,ND,NC,NV}
    SC = ApproxFun.Chebyshev(0..1)^2
    SU = ApproxFun.Ultraspherical(1, 0..1)^2

    Js = zeros(DT, NV)
    qσ = zeros(DT, ND, NV)
    qτ = zeros(DT, ND, NV)
    pσ = zeros(DT, ND, NV)
    pτ = zeros(DT, ND, NV)

    quote
        for j in axes(p,1)
            fq = Fun($SC, ApproxFun.transform($SC, q[j,:]))
            fp = Fun($SC, ApproxFun.transform($SC, p[j,:]))

            # compute derivatives of q and store values
            $qσ[j,:] .= compute_derivative1(fq)
            $qτ[j,:] .= compute_derivative2(fq)

            # compute derivatives of p and store values
            $pσ[j,:] .= compute_derivative1(fp)
            $pτ[j,:] .= compute_derivative2(fp)
        end

        # compute integrands of integral invariants at all points
        for k in 1:NV
            $Js[k] = $pσ[:,k] ⋅ $qτ[:,k] - $qσ[:,k] ⋅ $pτ[:,k]
        end

        # compute canonical integral invariant
        return compute_integral($Js)
    end
end

@generated function compute_noncanonical_invariant(pinv::PoincareInvariant2ndApproxFun{DT,ND,NC,NV}, t::DT, q::AbstractArray{DT,2}) where {DT,ND,NC,NV}
    SC = ApproxFun.Chebyshev(0..1)^2
    SU = ApproxFun.Ultraspherical(1, 0..1)^2
    CU = ApproxFun.Conversion(SC, SU)

    Is = zeros(DT, NV)
    qs = zeros(DT, ND, NV)
    qσ = zeros(DT, ND, NV)
    qτ = zeros(DT, ND, NV)
    Ω  = zeros(DT, ND, ND)

    quote
        # loop over dimensions of state vector
        for j in axes(q,1)
            # obtain Chebyshev coefficients
            fq = Fun($SC, ApproxFun.transform($SC, q[j,:]))

            # obtain values of q at the same points as the derivatives
            $qs[j,:] .= ApproxFun.itransform($SU, resize!(collect(($CU * fq).coefficients), NC))

            # compute derivatives of q and store values
            $qσ[j,:] .= compute_derivative1(fq)
            $qτ[j,:] .= compute_derivative2(fq)
        end

        # compute integrands of integral invariants at all points
        for k in 1:NV
            pinv.ω(t, $qs[:,k], $Ω)
            $Is[k] = vector_matrix_vector_product($qτ[:,k], $Ω, $qσ[:,k])
        end

        # compute noncanonical integral invariant
        return compute_integral($Is)
    end
end

@generated function compute_noncanonical_correction(pinv::PoincareInvariant2ndApproxFun{DT,ND,NC,NV}, t::DT, q::AbstractArray{DT,2}, λ::AbstractArray{DT,2}) where {DT,ND,NC,NV}
    SC = ApproxFun.Chebyshev(0..1)^2
    SU = ApproxFun.Ultraspherical(1, 0..1)^2
    CU = ApproxFun.Conversion(SC, SU)

    Ks = zeros(DT, NV)
    qs = zeros(DT, ND, NV)
    λs = zeros(DT, ND, NV)
    qσ = zeros(DT, ND, NV)
    qτ = zeros(DT, ND, NV)
    λσ = zeros(DT, ND, NV)
    λτ = zeros(DT, ND, NV)
    Ω  = zeros(DT, ND, ND)
    D²ϑ= zeros(DT, ND, ND)

    quote
        for j in axes(q,1)
            fq = Fun($SC, ApproxFun.transform($SC, q[j,:]))
            fλ = Fun($SC, ApproxFun.transform($SC, λ[j,:]))

            # obtain values of q at the same points as the derivatives
            $qs[j,:] .= ApproxFun.itransform($SU, resize!(($CU * fq).coefficients, NC))
            $λs[j,:] .= ApproxFun.itransform($SU, resize!(($CU * fλ).coefficients, NC))

            # compute derivatives of q and store values
            $qσ[j,:] .= compute_derivative1(fq)
            $qτ[j,:] .= compute_derivative2(fq)

            # compute derivatives of λ and store values
            $λσ[j,:] .= compute_derivative1(fλ)
            $λτ[j,:] .= compute_derivative2(fλ)
        end

        # compute integrands of integral invariants at all points
        for k in 1:NV
            pinv.ω(t, $qs[:,k], $Ω)
            $Ks[k] = vector_matrix_vector_product($λτ[:,k], $Ω, $λσ[:,k])

            for j in eachindex(pinv.D²ϑ)
                pinv.D²ϑ[j](t, $qs[:,k], $D²ϑ)
                $Ks[k] += 2 * $λs[j,k] * vector_matrix_vector_product($qτ[:,k], $D²ϑ, $λσ[:,k])
                $Ks[k] -= 2 * $λs[j,k] * vector_matrix_vector_product($qσ[:,k], $D²ϑ, $λτ[:,k])
            end
        end

        # compute correction to noncanonical integral invariant
        return compute_integral($Ks)
    end
end


text:
# GSA code v0.1.
# Coded by Jesús Mejía.
# Based on MATLAB code of Esmat Rashedi, 2010.
# Adopted from
# https://la.mathworks.com/matlcgsaentral/fileexchange/61116-gsa-+-chaotic-gravitational-constant
# " E. Rashedi, H. Nezamabadi-pour and S. Saryazdi,
# “GSA: A Gravitational Search Algorithm”, Information sciences, vol. 179,
# no. 13, pp. 2232-2248, 2009."
# ------------------------------------------
# Mirjalili, Seyedali, and Amir H. Gandomi.
# "Chaotic gravitational constants for the gravitational search algorithm."
# Applied Soft Computing 53 (2017): 407-419.
# ------------------------------------------

include("physics.jl")
include("chaos.jl")

mutable struct CGSA <: AbstractParameters
    N::Int
    chValueInitial::Real
    chaosIndex::Real
    ElitistCheck::Int
    Rpower::Int
    Rnorm::Int
    wMax::Real
    wMin::Real
    X::Matrix{Float64}
    V::Matrix{Float64}
    fitness::Vector{Float64}
end

"""

    CGSA(;
        N::Int    = 30,
        chValueInitial::Real   = 20,
        chaosIndex::Real   = 9,
        ElitistCheck::Int    = 1,
        Rpower::Int    = 1,
        Rnorm::Int    = 2,
        wMax::Real   = chValueInitial,
        wMin::Real   = 1e-10,
        information = Information(),
        options = Options()
    )

CGSA is an extension of the GSA algorithm but with Chaotic gravitational constants for
the gravitational search algorithm.

Ref. Chaotic gravitational constants for the gravitational search algorithm.
Applied Soft Computing 53 (2017): 407-419.

Parameters:

- N: Population size
- chValueInitial: Initial value for the chaos value
- chaosIndex: Integer 1 ≤ chaosIndex ≤ 10 is the function that model the chaos
- Rpower: power related to the distance norm(x)^Rpower
- Rnorm: is the value as in norm(x, Rnorm)

# Example

```jldoctest
julia> f(x) = sum(x.^2)
f (generic function with 1 method)

julia> optimize(f, [-1 -1 -1; 1 1 1.0], CGSA())
+=========== RESULT ==========+
  iteration: 500
    minimum: 8.63808e-08
  minimizer: [0.0002658098418993323, -1.140808975532608e-5, -0.00012488307670533095]
    f calls: 15000
 total time: 0.1556 s
+============================+

julia> optimize(f, [-1 -1 -1; 1 1 1.0], CGSA(N = 80, chaosIndex = 1))
+=========== RESULT ==========+
  iteration: 500
    minimum: 1.0153e-09
  minimizer: [-8.8507563788141e-6, -1.3050111801923072e-5, 2.7688577445980026e-5]
    f calls: 40000
 total time: 1.0323 s
+============================+
```

"""
function CGSA(;
        N::Int    = 30,
        chValueInitial::Real   = 20,
        chaosIndex::Real   = 8,
        ElitistCheck::Int    = 1,
        Rpower::Int    = 1,
        Rnorm::Int    = 2,
        wMax::Real   = chValueInitial,
        wMin::Real   = 1e-10,
        X::Matrix{Float64} = zeros(0,0),
        V::Matrix{Float64} = zeros(0,0),
        fitness::Vector{Float64} = zeros(0),
        information = Information(),
        options = Options()
    )

    parameters = CGSA(N, chValueInitial, chaosIndex, ElitistCheck, Rpower, Rnorm, wMax, wMin, X, V, fitness)

    Algorithm(
        parameters,
        information = information,
        options = options,
    )
end

function initialize!(
        status,
        parameters::CGSA,
        problem::AbstractProblem,
        information::Information,
        options::Options,
        args...;
        kargs...
       )

    Rnorm = parameters.Rnorm
    N = parameters.N
    D = size(problem.bounds, 2)
    fobj = problem.f

    # bounds vectors
    low, up = problem.bounds[1,:], problem.bounds[2,:]

    max_it = 500
    options.iterations = options.iterations == 0 ? max_it : options.iterations
    options.f_calls_limit = options.f_calls_limit == 0 ? options.iterations * N : options.f_calls_limit

    # random initialization for agents.
    P = generate_population(N, problem)

    # Current best
    theBest = get_best(P)

    status = State(theBest, P)
    status.f_calls = N

    # Velocity
    parameters.V = isempty(parameters.V) ? zeros(N,D) : parameters.V
    # Postions
    parameters.X = isempty(parameters.X) ? positions(status) : parameters.X
    # function values
    if isempty(parameters.fitness)

        parameters.fitness = fval.(status.population)
    end

    return status

end

function update_state!(
        status,
        parameters::CGSA,
        problem::AbstractProblem,
        information::Information,
        options::Options,
        args...;
        kargs...
        )

    wMax = parameters.wMax
    N = parameters.N
    wMin = parameters.wMin
    max_it = options.iterations
    iteration = status.iteration
    searchType = :minimize
    chaosIndex = parameters.chaosIndex
    Rnorm = parameters.Rnorm
    Rpower = parameters.Rpower
    ElitistCheck = parameters.ElitistCheck

    X = parameters.X
    V = parameters.V
    fitness = parameters.fitness

    P = status.population
    theBest = status.best_sol

    # iteration
    chValue = wMax-iteration*((wMax-wMin)/max_it)

    #Calculation of M. eq.14-20
    M = massCalculation(fitness,searchType)

    #Calculation of Gravitational constant. eq.13.
    G = Gconstant(iteration, max_it)

    if 1 <= chaosIndex <= 10
        G += chaos(chaosIndex,iteration,max_it,chValue)
    end

    #Calculation of accelaration in gravitational field. eq.7-10,21.
    a = Gfield(M,X,G,Rnorm,Rpower,ElitistCheck,iteration,max_it)

    #Agent movement. eq.11-12
    X, V = move(X,a,V)

    #
    # Checking allowable range.
    # X = correctPop(X, low, up)
    for i = 1:N
        x = reset_to_violated_bounds!(X[i,:], problem.bounds)
        X[i,:] = x
    end

    status.population = create_solutions(X, problem)

    parameters.X = X
    parameters.fitness = fvals(status.population)
    parameters.V = V

    #Evaluation of agents.
    currentBest = get_best(status.population)

    if is_better(currentBest, theBest)
        status.best_sol = currentBest
    end

    stop_criteria!(status, parameters, problem, information, options)

end

function final_stage!(
        status,
        parameters::CGSA,
        problem::AbstractProblem,
        information::Information,
        options::Options,
        args...;
        kargs...
    )
    status.final_time = time()
end


text:
struct LinkedList
    next :: Vector{Int32}
end

function run!(list :: LinkedList, currentcup :: Int, n :: Int)
    wrap(n) = mod1(n, length(list.next))
    for _ in 1:n
        pickedcup1 = list.next[currentcup]
        pickedcup2 = list.next[pickedcup1]
        pickedcup3 = list.next[pickedcup2]
        afterpickedup = list.next[pickedcup3]
        list.next[currentcup] = afterpickedup

        targetcup = wrap(currentcup - 1)
        while targetcup in (pickedcup1, pickedcup2, pickedcup3)
            targetcup = wrap(targetcup - 1)
        end
        aftertargetcup = list.next[targetcup]

        list.next[targetcup] = pickedcup1
        list.next[pickedcup3] = aftertargetcup
        currentcup = afterpickedup
    end
end

function cuporder(start, list, size)
    n = list.next[start]
    s = ""
    for _ in 2:size
        s = join([s, repr(n)])
        n = list.next[n]
    end
    return parse(Int64, s)
end

function cycle(numbers, cycles; part2=false, N=10^6)
    list = LinkedList(Int32[])
    if part2
        sizehint!(list.next, N)
    end
    indexof(n) = findfirst(i -> i == mod1(n, length(numbers)), numbers)
    for n in 1:length(numbers)
        push!(list.next, numbers[mod1(indexof(n) + 1, length(numbers))])
    end
    if part2
        for n in (length(numbers)+1):N
            push!(list.next, n+1)
        end
        list.next[N] = first(numbers)
        list.next[last(numbers)] = length(numbers)+1
    end
    run!(list, first(numbers), cycles)
    if part2
        return Int64(list.next[1]) * list.next[list.next[1]]
    else
        return cuporder(1, list, length(numbers))
    end
end

# test = [3,8,9,1,2,5,4,6,7] # 10:92658374 100:67384529; part2:149245887792
input = [5,3,8,9,1,4,7,6,2]

println("P1: ", cycle(input, 100))
println("P2: ", cycle(input, 10^7; part2 = true))

max_stars_repo_name: eduellery/adventofcode

max_issues_repo_name: eduellery/adventofcode

max_forks_repo_name: eduellery/adventofcode


text:
function MarshallSolarActivityFutureEstimation(arg0::JString, arg1::MarshallSolarActivityFutureEstimation_StrengthLevel)
    return MarshallSolarActivityFutureEstimation((JString, MarshallSolarActivityFutureEstimation_StrengthLevel), arg0, arg1)
end

function MarshallSolarActivityFutureEstimation(arg0::JString, arg1::MarshallSolarActivityFutureEstimation_StrengthLevel, arg2::DataProvidersManager, arg3::TimeScale)
    return MarshallSolarActivityFutureEstimation((JString, MarshallSolarActivityFutureEstimation_StrengthLevel, DataProvidersManager, TimeScale), arg0, arg1, arg2, arg3)
end

function get24_hours_kp(obj::MarshallSolarActivityFutureEstimation, arg0::AbsoluteDate)
    return jcall(obj, "get24HoursKp", jdouble, (AbsoluteDate,), arg0)
end

function get_ap(obj::MarshallSolarActivityFutureEstimation, arg0::AbsoluteDate)
    return jcall(obj, "getAp", Vector{jdouble}, (AbsoluteDate,), arg0)
end

function get_average_flux(obj::MarshallSolarActivityFutureEstimation, arg0::AbsoluteDate)
    return jcall(obj, "getAverageFlux", jdouble, (AbsoluteDate,), arg0)
end

function get_daily_flux(obj::MarshallSolarActivityFutureEstimation, arg0::AbsoluteDate)
    return jcall(obj, "getDailyFlux", jdouble, (AbsoluteDate,), arg0)
end

function get_file_date(obj::MarshallSolarActivityFutureEstimation, arg0::AbsoluteDate)
    return jcall(obj, "getFileDate", DateComponents, (AbsoluteDate,), arg0)
end

function get_instant_flux(obj::MarshallSolarActivityFutureEstimation, arg0::AbsoluteDate)
    return jcall(obj, "getInstantFlux", jdouble, (AbsoluteDate,), arg0)
end

function get_max_date(obj::MarshallSolarActivityFutureEstimation)
    return jcall(obj, "getMaxDate", AbsoluteDate, ())
end

function get_mean_flux(obj::MarshallSolarActivityFutureEstimation, arg0::AbsoluteDate)
    return jcall(obj, "getMeanFlux", jdouble, (AbsoluteDate,), arg0)
end

function get_min_date(obj::MarshallSolarActivityFutureEstimation)
    return jcall(obj, "getMinDate", AbsoluteDate, ())
end

function get_strength_level(obj::MarshallSolarActivityFutureEstimation)
    return jcall(obj, "getStrengthLevel", MarshallSolarActivityFutureEstimation_StrengthLevel, ())
end

function get_supported_names(obj::MarshallSolarActivityFutureEstimation)
    return jcall(obj, "getSupportedNames", JString, ())
end

function get_three_hourly_kp(obj::MarshallSolarActivityFutureEstimation, arg0::AbsoluteDate)
    return jcall(obj, "getThreeHourlyKP", jdouble, (AbsoluteDate,), arg0)
end

function load_data(obj::MarshallSolarActivityFutureEstimation, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end

function still_accepts_data(obj::MarshallSolarActivityFutureEstimation)
    return jcall(obj, "stillAcceptsData", jboolean, ())
end


text:
export generateSFDatabase, loadSFDatabase

function findSFFiles(d::BrukerDatasetStore)
  studies = readdir(d.path)
  bfiles = String[]

  for study in studies
    studypath = joinpath(d.path,study)
    if isdir(studypath)
      experiments = readdir(studypath)
      for exp in experiments
        path = joinpath(d.path,study,exp)
        if _iscalib(path)
          push!(bfiles, path)
        end
      end
    end
  end
  BrukerMDFSFs = readdir("/opt/data/MDF_SFs/")
  for BrukerMDFSF in BrukerMDFSFs
    prefix, ext = splitext(BrukerMDFSF)
     if ext == ".mdf"
      push!(bfiles,joinpath("/opt/data/MDF_SFs/",BrukerMDFSF))
    end
  end
  return bfiles
end

function findSFFiles(d::MDFDatasetStore)
  bfiles = String[]

  path = joinpath(d.path,"calibrations/")

  files = readdir(path)

  for file in files
    prefix, ext = splitext(file)
    if !isdir(file) && tryparse(Int64,prefix) != nothing &&
       (ext == ".mdf" || ext == ".hdf" || ext == ".h5") && !occursin("td.mdf",file)
      try
        push!(bfiles, joinpath(path,file))
      catch e
        @debug "" e
      end
    end
  end

  return bfiles
end

####

function generateSFDatabase(d::DatasetStore, filename::AbstractString)
  fileList = findSFFiles(d)
  A = generateSFDatabase(fileList)
  writedlm(filename, A, ',')
end

function generateSFDatabase(fileList::Vector)

  A = Array{Any}(undef,length(fileList)+1,17)

  # Headerrow
  A[1,1] = "Name"
  A[1,2] = "Gradient"
  A[1,3] = "DFx"
  A[1,4] = "DFy"
  A[1,5] = "DFz"
  A[1,6] = "Size x"
  A[1,7] = "Size y"
  A[1,8] = "Size z"
  A[1,9] = "Bandwidth"
  A[1,10] = "Tracer"
  A[1,11] = "TracerBatch"
  A[1,12] = "DeltaSampleConcentration"
  A[1,13] = "DeltaSampleVolume"
  A[1,14] = "Path"
  A[1,15] = "StartDate"
  A[1,16] = "MeasurementTime"
  A[1,17] = "ExperimentNumber"

  for (k,sf) in enumerate(fileList)
    i=k+1
    MPIFile(sf) do f
    _innerGenerateSFDatabase(A,i,sf,f)
    end
  end
  return A
end

function _innerGenerateSFDatabase(A,i,sf,b)
  A[i,1] = experimentName(b)
  A[i,2] = maximum(acqGradient(b))
  df = vec(dfStrength(b)).*1e3
  A[i,3:5] .= 0.0
  for l=1:min(length(df),3)
    A[i,l+2] = df[l]
  end
  N = calibSize(b)
  A[i,6] = N[1]
  A[i,7] = N[2]
  A[i,8] = N[3]
  A[i,9] = rxBandwidth(b) / 1e6
  A[i,10] = tracerName(b)[1]
  A[i,11] = tracerBatch(b)[1]
  A[i,12] = 0.0#deltaSampleConcentration(b)
  A[i,13] = 0.0#deltaSampleVolume(b)
  A[i,14] = sf #filepath(b)
  A[i,15] = string(acqStartTime(b))
  A[i,16] = 0.0#b["PVM_ScanTimeStr"]
  A[i,17] = experimentNumber(b)
end

function generateSFDatabase(d::MDFDatasetStore)
  oldfile = joinpath(d.path,"SF_DatabaseOld.csv")
  newfile = joinpath(d.path,"SF_Database.csv")
  generateSFDatabase_(d, oldfile, newfile)
end

# HAAACKKK
function generateSFDatabase(d::BrukerDatasetStore)
  oldfile = "/opt/data/SF_DatabaseOld.csv"
  newfile = "/opt/data/SF_Database.csv"
  generateSFDatabase_(d, oldfile, newfile)
end

function generateSFDatabase_(d::DatasetStore, oldfile, newfile)

  if isfile(newfile)
    if isfile(oldfile)
      cp(newfile, oldfile, force=true)
    else
      cp(newfile, oldfile, force=false)
    end
  end

  generateSFDatabase(d, newfile)
end

function loadSFDatabase(d::BrukerDatasetStore)
  if isfile("/opt/data/SF_Database.csv")
    A = readdlm("/opt/data/SF_Database.csv",',')
    if size(A,2) < 16
      A = readdlm("/opt/data/SF_Database.csv",'\t')
    end
    return A
  else
    return nothing
  end
end

function loadSFDatabase(d::MDFDatasetStore)
  files = readdir(calibdir(d))
  @debug "system function database" files
  mdffiles = files[endswith.(files,".mdf") .& .!endswith.(files,"td.mdf")]
  fileList = calibdir(d).*"/".*mdffiles
  A = generateSFDatabase(fileList)
  return A
end

####


text:
"""
    SplitLayer(xs...)

A layer that calls a number of sublayers/mappings with the same input and
returns a tuple of their outputs. Can be used in a regular Flux model:

```julia-repl
julia> m = Chain(Dense(2,3), SplitLayer(Dense(3,2), x->x .* 2))
julia> length(params(m)) == 4
julia> (x,y) = m(rand(2))
(Float32[-1.0541434, 1.1694773], Float32[-3.1472511, -0.86115724, -0.39665926])
```
Comes with a convenient constructor for a SplitLayer built from Dense layers
with given activation(s):
```julia-repl
julia> m = Chain(Dense(2,3), SplitLayer(3, [2,5], σ))
julia> (x,y) = m(rand(2))
(Float32[0.3069554, 0.3362006], Float32[0.437131, 0.4982477, 0.6465078, 0.4523438, 0.5068563])
```

You can also provide just a vector / scalar that should be trained but have the
same value for all inputs (like a lonely bias vector). This functionality is
provided by the `TrainableVector`/`TrainableScalar` types. For vector inputs
they simply return the array they are wrapping. For matrix (i.e. batch) inputs
they return appropriately repeated arrays:
```julia-repl
julia> m = SplitLayer(Dense(2,3), ones(Float32,3))
julia> length(params(m)) == 3
julia> (x,y) = m(rand(2,5))
julia> size(y) == (3,5)
julia> y
3×3 Array{Float32,2}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
```
"""
struct SplitLayer{T<:Tuple}
    layers::T
end

SplitLayer(layers...) = SplitLayer(map(maybe_trainable, layers))

function (m::SplitLayer)(x)
    Tuple(layer(x) for layer in m.layers)
end

@functor SplitLayer

# for use as e.g. shared variance
struct TrainableVector{T<:AbstractArray}
    v::T
end
(v::TrainableVector)(x::AbstractVector) = v.v
(v::TrainableVector)(x::AbstractMatrix) = v.v .* reshape(fillsimilar(v.v,size(x,ndims(x)),1),1,:)
(v::TrainableVector)() = v.v
@functor TrainableVector

# for use as e.g. shared variance
struct TrainableScalar{T<:Real}
    s::AbstractVector{T}
    TrainableScalar{T}(x::T) where T<:Real = new{T}([x])
end
TrainableScalar(x::T) where T<:Real = TrainableScalar{T}(x)
(s::TrainableScalar)(x::AbstractVector) = s.s[1]
(s::TrainableScalar)(x::AbstractMatrix) = fillsimilar(x,size(x,ndims(x)),1) .* s.s
@functor TrainableScalar

maybe_trainable(x) = x
maybe_trainable(x::AbstractArray) = TrainableVector(x)
maybe_trainable(x::Real) = TrainableScalar(x)

fillsimilar(x::AbstractArray, s::Tuple, value::Real) = fill!(similar(x, s...), value)
fillsimilar(x::AbstractArray, s, value::Real) = fill!(similar(x, s), value)
@non_differentiable fillsimilar(::Any, ::Any, ::Any)


text:
module MISC
using ..AGFFileReader: Glass, GlassID, AGF
export SEAWATER, POLYCARB, LAF3, BASF5, CAF2, BASF55, TEO2, TYPEA, KDP, WATER, F_SILICA, SAN, CR39, CDS, SILICA, PYREX, ACRYLIC, PMMA, QUARTZ, COC, POLYSTYR, N15, VACUUM

const SEAWATER = Glass(GlassID(AGF, 1705), 1, 1.78736713, -0.0165099148, 0.00147685351, 0.00140145252, -0.000162869896, 8.41096684e-6, 0.0, 0.0, 0.0, 0.0, 0.334, 2.325, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.339529, -1.0, -1.0, 0, 57.917652, 0, 1.0, 0)
const POLYCARB = Glass(GlassID(AGF, 1706), 1, 2.42838566, -3.86116645e-5, 0.0287574474, -0.000197897366, 0.000148358968, 1.38651935e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 1.014, -0.0002273, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 67.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.58547, -1.0, -1.0, 0, 29.909185, 0, 1.25, 0)
const LAF3 = Glass(GlassID(AGF, 1707), 2, 1.5376, 0.00776161, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.35, 0.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.604046, -1.0, -1.0, 0, 80.83154, 0, 1.0, 0)
const BASF5 = Glass(GlassID(AGF, 1708), 5, 1.57410327, 0.0139786149, 0.00083004511, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.334, 2.325, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.603233, -1.0, -1.0, 0, 42.504788, 0, 1.0, 0)
const CAF2 = Glass(GlassID(AGF, 1709), 2, 0.5675888, 0.00252643, 0.4710914, 0.010078333, 3.8484723, 1200.556, 0.0, 0.0, 0.0, 0.0, 0.23, 9.7, -2.662e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 18.9, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.433849, -1.0, -1.0, 0, 94.995854, 0, 3.181, 0)
const BASF55 = Glass(GlassID(AGF, 1710), 1, 2.8080853, -0.013076515, 0.024961324, 0.0019412734, -0.00015776742, 1.4562956e-5, 0.0, 0.0, 0.0, 0.0, 0.365, 1.06, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.69981, -1.0, -1.0, 0, 34.680447, 0, 3.95, 0)
const TEO2 = Glass(GlassID(AGF, 1711), 2, 2.584, 0.01800964, 1.157, 0.06959044, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 2.274935, -1.0, -1.0, 0, 16.22277, 0, 1.0, 0)
const TYPEA = Glass(GlassID(AGF, 1712), 5, 1.49224183, 0.0101876827, 0.000856657692, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.334, 2.325, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.51509, -1.0, -1.0, 0, 41.585023, 0, 1.0, 0)
const KDP = Glass(GlassID(AGF, 1713), 2, 1.256618, 0.0084478168, 33.89909, 1113.904, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 1.06, -9.718e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 22.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.509182, -1.0, -1.0, 0, 56.224997, 0, 2.339, 0)
const WATER = Glass(GlassID(AGF, 1714), 1, 1.75972105, -0.0112081602, 0.00793533773, -0.000413306965, 8.05068807e-5, -4.54289876e-6, 0.0, 0.0, 0.0, 0.0, 0.4, 0.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [], 1.333044, -1.0, -1.0, 0, 55.794322, 0, 1.0, 0)
const F_SILICA = Glass(GlassID(AGF, 1715), 2, 0.6961663, 0.004679148, 0.4079426, 0.013512063, 0.8974794, 97.9340025, 0.0, 0.0, 0.0, 0.0, 0.21, 3.71, 2.237e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.51, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.458464, -1.0, -1.0, 0, 67.821433, 0, 2.2, 0)
const SAN = Glass(GlassID(AGF, 1716), 1, 2.38687023, -0.00123063994, 0.0229467817, 0.000369810122, 2.67577106e-5, 2.84806182e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 1.014, -0.0002399, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 64.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.56744, -1.0, -1.0, 0, 34.812344, 0, 1.07, 0)
const CR39 = Glass(GlassID(AGF, 1717), 5, 1.485, 0.009, 0.00055, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.85, -0.00034513, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 114.0, -1.0, 3, -1.0, [(2.0, 0.9998, 2.0), (4.0, 0.9968, 2.0), (6.0, 0.9954, 2.0), (8.0, 0.9941, 2.0), (10.0, 0.9926, 2.0), (12.0, 0.956723, 2.0), (14.0, 0.9275, 2.0)], 1.0, -1.0, -1.0, 0, 0.0, 0, 0.0, 0)
const CDS = Glass(GlassID(AGF, 1718), 2, 3.9658282, 0.0558036869, 0.18113874, 0.233146066, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.51, 1.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.0, -1.0, -1.0, 0, 0.0, 0, 1.0, 0)
const SILICA = Glass(GlassID(AGF, 1719), 2, 0.6961663, 0.004679148, 0.4079426, 0.01351206, 0.8974794, 97.934, 0.0, 0.0, 0.0, 0.0, 0.21, 3.71, 2.237e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.51, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.458464, -1.0, -1.0, 0, 67.821443, 0, 2.2, 0)
const PYREX = Glass(GlassID(AGF, 1720), 5, 1.45217, 0.01258, 6.659e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.334, 2.325, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 3.25, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.474009, -1.0, -1.0, 0, 65.386423, 0, 2.23, 0)
const ACRYLIC = Glass(GlassID(AGF, 1721), 1, 2.16330492, 0.0136580764, 0.0256700975, -0.00280976174, 0.0002842279, -9.02124935e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 1.014, -0.0002592, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 60.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.491668, -1.0, -1.0, 0, 55.310192, 0, 1.18, 0)
const PMMA = Glass(GlassID(AGF, 1722), 1, 2.1864582, -0.00024475348, 0.014155787, -0.00044329781, 7.7664259e-5, -2.9936382e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 1.06, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.491756, -1.0, -1.0, 0, 57.440791, 0, 1.16, 0)
const QUARTZ = Glass(GlassID(AGF, 1723), 1, 2.35676495, -0.0113996924, 0.0108741656, 3.32066914e-5, 1.08609346e-5, -3.10123984e-7, 0.0, 0.0, 0.0, 0.0, 0.18, 0.71, -1.487e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.544296, -1.0, -1.0, 0, 70.132045, 0, 0.0, 0)
const COC = Glass(GlassID(AGF, 1724), 1, 2.28448546, 0.0102952211, 0.0373493703, -0.00928409653, 0.00173289808, -0.000115203047, 0.0, 0.0, 0.0, 0.0, 0.334, 2.325, -0.0002278, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 60.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.533732, -1.0, -1.0, 0, 56.227932, 0, 1.02, 0)
const POLYSTYR = Glass(GlassID(AGF, 1725), 1, 2.44598368, 2.21428933e-5, 0.0272988569, 0.000301210852, 8.88934888e-5, -1.75707929e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 1.014, -0.0002943, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 90.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.590481, -1.0, -1.0, 0, 30.866877, 0, 1.05, 0)
const N15 = Glass(GlassID(AGF, 1726), 5, 1.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.365, 1.014, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 1.5, -1.0, -1.0, 0, 0.0, 0, 1.0, 0)
const VACUUM = Glass(GlassID(AGF, 1727), 2, -9.983407e-5, -0.00531510156, -0.0031144725, 0.0103934642, 0.00267863389, 0.0107652513, 0.0, 0.0, 0.0, 0.0, 0.334, 2.325, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.334, 1.0, 25.0), (0.35, 1.0, 25.0), (0.365, 1.0, 25.0), (0.37, 1.0, 25.0), (0.38, 1.0, 25.0), (0.39, 1.0, 25.0), (0.4, 1.0, 25.0), (0.42, 1.0, 25.0), (0.46, 1.0, 25.0), (0.5, 1.0, 25.0), (0.66, 1.0, 25.0), (1.06, 1.0, 25.0), (1.529, 1.0, 25.0), (2.325, 1.0, 25.0)], 0.999728, -1.0, -1.0, 0, 89.195538, 0, 1.0, 0)
end # module

max_stars_repo_name: rambunctiousapple/GlassCat

max_issues_repo_name: rambunctiousapple/GlassCat

max_forks_repo_name: rambunctiousapple/GlassCat