text:
"""
Cache for all model results
"""
struct OptimizationOutputCaches
    data::Dict{OptimizationResultCacheKey, OptimizationOutputCache}
    max_size::Int
    min_flush_size::Int
end

function OptimizationOutputCaches()
    return OptimizationOutputCaches(
        Dict{OptimizationResultCacheKey, OptimizationOutputCache}(),
        0,
        0,
    )
end

function OptimizationOutputCaches(rules::CacheFlushRules)
    return OptimizationOutputCaches(
        Dict{OptimizationResultCacheKey, OptimizationOutputCache}(),
        rules.max_size,
        rules.min_flush_size,
    )
end

function Base.empty!(cache::OptimizationOutputCaches)
    for output_cache in values(cache.data)
        empty!(output_cache)
    end
end

# PERF: incremental improvement if we manually keep track of the current size.
# Could be prone to bugs if we miss a change.
# The number of containers is not expected to be more than 100.

#decrement_size!(cache::OptimizationOutputCaches, size) = cache.size -= size
#increment_size!(cache::OptimizationOutputCaches, size) = cache.size += size

get_max_size(cache::OptimizationOutputCaches) = cache.max_size
get_min_flush_size(cache::OptimizationOutputCaches) = cache.min_flush_size
get_size(cache::OptimizationOutputCaches) =
    reduce(+, (get_size(x) for x in values(cache.data)))

# Leave some buffer because we may slightly exceed the limit.
is_full(cache::OptimizationOutputCaches, cur_size) = cur_size >= cache.max_size * 0.95

function add_output_cache!(cache::OptimizationOutputCaches, model_name, key, flush_rule)
    cache_key = OptimizationResultCacheKey(model_name, key)
    cache.data[cache_key] = OptimizationOutputCache(cache_key, flush_rule)
    @debug "Added cache container for" LOG_GROUP_SIMULATION_STORE model_name key flush_rule
    return
end

"""
Return true if the cache has data that has not been flushed to storage.
"""
function has_dirty(cache::OptimizationOutputCaches)
    for output_cache in values(cache.data)
        if has_dirty(output_cache)
            return true
        end
    end

    return false
end

get_output_cache(cache::OptimizationOutputCaches, key::OptimizationResultCacheKey) =
    cache.data[key]

function get_output_cache(
    cache::OptimizationOutputCaches,
    model_name,
    key::OptimizationContainerKey,
)
    cache_key = OptimizationResultCacheKey(model_name, key)
    return get_output_cache(cache, cache_key)
end

"""
Return true if the data for `timestamp` is stored in cache.
"""
function is_cached(cache::OptimizationOutputCaches, model_name, key, index)
    cache_key = OptimizationResultCacheKey(model_name, key)
    return is_cached(cache, cache_key, index)
end

is_cached(cache::OptimizationOutputCaches, key, timestamp::Dates.DateTime) =
    has_timestamp(cache.data[key], timestamp::Dates.DateTime)

is_cached(cache::OptimizationOutputCaches, key, ::Int) = false

"""
Log the cache hit percentages for all caches.
"""
function log_cache_hit_percentages(cache::OptimizationOutputCaches)
    for key in keys(cache.data)
        output_cache = cache.data[key]
        cache_hit_pecentage = get_cache_hit_percentage(output_cache)
        @debug "Cache stats" LOG_GROUP_SIMULATION_STORE key cache_hit_pecentage
    end
    return
end

"""
Read the result from cache. Callers must first call [`is_cached`](@ref) to check if the
timestamp is present.
"""
function read_result(cache::OptimizationOutputCaches, model_name, key, timestamp)
    cache_key = OptimizationResultCacheKey(model_name, key)
    return read_result(cache, cache_key, timestamp)
end

read_result(cache::OptimizationOutputCaches, key, timestamp) =
    cache.data[key].data[timestamp]


text:
"""
NumDiffTools.gradient(f, x::AbstractArray, [abstol, reltol, h, t, maxiters])

Returns gradient `∇f` of `f` evaluted at `x`, assuming `f` is called as `f(x)`.

Optional Arguments:\n
`abstol`: Target absolute error in derivative estimate (default `1e-12`)\n
`reltol`: Target relative error in derivative estimate (default `1e-12`)\n
`h`: Initial step size (default `0.1`)\n
`t`: Ratio of step sizes for each Richardson iteration (default `2.0001`.  Should not be an exact integer for periodic functions)\n
`maxiters`: Maximum number of Richardson iterations to perform (default `6`)\n
"""
function gradient(f, x::AbstractArray; abstol::Real = 1e-12, reltol::Real = 1e-12, h::Real = 0.1, t::Real = 2.0001, maxiters::Integer = 6)
    fx = f(x)

    # Initialize Work Arrays
    A1 = Vector{typeof(fx)}(undef, maxiters)
    A2 = Vector{typeof(fx)}(undef, maxiters)

    # Vector of step sizes
    hh = [h/t^(i - 1) for i = 1:maxiters]
    # Vector of orders
    k = 1:maxiters

    result = Vector{typeof(fx)}(undef, length(x))
    dx = zero(x)
    for i = 1:length(x)
        dx[i] = 1

        result[i] = extrapolate(h->(f(x + h*dx) - fx)/h,          # Function to extrapolate
                                A1, A2,                           # Work Arrays
                                hh, k,                            # Sequence of step sizes (h)
                                abstol, reltol                    # tolerances
                                )
        dx[i] = 0
    end
    result
end

"""
NumDiffTools.gradient!(result::AbstractArray, f, x::AbstractArray, [abstol, reltol, h, t, maxiters])

Stores gradient `∇f` of `f` evaluted at `x` in `result`, assuming `f` is called as `f(x)`.

Optional Arguments:\n
`abstol`: Target absolute error in derivative estimate (default `1e-12`)\n
`reltol`: Target relative error in derivative estimate (default `1e-12`)\n
`h`: Initial step size (default `0.1`)\n
`t`: Ratio of step sizes for each Richardson iteration (default `2.0001`.  Should not be an exact integer for periodic functions)\n
`maxiters`: Maximum number of Richardson iterations to perform (default `6`)\n
"""
function gradient!(result::AbstractArray, f, x::AbstractArray; abstol::Real = 1e-12, reltol::Real = 1e-12, h::Real = 0.1, t::Real = 2.0001, maxiters::Integer = 6)
    fx = f(x)

    # Initialize Work Arrays
    A1 = Vector{typeof(fx)}(undef, maxiters)
    A2 = Vector{typeof(fx)}(undef, maxiters)

    # Vector of step sizes
    hh = [h/t^(i - 1) for i = 1:maxiters]
    # Vector of orders
    k = 1:maxiters

    dx = zero(x)
    for i = 1:length(x)
        dx[i] = 1

        result[i] = extrapolate(h->(f(x + h*dx) - fx)/h,          # Function to extrapolate
                                A1, A2,                           # Work Arrays
                                hh, k,                            # Sequence of step sizes (h)
                                abstol, reltol                    # tolerances
                                )
        dx[i] = 0
    end
    nothing
end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ a4f0a914-1ec9-11eb-3c22-db87a9776aaa
begin
	using JSON
	using DataFrames
end

# ╔═╡ 5f83510e-1fc5-11eb-3c40-83a844bf2907
using Turing

# ╔═╡ 2414104e-1d41-11eb-37cb-9b5cee78a98b
md"## Creating our conjectures"

# ╔═╡ 3b7f8a1e-1d42-11eb-1ab1-c11fda416e80
md"""If there is one thing we can all agree on, it is that the reality in which we live is complex. The explanation for the things we usually see, and have naturalized in our daily lives, is usually quite complex and requires abstraction from what "is simply seen".

In order to give an explanation and gain a deeper understanding of the things we see, we tend to generate models that seek to explain them in a simple and generalized way. In this way we can reduce the noise of our observations to general rules that "govern" them.

For example, it is obvious to everyone that if we push a glass it will move in the same direction as we did. We also know that if we keep pushing it and it goes beyond the limits of the table, it will fall to the floor. But one thing is to have the intuition of what's going to happen, and another is to have an understanding of the laws that govern that movement. In this case, they are the Newton´s Law´s of motion:

"""

# ╔═╡ 4e4385b4-203b-11eb-0c5b-cbe41156401b
md"$ \vec{F}^{\} = m*\vec{a}^{\} $"

# ╔═╡ 2dca2f42-203b-11eb-2cb3-0d6a3823cdcd
md"""In this way, and with only one formula, it is possible to gain an understanding that is generalizable to many aspects of reality.

#### Observable variables vs Latent variables

Now, it is worth noting that in this case all the variables that make up our model are observable. This means that they can be measured directly.

In the case of a glass, we could weigh it with a scale. Then, by pushing it, we could measure the acceleration it acquired and from these two measurements we could obtain the force we applied to it. So, every parameter of the model is fully defined.

However, as we try to advance in our understanding of reality, we arrive at more and more complex models and many times we are not so lucky to be able to define them with simple observable variables.

For example, this is very common in the economic sciences, where models are created with variables such as "quality of life". Economists will try to measure this latent variable with other variables that can be observed (such as quality of life, schooling rate, number of hospitals for a certain number of inhabitants, etc), but that do not have an obvious and direct relationship as if they had newton's equations.

This type of latent variables are used in different models to gain greater abstraction and to be able to obtain information that a priori is not found at first sight in the data. For example, in the case of economics, from concrete measures of a country's economy it is possible to generalize knowledge and be able to infer an abstract variable such as quality of life.

### Bayesian hierarchical models

The Bayesian framework allows us to build statistical models that can generalize the information obtained from the data and make inferences from latent variables.

A nice way to think about this kind of models is that they allow us to build our "story" about which are the variables that generate the data we are observing. Basically, they allow us to increase the "depth" of our model by indicating that the parameters of our prior distributions also follow other probability distributions.

This sure is sounding very strange. Don't worry, let's move on to an example to clarify it.

#### Football analysis

Let's imagine for a moment that we are brilliant statisticians. We find ourselves looking for new interesting challenges to solve and we come across a sports bookmaker. They tell us that they want to expand into football betting and that they would like us to be able to build a model that allows them to analyze the strengths and weaknesses of English Premier League teams. They are interested because they want to be able to predict possible outcomes and thus be able to price the bets.

The problem is that, as they have never worked in this sector before, they only have the results of the league matches. So what can we do? """

# ╔═╡ 0c0268dc-1ed3-11eb-2f73-b3706305b298
md"We have the data stored in a specific format called JSON, so the first thing to do is to parse and visualize it"

# ╔═╡ cbaaf72c-1eb4-11eb-3509-932b337f270b
begin
england_league = JSON.parsefile("matches_England.json")
matches_df = DataFrame(home = [], away = [], score_home = [], score_away = [])
end;

# ╔═╡ ab6fcb22-1ed2-11eb-2749-3bef16911972
begin
	matches = []
	for match in england_league
		push!(matches, split(match["label"], ","))
end
end

# ╔═╡ d2a81a14-1ed2-11eb-1307-bfb169aebbb4
begin
for match in matches
	home, away = split(match[1], " - ")
	score_home, score_away = split(match[2], " - ")

	push!(matches_df,[home, away, parse(Int,score_home), parse(Int,score_away)])
end
end

# ╔═╡ dd80fc38-1ed2-11eb-0e8c-49859d27c72c
matches_df

# ╔═╡ 7f0c72d4-1f9e-11eb-06f1-cb587d7e5436
teams = unique(collect(matches_df[1]))

# ╔═╡ 868e5c8c-1ed3-11eb-291c-fb26d512c103
md"""So, we have the data of the 380 matches that were played in the Premier League 2017/2018 and our challenge is to be able to analyze the characteristics of these teams.

A priori it may seem that we are missing data, that with the data we have we cannot infer "characteristics" specific to each team. At most, it might be possible to see who the teams that scored the most goals, the averages of goals per game or how the positions were after the tournament, but to obtain characteristics of the teams? how could we face this problem?

#### Creating our stories

Okay, let's see what information we have from our data:
On one hand we have specified the names of each team and which one is local. On the other hand, we have the number of goals scored.

A possible approach to this data is to realize that the goals scored by each team can be modeled with a poisson distribution.

Why? You have to remember that this distribution describes "arrivals" - discrete events - in a continuum. For example, it is widely used to describe customer arrivals to a location as time passes or failures in continuous industrial processes (e.g. failure in the production of a pipe).

In this particular case, we could propose that the goals scored by a team are the discrete events that occur in the time continuum that the game last:

"""

# ╔═╡ 04a8f16c-1ed5-11eb-3979-af35c6d9dea1
md"$Score \sim Poisson(θ)$"

# ╔═╡ 9f1e1618-1edb-11eb-02af-9fa5c8904163
md"""Well, we have an improvement. We've already told our model how to think about goals scored.

Now we can use the flexibility of Bayesianism to indicate what the "goal rate" of our Poisson depends on.  You can think of it literally as the number of goals a team scores per unit of time. And this is where we have to take advantage of all the information provided by the data set.

As expected, this rate has to be particular to each match the team plays and take into account the opponent. We can therefore propose that the scoring rate of each team (in each particular match) depends on the "attacking power" of the team on the one hand, and the "defensive power" of the opponent on the other:

"""

# ╔═╡ 42104900-1f80-11eb-019f-a7c561ca3e4b
md"$θ_{team1} \sim att_{team1} + def_{team2}$"

# ╔═╡ 9d0f3a1e-1f80-11eb-3422-fb7b4696cc13
md"""In this way we could be capturing, from the results of each game, the attack and defence strengths of each team.

Another latent variable that we could obtain, given the data, is if there is an effect that increases (or decreases) the goal rate related to whether the team is local or not. This would also help - in case there is indeed an effect - in not the attack and defence parameters be disrupted by having to "contain" that information."""

# ╔═╡ 0b1f9c68-1f83-11eb-11e7-e7614c021c05
md"$θ_{home} \sim home + att_{home} + def_{away}$
$θ_{away} \sim att_{away} + def_{home}$"

# ╔═╡ 043c438c-1f84-11eb-2cd7-6dd07a125910
md"""This leaves one attack and one defense parameter for each team, and a global league parameter that indicates the effect of being local on the scoring rate.

#### Letting the information flow

Okay, we are already getting much closer to the initial goal we set. As a last step, we must be able to make the information flow between the two independent poissons that we proposed to model the score of each of the two teams that are playing. We need to do that precisely because we have proposed that the poissons are independent, but we need that when making the inference of the parameters the model can access the information from both scores so it can catch the correlation between them. In other words, we have to find a way to interconnect our model.

And that is exactly what hierarchical Bayesian models allow us to do. How? By letting us choose probability distributions for the parameters that represent the characteristics of both equipment. With the addition that these parameters will share the same prior distributions. Let's see how:

The first thing to do, as we already know, is to assign the prior distributions of our attack and defense parameters. A reasonable idea would be to propose that they follow a normal distribution since it is consistent that there are some teams that have a very good defense, so the parameter would take negative values; or there may be others that have a very bad one, taking positive values (since they would "add up" to the goal rate of the opposing team). The normal distribution allows us to contemplate both cases.

Now, when choosing the parameters we are not going to stop and assign fixed numbers, but we will continue to deepen the model and add another layer of distributions:

"""

# ╔═╡ 4a31541e-1f97-11eb-038e-736446224c21
md"$att_{t} \sim Normal(μ_{att}, σ_{att})$
$def_{t} \sim Normal(μ_{def}, σ_{def})$"

# ╔═╡ e3a99868-1f97-11eb-1872-79c8f6abe1be
md"Where the t sub-index is indicating us that there are a couple of these parameters for each team.

Then, as a last step to have our model defined, we have to assign the priority distributions that follow the parameters of each normal distribution. We have to define our hyper priors.

"

# ╔═╡ c3ca1cc4-1f98-11eb-3c23-33acb95fbc03
md"$μ_{att}, μ_{def} \sim Normal(0, 0.1)$
$σ_{att}, σ_{def} \sim Exponential(1)$"

# ╔═╡ 4106e55a-1f99-11eb-110e-b79221774f12
md"We must not forget the parameter that represents the advantage of being local"

# ╔═╡ bff4fdae-1f99-11eb-2d94-aba7bf9d097f
md"$home \sim Normal(0,1)$"

# ╔═╡ 2b135202-1fa7-11eb-324d-93a5eda41e3e
md"""Now that our model is fully define, let's add one last restriction to the characteristics of the teams to make it easier to compare them: subtract the average of all the attack and defence powers from each one. In this way we will have the features centred on zero, with negative values for the teams that have less attacking power than the average and positive values for those that have more. As we already said, the opposite analysis applies to the defence, negative values are the ones that will indicate that a team has a strong defence as they will be "subtracting" from the scoring rate of the opponent. This is equivalent to introducing the restriction:
"""

# ╔═╡ ddc7b166-203a-11eb-0ac2-056c309bb590
md"$\sum att_{t} = 0$
$\sum def_{t} = 0$"

# ╔═╡ 8e164be0-203b-11eb-293a-6b32dba96133
md"Let's translate all this into Turing code:"

# ╔═╡ 51a309bc-2033-11eb-10c0-ed17545df33d
begin
	@model function football_matches(home_teams, away_teams, score_home, score_away, teams)
	#hyper priors
	σatt ~ Exponential(1)
	σdeff ~ Exponential(1)
	μatt ~ Normal(0,0.1)
	μdef ~ Normal(0,0.1)

	home ~ Normal(0,1)

	#Team-specific effects
	att ~ filldist(Normal(μatt, σatt), length(teams))
	def ~ filldist(Normal(μatt, σdeff), length(teams))

	dict = Dict{String, Int64}()
	for (i, team) in enumerate(teams)
		dict[team] = i
	end

	#Zero-sum constrains
	offset = mean(att) + mean(def)

	log_θ_home = Vector{Real}(undef, length(home_teams))
	log_θ_away = Vector{Real}(undef, length(home_teams))

	#Modeling score-rate and scores (as many as there were games in the league)
	for i in 1:length(home_teams)
		#score-rate
		log_θ_home[i] = home + att[dict[home_teams[i]]] + def[dict[away_teams[i]]] - offset
		log_θ_away[i] = att[dict[away_teams[i]]] + def[dict[home_teams[i]]] - offset
		#scores
		score_home[i] ~ LogPoisson(log_θ_home[i])
		score_away[i] ~ LogPoisson(log_θ_away[i])
	end

	end
end

# ╔═╡ f3f0b6f8-2033-11eb-0f9e-d951d791001d
md"""As you can see, the turing code is very clear and direct. In the first block we define our hyperpriors for the distributions of the characteristics of the equipment.

In the second one, we define the priors distributions that will encapsulate the information about the attack and defense powers of the teams. With the *filldist* function we are telling Turing that we need as many of these parameters as there are teams in the league *length(teams)*

Then, we calculate the average of the defense and attack parameters that we are going to use to centralize those variables, and we use the LogPoisson distribution to allow the theta to take some negative value in the inference process and give more sensitivity to the parameters that make it up.

As we said before, we will model the thetas for each game played in the league, that's why the *for* of the last block goes from 1 to *length(home_teams)*, which is the list that contains who was the local team of each game played.

So let´s run it and see if all of this effort was worth it:
"""

# ╔═╡ f3f4bfba-203f-11eb-142c-a187112744d2
model = football_matches(matches_df[1], matches_df[2], matches_df[3], matches_df[4], teams)

# ╔═╡ 2ce4435c-2040-11eb-1670-e39ad4cc690c
posterior = sample(model, NUTS(),3000);

# ╔═╡ 9e50e56e-5674-11eb-2919-4bde8964bc59
rand(LogPoisson(0.2), 100)

# ╔═╡ 8c1dd9de-2040-11eb-35a3-39e8196577a1
md"#### Analyzing the results
In order to compare and corroborate that the inference of our model makes sense, it is key to have the ranking table of how the teams actually performed in the 2017/2018 Premier League.
"

# ╔═╡ 0f9406ae-2054-11eb-3f0c-1d03ba779f18
begin
	table_positions =
	[11, 5, 9, 4, 13, 14, 1, 15, 12, 6, 2, 16, 10, 17, 20, 3, 7, 8, 19, 18]

	games_won =
	[32, 25, 23, 21, 21, 19, 14, 13, 12, 12, 11, 11, 10, 11, 9, 9, 7, 8, 7, 6]

	teams_ = []
	for i in table_positions
		push!(teams_, teams[i])
	end

	table_position_df = DataFrame(Table_of_positions = teams_, Wins = games_won)

end

# ╔═╡ 099e4188-2054-11eb-1e7e-69bdb2b0202c
md"Let's now explore a little bit the a posteriori values we obtained."

# ╔═╡ 16572944-2045-11eb-38d0-fb20d981e3e9
begin
	post_att = collect(get(posterior, :att)[1])
	post_def = collect(get(posterior, :def)[1])
	post_home = collect(get(posterior, :home)[1])
end;

# ╔═╡ 4ae7d968-206f-11eb-1925-c5b5a31e6940
begin
using Plots
histogram(post_home, legend=false, title="Posterior distribution of home parameter")
end

# ╔═╡ 93d3bf98-5684-11eb-361e-c7b0f1a3bb0e
get(posterior, :att)[:att]

# ╔═╡ 18f36e1a-5682-11eb-2826-ab4d4c457af3
posterior.value.data

# ╔═╡ 6f096138-5682-11eb-179f-f973803b9268
collect(get(posterior, :att)[1])

# ╔═╡ 8cec9318-206e-11eb-3caa-3bad9d788d3c
md"As a first measure to analyze, it is interesting to see and quantify (if any) the effect that being local has on the score rate:"

# ╔═╡ 8b77d12c-206f-11eb-19a5-f557157ac05e
mean(post_home)

# ╔═╡ b56df7ea-206f-11eb-18b9-0544cc596ca5
md"So, to include in the model the parameter home was a good idea. indeed being local provides a very big advantage.

Beyond the fact that it is interesting to be able to quantify how much the location influences the scoring rate of the teams, including it in the analysis allow us to have better estimates of the defense and attack parameters of the teams. This is true because if it had not been included, this positive effect would have manifested itself in the only parameters it would have found, the attack and defense parameters, deforming the real measure of these.

So, being confident that we are on the right track, let´s find the attack and defence parameters of each team."

# ╔═╡ 554b4ef0-2045-11eb-3218-050ec936f1aa
begin
	teams_att = []
	teams_def = []
	for i in 1:length(post_att)
		push!(teams_att, post_att[i])
		push!(teams_def, post_def[i])
	end
end

# ╔═╡ 8a4df438-2045-11eb-2798-af592dbbffb5
md"This way we obtain all the samples of the posterior distributions for each one of the parameters of each equipment. Scroll right to explore the entire array."

# ╔═╡ 6b1a1520-2045-11eb-1302-0778b4e7a836
teams_att

# ╔═╡ da93bf52-2045-11eb-2b6f-3f31f57d8653
md"For example, if we would like to se the posterior distribution of the attack parameter for Burnley:"

# ╔═╡ 1fa92570-2046-11eb-3541-e5fccf7afabd
teams[1]

# ╔═╡ 255c4d08-2046-11eb-1cf0-b32d2cd711bf
histogram(teams_att[1], legend=false, title="Posterior distribution of Burnley´s attack power")

# ╔═╡ 473e6a46-205b-11eb-2731-7b47b99d86a5
mean(teams_att[1])

# ╔═╡ 917d4550-2058-11eb-174b-65436a4e6cc8
md"Comparing it to the attacking power of Manchester City, champion of the Premier league:"

# ╔═╡ c2de4714-2058-11eb-1a0d-f1bf96649844
teams[11]

# ╔═╡ 0103406e-2059-11eb-333f-7b5814f21a54
begin
	histogram(teams_att[11], legend=false, title="Posterior distribution of Manchester City´s attack power")
end

# ╔═╡ 579a4d9e-205b-11eb-2352-83beca145211
mean(teams_att[11])

# ╔═╡ 6924f372-2059-11eb-2b98-6f4b65b777f4
md"Okay, when comparing the league champion against a mid-table team, we can clearly see the superiority in attack. For now, it seems that the inference comes in handy.

Let's try now to have an overview of the attacking powers of each team. To do this, just take the average of each and plot it next to the standard deviation
"

# ╔═╡ ece1d052-2060-11eb-2b50-ef46663ed88f
begin
	teams_att_μ = mean.(teams_att)
	teams_def_μ = mean.(teams_def)
	teams_att_σ = std.(teams_att)
	teams_def_σ = std.(teams_def)
end;

# ╔═╡ 43a39b30-2061-11eb-28c5-37cbc4aed84b
md"""Remember that the "." operator is used for broadcasting. This means that it will apply the function to each component of the array"""

# ╔═╡ 41ca2316-2062-11eb-3fe6-f33a789fd578
begin
	teams_att_μ
	sorted_att = sortperm(teams_att_μ)
	abbr_names = [t[1:3] for t in teams]
end;

# ╔═╡ 0ef8a16e-2063-11eb-1391-6b3b959de541
begin
	abbr_names[5] = "Mun"
	abbr_names[10] = "Whu"
	abbr_names[11] = "Mci"
	abbr_names[16] = "Bou"
	abbr_names[18] = "Wba"
	abbr_names[19] = "Stk"
	abbr_names[20] = "Bha"
end;

# ╔═╡ 1eb73246-2063-11eb-0f20-416e433b5144
sorted_names = abbr_names[sorted_att]

# ╔═╡ abace764-2062-11eb-2949-fd9e86b50f17
begin
	scatter(1:20, teams_att_μ[sorted_att], grid=false, legend=false, yerror=teams_att_σ[sorted_att], color=:blue, title="Premier league 17/18 teams attack power")
	annotate!(collect(1:20), teams_att_μ[sorted_att] .+ 0.238, text.(sorted_names, :black, :center, 8))
	ylabel!("Mean team attack")
end

# ╔═╡ ebd267a4-2064-11eb-02a7-b93465de6319
md"""Although there is a high correlation between the attacking power of each team and its position on the table after the league ends, it is clear that this is not enough to explain the results. For example, Manchester City was the league's sub-champion, but only appeared in fifth place.

Let's explore what happens to the defence power: """

# ╔═╡ f6ad5f52-2065-11eb-0be7-f159941e0d89
begin
	sorted_def = sortperm(teams_def_μ)
	sorted_names_def = abbr_names[sorted_def]
end

# ╔═╡ 539e68b4-2066-11eb-09c2-b337241c36bc
begin
	scatter(1:20, teams_def_μ[sorted_def], grid=false, legend=false, yerror=teams_def_σ[sorted_def], color=:blue, title="Premier league 17/18 teams defence power")
	annotate!(collect(1:20), teams_def_μ[sorted_def] .+ 0.2, text.(sorted_names_def, :black, :center, 8))
	ylabel!("Mean team defence")
end

# ╔═╡ 69e7e54c-2067-11eb-21a1-cb72f4b423cd
md"To read this graph we have to remember that the defense effect is better the more negative it is, since it is representing the scoring rate that takes away from the opponent team. As we already said:"

# ╔═╡ f59d488a-2067-11eb-1d4b-1d9dd97ea39e
md"$θ_{team1} \sim att_{team1} + def_{team2}$"

# ╔═╡ 27d5a290-2068-11eb-06f5-0b16f6ec6cc5
md"As the $def_{team2}$ is adding up in the equation, if it take negative values, it is going to start substracting the scoring rate of the oponent.

Things, then, begin to make a little more sense. Now we can see that Manchester United is the team with the strongest defence, so being second in the overall is not extrange.

To gain a deeper understanding of what´s going on here, let's chart both characteristics together. This is going to let us see the combined effect they have. Also i´m going to add the final position of each team to improve the interpretability."

# ╔═╡ b8f9bae8-206a-11eb-1426-031f6fd05fd6
begin
	table_position =
	[11, 5, 9, 4, 13, 14, 1, 15, 12, 6, 2, 16, 10, 17, 20, 3, 7, 8, 19, 18]
	position = sortperm(table_position)
end

# ╔═╡ ee45d48e-206a-11eb-0edf-2b8b893bb583
begin
scatter(teams_att_μ, teams_def_μ, legend=false)
annotate!(teams_att_μ, teams_def_μ.+ 0.016, text.(abbr_names, :black, :center, 6))
annotate!(teams_att_μ, teams_def_μ.- 0.016, text.(position, :left, :center, 5))

xlabel!("Mean team attack")
ylabel!("Mean team defence")
end

# ╔═╡ c794f45e-206b-11eb-33c4-05bc429ba846
md"Well, great! Now we have some interesting information to analyze the teams and the league in general. It´s easier now to perceive how the two features interact with each other, comparing between teams and being able to see how that affects the final position.

For example, looking at the cases of Liverpool and Tottenham, or Leicester City and Everton; one could say (against general common sense) that the power of defense has a greater effect on the performance of each team than the attack. But we leave you to do those analysis for the betting house.

Well, we went from having a problem that seemed almost impossible to have a solid solution, with a quantitative analysis of the characteristics of each team. We even know how much the localization of the teams increases the scoring rate. We were able to achieve this thanks to the hierarchical framework that Bayesianism provides us. Using this tool allows us to create models proposing latent variables that cannot be observed, to infer them and to gain a much deeper and more generalized knowledge than we had at first. You just have to imagine a good story."

# ╔═╡ fcb7e8d2-2071-11eb-020d-7bb1d53f8a6d
md"### Bibliography

- [Paper of Gianluca Baio and Marta A. Blangiardo](https://discovery.ucl.ac.uk/id/eprint/16040/1/16040.pdf)
- [Post of Daniel Weitzenfeld](http://danielweitzenfeld.github.io/passtheroc/blog/2014/10/28/bayes-premier-league/)

"

# ╔═╡ fac413ec-5684-11eb-3167-d56d967247fc
mun_att_post = collect(get(posterior, :att)[:att])[5][:,1];

# ╔═╡ cc6681a2-5688-11eb-1d8e-e1edcb80a20d
mun_def_post = collect(get(posterior, :def)[:def])[5][:,1];

# ╔═╡ f01e7438-5688-11eb-2900-8df62e9feb18
liv_att_post = collect(get(posterior, :att)[:att])[4][:,1];

# ╔═╡ fae4e938-5688-11eb-3ab8-e3a8a4a61692
liv_def_post = collect(get(posterior, :def)[:def])[4][:,1];

# ╔═╡ 84d5e7a2-56dc-11eb-0d54-b790c7c714fb

# ╔═╡ 7704b9e6-56a0-11eb-0ade-2194a0845f61
post_home;

# ╔═╡ 03f063cc-5689-11eb-3d90-b7a5991332ba
# This function simulates matches given the attach, defense and home parameters.
# The first pair of parameters alwas correspond to the home team.
function simulate_matches_(att₁, def₁, att₂, def₂, home, n_matches, home_team=1)
	if home_team == 1
		logθ₁ = home + att₁ + def₂
		logθ₂ = att₂ + def₁

	elseif home_team == 2
		logθ₁ = att₁ + def₂
		logθ₂ = home + att₂ + def₁
	else
		return DomainError(home_team, "Invalid home_team value")
	end

	scores₁ = rand(LogPoisson(logθ₁), n_matches)
	scores₂ = rand(LogPoisson(logθ₂), n_matches)

	results = [(s₁, s₂) for (s₁, s₂) in zip(scores₁, scores₂)]

	return results
end

# ╔═╡ a8a9179c-5693-11eb-1758-e7d34722c7fa
mat = simulate_matches_(0.2, -0.1, 0.4, -0.3, 0.2, 1000, 2)

# ╔═╡ 7e82f15c-5695-11eb-20f1-85763cbfa572
max_h = maximum(map(x -> x[1], mat))

# ╔═╡ bd11306c-5697-11eb-22ed-99ae7703be59
max_a = maximum(map(x -> x[2], mat))

# ╔═╡ d99e52fa-5697-11eb-12fe-ab4ae36ab362
function simulate_matches(team1_att_post, team1_def_post, team2_att_post,				team2_def_post, home_post, n_matches)

	team1_as_home_results = Tuple{Int64,Int64}[]
	team2_as_home_results = Tuple{Int64,Int64}[]

	for (t1_att, t1_def, t2_att, t2_def, home) in zip(team1_att_post, team1_def_post, team2_att_post, team2_def_post, home_post)

		team1_as_home_results = vcat(team1_as_home_results, simulate_matches_(t1_att, t1_def, t2_att, t2_def, home, n_matches, 1))

		team2_as_home_results = vcat(team2_as_home_results, simulate_matches_(t1_att, t1_def, t2_att, t2_def, home, n_matches, 2))
	end

	return team1_as_home_results, team2_as_home_results
end

# ╔═╡ 5a3d6966-56a0-11eb-0db1-c5a06ea2dd83
kjj, kjjj = simulate_matches(mun_att_post, mun_def_post, liv_att_post, liv_def_post, post_home, 1000)

# ╔═╡ 9d1f1e7e-56a1-11eb-2539-f5da33710a32
function match_heatmaps(team1_as_home_results, team2_as_home_results)

	max_t1_as_home = maximum(map(x -> x[1], team1_as_home_results))
	max_t2_as_away = maximum(map(x -> x[2], team1_as_home_results))

	max_t1_as_away = maximum(map(x -> x[1], team2_as_home_results))
	max_t2_as_home = maximum(map(x -> x[2], team2_as_home_results))

	matrix_t1_as_home = zeros(Int64, (max_t1_as_home + 1, max_t2_as_away + 1))
	matrix_t2_as_home = zeros(Int64, (max_t1_as_away + 1, max_t2_as_home + 1))

	for match in team1_as_home_results
		matrix_t1_as_home[match[1] + 1, match[2] + 1] += 1
	end

	for match in team2_as_home_results
		matrix_t2_as_home[match[1] + 1, match[2] + 1] += 1
	end

	gr()
	#heat_t1_home = heatmap(0:max_t1_as_home,
	heatmap(0:max_t1_as_home,
								0:max_t2_as_away,
								matrix_t1_as_home,
								c=cgrad([:blue, :white,:red, :yellow]),
								xlabel="Team1 score", ylabel="Team2 score",
								title="Team1 as home")
	savefig("./t1_as_home")

	#heat_t2_home = heatmap(0:max_t1_as_away,
	heatmap(0:max_t1_as_away,
								0:max_t2_as_home,
								matrix_t2_as_home,
								c=cgrad([:blue, :white,:red, :yellow]),
								xlabel="Team1 Score", ylabel="Team2 Score",
								title="Team2 as home")
	savefig("./t2_as_home")
	#plot(heat_t1_home, heat_t2_home, layout=(2,1), size=(900, 400))
	#current()

	return matrix_t1_as_home, matrix_t2_as_home
end

# ╔═╡ d7bfc106-56d0-11eb-21b3-3b9036905041
m_1_home, m_t2_home = match_heatmaps(kjj, kjjj)

# ╔═╡ 9ceba0fa-56d4-11eb-2696-21f4babf699c
begin
	gr()
	heatmap(collect(0:(size(m_t2_home, 1)-1)),
			collect(0:(size(m_t2_home, 2)-1)),
			m_t2_home,
			c=cgrad([:blue, :white,:red, :yellow]),
			xlabel="Team1 score", ylabel="Team2 score",
			title="Team1 as home")
end

# ╔═╡ 0da9d3a0-56d5-11eb-08a9-c9e83106dfc3
m_t1_home

# ╔═╡ db8ed140-5735-11eb-2078-2f90d15d6325
"hola"

# ╔═╡ Cell order:
# ╟─2414104e-1d41-11eb-37cb-9b5cee78a98b
# ╟─3b7f8a1e-1d42-11eb-1ab1-c11fda416e80
# ╟─4e4385b4-203b-11eb-0c5b-cbe41156401b
# ╟─2dca2f42-203b-11eb-2cb3-0d6a3823cdcd
# ╠═a4f0a914-1ec9-11eb-3c22-db87a9776aaa
# ╠═0c0268dc-1ed3-11eb-2f73-b3706305b298
# ╠═cbaaf72c-1eb4-11eb-3509-932b337f270b
# ╠═ab6fcb22-1ed2-11eb-2749-3bef16911972
# ╠═d2a81a14-1ed2-11eb-1307-bfb169aebbb4
# ╠═dd80fc38-1ed2-11eb-0e8c-49859d27c72c
# ╠═7f0c72d4-1f9e-11eb-06f1-cb587d7e5436
# ╟─868e5c8c-1ed3-11eb-291c-fb26d512c103
# ╟─04a8f16c-1ed5-11eb-3979-af35c6d9dea1
# ╟─9f1e1618-1edb-11eb-02af-9fa5c8904163
# ╟─42104900-1f80-11eb-019f-a7c561ca3e4b
# ╟─9d0f3a1e-1f80-11eb-3422-fb7b4696cc13
# ╟─0b1f9c68-1f83-11eb-11e7-e7614c021c05
# ╟─043c438c-1f84-11eb-2cd7-6dd07a125910
# ╟─4a31541e-1f97-11eb-038e-736446224c21
# ╟─e3a99868-1f97-11eb-1872-79c8f6abe1be
# ╟─c3ca1cc4-1f98-11eb-3c23-33acb95fbc03
# ╟─4106e55a-1f99-11eb-110e-b79221774f12
# ╟─bff4fdae-1f99-11eb-2d94-aba7bf9d097f
# ╟─2b135202-1fa7-11eb-324d-93a5eda41e3e
# ╟─ddc7b166-203a-11eb-0ac2-056c309bb590
# ╟─8e164be0-203b-11eb-293a-6b32dba96133
# ╠═5f83510e-1fc5-11eb-3c40-83a844bf2907
# ╠═51a309bc-2033-11eb-10c0-ed17545df33d
# ╟─f3f0b6f8-2033-11eb-0f9e-d951d791001d
# ╠═f3f4bfba-203f-11eb-142c-a187112744d2
# ╠═2ce4435c-2040-11eb-1670-e39ad4cc690c
# ╠═9e50e56e-5674-11eb-2919-4bde8964bc59
# ╟─8c1dd9de-2040-11eb-35a3-39e8196577a1
# ╟─0f9406ae-2054-11eb-3f0c-1d03ba779f18
# ╟─099e4188-2054-11eb-1e7e-69bdb2b0202c
# ╠═16572944-2045-11eb-38d0-fb20d981e3e9
# ╠═93d3bf98-5684-11eb-361e-c7b0f1a3bb0e
# ╠═18f36e1a-5682-11eb-2826-ab4d4c457af3
# ╠═6f096138-5682-11eb-179f-f973803b9268
# ╟─8cec9318-206e-11eb-3caa-3bad9d788d3c
# ╠═4ae7d968-206f-11eb-1925-c5b5a31e6940
# ╠═8b77d12c-206f-11eb-19a5-f557157ac05e
# ╟─b56df7ea-206f-11eb-18b9-0544cc596ca5
# ╠═554b4ef0-2045-11eb-3218-050ec936f1aa
# ╟─8a4df438-2045-11eb-2798-af592dbbffb5
# ╠═6b1a1520-2045-11eb-1302-0778b4e7a836
# ╟─da93bf52-2045-11eb-2b6f-3f31f57d8653
# ╠═1fa92570-2046-11eb-3541-e5fccf7afabd
# ╠═255c4d08-2046-11eb-1cf0-b32d2cd711bf
# ╠═473e6a46-205b-11eb-2731-7b47b99d86a5
# ╟─917d4550-2058-11eb-174b-65436a4e6cc8
# ╠═c2de4714-2058-11eb-1a0d-f1bf96649844
# ╠═0103406e-2059-11eb-333f-7b5814f21a54
# ╠═579a4d9e-205b-11eb-2352-83beca145211
# ╟─6924f372-2059-11eb-2b98-6f4b65b777f4
# ╠═ece1d052-2060-11eb-2b50-ef46663ed88f
# ╟─43a39b30-2061-11eb-28c5-37cbc4aed84b
# ╠═41ca2316-2062-11eb-3fe6-f33a789fd578
# ╟─0ef8a16e-2063-11eb-1391-6b3b959de541
# ╠═1eb73246-2063-11eb-0f20-416e433b5144
# ╠═abace764-2062-11eb-2949-fd9e86b50f17
# ╟─ebd267a4-2064-11eb-02a7-b93465de6319
# ╠═f6ad5f52-2065-11eb-0be7-f159941e0d89
# ╠═539e68b4-2066-11eb-09c2-b337241c36bc
# ╟─69e7e54c-2067-11eb-21a1-cb72f4b423cd
# ╟─f59d488a-2067-11eb-1d4b-1d9dd97ea39e
# ╟─27d5a290-2068-11eb-06f5-0b16f6ec6cc5
# ╠═b8f9bae8-206a-11eb-1426-031f6fd05fd6
# ╠═ee45d48e-206a-11eb-0edf-2b8b893bb583
# ╟─c794f45e-206b-11eb-33c4-05bc429ba846
# ╟─fcb7e8d2-2071-11eb-020d-7bb1d53f8a6d
# ╠═fac413ec-5684-11eb-3167-d56d967247fc
# ╠═cc6681a2-5688-11eb-1d8e-e1edcb80a20d
# ╠═f01e7438-5688-11eb-2900-8df62e9feb18
# ╠═fae4e938-5688-11eb-3ab8-e3a8a4a61692
# ╠═84d5e7a2-56dc-11eb-0d54-b790c7c714fb
# ╠═7704b9e6-56a0-11eb-0ade-2194a0845f61
# ╠═03f063cc-5689-11eb-3d90-b7a5991332ba
# ╠═a8a9179c-5693-11eb-1758-e7d34722c7fa
# ╠═7e82f15c-5695-11eb-20f1-85763cbfa572
# ╠═bd11306c-5697-11eb-22ed-99ae7703be59
# ╠═d99e52fa-5697-11eb-12fe-ab4ae36ab362
# ╠═5a3d6966-56a0-11eb-0db1-c5a06ea2dd83
# ╠═9d1f1e7e-56a1-11eb-2539-f5da33710a32
# ╠═d7bfc106-56d0-11eb-21b3-3b9036905041
# ╠═9ceba0fa-56d4-11eb-2696-21f4babf699c
# ╠═0da9d3a0-56d5-11eb-08a9-c9e83106dfc3
# ╠═db8ed140-5735-11eb-2078-2f90d15d6325


text:
using BandedMatrices
using BenchmarkTools
using LinearAlgebra
using LinearRationalExpectations
using Test
function make_orthogonal_matrix(n)
    W = randn(n, n)
    w = view(W, :,1)
    A = w/norm(w)
    for i=2:n
        w = view(W, :, i)
        v = (I - A*transpose(A))*w
        v = v/norm(v)
        A = hcat(A, v)
    end
    return A
end

n = 10
A = make_orthogonal_matrix(n)
@test A'*A ≈ I(n)

function make_complex2x2!(V::AbstractMatrix{Float64}, W::AbstractVector{Float64})
    a = W[1]
    b = W[2]
    c = W[3]
    V[1, 1] = V[2, 2] = a
    V[2, 1] = b
    x = abs(c)
    V[1, 2] = (b < 0) ? x : -x
end

function make_stationary_matrix_complex_roots(n::Int64; p::Float64 = 1/5, m::Float64 = 0.99)
    Q = make_orthogonal_matrix(n)
    values = rand(-m:0.001:m, n)
    V = BandedMatrix(Pair(-1, zeros(n - 1)), Pair(0, values), Pair(1, zeros(n - 1)))
    pn = Integer(floor(p*n))
    KK = rand(1:n, pn)
    W = randn(pn, 3)
    mc = maximum(abs.(view(W, :, 3)))
    @views W[:, 3] ./= W[:, 3]/mc
    @views for (i, k) in enumerate(KK)
        if (k == 1 ||
            (k < n && V[k, k-1] == 0 && V[k, k+1] == 0))
            make_complex2x2!(V[k:k+1, k:k+1], W[i, :])
        end
    end
    return Q*V
end

function make_ABΣe!(n, ns, nns, nx)
    state_nbr = ns + nns
    nonstate_nbr = n - nns - ns
    A1a = triu(randn(nns, nns), 1) + I
    A1b = randn(nns, ns)
    A2a = randn(nonstate_nbr, ns)
    A2b = make_stationary_matrix_complex_roots(ns)
    A = hcat(zeros(n, nonstate_nbr),
             vcat(hcat(zeros(nonstate_nbr, nns),
                       A2a),
                  hcat(A1a, A1b),
                  hcat(zeros(ns, nns),
                       A2b)))
    B = rand(n, nx)
    Σe = randn(nx, nx)
    Σe = Σe*transpose(Σe)
    return (A, B, Σe)
end

n = 100
ns = 50
nns = 20
nx = 30
state_nbr = ns + nns
nonstate_nbr = n - nns - ns
backward_indices = nonstate_nbr + 1:n
non_backward_indices = 1:nonstate_nbr

(A, B, Σe) = make_ABΣe!(n, ns, nns, nx)

A1 = A[backward_indices, backward_indices]
A2 = A[non_backward_indices, backward_indices]

B1 = B[backward_indices, :]
B2 = B[non_backward_indices, :]

Σyss = Matrix{Float64}(undef, state_nbr, state_nbr)

@btime wsl = LinearRationalExpectations.LyapdWs(state_nbr)
wsl = LinearRationalExpectations.LyapdWs(state_nbr)
@btime LinearRationalExpectations.extended_lyapd!(Σyss, A1, B1*Σe*transpose(B1), wsl)

Σy = zeros(n, n)

algo = "GS"
forward_indices = collect(1:nonstate_nbr - 1)
current_indices = collect(1:n)
both_indices = Vector{Int64}(undef, 0)
static_indices = [nonstate_nbr]
lre_ws = LinearRationalExpectationsWs(algo,
                                      n,
                                      nx,
                                      0,
                                      forward_indices,
                                      current_indices,
                                      collect(backward_indices),
                                      both_indices,
                                      static_indices)

@btime ws = LinearRationalExpectations.VarianceWs(n,
                                                  state_nbr,
                                                  nx,
                                                  lre_ws)
ws = LinearRationalExpectations.VarianceWs(n,
                                           state_nbr,
                                           nx,
                                           lre_ws)
@btime LinearRationalExpectations.compute_variance!(Σy, A1, A2, B1, B2, Σe, ws)


text:
using Test
using QuantumControlBase
import QuantumControlBase: optimize

@testset "optimize-kwargs" begin

    # test that we can call optimize with kwargs that override the kwargs of
    # the `problem` without permantly changing `problem`

    struct Result
        iter_stop :: Int
        flag :: Bool
    end

    struct DummyObjective <: QuantumControlBase.AbstractControlObjective
        initial_state
        generator
    end

    function optimize_kwargstest(problem)
        return Result(problem.kwargs[:iter_stop], problem.kwargs[:flag])
    end

    QuantumControlBase.optimize(problem, method::Val{:kwargstest}) = optimize_kwargstest(problem)

    problem = ControlProblem(
        objectives=[DummyObjective(nothing, nothing)],
        pulse_options=Dict(),
        tlist=[0.0, 10.0],
        iter_stop=2,
        flag=false
    )

    res = QuantumControlBase.optimize(problem; method=:kwargstest)
    @test res.iter_stop == 2
    @test !res.flag

    res2 = QuantumControlBase.optimize(problem; method=:kwargstest, iter_stop=10, flag=true)
    @test res2.iter_stop == 10
    @test res2.flag
    @test problem.kwargs[:iter_stop] == 2
    @test !problem.kwargs[:flag]

end


text:
using Test,
      RandomQuantum,
      QuantumInfo,
      Cliffords,
      Distributions,
      QuantumTomography,
      Qlab

import SchattenNorms.trnorm

"""
Basic trace-norm comparison of the true state ρ and the estimate ρest
"""
function check_tomo_result(ρ::Array{<:Number,2},
                           ρest::Array{<:Number,2})
    @test isapprox(trnorm(ρ-ρest), 0, atol=3e-2)
end

################################################################################
# Setup tomography problem #####################################################
################################################################################

"""
Create the quantum state tomography problem for a given number of qubits

This setup defaults to the obervables expected in the python QGL package and
assumed in tht rest of the code base.

Returns a random density matrix, the set of observables, the mean values of
those observables and the variance in the case of least squares mean values.

If asymptotic is true, return the true predicted means.  If false, add noise.
"""
function qlab_qst_setup(num_Qubits;num_obs = 6,
                                   asymptotic=false,
                                   n = 10_000)

    if num_obs == 4
        # observe along the Z, X, -Y, -Z basis created in QGL
        obs_ = [(complex(Pauli(2)) + QuantumInfo.eye(2))/2,
               (complex(Pauli(3)) + QuantumInfo.eye(2))/2,
               (complex(-Pauli(1)) + QuantumInfo.eye(2))/2,
               (complex(-Pauli(2)) + QuantumInfo.eye(2))/2]
    elseif num_obs == 6
        # observe along the Z, X, -X, Y, -Y, -Z basis created in QGL
        obs_ = [(complex(Pauli(2)) + QuantumInfo.eye(2))/2,
               (complex(Pauli(3)) + QuantumInfo.eye(2))/2,
               (complex(-Pauli(3)) + QuantumInfo.eye(2))/2,
               (complex(-Pauli(1)) + QuantumInfo.eye(2))/2,
               (complex(Pauli(1)) + QuantumInfo.eye(2))/2,
               (complex(-Pauli(2)) + QuantumInfo.eye(2))/2]
    end

    # expand the state space based on the number of qubits
    obs = Array{ComplexF64,2}[]
    for gates in Base.product([obs_ for _ in 1:num_Qubits]...)
        push!(obs, kron(1, reverse(gates)...))
    end

    # create a random num_Qubit state
    ρ = .98*projector(rand(FubiniStudyPureState(2^num_Qubits))) +
                0.02*rand(HilbertSchmidtMixedState(2^num_Qubits))

    # restrict to least-squares if we only have four observables
    if num_obs == 4
        LSQtomo = LSStateTomo(obs)

        LSQ_means = real(predict(LSQtomo, ρ))

        # if asymptotic is true, add Binomial noise
        LSQ_samples = asymptotic ? LSQ_means[1:num_obs^num_Qubits] : Float64[rand(Distributions.Binomial(n,μ))/n for μ in LSQ_means[1:num_obs^num_Qubits]]
        LSQ_var = asymptotic ? ones(length(LSQ_samples)) : n*(LSQ_samples - LSQ_samples.^2)/(n-1)

        return ρ, obs, LSQ_samples, LSQ_var, []

    elseif num_obs == 6
        LSQtomo = LSStateTomo(obs)
        MLtomo  = MLStateTomo(obs)

        LSQ_means = real(predict(LSQtomo, ρ))
        ML_means = real(predict(MLtomo, ρ))

        # if asymptotic is true, add Binomial noise
        LSQ_samples = asymptotic ? LSQ_means[1:num_obs^num_Qubits] : Float64[rand(Distributions.Binomial(n,μ))/n for μ in LSQ_means[1:num_obs^num_Qubits]]
        LSQ_var = asymptotic ? ones(length(LSQ_samples)) : n*(LSQ_samples - LSQ_samples.^2)/(n-1)

        ML_samples = asymptotic ? ML_means[1:num_obs^num_Qubits] : Float64[rand(Distributions.Binomial(n,μ))/n for μ in ML_means[1:num_obs^num_Qubits]]

        return ρ, obs, LSQ_samples, LSQ_var, ML_samples
    end
end

################################################################################
# Test single-qubit LSQ_QST and ML_QST #########################################
################################################################################
"""
Test single qubit least squares tomography with the manual version of the
analyzeStateTomo function
"""
function do_1q_LSStateTomo(num_obs)

    ρ, obs, LSQ_data, LSQ_var, _ = qlab_qst_setup(1,
                                                  num_obs = num_obs,
                                                  asymptotic=false);

    # did reconstruction work without Qlab.jl??
    LSQ_tomo = LSStateTomo(obs);
    ρest, obj, status = fit(LSQ_tomo, LSQ_data, LSQ_var);
    @test status == :Optimal
    @test trnorm(ρ-ρest) < 3e-2

    # pack the data as the analyzeStateTomo function expects
    data = Dict{String,Dict{String,Array{Any,N} where N}}("q1-main"=>Dict("data"=>LSQ_data, "variance"=>LSQ_var))

    rhoLSQ, _  = Qlab.analyzeStateTomo(data,1,num_obs,nbrCalRepeats=0)
    check_tomo_result(ρ, rhoLSQ)
end

"""
Test single qubit, maximum likelihood by manually creating the set of
observables, data and projector using the QST_ML function
"""
function do_1q_MLStateTomo()
    # Note here we only test with num_obs = 6 because the num_obs = 4 case does
    # not form a POVM
    ρ, obs, _, _, ML_data = qlab_qst_setup(1,
                                           num_obs = 6,
                                           asymptotic=false);

    numQubits = 1
    numAxes = 6
    numCals = 0
    numMeas = 1
    proj = Qlab._create_ml_POVM(numQubits)

    # did reconstruction work without Qlab.jl??
    ML_tomo = MLStateTomo(obs);
    ρest, obj, status = fit(ML_tomo, ML_data);
    @test status == :Optimal
    @test trnorm(ρ-ρest) < 3e-2

    # Map each experiment to the appropriate readout pulse
    # These are integers 1:length(data), each repeated numAxes^nbrQubits times
    measOpMap = repeat(1:numMeas, inner=numAxes^numQubits)
    # These are integers 1:nbrAxes^nbrQubits, unrolled length(data) times
    measPulseMap = repeat(1:numAxes^numQubits, outer=numMeas)
    # Use a helper to get the measurement unitaries.
    measPulseUs = Qlab.tomo_gate_set(numQubits, numAxes);

    rhoML = Qlab.QST_ML(ML_data,
                     measPulseMap,
                     measOpMap,
                     measPulseUs,
                     proj)
    check_tomo_result(ρ, rhoML)
end

"""
Test the creation of the StateTomo object with single qubit least squares.
This test reads in data, builds the object and passes it to the analyzeStateTomo
function.
"""
function test_tomo_obj(num_obs)
    ρ, obs, LSQ_data, LSQ_var, _ = qlab_qst_setup(1,
                                                  num_obs = num_obs,
                                                  asymptotic=false);

    # did reconstruction work without Qlab.jl??
    LSQ_tomo = LSStateTomo(obs);
    ρest, obj, status = fit(LSQ_tomo, LSQ_data, LSQ_var);
    @test status == :Optimal
    @test trnorm(ρ-ρest) < 5e-2

    # pack the data as the analyzeStateTomo function expects
    data = Dict{String,Dict{String,Array{Any,N} where N}}("q1-main"=>Dict("data"=>LSQ_data, "variance"=>LSQ_var))
    desc = Dict{String,Any}("q1-main"=>Any[0])

    tomo = Qlab.StateTomo(data, desc)
    rhoLSQ, _  = Qlab.analyzeStateTomo(tomo)
    check_tomo_result(ρ, rhoLSQ)
end

################################################################################
# Test two-qubit LSQ_QST and ML_QST #########################################
################################################################################

# Note: only LSQ is supported`

function do_2q_LSStateTomo(num_obs)

    ρ, obs, LSQ_data, LSQ_var, _ = qlab_qst_setup(2,num_obs = num_obs,
                                                    asymptotic=false);

    # did reconstruction work without Qlab.jl??
    LSQ_tomo = LSStateTomo(obs);
    ρest, obj, status = fit(LSQ_tomo, LSQ_data, LSQ_var);
    @test status == :Optimal
    @test trnorm(ρ-ρest) < 3e-2

    # pack the data as the analyzeStateTomo function expects
    data = Dict{String,Dict{String,Array{Any,N} where N}}("q1-main"=>Dict("data"=>LSQ_data, "variance"=>LSQ_var))
    desc = Dict{String,Any}("q1-main"=>Any[0])

    rhoLSQ, _  = Qlab.analyzeStateTomo(data,2,num_obs,nbrCalRepeats=0)
    check_tomo_result(ρ, rhoLSQ)
end

# execute the tests
do_1q_LSStateTomo(4)
do_1q_LSStateTomo(6)
do_1q_MLStateTomo()
test_tomo_obj(4)
test_tomo_obj(6)
do_2q_LSStateTomo(6)


text:
using Test
elm = ExtremeLearningMachine(100)

# Testing XOR
x = [1.0 1.0; 0.0 1.0; 0.0 0.0; 1.0 0.0]
y = [0.0, 1.0, 0.0, 1.0]
println("Fitting")
fit!(elm, x, y)
println("predicting")
y_pred = predict(elm, [1 1.; 0 1; 1 1])
@testset "Normal test predictions" begin
    @test y_pred[1] < 0.2
    @test y_pred[2] > 0.8
    @test y_pred[3] < 0.2
end
println("Done")
fitSVD!(elm, x, y)
println("predicting")
y_pred = predict(elm, [1 1.; 0 1; 1 1])
@testset "SVD test predictions" begin
    @test y_pred[1] < 0.2
    @test y_pred[2] > 0.8
    @test y_pred[3] < 0.2
end
println("Done")


text:
function vKITTI_Conv_RGBD(w,x)
    #RGB
    x1 = KnetArray(Array(x)[:,:,1:3,:])
    #Colorized Depth
    x2 = KnetArray(Array(x)[:,:,4:6,:])
    #RGB
    x1 = relu(conv4(w[1],x1;padding=1).+w[2])
    x1 = relu(conv4(w[3],x1;padding=1).+w[4])
    x1 = pool(x1)
    x1 = relu(conv4(w[5],x1;padding=1).+w[6])
    x1 = relu(conv4(w[7],x1;padding=1).+w[8])
    x1 = pool(x1)
    x1 = relu(conv4(w[9],x1;padding=1).+w[10])
    x1 = relu(conv4(w[11],x1;padding=1).+w[12])
    x1 = relu(conv4(w[13],x1;padding=1).+w[14])
    x1 = relu(conv4(w[15],x1;padding=1).+w[16])
    x1 = pool(x1)
    x1 = relu(conv4(w[17],x1;padding=1).+w[18])
    x1 = relu(conv4(w[19],x1;padding=1).+w[20])
    x1 = relu(conv4(w[21],x1;padding=1).+w[22])
    x1 = relu(conv4(w[23],x1;padding=1).+w[24])
    #Colorized Depth
    x2 = relu(conv4(w[25],x2;padding=1).+w[26])
    x2 = relu(conv4(w[27],x2;padding=1).+w[28])
    x2 = pool(x2)
    x2 = relu(conv4(w[29],x2;padding=1).+w[30])
    x2 = relu(conv4(w[31],x2;padding=1).+w[32])
    x2 = pool(x2)
    x2 = relu(conv4(w[33],x2;padding=1).+w[34])
    x2 = relu(conv4(w[35],x2;padding=1).+w[36])
    x2 = relu(conv4(w[37],x2;padding=1).+w[38])
    x2 = relu(conv4(w[39],x2;padding=1).+w[40])
    x2 = pool(x2)
    x2 = relu(conv4(w[41],x2;padding=1).+w[42])
    x2 = relu(conv4(w[43],x2;padding=1).+w[44])
    x2 = relu(conv4(w[45],x2;padding=1).+w[46])
    x2 = relu(conv4(w[47],x2;padding=1).+w[48])
    #Scores
    x3 = x1.+x2
    x3 = relu(conv4(w[49],x3).+w[50])
    x3 = deconv4(w[51],x3;padding=1,stride=6).+w[52]
end

function vKITTI_Conv_RGBD_w()
    w = Array[
        #RGB
        #Conv64
        xavier(Float32,3,3,3,64),        zeros(Float32,1,1,64,1),
        xavier(Float32,3,3,64,64),       zeros(Float32,1,1,64,1),
        #Conv128
        xavier(Float32,3,3,64,128),      zeros(Float32,1,1,128,1),
        xavier(Float32,3,3,128,128),     zeros(Float32,1,1,128,1),
        #Conv256
        xavier(Float32,3,3,128,256),     zeros(Float32,1,1,256,1),
        xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
        xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
        xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
        #Conv512
        xavier(Float32,3,3,256,512),     zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
        #Colorized Depth
        #Conv64
        xavier(Float32,3,3,3,64),        zeros(Float32,1,1,64,1),
        xavier(Float32,3,3,64,64),       zeros(Float32,1,1,64,1),
        #Conv128
        xavier(Float32,3,3,64,128),      zeros(Float32,1,1,128,1),
        xavier(Float32,3,3,128,128),     zeros(Float32,1,1,128,1),
        #Conv256
        xavier(Float32,3,3,128,256),     zeros(Float32,1,1,256,1),
        xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
        xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
        xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
        #Conv512
        xavier(Float32,3,3,256,512),     zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
        #Scores
        xavier(Float32,1,1,512,13),      zeros(Float32,1,1,13,1),
        #Deconv
        bilinear(Float32,32,32,13,13),   zeros(Float32,1,1,13,1),
    ]
    rgb_w = load(pretrainedPath("vKITTI_Conv_RGB.jld"),"w")
    d_w = load(pretrainedPath("vKITTI_Conv_D.jld"),"w")
    w[1:24] = rgb_w[1:24]
    w[25:48] = d_w[1:24]
    return w
end

function vKITTI_Conv_RGBD_loss(w,x,y)
    pixelwiseSoftloss(Array(vKITTI_Conv_RGBD(w,x)),y)
end

max_stars_repo_name: eeyrdkl/MScThesis

max_issues_repo_name: eeyrdkl/MScThesis

max_forks_repo_name: eeyrdkl/MScThesis


text:
module BSONLite

import Base64: base64encode, base64decode
import Base: Vector, show, string, getindex
import Dates: DateTime, UTC, UTM, UNIXEPOCH, value, now
import Random: rand

export AbstractBSONCodec, BSONCodec, JSONCodec
export BSONType, Element, ObjectId, Document, Binary, BSONArray
export read_bson, write_bson

const _oid_context = Dict{Symbol,Any}()

function __init__()
    global _oid_context

    io = IOBuffer()
    write(io, rand(UInt64))
    _oid_context[:rand] = take!(io)[1:5]
    _oid_context[:counter] = rand(UInt32)
end

include("type.jl")
include("decode.jl")
include("encode.jl")
include("read.jl")
include("write.jl")

end


text:
export Circle

##  Circle

# T Must be in an Algebra
immutable Circle{T,V<:Real,TT} <: PeriodicDomain{TT}
	center::T
	radius::V
	orientation::Bool

end

Circle(c::Number,r::Real,o::Bool) = Circle{typeof(c),typeof(r),Complex{typeof(r)}}(c,r,o)
Circle(c::Vec,r::Real,o::Bool) = Circle{typeof(c),typeof(r),typeof(c)}(c,r,o)

Circle{T1,T2,V<:Real}(::Type{T1},c::T2,r::V,o::Bool) = Circle(convert(promote_type(T1,T2,V),c),
															  convert(promote_type(real(T1),real(T2),V),r),o)
Circle{T1<:Number}(::Type{T1},c,r::Bool) = Circle(T1,c,r)
Circle{T1<:Number}(::Type{T1},c,r::Real) = Circle(T1,c,r)
Circle(c,r::Real) = Circle(c,r,true)
Circle(r::Real) = Circle(zero(r),r)
Circle(r::Int) = Circle(Float64,0.,r)
Circle(a::Tuple,r::Real) = Circle(Vec(a...),r)

Circle{V<:Real}(::Type{V}) = Circle(one(V))
Circle() = Circle(1.0)

isambiguous{T<:Number}(d::Circle{T}) = isnan(d.center) && isnan(d.radius)
isambiguous{T<:Vec}(d::Circle{T}) = all(isnan,d.center) && isnan(d.radius)
Base.convert{T<:Number,V<:Real}(::Type{Circle{T,V}},::AnyDomain) = Circle{T,V}(NaN,NaN)
Base.convert{IT<:Circle}(::Type{IT},::AnyDomain) = Circle(NaN,NaN)

function tocanonical{T<:Number}(d::Circle{T},ζ)
    v=mappoint(d,Circle(),ζ)
    atan2(imag(v)-0.0,real(v))  # -0.0 to get branch cut right
end

function tocanonical{T<:Vec}(d::Circle{T},ζ)
    v=mappoint(d,Circle((0.0,0.0),1.0),ζ)
    atan2(v[2]-0.0,v[1])  # -0.0 to get branch cut right
end

fromcanonical{T<:Number,V<:Real}(d::Circle{T,V,Complex{V}},v::AbstractArray) =
	eltype(d)[fromcanonical(d,vk) for vk in v]
fromcanonical{T<:Number,V<:Real}(d::Circle{T,V,Complex{V}},θ) =
	d.radius*exp((d.orientation?1:-1)*1.0im*θ) + d.center
fromcanonicalD{T<:Number}(d::Circle{T},θ) =
	(d.orientation?1:-1)*d.radius*1.0im*exp((d.orientation?1:-1)*1.0im*θ)

fromcanonical{T<:Vec}(d::Circle{T},θ::Number) =
	d.radius*Vec(cos((d.orientation?1:-1)*θ),sin((d.orientation?1:-1)*θ)) + d.center
fromcanonicalD{T<:Vec}(d::Circle{T},θ::Number) =
	d.radius*(d.orientation?1:-1)*Vec(-sin((d.orientation?1:-1)*θ),cos((d.orientation?1:-1)*θ))

Base.in(z,d::Circle) = isapprox(abs(z-d.center),d.radius)

arclength(d::Circle) = 2π*d.radius
complexlength(d::Circle) = (d.orientation?1:-1)*im*arclength(d)  #TODO: why?

==(d::Circle,m::Circle) = d.center == m.center && d.radius == m.radius && d.orientation == m.orientation

mappoint{T<:Vec,V<:Number}(d1::Circle{T},d2::Circle{V},z) =
	mappoint(Circle(complex(d1.center...),d1.radius),d2,z[1]+im*z[2])

mappoint{T<:Number,V<:Vec}(d1::Circle{T},d2::Circle{V},z) =
	mappoint(Circle(Vec(d1.center...),d1.radius),d2,Vec(real(z),imag(z)))

function mappoint(d1::Circle,d2::Circle,z)
   v=(z-d1.center)/d1.radius
   d1.orientation != d2.orientation && (v=1./v)
   v*d2.radius+d2.center
end

Base.reverse(d::Circle) = Circle(d.center,d.radius,!d.orientation)
Base.conj(d::Circle) = Circle(conj(d.center),d.radius,!d.orientation)

for op in (:+,:-,:.+,:.-)
    @eval begin
        $op(c::Number,d::Circle) = Circle($op(c,d.center),d.radius,d.orientation)
        $op(d::Circle,c::Number) = Circle($op(d.center,c),d.radius,d.orientation)
    end
end

for op in (:*,:.*)
    @eval begin
        $op(c::Real,d::Circle) = Circle($op(c,d.center),$op(abs(c),d.radius),sign(c)<0?!d.orientation:d.orientation)
        $op(d::Circle,c::Real) = Circle($op(c,d.center),$op(abs(c),d.radius),sign(c)<0?!d.orientation:d.orientation)
    end
end

/(c::Number,d::Circle) =
	c==1?(d.center==0?Circle(d.center,1/d.radius,!d.orientation):
				Circle(1/d.center,abs(1/(d.center+d.radius)-1/(d.center)),!d.orientation)):
				c*(1/d)


text:
export quantize_norms, splitarray, K2vec, quantize_norms, get_norms_codebook

"Get the codebook of the norms with k-means"
function get_norms_codebook(
  B::Matrix{T1},
  C::Vector{Matrix{T2}}) where {T1<:Integer, T2<:AbstractFloat}

  m, n = size(B)
  d, h = size(C[1])

  # Make sure there are m codebooks
  @assert m == length(C)

  # Reconstruct the approximation and compute its norms
  CB      = reconstruct(B, C)
  dbnorms = sum(CB.^2, dims=1)

  # Quantize the norms with k-means
  dbnormsq = Clustering.kmeans(dbnorms, h)

  norms_codes     = vec(reshape(dbnormsq.assignments, 1, n)[:])
  norms_codebook  = vec(dbnormsq.centers)

  # Add the dbnorms to the codes
  return norms_codes, norms_codebook
end

"Quantize the norms of an encoding"
function quantize_norms(
  B::Matrix{T1},         # m-by-n. Codes
  C::Vector{Matrix{T2}}, # m-long h-by-d each. Codebooks
  cbnorms::Vector{T2}) where {T1<:Integer, T2<:AbstractFloat} # h-long. norms codebooks

  CB = reconstruct(B,C)
  d, n = size( CB )
  _, h = size( C[1] )

  dbnormsB   = Vector{T1}(undef, n)
  dbnormsX   = Vector{T2}(undef, n)
  dists2norm = Vector{T2}(undef, h)

  @inbounds for i = 1:n
    ithnorm::T2 = zero(T2)

    @simd for j = 1:d
      ithnorm += CB[j,i].^2
    end
    dbnormsX[i] = ithnorm

    for j = 1:h
      dists2norm[j] = (ithnorm - cbnorms[j]).^2
    end

    # Store of the index of the minimum distance
    _, dbnormsB[i] = findmin( dists2norm )
  end

  return dbnormsB, dbnormsX
end

# Compute dot products between codebooks
function comp_prod_C(C::Vector{Matrix{T}}) where T <: AbstractFloat  # Codebooks
  m   = length(C)
  dps = Matrix{Matrix{T}}(m, m) # dot products
  for i = 1:m, j = (i+1):m
    dps[i,j] = 2 * C[i]' * C[j]
  end
  return dps
end

############################################
## Functions that support codebook update ##
############################################

# Creates a sparse matrix out of codes
function sparsify_codes(
  B::Matrix{T}, # m-by-n matrix. Codes to sparsify
  h::Integer,
  T2=Float32) where T <: Integer # Number of entries per codebook

  m, n   = size(B)
  ncodes = length(B)

  # Storing the indices of the sparse matrix.
  I = zeros(Int32, ncodes) # row indices
  J = zeros(Int32, ncodes) # column indices

  for i = 1:m
    I[ (i-1)*n+1 : i*n ] = 1:n
    J[ (i-1)*n+1 : i*n ] = vec(B[i,:]) .+ (i-1) * h
  end

  C = sparse(I, J, ones(T2, ncodes), n, m * h)

  return C
end

# Transform the output of LSQR/SPGL1 into a vector of matrices again
function K2vec(
  K::Union{SharedMatrix{T}, Matrix{T}}, # d-by-(h*m) matrix with the new codebooks
  m::Integer, # Number of codebooks
  h::Integer) where T <: AbstractFloat # Number of elements in each codebook

  @assert size(K,2) == m*h

  C = Vector{Matrix{T}}(undef, m)

  subdims = splitarray(1:(h*m), m)
  for i = 1:m
    C[i] = sdata(K[ :, subdims[i] ])
  end

  return C
end

###########################################
## Functions to get unaries and binaries ##
###########################################

# Get unaries terms
function get_unaries(
  X::Matrix{T},         # data to get unaries from
  C::Vector{Matrix{T}}, # codebooks
  V::Bool = false) where T <: AbstractFloat

  d, n = size( X )
  m    = length(C)
  _, h = size( C[1] )

  unaries = Vector{Matrix{T}}(undef, m)
  #if V print("Computing unaries... "); st=time(); end

  @inbounds for i = 1:m
    unaries[i] = -2 * C[i]' * X
    sci        = diag( C[i]' * C[i] )
    #unaries[i] = broadcast(+, sci, unaries[i])

    ui = unaries[i]
    for j = 1:n
      @simd for k = 1:h
        ui[k, j] += sci[k]
      end
    end
  end

  #if V @printf("done in %.3f seconds.\n", time()-st); end
  return unaries
end

# Computes all the binaries for a set of codebooks
function get_binaries(
  C::Vector{Matrix{T}}) where T <: AbstractFloat # codebooks

  m = length( C )

  ncbi     = sum(1:(m-1)) # The number of pairwise interactions
  binaries = Vector{Matrix{T}}(undef, ncbi) # We'll store the binaries here
  cbi      = zeros(Int32, 2, ncbi) # Store the pairwise indices here

  idx=1
  for i = 1:m
    for j = (i+1):m
      binaries[idx] = 2 * C[i]' * C[j]
      cbi[:,idx]    = [i;j]
      idx = idx + 1
    end
  end

  return binaries, cbi
end

# Split an array in equal parts
# if the number of x is a multiple of nparts, then each parts has the same number of
# elements.
# if that's not the case, then the n rest elements, e.t. the rest of the division, will be
# added to the first n parts. So the first n parts contain each 1 item more than the the
# parts from n + 1 to nparts
function splitarray(
  x::Union{UnitRange, Vector}, # Vector to split in equal parts
  nparts::Integer)             # Number of parts to split the vector on

  n       = length(x)
  perpart = div(n, nparts)
  xtra    = mod(n, nparts)
  out     = Array{Any}(undef, nparts)

  #fills the parts, which have 1 item more than the other parts,
  #so these parts have size perpart + 1
  glidx = 1
  for i = 1:xtra
    out[i] = x[glidx : (glidx+perpart)]
    glidx  = glidx+perpart+1
  end

  #fills the parts with the normal size
  for i = (xtra+1):nparts
    out[i] = x[glidx : (glidx+perpart-1)]
    glidx  = glidx+perpart
  end

  return out
end


text:
using Pkg
println("Activating environment in $(pwd())...")
pkg"activate ."

println("Installing packages..."); flush(stdout);
pkg"instantiate"
#pkg"precompile"

if v"1.5" <= VERSION < v"1.6"
    println("Installing IJulia.jl into global environment...")
    pkg"activate @v1.5"
    pkg"add IJulia"
else
    @warn "You aren't using Julia 1.5! Didn't install IJulia into global environment."
end

println("Installing WebIOs jupyter extension...")
pkg"activate ."
using WebIO
WebIO.install_jupyter_nbextension()

println("Done!")


text:
module Baysor

using Distributions
using LinearAlgebra
using ProgressMeter
using Statistics

import Distributions.pdf
import Distributions.logpdf
import Statistics.rand

export
    BmmData, bmm!,
    build_molecule_graph, load_df, initial_distribution_arr,
    boundary_polygons,
    neighborhood_count_matrix, gene_composition_transformation, gene_composition_colors,
    append_confidence!,
    plot_molecules, plot_molecules!, plot_expression_vectors,
    cluster_molecules_on_mrf

include("utils/utils.jl")
include("utils/logging.jl")
include("utils/convex_hull.jl")
include("utils/spatial.jl")

include("distributions/MvNormal.jl")
include("distributions/CategoricalSmoothed.jl")

include("models/InitialParams.jl")
include("models/Component.jl")
include("models/BmmData.jl")

include("data_processing/prior_segmentation.jl")
include("data_processing/triangulation.jl")
include("data_processing/umap_wrappers.jl")
include("data_processing/expression_local_vectors.jl")
include("data_processing/noise_estimation.jl")
include("data_processing/initialization.jl")
include("data_processing/plots.jl")
include("data_processing/boundary_estimation.jl")
include("data_processing/preview_diagnostics.jl")

include("bmm_algorithm/molecule_clustering.jl")
include("bmm_algorithm/tracing.jl")
include("bmm_algorithm/distribution_samplers.jl")
include("bmm_algorithm/history_analysis.jl")
include("bmm_algorithm/bmm_algorithm.jl")

include("cli_wrappers.jl")

end # module

max_stars_repo_name: YinuoJin/Baysor

max_issues_repo_name: YinuoJin/Baysor

max_forks_repo_name: YinuoJin/Baysor


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Deterministic equivalent evaluation #
# =================================== #
function evaluate_decision(structure::DeterministicEquivalent{2}, decision::AbstractVector)
    # Update decisions (checks handled by structure model)
    take_decisions!(structure.model, all_decision_variables(structure.model, 1), decision)
    # Optimize model
    optimize!(structure.model)
    # Switch on return status
    status = termination_status(structure.model)
    result = if status in AcceptableTermination
        result = objective_value(structure.model)
    else
        result = if status == MOI.INFEASIBLE
            result = objective_sense(structure.model) == MOI.MAX_SENSE ? -Inf : Inf
        elseif status == MOI.DUAL_INFEASIBLE
            result = objective_sense(structure.model) == MOI.MAX_SENSE ? Inf : -Inf
        else
            error("Deterministically equivalent model could not be solved, returned status: $status")
        end
    end
    # Revert back to untaken decisions
    untake_decisions!(structure.model, all_decision_variables(structure.model, 1))
    # Return evaluation result
    return result
end

function statistically_evaluate_decision(structure::DeterministicEquivalent{2}, decision::AbstractVector)
    # Update decisions (checks handled by structure model)
    take_decisions!(structure.model, all_decision_variables(structure.model, 1), decision)
    # Optimize model
    optimize!(structure.model)
    # Get sense-correted objective value
    status = termination_status(structure.model)
    Q̂ = if status in AcceptableTermination
        Q̂ = objective_value(structure.model)
    else
        Q̂ = if status == MOI.INFEASIBLE
            return objective_sense(structure.model) == MOI.MAX_SENSE ? (-Inf, 0) : (Inf, 0)
        elseif status == MOI.DUAL_INFEASIBLE
            return objective_sense(structure.model) == MOI.MAX_SENSE ? (Inf, 0) : (-Inf, 0)
        else
            error("Deterministically equivalent model could not be solved, returned status: $status")
        end
    end
    # Calculate subobjectives
    N = num_scenarios(structure)
    Q = Vector{Float64}(undef, N)
    obj_sense = objective_sense(structure.model)
    for (i, sub_objective) in enumerate(structure.decisions.stage_objectives[2])
        (sub_sense, sub_obj) = sub_objective
        Qᵢ = MOIU.eval_variables(sub_obj) do idx
            return MOI.get(backend(structure.model), MOI.VariablePrimal(), idx)
        end
        if obj_sense == sub_sense
            Q[i] = Qᵢ
        else
            Q[i] = -Qᵢ
        end
    end
    probabilities = map(1:num_scenarios(structure, 2)) do i
        probability(structure, 2, i)
    end
    weights = ProbabilityWeights(probabilities)
    σ = std(Q, weights, corrected = true)
    return Q̂, σ
end


text:
using JuMP
using LightGraphs, SimpleWeightedGraphs

function multi_commodity_flow!(model,
                                graph :: SimpleWeightedGraph,
                                k :: Int;
                                generate_timeout_sec::Float64 = Inf64,
                                isdebug::Bool=false,
                                total_memory::Union{Int,UInt,Float64}=Base.Sys.total_memory(),
                                memory_limit_ratio::Union{Rational{Int},Rational{UInt},Float64}=1//20)
    n::Int = nv(graph)
    nes::Int = ne(graph)
    es::Vector{SimpleWeightedEdge} = collect(edges(graph))
    sort!(es, by=dst)
    sort!(es, by=src, alg=MergeSort)

    progress::Progress = Progress(
        (n-1) * 4               # z[2:n], Bin, >= 0, <= 1
        + 2 * 3 * (n-1) * nes   # f[m,i,j] f[m,j,i], <= 0, >= 1
        + (n-1)                 # fix(y[i,alpha])
        + (n-1) * (n-1)         # fix(f[m,i,alpha])
        + 1                     # one_real_root
        + nes * (n-1) * 2       # f[m,i,j] <= y[i,j], f[m,j,i] <= y[j,i]
        + n-1                   # sum(f[m,m,j])- sum(f[m,j,m]) == z[m]
        + n-1                   # sum(f[m,alpha,j]) == z[m]
        + (n-1) * (n-1) - n     # sum(f[m,i,v]) - sum(f[m,v,j]) == 0
        + 1                     # treesize
        ,
                                    dt=0.1,
                                    barglyphs=BarGlyphs('[', '#', ['-','~','+','*','=','>','#'], ' ', ']'),
                                    enabled=isdebug)
    generate_time::Float64 = 0

    generate_time += @elapsed begin
        @variable(model, z[2:n], Bin, lower_bound=0, upper_bound=1)
    end
    ProgressMeter.update!(progress, (n-1) * 4)
    f::Dict{Tuple{Int, Int, Int},VariableRef} = Dict()
    for m in 2:n
        for e in es
            i::Int = src(e)
            j::Int = dst(e)
            generate_time += @elapsed begin
                f[m, i, j] = @variable(model, base_name="f[$m,$i,$j]", lower_bound=0, upper_bound=1)
                ProgressMeter.update!(progress, 3)
                f[m, j, i] = @variable(model, base_name="f[$m,$j,$i]", lower_bound=0, upper_bound=1)
                ProgressMeter.update!(progress, 3)
            end
            if (generate_timeout_sec != Inf64 && generate_time > generate_timeout_sec) || (Base.Sys.free_memory()//total_memory <= memory_limit_ratio)
                cancel(progress)
                return
            end
        end
    end

    alpha::Int = 1

    if (generate_timeout_sec != Inf64 && generate_time > generate_timeout_sec) || (Base.Sys.free_memory()//total_memory <= memory_limit_ratio)
        cancel(progress)
        return
    end

    for i in 2:n
        if has_edge(graph, i, alpha)
            generate_time += @elapsed begin
                fix(variable_by_name(model, "y[$i,$alpha]"), 0, force=true) # Constraint (24)
            end
            next!(progress)
            for m in 2:n
                generate_time += @elapsed begin
                    fix(f[m, i, alpha], 0, force=true) # Constraint (25)
                end
                next!(progress)
                if (generate_timeout_sec != Inf64 && generate_time > generate_timeout_sec) || (Base.Sys.free_memory()//total_memory <= memory_limit_ratio)
                    cancel(progress)
                    return
                end
            end
        end
    end

    #@constraint(model, source_flow, sum(sum(f[m, alpha, j] for j in 2:n if has_edge(graph, alpha, j)) for m in 2:n) == k)
    generate_time += @elapsed @constraint(model, one_real_root, sum(variable_by_name(model, "y[$alpha,$j]") for j in 2:n if has_edge(graph, alpha, j)) == 1) # Constraint (26)
    next!(progress)
    for e in es
        i :: Int = src(e)
        j :: Int = dst(e)
        for m in 2:n
            generate_time += @elapsed begin
                @constraint(model, f[m, i, j] <= variable_by_name(model, "y[$i,$j]")) # Constraint (27)
            end
            next!(progress)
            generate_time += @elapsed begin
                @constraint(model, f[m, j, i] <= variable_by_name(model, "y[$j,$i]")) # Constraint (27)
            end
            next!(progress)
            if (generate_timeout_sec != Inf64 && generate_time > generate_timeout_sec) || (Base.Sys.free_memory()//total_memory <= memory_limit_ratio)
                cancel(progress)
                return
            end
        end
    end
    for m in 2:n
        generate_time += @elapsed begin
            @constraint(model, sum(f[m, i, m] for i in 1:n if has_edge(graph, i, m))
                            - sum(f[m, m, j] for j in 2:n if has_edge(graph, m, j)) == z[m]) # Constraint (28)
        end
        next!(progress)
        if (generate_timeout_sec != Inf64 && generate_time > generate_timeout_sec) || (Base.Sys.free_memory()//total_memory <= memory_limit_ratio)
            cancel(progress)
            return
        end
        generate_time += @elapsed begin
            @constraint(model, sum(f[m, alpha, j] for j in 2:n if has_edge(graph, alpha, j)) == z[m]) # Constraint (29)
        end
        next!(progress)
        if (generate_timeout_sec != Inf64 && generate_time > generate_timeout_sec) || (Base.Sys.free_memory()//total_memory <= memory_limit_ratio)
            cancel(progress)
            return
        end
        for v in 2:n
            if m != v
                generate_time += @elapsed begin
                    @constraint(model, sum(f[m, i, v] for i in 1:n if has_edge(graph, i, v))
                                    - sum(f[m, v, j] for j in 2:n if has_edge(graph, v, j)) == 0) # Constraint (30)
                end
                next!(progress)
            end
            if (generate_timeout_sec != Inf64 && generate_time > generate_timeout_sec) || (Base.Sys.free_memory()//total_memory <= memory_limit_ratio)
                cancel(progress)
                return
            end
        end
    end

    @constraint(model, treesize, sum(z[i] for i in 2:n) == k) # Constraint (31)
    next!(progress)
    finish!(progress)
end

max_stars_repo_name: Entze/jkMST

max_issues_repo_name: Entze/jkMST

max_forks_repo_name: Entze/jkMST


text:
module TestUtils

using CSV
using DataFrames
using HTTP
using OrderedCollections

financial_url = "https://raw.githubusercontent.com/scikit-learn/examples-data/master/financial-data/"

financial_dict = OrderedDict(
    "TOT" => "Total",
    "XOM" => "Exxon",
    "CVX" => "Chevron",
    "COP" => "ConocoPhillips",
    "VLO" => "Valero Energy",
    "MSFT" => "Microsoft",
    "IBM" => "IBM",
    "TWX" => "Time Warner",
    "CMCSA" => "Comcast",
    "CVC" => "Cablevision",
    "YHOO" => "Yahoo",
    "DELL" => "Dell",
    "HPQ" => "HP",
    "AMZN" => "Amazon",
    "TM" => "Toyota",
    "CAJ" => "Canon",
    "SNE" => "Sony",
    "F" => "Ford",
    "HMC" => "Honda",
    "NAV" => "Navistar",
    "NOC" => "Northrop Grumman",
    "BA" => "Boeing",
    "KO" => "Coca Cola",
    "MMM" => "3M",
    "MCD" => "McDonald\"s",
    "PEP" => "Pepsi",
    "K" => "Kellogg",
    "UN" => "Unilever",
    "MAR" => "Marriott",
    "PG" => "Procter Gamble",
    "CL" => "Colgate-Palmolive",
    "GE" => "General Electrics",
    "WFC" => "Wells Fargo",
    "JPM" => "JPMorgan Chase",
    "AIG" => "AIG",
    "AXP" => "American express",
    "BAC" => "Bank of America",
    "GS" => "Goldman Sachs",
    "AAPL" => "Apple",
    "SAP" => "SAP",
    "CSCO" => "Cisco",
    "TXN" => "Texas Instruments",
    "XRX" => "Xerox",
    "WMT" => "Wal-Mart",
    "HD" => "Home Depot",
    "GSK" => "GlaxoSmithKline",
    "PFE" => "Pfizer",
    "SNY" => "Sanofi-Aventis",
    "NVS" => "Novartis",
    "KMB" => "Kimberly-Clark",
    "R" => "Ryder",
    "GD" => "General Dynamics",
    "RTN" => "Raytheon",
    "CVS" => "CVS",
    "CAT" => "Caterpillar",
    "DD" => "DuPont de Nemours"
)

function get_financial_data(; D::Int=56, N::Int=1258)
    financial_data = DataFrame()
    for key in collect(keys(financial_dict))[1:D]
        financial_data[!, key] = DataFrame(CSV.File(HTTP.get(financial_url*key*".csv").body))[!, :close]
    end
    return financial_data[1:N, :]
end

end

max_stars_repo_name: molet/DiffDRO

max_issues_repo_name: molet/DiffDRO

max_forks_repo_name: molet/DiffDRO


text:
# using LazySets

include("necessary_check.jl")
include("sufficient_check.jl")
include("refinement.jl")

function analyze_piece(dyn1, dyn2, x0::AbstractVector{N}, P0::LazySet,
                       ε::Number, time, m::Integer, ζ::Number, compute_ua::Bool;
                       log::Bool=false) where {N}
    dim = length(x0)
    T = time[2] - time[1]

    # necessary check
    P1_oa = check_necessary(dyn1, dyn2, x0, P0, ε, T, m; log=log)
    if !compute_ua
        return P1_oa
    elseif P1_oa isa EmptySet
        if log
            println("no solution exists")
        end
        return EmptySet{N}(dim), EmptySet{N}(dim), Vector{Vector{Float64}}(), x0
    end

    # invert system
    dyn1⁻ = -dyn1
    dyn2⁻ = -dyn2
    x1 = successor(dyn1, T, x0)

    # check sufficient criterion for all corners
    vertices = vertices_list(P1_oa)
    vertices_sat, vertices_unsat =
        check_sufficient_points(dyn1⁻, dyn2⁻, x1, vertices, time, ε; log=log)

    # refinement
    refine!(vertices_sat, P1_oa, vertices_unsat, ζ, dyn1⁻, dyn2⁻, x1, time, ε;
            log=log)

    if isempty(vertices_sat)
        P1_ua = EmptySet{N}(dim)
    else
        P1_ua = VPolytope(vertices_sat)
    end
    return P1_ua, P1_oa, vertices_unsat, x1
end

max_stars_repo_name: HySynth/HySynth

max_issues_repo_name: HySynth/HySynth

max_forks_repo_name: HySynth/HySynth


text:
struct WilsonParam <: EoSParam
    g::PairParam{Float64}
    Tc::SingleParam{Float64}
    Pc::SingleParam{Float64}
    ZRA::SingleParam{Float64}
    Mw::SingleParam{Float64}
end

abstract type WilsonModel <: ActivityModel end

struct Wilson{c<:EoSModel} <: WilsonModel
    components::Array{String,1}
    icomponents::UnitRange{Int}
    params::WilsonParam
    puremodel::Vector{c}
    absolutetolerance::Float64
    references::Array{String,1}
end

has_sites(::Type{<:WilsonModel}) = false
has_groups(::Type{<:WilsonModel}) = false
built_by_macro(::Type{<:WilsonModel}) = false

function Base.show(io::IO, mime::MIME"text/plain", model::Wilson)
    return eosshow(io, mime, model)
end

function Base.show(io::IO, model::Wilson)
    return eosshow(io, model)
end

Base.length(model::Wilson) = Base.length(model.icomponents)

molecular_weight(model::Wilson,z=SA[1.0]) = comp_molecular_weight(mw(model),z)

export Wilson

function Wilson(components::Vector{String}; puremodel=PR,
    userlocations=String[],
     verbose=false)
    params = getparams(components, ["properties/critical.csv", "properties/molarmass.csv","Activity/Wilson/Wilson_unlike.csv"]; userlocations=userlocations, asymmetricparams=["g"], ignore_missing_singleparams=["g"], verbose=verbose)
    g  = params["g"]
    Tc        = params["Tc"]
    pc        = params["pc"]
    Mw        = params["Mw"]
    ZRA       = SingleParam(params["w"],"acentric factor")
    ZRA.values .*= -0.08775
    ZRA.values .+= 0.29056
    icomponents = 1:length(components)

    init_puremodel = [puremodel([components[i]]) for i in icomponents]
    packagedparams = WilsonParam(g,Tc,pc,ZRA,Mw)
    references = String[]
    model = Wilson(components,icomponents,packagedparams,init_puremodel,1e-12,references)
    return model
end

function activity_coefficient(model::WilsonModel,p,T,z)
    ZRA = model.params.ZRA.values
    Tc  = model.params.Tc.values
    Pc  = model.params.Pc.values

    Tr  = T ./ Tc
    V =  @. (R̄ *Tc/Pc)*ZRA^(1 + (1-Tr)^2/7)
    Λ = (V' ./ V) .*exp.(-model.params.g.values/R̄/T)
    x = z ./ sum(z)
    lnγ = 1 .- log.(sum(x[i]*Λ[:,i] for i ∈ @comps)) .-sum(x[j] .*Λ[j,:] ./(sum(x[i]*Λ[j,i] for i ∈ @comps)) for j ∈ @comps)
    return exp.(lnγ)
end

is_splittable(::Wilson) = true


text:
module LyceumDocs

using Lyceum, Pkg, Documenter, Literate, DevTools, Markdown, YAML
using LyceumBase
using LyceumMuJoCoViz
using Shapes

const AbsStr = AbstractString
const TupleN{T,N} = NTuple{N,T}

const REPO_DIR = normpath(joinpath(@__DIR__, ".."))
const DOCS_DIR = joinpath(REPO_DIR, "docs")

const SRC_DIR = joinpath(DOCS_DIR, "src")
const STAGING_DIR = joinpath(REPO_DIR, "staging")
const BUILD_DIR = joinpath(REPO_DIR, "build")

const ASSETS_DIR = joinpath(DOCS_DIR, "assets")
const EXAMPLE_DIR = joinpath(DOCS_DIR, "assets/LyceumExamples")

const PATHS = begin
    src = "src"
    static = "static"
    (
        src = src,
        static = static,
        script = joinpath(static, "scripts"),
        notebook = joinpath(static, "notebooks"),
        examples_tarfile = joinpath(static, "examples.tar.gz"),
    )
end

include("document.jl")
include("utils.jl")
include("processors.jl")
include("examples.jl")
include("package_definition.jl")

function make(; clean::Bool = false, skipliterate::Bool = false, config = Dict())
    if isdir(STAGING_DIR) && clean
        @info "Cleaning staging dir ($(STAGING_DIR))"
        println()
        rm(STAGING_DIR, recursive = true, force = true)
    elseif isdir(STAGING_DIR)
        error("$STAGING_DIR exists but clean was false")
    end

    # TODO bundle_examples() fails unless these exist, which only happens if there is
    # at least one example
    mkpath(joinpath(STAGING_DIR, PATHS.script))
    mkpath(joinpath(STAGING_DIR, PATHS.notebook))

    @info "Building Source Tree"
    rootgroup = group(SRC_DIR, skipliterate = skipliterate)
    isempty(rootgroup.children) && error("No source files found in $(SRC_DIR)")

    @info "Processing Files"
    println()
    # disable execution of Jupyter cells for Literate files
    config = Dict("execute" => false)
    process(rootgroup, config = config)

    pages = build_pages(rootgroup)
    isempty(pages) && error("No generated source files found in $(STAGING_DIR)")
    @info "Pages Index:"
    println()
    print_pages(pages)

    @info "Bundling examples"
    println()
    r, status, _, str = Literate.Documenter.withoutput() do
        bundle_examples()
    end
    status || throw(r)

    cp(joinpath(DOCS_DIR, "assets"), joinpath(STAGING_DIR, "assets"))

    @info "Generating Docs"
    println()
    makedocs(;
        # modules = [Lyceum, Lyceum.LYCEUM_PACKAGES...], # TODO
        modules = [LyceumBase, LyceumAI, LyceumMuJoCo, LyceumMuJoCoViz, Shapes],
        format = Documenter.HTML(
            canonical = "https://docs.lyceum.ml/dev/",
            prettyurls = !islocalbuild(),
            assets = ["assets/custom.css"],
        ),
        pages = pages,
        sitename = "Lyceum",
        authors = "Colin Summers",
        strict = false, #!islocalbuild(),  # TODO

        # source/build are specified relative to root
        root = DOCS_DIR,
        build = relpath(BUILD_DIR, DOCS_DIR),
        source = relpath(STAGING_DIR, DOCS_DIR),
    )

    if !islocalbuild()
        @info "Deploying docs"
        deploydocs(
            repo = "github.com/Lyceum/LyceumDocs.jl.git",
            push_preview = true,
            root = DOCS_DIR,
            target = relpath(BUILD_DIR, DOCS_DIR),
        )
    else
        @info "Not deploying docs (local build detected)"
    end
end

end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.Frames (ModiaMath/Frames/_module.jl)
#

"""
    const ModiaMath.Quaternion = SVector{4,Float64}

Describes the rotation from a frame 1 into a frame 2 with a quaternion vector.
If `e` is the (normalized) axis of rotation to rotate frame 1 into frame 2
(either resolved in frame 1 or frame 2) and `angle`
is the rotation angle for this rotation then the quaternion vector
`q::ModiaMath.Quaternions` is defined as:

```julia
q = [e*sin(angle/2),
       cos(angle/2]
```
"""
const Quaternion = SVector{4,Float64}

"""
    const ModiaMath.NullQuaternion = Quaternion(0,0,0,1)

Constant Quaternion vector of a null rotation (= no rotation from frame 1 to frame 2)
"""
NullQuaternion = Quaternion(0.0, 0.0, 0.0, 1.0)

"""
    ModiaMath.assertQuaternion(q::AbstractVector)

Assert that vector `q` has the properties of a `Quaternion` vector
(has 4 elements, `norm(q) = 1`)
"""
function assertQuaternion(q::AbstractVector)
    @assert(length(q) == 4)
    @assert(abs(norm(q) - 1.0) <= 1e-10)
end

"""
    R = ModiaMath.from_q(q::ModiaMath.Quaternion)

Return RotationMatrix `R` from Quaternion `q`.
"""
from_q(q::Quaternion) = RotationMatrix(2.0 * (q[1] * q[1] + q[4] * q[4]) - 1.0,
                                       2.0 * (q[2] * q[1] - q[3] * q[4]),
                                       2.0 * (q[3] * q[1] + q[2] * q[4]),
                                       2.0 * (q[1] * q[2] + q[3] * q[4]),
                                       2.0 * (q[2] * q[2] + q[4] * q[4]) - 1.0,
                                       2.0 * (q[3] * q[2] - q[1] * q[4]),
                                       2.0 * (q[1] * q[3] - q[2] * q[4]),
                                       2.0 * (q[2] * q[3] + q[1] * q[4]),
                                       2.0 * (q[3] * q[3] + q[4] * q[4]) - 1.0)

const p4limit = 0.1
const c4limit = 4.0 * p4limit * p4limit

"""
    q = ModiaMath.from_R(R::ModiaMath.RotationMatrix;
                         q_guess = NullQuaternion)

Return `Quaternion q` from `RotationMatrix R`.

From the two possible solutions `q` the one is returned that is closer
to `q_guess` (note, `q` and `-q` define the same rotation).
"""
function from_R(R::RotationMatrix;
                q_guess::Quaternion=NullQuaternion)::Quaternion
   #=
   Below it is guaranteed that c1>=0, c2>=0, c3>=0, c4>=0 and
   that not all of them can be zero at the same time
   (e.g., if 3 of them are zero, the 4th variable is 1).
   Since the sqrt(..) has to be performed on one of these variables,
   it is applied on a variable which is far enough from zero.
   This guarantees that the sqrt(..) is never taken near zero
   and therefore the derivative of sqrt(..) can never be infinity.
   There is an ambiguity for quaternions, since q and -q
   lead to the same RotationMatrix. This ambiguity
   is resolved here by selecting the q that is closer to
   the input argument q_guess.
   =#
    c1::Float64 = 1.0 + R[1,1] - R[2,2] - R[3,3]
    c2::Float64 = 1.0 + R[2,2] - R[1,1] - R[3,3]
    c3::Float64 = 1.0 + R[3,3] - R[1,1] - R[2,2]
    c4::Float64 = 1.0 + R[1,1] + R[2,2] + R[3,3]
    paux::Float64  = 0.0
    paux4::Float64 = 0.0

    if c4 > c4limit || (c4 > c1 && c4 > c2 && c4 > c3)
        paux  = sqrt(c4) / 2
        paux4 = 4 * paux
        q = Quaternion((R[2,3] - R[3,2]) / paux4,
                       (R[3,1] - R[1,3]) / paux4,
                       (R[1,2] - R[2,1]) / paux4,
                       paux)

    elseif c1 > c2 && c1 > c3 && c1 > c4
        paux  = sqrt(c1) / 2
        paux4 = 4 * paux
        q = Quaternion(paux,
                    (R[1,2] + R[2,1]) / paux4,
                    (R[1,3] + R[3,1]) / paux4,
                    (R[2,3] - R[3,2]) / paux4)

    elseif c2 > c1 && c2 > c3 && c2 > c4
        paux  = sqrt(c2) / 2
        paux4 = 4 * paux
        q = Quaternion((R[1,2] + R[2,1]) / paux4,
                      paux,
                     (R[2,3] + R[3,2]) / paux4,
                     (R[3,1] - R[1,3]) / paux4)

    else
        paux  = sqrt(c3) / 2
        paux4 = 4 * paux
        q = Quaternion((R[1,3] + R[3,1]) / paux4,
                     (R[2,3] + R[3,2]) / paux4,
                      paux,
                     (R[1,2] - R[2,1]) / paux4)
    end

    return dot(q, q_guess) >= 0 ? q : -q
end

from_R(R::AbstractMatrix; q_guess::AbstractVector=NullQuaternion)::Quaternion =
    from_R(RotationMatrix(R); q_guess=Quaternion(q_guess))

"""
    q = ModiaMath.qrot1(angle; q_guess = NullQuaternion)

Return Quaternion `q` that rotates with angle `angle` along the x-axis of frame 1.

From the two possible solutions `q` the one is returned that is closer
to `q_guess` (note, `q` and `-q` define the same rotation).
"""
@inline function qrot1(angle::Number; q_guess::Quaternion=NullQuaternion)::Quaternion
    q = Quaternion(sin(angle / 2), 0.0, 0.0, cos(angle / 2))

    return dot(q, q_guess) >= 0 ? q : -q
end

"""
    q = ModiaMath.qrot2(angle; q_guess = NullQuaternion)

Return Quaternion `q` that rotates with angle `angle` along the y-axis of frame 1.

From the two possible solutions `q` the one is returned that is closer
to `q_guess` (note, `q` and `-q` define the same rotation).
"""
@inline function qrot2(angle::Number; q_guess::Quaternion=NullQuaternion)::Quaternion
    q = Quaternion(0.0, sin(angle / 2), 0.0, cos(angle / 2))

    return dot(q, q_guess) >= 0 ? q : -q
end

"""
    q = ModiaMath.qrot3(angle; q_guess = NullQuaternion)

Return Quaternion `q` that rotates with angle `angle` along the z-axis of frame 1.

From the two possible solutions `q` the one is returned that is closer
to `q_guess` (note, `q` and `-q` define the same rotation).
"""
@inline function qrot3(angle::Number; q_guess::Quaternion=NullQuaternion)::Quaternion
    q = Quaternion(0.0, 0.0, sin(angle / 2), cos(angle / 2))

    return dot(q, q_guess) >= 0 ? q : -q
end

absoluteRotation(q1::Quaternion, q_rel::Quaternion)::Quaternion =
     (@SMatrix [ q_rel[4]  q_rel[3] -q_rel[2] q_rel[1];
                -q_rel[3]  q_rel[4]  q_rel[1] q_rel[2];
                 q_rel[2] -q_rel[1]  q_rel[4] q_rel[3];
                -q_rel[1] -q_rel[2] -q_rel[3] q_rel[4]]) * q1

"""
    q = ModiaMath.qrot123(angle1, angle2, angle3)

Return Quaternion `q` by rotating with angle1 along the x-axis of frame 1,
then with angle2 along the y-axis of this frame and then with angle3 along
the z-axis of this frame.

From the two possible solutions `q` the one is returned that is closer
to `q_guess` (note, `q` and `-q` define the same rotation).
"""
qrot123(angle1::Number, angle2::Number, angle3::Number)::Quaternion = absoluteRotation(absoluteRotation(qrot1(angle1), qrot2(angle2)), qrot3(angle3))

"""
    q = ModiaMath.qrot_e(e, angle; q_guess = NullQuaternion)

Return Quaternion `q` that rotates with angle `angle` along unit axis `e`.
This function assumes that `norm(e) == 1`.

From the two possible solutions `q` the one is returned that is closer
to `q_guess` (note, `q` and `-q` define the same rotation).
"""
@inline function qrot_e(e::Vector3D, angle::Number; q_guess::Quaternion=NullQuaternion)::Quaternion
    sa = sin(angle / 2)
    q = Quaternion(e[1] * sa, e[2] * sa, e[3] * sa, cos(angle / 2))

    return dot(q, q_guess) >= 0 ? q : -q
end
qrot_e(e::AbstractVector, angle::Number)::Quaternion = qrot_e(Vector3D(e), convert(Float64, angle))

"""
    q = ModiaMath.qrot_nxy(nx, ny)

It is assumed that the two input vectors `nx` and `ny` are resolved in frame 1 and
are directed along the x and y axis of frame 2.
The function returns the Quaternion `q` to rotate from frame 1 to frame 2.

The function is robust in the sense that it returns always a Quaternion `q`,
even if `ny` is not orthogonal to `nx` or if one or both vectors have zero length.
This is performed in the following way:
If `nx` and `ny` are not orthogonal to each other, first a unit vector `ey` is
determined that is orthogonal to `nx` and is lying in the plane spanned by
`nx` and `ny`. If `nx` and `ny` are parallel or nearly parallel to each other
or `ny` is a vector with zero or nearly zero length, a vector `ey` is selected
arbitrarily such that `ex` and `ey` are orthogonal to each other.
If both `nx` and `ny` are vectors with zero or nearly zero length, an
arbitrary Quaternion `q` is returned.

# Example

```julia
using Unitful
import ModiaMath

q1 = ModiaMath.qrot1(90u"°")
q2 = ModiaMath.qrot_nxy([1  , 0, 0], [0  , 0, 1  ])
q3 = ModiaMath.qrot_nxy([0.9, 0, 0], [1.1, 0, 1.1])
isapprox(q1,q2)   # returns true
isapprox(q1,q3)   # returns true
```
"""
qrot_nxy(nx, ny)::Quaternion = from_R(rot_nxy(nx, ny))

resolve1(q::Quaternion, v2::Vector3D)::Vector3D =
    2 * ((q[4] * q[4] - 0.5) * v2 + dot(q[1:3], v2) * q[1:3] + q[4] * cross(q[1:3], v2))

resolve2(q::Quaternion, v1::Vector3D)::Vector3D =
    2 * ((q[4] * q[4] - 0.5) * v1 + dot(q[1:3], v1) * q[1:3] - q[4] * cross(q[1:3], v1))

relativeRotation(q1::Quaternion, q2::Quaternion)::Quaternion =
     (@SMatrix [ q1[4]  q1[3] -q1[2] -q1[1];
                -q1[3]  q1[4]  q1[1] -q1[2];
                 q1[2] -q1[1]  q1[4] -q1[3];
                 q1[1]  q1[2]  q1[3]  q1[4]]) * q2

inverseRotation(q::Quaternion)::Quaternion = Quaternion(-q[1], -q[2], -q[3], q[4])