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

using Markdown
using InteractiveUtils

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

# ╔═╡ ea9b4e8b-cfc1-4177-bd3e-45608aca9398
function skip_as_script(m::Module)
	if isdefined(m, :PlutoForceDisplay)
		return m.PlutoForceDisplay
	else
		isdefined(m, :PlutoRunner) && parentmodule(m) == Main
	end
end

# ╔═╡ f85cb5fc-246d-46a1-9875-9c06b7f102d6
if skip_as_script(@__MODULE__)
	import Pkg
	Pkg.activate(Base.current_project(@__DIR__))
	Text("Project env active")
end

# ╔═╡ eab272db-413a-44e4-9c07-b1c1b96e9a5c
if skip_as_script(@__MODULE__)
	using PlutoUI
end

# ╔═╡ 4d94fbcd-f53f-4a90-b01c-f4dbe83313d1
"""
	@skip_as_script expression

Marks a expression as Pluto-only, which means that it won't be executed when running outside Pluto. Do not use this for your own projects.
"""
macro skip_as_script(ex) skip_as_script(__module__) ? esc(ex) : nothing end

# ╔═╡ 460e57a6-d40e-4892-a945-55b03141b1d7
macro only_as_script(ex) skip_as_script(__module__) ? nothing : esc(ex) end

# ╔═╡ b88673ff-30b2-4aca-b0d1-b55a1cd393e7
import AbstractPlutoDingetjes.Bonds

# ╔═╡ b230e180-8e51-4490-8c3e-66a8829a3b7e
import AbstractPlutoDingetjes

# ╔═╡ 80906b18-59b2-4fdf-ba8b-e8080da7066e
md"""
# `transformed_value`

A high-level widget to transform a bond value.
"""

# ╔═╡ 36ff1ff5-cd12-4d66-835b-2b1a6aa9020d
md"""
## Examples
"""

# ╔═╡ 906a7a0c-3ae6-42e6-be62-89d504c348ba
md"""
## The magic
"""

# ╔═╡ 70f18bbf-c933-4017-a1e4-e4dc84f6e696
const compat_error = HTML("<span>❌ You need to update Pluto to use this PlutoUI element.</span>")

# ╔═╡ 666c938a-87fd-4db2-b4f6-992e3c7ef0d9
begin
	struct TransformedWidget{T}
		x::T
		transform::Function
		get_initial_value::Union{Nothing,Function}
	end

	function Base.show(io::IO, m::MIME"text/html", tw::TransformedWidget)
		supported = AbstractPlutoDingetjes.is_supported_by_display(io, Bonds.transform_value) && AbstractPlutoDingetjes.is_supported_by_display(io, Bonds.initial_value)

		return Base.show(io, m, supported ? tw.x : compat_error)
	end

	# AbstractPlutoDingetjes.jl

	function Bonds.transform_value(tw::TransformedWidget, from_js)
		tw.transform(Bonds.transform_value(tw.x, from_js))
	end

	function Bonds.initial_value(tw::TransformedWidget)
		if tw.get_initial_value !== nothing
			tw.get_initial_value()
		else
			child_initial_value = Bonds.initial_value(tw.x)
			if child_initial_value === missing
				missing
			else
				try
					tw.transform(child_initial_value)
				catch e
					@warn "`PlutoUI.transformed_value`: Failed to apply my transform function to the initial value of my contained element." child_initial_value exception=(e,catch_backtrace())
					missing
				end
			end
		end
	end

	# These next two methods are about the value *before* transformation
	# so the user does not need to define those. Yay!
	Bonds.possible_values(tw::TransformedWidget) =
		Bonds.possible_values(tw.x)

	Bonds.validate_value(tw::TransformedWidget, from_browser) =
		Bonds.validate_value(tw.x, from_browser)

	TransformedWidget
end

# ╔═╡ 78795f88-46b9-40af-8100-4e05cfaf3b85
"""
```julia
transformed_value(transform::Function, widget::Any)
```

Create a new widget that wraps around an existing one, with a **value transformation**.

This function creates a so-called *high-level widget*: it returns your existing widget, but with additional functionality. You can use it in your package

# Example
A simple example to get the point accross:
```julia
function RepeatedTextSlider(text::String)
	old_widget = PlutoUI.Slider(1:10)

	# our transformation function
	transform = input -> repeat(text, input)

	# use `transformed_value` to add the value tranformation to our widget
	new_widget = transformed_value(transform, old_widget)
	return new_widget
end

@bind greeting RepeatedTextSlider("hello")

# moving the slider to the right...

greeting == "hellohellohello"
```

![screenshot of the above code in action](https://user-images.githubusercontent.com/6933510/146782076-a993f50c-de27-4a6b-956d-264a5002bfba.gif)

---

This function is very useful in combination with `PlutoUI.combine`. Let's enhance our previous example by **adding a text box** where the repeated text can be entered. If you have not used `PlutoUI.combine` yet, you should read about that first.

```julia
function RepeatedTextSlider()
	old_widget = PlutoUI.combine() do Child
		md""\" \$(Child(PlutoUI.TextField())) \$(Child(PlutoUI.Slider(1:10)))""\"
	end

	# Note that the input to `transform` is now a Tuple!
	# (This is the output of `PlutoUI.combine`)
	transform = input -> repeat(input[1], input[2])

	# use `transformed_value` to add the value tranformation to our widget
	new_widget = transformed_value(transform, old_widget)
	return new_widget
end
```

![screenshot of the above code in action](https://user-images.githubusercontent.com/6933510/146782947-45d67770-03fe-4cf7-82ce-0b9f877688f4.gif)

"""
function transformed_value(f::Function, x::Any; get_initial_value::Union{Nothing,Function}=nothing)
	TransformedWidget(
		x,
		f,
		get_initial_value
	)
end

# ╔═╡ 22109e0c-1815-4b70-9f8c-182b8fe186ea
@skip_as_script function RepeatedTextSlider(text::String)
	old_widget = PlutoUI.Slider(1:10)

	# our transformation function
	transform = input -> repeat(text, input)

	# use `transformed_value` to add the value tranformation to our widget
	new_widget = transformed_value(transform, old_widget)
	return new_widget
end

# ╔═╡ 99df4d14-4fa3-47f1-a002-047d1be18e27
@skip_as_script function RepeatedTextSlider()
	old_widget = PlutoUI.combine() do Child
		md""" $(Child(PlutoUI.TextField())) $(Child(PlutoUI.Slider(1:10)))"""
	end

	# Note that the input to `transform` is now a Tuple!
	# (This is the output of `PlutoUI.combine`)
	transform = input -> repeat(input[1], input[2])

	# use `transformed_value` to add the value tranformation to our widget
	new_widget = transformed_value(transform, old_widget)
	return new_widget
end

# ╔═╡ 7d449272-2e96-418a-80a8-50479cefea7c
@skip_as_script @bind greeting RepeatedTextSlider("hello")

# ╔═╡ 3ab47131-fa65-4e6b-8205-b11620795bfc
@skip_as_script greeting

# ╔═╡ 22004433-904a-42ae-8407-db87a6417389
@skip_as_script @bind custom_greeting RepeatedTextSlider()

# ╔═╡ aecc84f8-3024-4441-ae00-ca1c00db72c0
@skip_as_script custom_greeting

# ╔═╡ Cell order:
# ╟─ea9b4e8b-cfc1-4177-bd3e-45608aca9398
# ╟─4d94fbcd-f53f-4a90-b01c-f4dbe83313d1
# ╟─460e57a6-d40e-4892-a945-55b03141b1d7
# ╟─f85cb5fc-246d-46a1-9875-9c06b7f102d6
# ╠═b88673ff-30b2-4aca-b0d1-b55a1cd393e7
# ╠═b230e180-8e51-4490-8c3e-66a8829a3b7e
# ╠═eab272db-413a-44e4-9c07-b1c1b96e9a5c
# ╟─80906b18-59b2-4fdf-ba8b-e8080da7066e
# ╟─36ff1ff5-cd12-4d66-835b-2b1a6aa9020d
# ╠═22109e0c-1815-4b70-9f8c-182b8fe186ea
# ╠═7d449272-2e96-418a-80a8-50479cefea7c
# ╠═3ab47131-fa65-4e6b-8205-b11620795bfc
# ╠═99df4d14-4fa3-47f1-a002-047d1be18e27
# ╠═22004433-904a-42ae-8407-db87a6417389
# ╠═aecc84f8-3024-4441-ae00-ca1c00db72c0
# ╟─906a7a0c-3ae6-42e6-be62-89d504c348ba
# ╠═78795f88-46b9-40af-8100-4e05cfaf3b85
# ╠═666c938a-87fd-4db2-b4f6-992e3c7ef0d9
# ╟─70f18bbf-c933-4017-a1e4-e4dc84f6e696


text:
using SafeTestsets
@safetestset "snowflakes test" begin include("testscripts.jl") end

ext: jl


text:
abstract type Velocity end

struct NoVelocity<:Velocity end

struct GaussianVelocity<:Velocity
    T::Float64
    seed::Int64
    function GaussianVelocity(  T::Float64,
                                seed::Int64)
        new(T::Float64,
            seed::Int64)
    end
end

struct ZeroVelocity<:Velocity end


text:
module DescentMethods

using Compat, Printf
using NLPModels
using Stopping, LineSearch
using LinearAlgebra
using LinearOperators
using Krylov
using HSL
using LDLFactorizations
using SparseArrays
using SolverTools, Logging

include("HessianDense.jl")
include("HessianOp.jl")
include("HessianSparse.jl")
include("ldlt_symm.jl")
include("Newton.jl")
include("cgTN.jl")
include("NewtonCG.jl")
include("NewtonLDLT.jl")

include("formulae.jl")
include("CG_generic.jl")
include("CG_FR.jl")
include("CG_PR.jl")
include("CG_HS.jl")
include("CG_HZ.jl")

include("hess_approx.jl")
include("shamanskii-LS.jl")
include("shamanskii-LS-MA57.jl")

include("Cheb-dir-AD.jl")
include("Chebyshev.jl")

include("Halley-dir-AD.jl")
include("Halley.jl")

include("SHalley-dir-AD.jl")
include("SuperHalley.jl")

include("NewtonLDLtAbs.jl")
include("NewtonSpectralAbs.jl")

include("lbfgs.jl")

include("ShamanskiiDir.jl")
include("shamanskii-m.jl")

include("autodiff_high_order_model.jl")

end # module


text:
export HONData,
    SpIntMat,
    SpFltMat,
    example_dataset,
    NbrSetMap,
    common_neighbors_map,
    neighbor_pairs,
    num_open_closed_triangles,
    sorted_tuple,
    enum_open_triangles,
    new_closures,
    make_sparse_ones,
    basic_matrices,
    simplex_degree_order,
    proj_graph_degree_order,
    nz_row_inds,
    remove_diagonal,
    split_data,
    triangle_closed,
    tetrahedron_closed

"""
example_dataset
---------------

Returns one of the example datasets.

example_dataset(dataset::String)

Input parameter:
- dataset::String: one of "email-Enron", "contact-primary-school", "contact-high-school", "example1", or "example2"

"""
function example_dataset(dataset::String)
    if !(dataset in ["email-Enron", "contact-primary-school", "contact-high-school", "example1", "example2"])
        error("Unknown dataset $dataset")
    end
    dir = joinpath(dirname(dirname(@__FILE__)),"data")
    read(filename::String) = convert(Vector{Int64}, readdlm(filename, Int64)[:, 1])
    return HONData(read("$(dir)/$(dataset)/$(dataset)-simplices.txt"),
                   read("$(dir)/$(dataset)/$(dataset)-nverts.txt"),
                   read("$(dir)/$(dataset)/$(dataset)-times.txt"),
                   dataset)
end

"""
SpIntMat
--------

const SpIntMat = SparseMatrixCSC{Float64,Int64}
"""
const SpIntMat = SparseMatrixCSC{Int64,Int64}

"""
SpFltMat
--------

const SpFltMat = SparseMatrixCSC{Float64,Int64}
"""
const SpFltMat = SparseMatrixCSC{Float64,Int64}

"""
SpMat
--------

const SpMat = Union{SpIntMat,SpFltMat}
"""
const SpMat = Union{SpIntMat,SpFltMat}

"""
NbrSetMap
--------

const NbrSetMap = Dict{NTuple{2, Int64}, Set{Int64}}
"""
const NbrSetMap = Dict{NTuple{2, Int64}, Set{Int64}}

"""
HONData
-------

Data structure for a temporal higher-order network.

Each dataset consists of three integer vectors: simplices, nverts, and times.
- The simplices is a contiguous vector of nodes comprising the simplices.
- The nverts vector contains the number of vertices within each simplex.
- The times vector contains the timestamps of the simplices (same length as nverts).

For example, consider a dataset consisting of three simplices:

    1. {1, 2, 3} at time 10
    2. {2, 4} at time 15.
    3. {1, 3, 4, 5} at time 21.

Then the data structure would be
- simplices = [1, 2, 3, 2, 4, 1, 3, 4, 5]
- nverts = [3, 2, 4]
- times = [10, 15, 21]
There is an additional name variable attached to the dataset.
"""
struct HONData
    simplices::Vector{Int64}
    nverts::Vector{Int64}
    times::Vector{Int64}
    name::String
end

"""
common_neighbors_map
-------------------

Construct a map where a key is an edge in the graph B participating in at least
one triangle of interest and a value is the vector of common neighbors of the
end points of the edge. The graph B is assumed to be undirected, and the keys
are ordered by the pair with smallest ID first.

common_neighbors_map(B::SpIntMat, triangles::Vector{NTuple{3,Int64}})

Input parameters:
- B::SpIntMat: the graph
- triangles::Vector{NTuple{3,Int64}}: the triangles of interest
"""
function common_neighbors_map(B::SpIntMat, triangles::Vector{NTuple{3,Int64}})
    I = zeros(Int64, 3 * length(triangles))
    J = zeros(Int64, 3 * length(triangles))
    Threads.@threads for ind = 1:length(triangles)
        (i, j, k) = triangles[ind]
        ran = ((ind - 1) * 3 + 1):(ind * 3)
        a, b, c = sort([i, j, k], alg=InsertionSort)
        I[ran] = [a, a, b]
        J[ran] = [b, c, c]
    end
    n = size(B, 2)
    T = sparse(I, J, ones(length(I)), n, n)

    nthreads = Threads.nthreads()
    common_nbrs_vec = Vector{NbrSetMap}(undef, nthreads)
    Threads.@threads for tid = 1:nthreads
        common_nbrs_vec[tid] = NbrSetMap()
    end

    Threads.@threads for j = 1:n
        tid = Threads.threadid()
        if tid == 1
            print(stdout, "$j of $n \r")
            flush(stdout)
        end
        Bj = Set{Int64}(nz_row_inds(B, j))
        # only collect data on edges that appear in triangles
        for i in filter(v -> (v < j) && T[v, j] > 0, Bj)
            Bi = Set{Int64}(nz_row_inds(B, i))
            common_nbrs_vec[tid][(i, j)] = intersect(Bi, Bj)
        end
    end

    # Combine the maps
    common_nbrs = common_nbrs_vec[1]
    for tid = 2:nthreads
        for (key, val) in common_nbrs_vec[tid]
            common_nbrs[key] = val
        end
    end
    return common_nbrs
end

function new_closures(old_simplices::Vector{Int64}, old_nverts::Vector{Int64},
                      new_simplices::Vector{Int64}, new_nverts::Vector{Int64})
    A_old, A_old_t, B_old = basic_matrices(old_simplices, old_nverts)
    simp_order = simplex_degree_order(A_old_t)
    n = size(B_old, 1)
    is_new_node = sum(A_old_t, dims=1) .== 0
    new_triangles = Set{NTuple{3,Int64}}()

    curr_ind = 1
    for nvert in new_nverts
        simp = new_simplices[curr_ind:(curr_ind + nvert - 1)]
        curr_ind += nvert
        for (i, j, k) in combinations(simp, 3)
            # SKIP if any duplicates
            if i == j || i == k || j == k; continue; end
            # SKIP if 0-simplex is missing
            if max(i, j, k) > n || any(is_new_node[[i, j, k]]); continue; end
            # SKIP if already closed
            if triangle_closed(A_old, A_old_t, simp_order, i, j, k); continue; end
            # SKIP if we have already looked at this triangle
            tri_key = sorted_tuple(i, j, k)
            if tri_key in new_triangles; continue; end
            # Newly closed triangle
            push!(new_triangles, tri_key)
        end
    end

    return new_triangles
end

function enum_open_triangles(simplices::Vector{Int64}, nverts::Vector{Int64})
    A, At, B = basic_matrices(simplices, nverts)
    simp_order = simplex_degree_order(At)
    tri_order = proj_graph_degree_order(B)
    n = size(B, 2)

    nthreads = Threads.nthreads()
    triangles = Vector{Vector{NTuple{3,Int64}}}(undef, nthreads)
    for i in 1:nthreads; triangles[i] = Vector{NTuple{3,Int64}}(); end

    # Shuffle so that data better distributes over threads
    shuffled_inds = shuffle(collect(1:n))
    Threads.@threads for ii = 1:length(shuffled_inds)
        tid = Threads.threadid()
        i = shuffled_inds[ii]
        for (j, k) in neighbor_pairs(B, tri_order, i)
            if B[j, k] > 0 && !triangle_closed(A, At, simp_order, i, j, k)
                push!(triangles[tid], (i, j, k))
            end
        end
    end

    # Combine arrays
    total = sum([length(ti) for ti in triangles])
    combined_triangles = Vector{NTuple{3,Int64}}(undef, total)
    curr_ind = 1
    for i in 1:nthreads
        size = length(triangles[i])
        combined_triangles[curr_ind:(curr_ind + size - 1)] = triangles[i][:]
        curr_ind += size
    end
    return combined_triangles
end

""" Turns 3 integers into a sorted tuple. """
sorted_tuple(a::Int64, b::Int64, c::Int64) =
    NTuple{3, Int64}(sort([a, b, c], alg=InsertionSort))
""" Turns 4 integers into a sorted tuple. """
sorted_tuple(a::Int64, b::Int64, c::Int64, d::Int64) =
    NTuple{4, Int64}(sort([a, b, c, d], alg=InsertionSort))

function bipartite_graph(simplices::Vector{Int64}, nverts::Vector{Int64})
    if length(simplices) == 0
        return convert(SpIntMat, sparse([], [], []))
    end
    I, J = Int64[], Int64[]
    curr_ind = 1
    for (simplex_ind, nv) in enumerate(nverts)
        for vert in simplices[curr_ind:(curr_ind + nv - 1)]
            push!(I, vert)
            push!(J, simplex_ind)
        end
        curr_ind += nv
    end
    return convert(SpIntMat, sparse(I, J, 1, maximum(simplices), length(nverts)))
end

"""
basic_matrices
--------------

Computes some simple matrices associated with a dataset.

basic_matrices(simplices::Vector{Int64}, nverts::Vector{Int64})

Input parameters:
- simplices::Vector{Int64}: the contiguous vector of simplices
- nverts::Vector{Int64}: the vector of sizes of simplices

Outputs tuple (A, At, B):
- A::SpIntMat: (# nodes) x (# simplices) adjacency matrix
- At::SpIntMat: the transpose of A
- B::SpIntMat: Projected graph as a Sparse integer matrix, where B[i, j] is the number of times that i and j co-appear in a simplex.
"""
function basic_matrices(simplices::Vector{Int64}, nverts::Vector{Int64})
    A = bipartite_graph(simplices, nverts)
    At = A'
    B = A * At
    B -= sparse(Diagonal(B))
    dropzeros!(B)
    return (A, convert(SpIntMat, At), B)
end

"""
basic_matrices
--------------

Computes some simple matrices associated with a dataset.

basic_matrices(dataset::HONData)

Input parameter:
- dataset::HONData: the dataset

Outputs tuple (A, At, B):
- A::SpIntMat: (# nodes) x (# simplices) adjacency matrix
- At::SpIntMat: the transpose of A
- B::SpIntMat: Projected graph as a Sparse integer matrix, where B[i, j] is the number of times that i and j co-appear in a simplex.
"""
basic_matrices(dataset::HONData) =
    basic_matrices(dataset.simplices, dataset.nverts)

nz_row_inds(A::SpIntMat, ind::Int64) = A.rowval[A.colptr[ind]:(A.colptr[ind + 1] - 1)]
nz_row_inds(A::SpFltMat, ind::Int64) = A.rowval[A.colptr[ind]:(A.colptr[ind + 1] - 1)]
nz_row_vals(A::SpIntMat, ind::Int64) = A.nzval[A.colptr[ind]:(A.colptr[ind + 1] - 1)]
nz_row_vals(A::SpFltMat, ind::Int64) = A.nzval[A.colptr[ind]:(A.colptr[ind + 1] - 1)]

"""
triangle_closed
--------------

Checks if a triangle is closed.

triangle_closed(A::SpIntMat, At::SpIntMat, order::Vector{Int64},
                i::Int64, j::Int64, k::Int64)

Input parameter:
- A::SpIntMat: (# nodes) x (# simplices) adjacency matrix
- At::SpIntMat: the transpose of A
- order::Vector{Int64}: ordering of the nodes
- i::Int64: first node
- j::Int64: second node
- k::Int64: third node

Outputs true or false for closed or open triangle
"""
function triangle_closed(A::SpIntMat, At::SpIntMat, order::Vector{Int64},
                         i::Int64, j::Int64, k::Int64)
    ind, nbr1, nbr2 = sort([i, j, k], by=(v -> order[v]), alg=InsertionSort)
    # Search all simplices of least common vertex
    for simplex_id in nz_row_inds(At, ind)
        if A[nbr1, simplex_id] > 0 && A[nbr2, simplex_id] > 0
            return true
        end
    end
    return false
end

"""
tetrahedron_closed
--------------

Checks if a tetrahedron is closed.

tetrahedron_closed(A::SpIntMat, At::SpIntMat, order::Vector{Int64},
                   i::Int64, j::Int64, k::Int64, l::Int64)

Input parameter:
- A::SpIntMat: (# nodes) x (# simplices) adjacency matrix
- At::SpIntMat: the transpose of A
- order::Vector{Int64}: ordering of the nodes
- i::Int64: first node
- j::Int64: second node
- k::Int64: third node
- l::Int64: fourth node

Outputs true or false for closed or open tetrahedron
"""
function tetrahedron_closed(A::SpIntMat, At::SpIntMat, order::Vector{Int64},
                            i::Int64, j::Int64, k::Int64, l::Int64)
    ind, nbr1, nbr2, nbr3 = sort([i, j, k, l], by=(v -> order[v]), alg=InsertionSort)
    # Search all simplices of least common vertex
    for simplex_id in nz_row_inds(At, ind)
        if A[nbr1, simplex_id] > 0 && A[nbr2, simplex_id] > 0 && A[nbr3, simplex_id] > 0
            return true
        end
    end
    return false
end

"""
make_sparse_ones
---------------

Returns a new sparse matrix with the same non-zero pattern as the input but
where all non-zeros are set to 1.

make_sparse_ones(A::SpIntMat)

Input parameter:
- A::SpIntMat: a sparse matrix
"""
function make_sparse_ones(A::SpMat)
    C = copy(A)
    LinearAlgebra.fillstored!(C, 1)
    return C
end

function neighbors(B::SpIntMat, order::Vector{Int64}, node::Int64)
    node_order = order[node]
    return filter(nbr -> order[nbr] > node_order, nz_row_inds(B, node))
end

neighbor_pairs(B::SpIntMat, order::Vector{Int64}, node::Int64) =
    combinations(neighbors(B, order, node), 2)

""" Ordering of nodes by the number of simplices in which they appear """
function simplex_degree_order(At::SpIntMat)
    n = size(At, 2)
    simplex_order = zeros(Int64, n)
    simplex_order[sortperm(vec(sum(At, dims=1)))] = collect(1:n)
    return simplex_order
end

""" Ordering of nodes by their degree """
function proj_graph_degree_order(B::SpIntMat)
    n = size(B, 1)
    triangle_order = zeros(Int64, n)
    triangle_order[sortperm(vec(sum(make_sparse_ones(B), dims=1)))] = collect(1:n)
    return triangle_order
end

function num_open_closed_triangles(A::SpIntMat, At::SpIntMat, B::SpIntMat)
    simp_order = simplex_degree_order(At)
    tri_order = proj_graph_degree_order(B)
    n = size(B, 2)
    counts = zeros(Int64, 2, Threads.nthreads())
    Threads.@threads for i = 1:n
        for (j, k) in neighbor_pairs(B, tri_order, i)
            if B[j, k] > 0
                tid = Threads.threadid()
                closed = triangle_closed(A, At, simp_order, i, j, k)
                counts[1 + closed, tid] += 1
            end
        end
    end
    return tuple(vec(sum(counts, dims=2))...)
end

"""
num_open_closed_triangles
-------------------------

Computes the number of open and closed triangles in a dataset.

num_open_closed_triangles(data::HONData)

Input parameter:
- data::HONData: the dataset

Outputs tuple (no, nc):
- no: number of open triangles
- nc: number of closed triangles
"""
num_open_closed_triangles(data::HONData) =
    num_open_closed_triangles(basic_matrices(data)...)

# Get events in the [start_time, end_time] interval
function window_data(start_time::Int64, end_time::Int64, simplices::Vector{Int64},
                     nverts::Vector{Int64}, times::Vector{Int64})
    curr_ind = 1
    window_simplices, window_nverts, window_times = Int64[], Int64[], Int64[]
    for (nv, time) in zip(nverts, times)
        end_ind = curr_ind + nv - 1
        if time >= start_time && time <= end_time
            push!(window_nverts, nv)
            push!(window_times, time)
            append!(window_simplices, simplices[curr_ind:end_ind])
        end
        curr_ind += nv
    end
    return (window_simplices, window_nverts, window_times)
end

"""
split_data
----------

Split data by timestamps into quantiles specified by quantile1 and
quantile2. Returns a 4-tuple (old_simps, old_nverts, new_simps, new_nverts),
where (old_simps, old_nverts) are the data in the quantile [0, quantile1]
and (new_simps, new_nverts) are the data in the quantile (quantile1, quantile2].

split_data(simplices::Vector{Int64}, nverts::Vector{Int64},
           times::Vector{Int64}, quantile1::Int64,
           quantile2::Int64)

Input parameters:
- simplices::Vector{Int64}: contiguous vector of simplices from dataset
- nverts::Vector{Int64}: vector of simplex sizes from dataset
- times::Vector{Int64}: vector of timestamps of simplices
- quantile1::Int64: first quantile
- quantile2::Int64: second quantile

Returns a tuple (old_simps, old_nverts, new_simps, new_nverts):
- old_simplices::Vector{Int64}: simplices in quantile1
- old_nverts::Vector{Int64}: sizes of simplices in old_simplices
- new_simplices::Vector{Int64}: simplices between quantile1 and quantile 2
- new_nverts::Vector{Int64}: sizes of simplices in new_simplices
"""
function split_data(simplices::Vector{Int64}, nverts::Vector{Int64},
                    times::Vector{Int64}, quantile1::Int64,
                    quantile2::Int64)
    if quantile1 > quantile2
        error("First quantile ($quantile1) needs to be <= second quantile ($quantile2)")
    end
    cutoff(prcntl::Int64) = convert(Int64, round(percentile(times, prcntl)))

    cutoff1 = cutoff(quantile1)
    old_simps, old_nverts =
        window_data(minimum(times), cutoff1, simplices, nverts, times)[1:2]

    cutoff2 = (quantile2 == 100) ? (maximum(times) + 1) : cutoff(quantile2)
    new_simps, new_nverts =
        window_data(cutoff1 + 1, cutoff2, simplices, nverts, times)[1:2]

    return old_simps, old_nverts, new_simps, new_nverts
end

function backbone(simplices::Vector{Int64}, nverts::Vector{Int64},
                  times::Vector{Int64})
    # backbone data
    bb_simplices, bb_nverts, bb_times = Int64[], Int64[], Int64[]

    # contains for all simplices
    max_size = maximum(nverts)
    all_simplices = Vector{Set}(undef, max_size)
    for i in 1:max_size; all_simplices[i] = Set{Any}(); end

    curr_ind = 1
    for (nvert, time) in zip(nverts, times)
        simplex = simplices[curr_ind:(curr_ind + nvert - 1)]
        nvert = length(simplex)
        # Add to data if we have not seen it yet
        if !(simplex in all_simplices[nvert])
            push!(all_simplices[nvert], simplex)
            # add to backbone
            append!(bb_simplices, simplex)
            push!(bb_nverts, nvert)
            push!(bb_times, time)
        end
        curr_ind += nvert
    end

    return (bb_simplices, bb_nverts, bb_times)
end

# Get a configuration that preserves the number of k-vertex simplices that every
# vertex participates in.
function configuration_sizes_preserved(simplices::Vector{Int64},
                                       nverts::Vector{Int64})
    config = zeros(Int64, length(simplices))
    app = copy(nverts)
    pushfirst!(app, 1)
    capp = cumsum(app)
    for val in unique(nverts)
        # Get the simplices with this number of vertices
        inds = Int64[]
        cnt = 0
        for ind in findall(nverts .== val)
            append!(inds, capp[ind]:(capp[ind] + val  - 1))
            cnt += 1
        end
        # TODO(arb): this isn't necessarily creating a simple graph
        config[inds] = shuffle(simplices[inds])
    end
    return config
end


text:
function initTxUHD(pointerUSRP)
	# ----------------------------------------------------
	# --- Tx Streamer
	# ----------------------------------------------------
	# --- Create a pointer related to the Tx streamer
	addressStream = Ref{Ptr{uhd_tx_streamer}}();
	# --- Cal the init
	@assert_uhd ccall((:uhd_tx_streamer_make, libUHD), uhd_error, (Ptr{Ptr{uhd_tx_streamer}},), addressStream);
	streamerPointer = addressStream[];
	# ----------------------------------------------------
	# --- Tx Metadata
	# ----------------------------------------------------
	# --- Create a pointer related to Metadata
	addressMD = Ref{Ptr{uhd_tx_metadata}}();
	# --- Cal the init
	@assert_uhd ccall((:uhd_tx_metadata_make, libUHD), uhd_error, (Ptr{Ptr{uhd_tx_metadata}}, Cuchar, FORMAT_LONG, Cdouble, Cuchar, Cuchar), addressMD, false, 0, 0.1, true, false);
	# --- Get the usable object
	metadataPointer = addressMD[];
	# ----------------------------------------------------
	# --- Create the USRP wrapper object
	# ----------------------------------------------------
	uhd  = UHDTxWrapper(true, pointerUSRP, streamerPointer, metadataPointer, addressStream, addressMD);
	# @infotx("Done init \n");
	return uhd;
end

function openUHDTx(pointerUSRP,carrierFreq, samplingRate, gain, antenna = "TX-RX";args="")
	# ----------------------------------------------------
	# --- Init  UHD object
	# ----------------------------------------------------
	uhd	  = initTxUHD(pointerUSRP);
	# ----------------------------------------------------
	# --- Creating Runtime structures
	# ----------------------------------------------------
	# --- Create structure for UHD argument
	# TODO Adding custom levels here for user API
	channel		   =  Ref{Csize_t}(0);
	a1			   =  Base.unsafe_convert(Cstring, "fc32");
	a2			   =  Base.unsafe_convert(Cstring, "sc16");
	a3			   =  Base.unsafe_convert(Cstring, "");
	uhdArgs		   = uhd_stream_args_t(a1, a2, a3, channel, 1);
	# ----------------------------------------------------
	# --- Sampling rate configuration
	# ----------------------------------------------------
	# --- Update the Tx sampling rate
	ccall((:uhd_usrp_set_tx_rate, libUHD), Cvoid, (Ptr{uhd_usrp}, Cdouble, Csize_t), pointerUSRP, samplingRate, 0);
	# --- Get the Tx rate from the radio
	pointerRate  = Ref{Cdouble}(0);
	ccall((:uhd_usrp_get_tx_rate, libUHD), Cvoid, (Ptr{uhd_usrp}, Csize_t, Ref{Cdouble}), pointerUSRP, 0, pointerRate);
	updateRate  = pointerRate[];
	# --- Print a flag
	if updateRate != samplingRate
		@warntx "Effective Rate is $(updateRate / 1e6) MHz and not $(samplingRate / 1e6) MHz\n"
	else
		# @infotx "Effective Rate is $(updateRate / 1e6) MHz\n";
	end
	# ----------------------------------------------------
	# --- Carrier Frequency configuration
	# ----------------------------------------------------
	# --- Create structure for request
	a5			   =  Base.unsafe_convert(Cstring,"");
	tuneRequest	   = uhd_tune_request_t(carrierFreq,UHD_TUNE_REQUEST_POLICY_AUTO,carrierFreq,UHD_TUNE_REQUEST_POLICY_AUTO,200e6,a5)
	# tuneRequest	   = uhd_tune_request_t(carrierFreq, UHD_TUNE_REQUEST_POLICY_AUTO, UHD_TUNE_REQUEST_POLICY_AUTO);
	tunePointer	  = Ref{uhd_tune_request_t}(tuneRequest);
	pointerTuneResult	  = Ref{uhd_tune_result}();
	ccall((:uhd_usrp_set_tx_freq, libUHD), Cvoid, (Ptr{uhd_usrp}, Ptr{uhd_tune_request_t}, Csize_t, Ptr{uhd_tune_result}), pointerUSRP, tunePointer, 0, pointerTuneResult);
	pointerCarrierFreq = Ref{Cdouble}(0);
	ccall((:uhd_usrp_get_tx_freq, libUHD), Cvoid, (Ptr{uhd_usrp}, Csize_t, Ref{Cdouble}), pointerUSRP, 0, pointerCarrierFreq);
	updateCarrierFreq	= pointerCarrierFreq[];
	if updateCarrierFreq != carrierFreq
		@warntx "Effective carrier frequency is $(updateCarrierFreq / 1e6) MHz and not $(carrierFreq / 1e6) Hz\n"
	else
		# @infotx "Effective carrier frequency is $(updateCarrierFreq / 1e6) MHz\n";
	end
	# ----------------------------------------------------
	# --- Gain configuration
	# ----------------------------------------------------
	# Update the UHD sampling rate
	ccall((:uhd_usrp_set_tx_gain, libUHD), Cvoid, (Ptr{uhd_usrp}, Cdouble, Csize_t, Cstring), pointerUSRP, gain, 0, "");
	# Get the updated gain from UHD
	pointerGain	  = Ref{Cdouble}(0);
	ccall((:uhd_usrp_get_tx_gain, libUHD), Cvoid, (Ptr{uhd_usrp}, Csize_t, Cstring, Ref{Cdouble}), pointerUSRP, 0, "", pointerGain);
	updateGain	  = pointerGain[];
	# --- Print a flag
	if updateGain != gain
		@warntx "Effective gain is $(updateGain) dB and not $(gain) dB\n"
	else
		# @infotx "Effective gain is $(updateGain) dB\n";
	end
	# ----------------------------------------------------
	# --- Antenna configuration
	# ----------------------------------------------------
	# ccall((:uhd_usrp_set_tx_antenna, libUHD), Cvoid, (Ptr{uhd_usrp}, Cstring, Csize_t), uhd.pointerUSRP, antenna, 0);
	# ----------------------------------------------------
	# --- Setting up streamer
	# ----------------------------------------------------
	# --- Setting up arguments
	pointerArgs	  = Ref{uhd_stream_args_t}(uhdArgs);
	ccall((:uhd_usrp_get_tx_stream, libUHD), Cvoid, (Ptr{uhd_usrp}, Ptr{uhd_stream_args_t}, Ptr{uhd_tx_streamer}), pointerUSRP, pointerArgs, uhd.pointerStreamer);
	# --- Getting number of samples ber buffer
	pointerSamples	  = Ref{Csize_t}(0);
	ccall((:uhd_tx_streamer_max_num_samps, libUHD), Cvoid, (Ptr{uhd_stream_args_t}, Ref{Csize_t}), uhd.pointerStreamer, pointerSamples);
	nbSamples		  = pointerSamples[];
	# ----------------------------------------------------
	# --- Create object and return
	# ----------------------------------------------------
	# --- Return
	return UHDTx(uhd, updateCarrierFreq, updateRate, updateGain, antenna, nbSamples, 0);
end

function updateSamplingRate!(radio::UHDTx, samplingRate)
	# ----------------------------------------------------
	# --- Sampling rate configuration
	# ----------------------------------------------------
	# @infotx  "Try to change rate from $(radio.samplingRate / 1e6) MHz to $(samplingRate / 1e6) MHz";
	# --- Update the Tx sampling rate
	ccall((:uhd_usrp_set_tx_rate, libUHD), Cvoid, (Ptr{uhd_usrp}, Cdouble, Csize_t), radio.uhd.pointerUSRP, samplingRate, 0);
	# --- Get the Tx rate from the radio
	pointerRate  = Ref{Cdouble}(0);
	ccall((:uhd_usrp_get_tx_rate, libUHD), Cvoid, (Ptr{uhd_usrp}, Csize_t, Ref{Cdouble}), radio.uhd.pointerUSRP, 0, pointerRate);
	updateRate  = pointerRate[];
	# --- Print a flag
	if updateRate != samplingRate
		@warntx "Effective Rate is $(updateRate / 1e6) MHz and not $(samplingRate / 1e6) MHz\n"
	else
		# @infotx "Effective Rate is $(updateRate / 1e6) MHz\n";
	end
	radio.samplingRate = updateRate;
end

function updateGain!(radio::UHDTx, gain)
	# ----------------------------------------------------
	# --- Sampling rate configuration
	# ----------------------------------------------------
	# @infotx  "Try to change gain from $(radio.gain) dB to $(gain) dB";
	# Update the UHD sampling rate
	ccall((:uhd_usrp_set_tx_gain, libUHD), Cvoid, (Ptr{uhd_usrp}, Cdouble, Csize_t, Cstring), radio.uhd.pointerUSRP, gain, 0, "");
	# Get the updated gain from UHD
	pointerGain	  = Ref{Cdouble}(0);
	ccall((:uhd_usrp_get_tx_gain, libUHD), Cvoid, (Ptr{uhd_usrp}, Csize_t, Cstring, Ref{Cdouble}), radio.uhd.pointerUSRP, 0, "", pointerGain);
	updateGain	  = pointerGain[];
	# --- Print a flag
	if updateGain != gain
		@warntx "Effective gain is $(updateGain) dB and not $(gain) dB\n"
	else
		# @infotx "Effective gain is $(updateGain) dB\n";
	end
	radio.gain = updateGain;
	return updateGain;
end

function updateCarrierFreq!(radio::UHDTx, carrierFreq)
	# ----------------------------------------------------
	# --- Carrier Frequency configuration
	# ----------------------------------------------------
	# @infotx  "Try to change carrier frequency from $(radio.carrierFreq / 1e6) MHz to $(carrierFreq / 1e6) MHz";
	a5			   =  Base.unsafe_convert(Cstring,"");
	tuneRequest	   = uhd_tune_request_t(carrierFreq,UHD_TUNE_REQUEST_POLICY_AUTO,carrierFreq,UHD_TUNE_REQUEST_POLICY_AUTO,200e6,a5)
	# tuneRequest   = uhd_tune_request_t(carrierFreq, UHD_TUNE_REQUEST_POLICY_AUTO, UHD_TUNE_REQUEST_POLICY_AUTO);
	tunePointer	  = Ref{uhd_tune_request_t}(tuneRequest);
	pointerTuneResult	  = Ref{uhd_tune_result}();
	ccall((:uhd_usrp_set_tx_freq, libUHD), Cvoid, (Ptr{uhd_usrp}, Ptr{uhd_tune_request_t}, Csize_t, Ptr{uhd_tune_result}), radio.uhd.pointerUSRP, tunePointer, 0, pointerTuneResult);
	pointerCarrierFreq = Ref{Cdouble}(0);
	ccall((:uhd_usrp_get_tx_freq, libUHD), Cvoid, (Ptr{uhd_usrp}, Csize_t, Ref{Cdouble}), radio.uhd.pointerUSRP, 0, pointerCarrierFreq);
	updateCarrierFreq	= pointerCarrierFreq[];
	if updateCarrierFreq != carrierFreq
		@warntx "Effective carrier frequency is $(updateCarrierFreq / 1e6) MHz and not $(carrierFreq / 1e6) MHz\n"
	else
		# @infotx "Effective carrier frequency is $(updateCarrierFreq / 1e6) MHz\n";
	end
	radio.carrierFreq = updateCarrierFreq;
	return updateCarrierFreq;
end

"""
Send a buffer though the radio device. It is possible to force a cyclic buffer send (the radio uninterruply send the same buffer) by setting the cyclic parameter to true

# --- Syntax

send(radio,buffer,cyclic=false)
# --- Input parameters
- radio	  	: UHD device [UHDRx]
- buffer 	: Buffer to be send [Union{Array{Complex{Cfloat}},Array{Cfloat}}]
- cyclic 	: Send same buffer multiple times (default false) [Bool]
# --- Output parameters
- nbEch 	: Number of samples effectively send [Csize_t]. It corresponds to the number of complex samples sent.
"""
function send(radio::UHDTx, buffer::Union{Array{Complex{Cfloat}},Array{Cfloat}}, cyclic::Bool = false)
	# --- Pointer to number of samples transmitted
	pointerSamples 	= Ref{Csize_t}(0);
	# --- Size of buffer
	# We should handle Complex and non complex inputs
	# If Complex, we have length(buffer) elements. If input is not complex, half element are effectively transmitted (real, imag, real, imag...) and we have length(buffer)÷2.
	# This can done in one instruction with sizeof(Cfloat) and sizeof(buffer)
	nbSamples 				= sizeof(buffer) ÷ sizeof(Cfloat) ÷ 2;
	# --- Fragmentation handling
	posT 	= Csize_t(0);
	filled 	= false;
	# --- Accumulated value of number of samples transmitted
	nbEch 			= Csize_t(0);
	try
		# --- First while loop is to handle cyclic transmission
		# It turns to false in case of interruption or cyclic to false
		while (true)
			# --- Getting pointer position
			(posT+radio.packetSize> nbSamples) ? n = nbSamples - posT : n = radio.packetSize;
			# --- Create pointer pointing to the fragment
			ptr=Ref(Ptr{Cvoid}(pointer(buffer,1+posT)));
			# --- Effectively transmit data
			ccall((:uhd_tx_streamer_send, libUHD), uhd_error, (Ptr{uhd_tx_streamer}, Ptr{Ptr{Cvoid}}, Csize_t, Ptr{Ptr{uhd_tx_metadata}}, Cfloat, Ref{Csize_t}), radio.uhd.pointerStreamer, ptr, n, radio.uhd.addressMD, 0.1, pointerSamples);
			# --- Getting number of transmitted samples
			nbEch 		+= pointerSamples[];
			posT 		+= pointerSamples[];
			(posT == nbSamples) ? filled = true : filled = false;
			# --- Detection of cyclic mode
			(cyclic == false && filled == true) && break
			# --- We are in cyclic mode here
			if filled == true
				# --- We are not out and filled is true => cyclic is true and we have done one complete Tx buffer transmit. We shall reinit all the counters here
				posT 	= 0;
				# --- Filled to false
				filled = false;
			end
			# --- Forcing refresh
			yield();
		end
	catch e;
		# --- Interruption handling
		print(e);
		print("\n");
		@infotx "Interruption detected";
        return Csize_t(0);
	end
	# --- Return number of complex samples transmitted
	# We accumulate number of samples transmitted, we should divide by 2 to get number of Complex samples
	return (nbEch ÷ 2);
end

function Base.close(radio::UHDTx)
	# --- Checking realease nature
	# There is one flag to avoid double free (that leads to seg fault)
	if radio.released == 0
		# C Wrapper to ressource release
		@assert_uhd ccall((:uhd_tx_streamer_free, libUHD), uhd_error, (Ptr{Ptr{uhd_tx_streamer}},), radio.uhd.addressStream);
		@assert_uhd ccall((:uhd_tx_metadata_free, libUHD), uhd_error, (Ptr{Ptr{uhd_tx_metadata}},), radio.uhd.addressMD);
	else
		# print a warning
		@warn "UHD ressource was already released, abort call";
	end
	# --- Force flag value
	radio.released = 1;
end

function Base.print(radio::UHDTx)
	# Get the gain from UHD
	pointerGain	  = Ref{Cdouble}(0);
	ccall((:uhd_usrp_get_tx_gain, libUHD), Cvoid, (Ptr{Cvoid}, Csize_t, Cstring, Ref{Cdouble}), radio.uhd.pointerUSRP, 0, "", pointerGain);
	updateGain	  = pointerGain[];
	# Get the rate from UHD
	pointerRate	  = Ref{Cdouble}(0);
	ccall((:uhd_usrp_get_tx_rate, libUHD), Cvoid, (Ptr{Cvoid}, Csize_t, Ref{Cdouble}), radio.uhd.pointerUSRP, 0, pointerRate);
	updateRate	  = pointerRate[];
	# Get the freq from UHD
	pointerFreq	  = Ref{Cdouble}(0);
	ccall((:uhd_usrp_get_tx_freq, libUHD), Cvoid, (Ptr{Cvoid}, Csize_t, Ref{Cdouble}), radio.uhd.pointerUSRP, 0, pointerFreq);
	updateFreq	  = pointerFreq[];
	# Print message
	strF  = @sprintf(" Carrier Frequency: %2.3f MHz\n Sampling Frequency: %2.3f MHz\n Tx Gain: %2.2f dB\n",updateFreq / 1e6,updateRate / 1e6,updateGain);
	@infotx "Current UHD Configuration in Tx mode\n$strF";
end


text:
struct Quadrotor2DArm{I, T} <: Model{I, T}
    n::Int
    m::Int
    d::Int

    # body
    lb    # length
    mb # mass
    Jb    # inertia

    # link 1
    l1
    lc1
    m1
    J1

    # link
    l2
    lc2
    m2
    J2

    g    # gravity
end

function kinematics(model::Quadrotor2DArm, q)
    @SVector [q[1] + model.l1 * sin(q[3] + q[4]) + model.l2 * sin(q[3] + q[4] + q[5]),
              q[2] - model.l1 * cos(q[3] + q[4]) - model.l2 * cos(q[3] + q[4] + q[5])]
end

function jacobian(model::Quadrotor2DArm, q)
    a = model.l1 * cos(q[3] + q[4]) + model.l2 * cos(q[3] + q[4] + q[5])
    b = model.l2 * cos(q[3] + q[4] + q[5])
    c = model.l1 * sin(q[3] + q[4]) + model.l2 * sin(q[3] + q[4] + q[5])
    d = model.l2 * sin(q[3] + q[4] + q[5])
    @SMatrix [1.0 0.0 a a b;
              0.0 1.0 c c d]
end

function B_func(model::Quadrotor2DArm, q)
    a = model.l1 * cos(q[3] + q[4]) + model.l2 * cos(q[3] + q[4] + q[5])
    b = model.l2 * cos(q[3] + q[4] + q[5])
    c = model.l1 * sin(q[3] + q[4]) + model.l2 * sin(q[3] + q[4] + q[5])
    d = model.l2 * sin(q[3] + q[4] + q[5])

    @SMatrix [-0.5 * model.lb * sin(q[3]) -0.5 * model.lb * sin(q[3]) 0.0 0.0 1.0 0.0;
              0.5 * model.lb * cos(q[3]) 0.5 * model.lb * cos(q[3]) 0.0 0.0 0.0 1.0;
              0.0 0.0 0.0 0.0 a c;
              0.0 0.0 1.0 0.0 a c;
              0.0 0.0 0.0 1.0 b d]
end


text:
##################################################
# Compute a map of the dF/dt in action space
# Linear sampling
##################################################
using HDF5
########################################
include("../sources/julia/Main.jl") # Loading the main code
#include("../sources/julia/test/OldFunctions.jl")
########################################

m_field  = _M/nbGlobularCluster # Mass of field star
epsRef = _L0*10^(-5)

########################################
# Action space parameter
########################################

JrminMeasure, JrmaxMeasure = _L0*0.0001, _L0*0.5 # Jr range
LminMeasure, LmaxMeasure = _L0*0.0001, _L0*1.0 # L range
nbJrMeasure = 100 # Number of Jr sampling points
nbLMeasure = 500 # Number of L sampling points

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

tabJrMeasure = collect(range(JrminMeasure,length=nbJrMeasure,JrmaxMeasure))
tabLMeasure = collect(range(LminMeasure,length=nbLMeasure,LmaxMeasure))
nbJrLGrid = nbJrMeasure*nbLMeasure

const tabLJrGrid  = zeros(Float64,2,nbJrLGrid) # Location (L,Jr) of the grid points where the diffusion coefficients are computed
const tabdFdt  = zeros(Float64,nbJrLGrid) # Values of dF/dt on the (L,Jr)-grid

########################################
# Functions to fill the arrays of coefficients
########################################

function tabLJrGrid!()
    index = 1
    for iJr=1:nbJrMeasure
    JrMeasure = tabJrMeasure[iJr]
        for iL=1:nbLMeasure
            LMeasure = tabLMeasure[iL]
            tabLJrGrid[1,index], tabLJrGrid[2,index] = LMeasure, JrMeasure
            index += 1
        end
    end
end

function tabdFdt!()
    if (PARALLEL == "yes")
        println("parallel")
        Threads.@threads for iGrid=1:nbJrLGrid
            LMeasure, JrMeasure = tabLJrGrid[1,iGrid], tabLJrGrid[2,iGrid]
            eps = min(epsRef,JrMeasure*10^(-3),LMeasure*10^(-3))

            #dfdt = -divflux_NR_num_expand(JrMeasure,LMeasure,m_field,eps)
            dfdt = -divflux_NR_num(JrMeasure,LMeasure,m_field,eps)
            tabdFdt[iGrid] = dfdt
        end
    else # Computation is not made in parallel
        for iGrid=1:nbJrLGrid
            LMeasure, JrMeasure = tabLJrGrid[1,iGrid], tabLJrGrid[2,iGrid]
            eps = min(epsRef,JrMeasure*10^(-3),LMeasure*10^(-3))

            #dfdt = -divflux_NR_num_expand(JrMeasure,LMeasure,m_field,eps)
            dfdt = -divflux_NR_num(JrMeasure,LMeasure,m_field,eps)
            tabdFdt[iGrid] = dfdt
        end
    end
end

########################################
namefile = "../data/Dump_dFdt_Map_q_"*string(qCalc)*".hf5"
########################################
# Function that writes to .hf5 files
########################################
function writedump!(namefile)
    file = h5open(namefile,"w") # Opening the file
    write(file,"q",qCalc)
    write(file,"tabJr",tabJrMeasure)
    write(file,"tabL",tabLMeasure)
    write(file,"tabLJr",tabLJrGrid)
    write(file,"tabdFdt",tabdFdt)

    write(file,"nbJrMeasure",nbJrMeasure)
    write(file,"nbLMeasure",nbLMeasure)
    write(file,"Mtot",_M)
    write(file,"Npart",nbGlobularCluster)

    close(file) # Closing the file
end
########################################

@time tabLJrGrid!()
@time tabdFdt!()

########################################
writedump!(namefile) # Dumping the computed table

max_stars_repo_name: KerwannTEP/CAT

max_issues_repo_name: KerwannTEP/CAT

max_forks_repo_name: KerwannTEP/CAT


text:
module Language
export Lexer, advanceLexer

# Each kind of token
const SOF = Symbol("<SOF>")
const EOF = Symbol("<EOF>")
const BANG = Symbol("!")
const DOLLAR = Symbol("\$")
const PAREN_L = Symbol("(")
const PAREN_R = Symbol(")")
const SPREAD = Symbol("...")
const COLON = Symbol(":")
const EQUALS = Symbol("=")
const AT = Symbol("@")
const BRACKET_L = Symbol("[")
const BRACKET_R = Symbol("]")
const BRACE_L = Symbol("{")
const PIPE = Symbol("|")
const BRACE_R = Symbol("}")
const NAME = Symbol("Name")
const INT = Symbol("Int")
const FLOAT = Symbol("Float")
const STRING = Symbol("String")
const COMMENT = Symbol("Comment")

struct Token
  kind::Symbol
  startPosition::Int64
  endPosition::Int64
  line::Int64
  column::Int64
  prev::Union{Symbol,Nothing}
  value::Union{String,Nothing}
end

mutable struct Lexer
  body::String
  lastToken::Token
  currentToken::Token
  line::Int64
  column::Int64

  Lexer(body::String) = new(body,Token(SOF,0,0,0,0,nothing,nothing),Token(SOF,0,0,0,0,nothing,nothing),0,0)
end

function advanceLexer(lexer::Lexer)

end

end


text:
module Test_rtf

using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

println("Test_rtf")

@testset "test_rtf" begin

# two polynomial entries
a = 1; b = 2; c = 3; d = 4;
@time sys = rtf(Polynomial([b, a],:s), Polynomial([d, c],:s), Ts = 0)
s = Polynomial([0,1],:s)
@time sys1 = rtf(a*s+b, c*s+d)
@time sys2 = rtf(Polynomial([b, a]), Polynomial([d, c]), Ts = 0, var = :s)

@test sys == sys1 && sys == sys2
@test sys ≈ sys1  && sys ≈ sys2

# some basic tests
@test propertynames(sys) == (:Ts, :var, :zeros, :poles, :gain, :num, :den)
@test poles(sys) ≈ [-d/c] && gpole(sys) ≈ [-d/c]
@test gzero(sys) ≈ [-b/a] && gain(sys) ≈ a/c #&& length(sys) == 1
@test poles(sys.num) ≈ Int[] && gain(sys.den) ≈ c
@test !isconstant(sys) && !isconstant(sys.num) && isconstant(a)
@test variable(sys) == variable(sys.num)
@test convert(RationalFunction,sys) == RationalFunction(sys.num,sys.den)
@test all(zpk(sys) .≈ ([-2.0], [-1.3333333333333333], 0.3333333333333333))
@test all(zpk(sys.num) .≈ ([-2.0], Int64[], 1))
@test sys' == rtf(-a*s+b, -c*s+d)
@test numpoly(5) == Polynomial(5) && denpoly(5) == Polynomial(1)
@test order(sys) == 1
#@test convert(RationalTransferFunction{Float64},5) == rtf(5.)
@test DescriptorSystems.promote_var(a*s+b,Polynomial(1)) == :s &&
      DescriptorSystems.promote_var(Polynomial(1),a*s+b) == :s &&
      DescriptorSystems.promote_var(Polynomial(1),Polynomial(1)) == :s
@test zero(sys) == rtf(0.) && zero(sys) == 0 && zero(sys) == Polynomial(0) &&
      one(sys) == rtf(1.)
@test isnothing(sampling_time(3.)) && isnothing(sampling_time(Int))
@test DescriptorSystems.promote_Ts(sampling_time.((sys, 3))...) == 0
@test DescriptorSystems._eltype([sys.num;sys.den]) == Int

# check some standard constructors
@time r = RationalTransferFunction{Complex{Float64},:s,Polynomial{Complex{Float64},:s},1}(2,Polynomial([1,1. *im],:s),1)
@time r = RationalTransferFunction{Float64,:s,Polynomial{Float64,:s},1}(Polynomial([1,1.],:s),2//3,1)
@time t = convert(RationalFunction,r)
@time r = RationalTransferFunction(Polynomial([1,1.],:s),Polynomial(2//3*im,:s),1)
@time r = RationalTransferFunction(Polynomial([1,1.],:s),2//3*im,1)
@time r = RationalTransferFunction(2//3*im,Polynomial([1,1.],:s),1)
@time r = RationalTransferFunction(2.1*im,Polynomial([1,1.],:s),1)
@time r = RationalTransferFunction(Polynomial([1,1.],:s),2.1*im,1)
@time r = rtf(2//3,2.1*im,Ts=1)
@time convert(typeof(r),Polynomial([0,1],:s))
rr = RationalFunction(Polynomial(2,:s),Polynomial([1,2],:s));
@test promote_rule(typeof(r),typeof(rr)) == RationalTransferFunction{ComplexF64, :s, Polynomial{ComplexF64, :s}, 1.0}
pol = Polynomial(2,:s);
@test promote_rule(typeof(r),typeof(pol)) == RationalTransferFunction{ComplexF64, :s, Polynomial{ComplexF64, :s}, 1.0}

a = 1; b = 2; c = 3; d = 4;
sysd = rtf(Polynomial([b, a],:z), Polynomial([d, c],:z), Ts = 1)
z = Polynomial([0,1],:z)
sysd1 = rtf(a*z+b, c*z+d, Ts = 1)
sysd2 = rtf(Polynomial([b, a]), Polynomial([d, c]), Ts = 1, var = :z)
@test sysd == sysd1 && sysd == sysd2
@test sysd ≈ sysd1  && sysd == sysd2
@test sysd' == rtf(a+b*z, c+d*z, Ts = 1)
@test rtf(a*z+b, 1, Ts = 1)' == rtf(a+b*z, z, Ts = 1)
@test rtf(1, a*z+b, Ts = 1)' == rtf(z, a+b*z, Ts = 1)
#@test promote_type(typeof(sysd),Float64) == RationalTransferFunction{Float64,:z}

a = 1; b = 2; c = 3; d = 4;
sys = rtf(Polynomial([b, a],:s), d, Ts = 0)
s = Polynomial([0,1],:s)
sys1 = rtf(a*s+b, d)
sys2 = rtf(Polynomial([b, a]), Polynomial([d]), Ts = 0, var = :s)

@test sys == sys1 && sys == sys2
@test sys ≈ sys1  && sys ≈ sys2

a = 1; b = 2; c = 3; d = 4;
sysd = rtf(Polynomial([b],:z), Polynomial([d, c],:z), Ts = 1)
z = Polynomial([0,1],:z)
sysd1 = rtf(b, c*z+d, Ts = 1)
sysd2 = rtf(Polynomial([b]), Polynomial([d, c]), Ts = 1, var = :z)
@test sysd == sysd1 && sysd == sysd2
@test sysd ≈ sysd1  && sysd == sysd2

a = 1; b = 2; c = 3; d = 4;
sys = rtf(Polynomial(b,:s), d, Ts = 0)
s = Polynomial([0,1],:s)
sys1 = rtf(b, d, var=:s)
sys2 = rtf(Polynomial([b]), Polynomial([d]), Ts = 0, var = :s)

@test sys == sys1 && sys == sys2
@test sys ≈ sys1  && sys ≈ sys2

# one polynomial
a = 1; b = 2; c = 3; d = 4;
sys = rtf(Polynomial([b, a],:s), Ts = 0)
s = Polynomial([0,1],:s)
sys1 = rtf(a*s+b,Ts=0)
sys2 = rtf(Polynomial([b, a]), Ts = 0, var = :s)

@test sys == sys1 && sys == sys2
@test sys ≈ sys1  && sys ≈ sys2

a = 1; b = 2; c = 3; d = 4;
sys = rtf(Polynomial([b],:s), Ts = 0)
s = Polynomial([0,1],:s)
sys1 = rtf(b,Ts=0, var=:s)
sys2 = rtf(Polynomial([b]), Ts = 0, var = :s)

@test sys == sys1 && sys == sys2
@test sys ≈ sys1  && sys ≈ sys2

# one rational transfer function
a = 1; b = 2; c = 3; d = 4;
sys = rtf(Polynomial([b, a],:s), Polynomial([d, c],:s), Ts = 0)
z = Polynomial([0,1],:z)
sys1 = rtf((a*z+b)/(c*z+d),Ts=1)

@test rtf(sys,Ts=1,var=:z) == sys1

# Special constructors

@time t = rtf('s')
@test t == rtf("s") &&  t == rtf(:s) && t.var == :s && t.Ts == 0
@time t = rtf('z')
@test t == rtf("z") &&  t == rtf(:z) && t.var == :z && t.Ts == -1
@time t = rtf('z',Ts=2)
@test t == rtf("z",Ts=2) &&  t == rtf(:z,Ts=2) && t.var == :z && t.Ts == 2
@time t = rtf('λ')
@test t == rtf("λ") &&  t == rtf(:λ) && t.var == :λ

a = 1; b = 2; c = 3; d = 4;
sys = rtf(Polynomial([b, a],:s), Polynomial([d, c],:s), Ts = 0)
s = rtf('s')
@test sys == (a*s+b)/(c*s+d)

a = 1; b = 2; c = 3; d = 4;
sysd = rtf(Polynomial([b, a],:z), Polynomial([d, c],:z), Ts = 1)
z = rtf('z',Ts = 1)
@test sysd == (a*z+b)/(c*z+d)

# zpk
a = 1; b = 2; c = 3; d = 4;
@time sys = rtf(Polynomial([b, a],:s), Polynomial([d, c],:s), Ts = 0)
@time sys1 = rtf([-b/a],[-d/c],a/c,Ts = 0, var = :s)
@test sys ≈ sys1 && sys1.zeros ≈ [-b/a] && sys1.poles ≈ [-d/c] &&
      sys1.gain ≈ a/c && sys1.Ts == 0 && sys.var == :s

zer = eigvals(rand(5,5)) #fail
pol = eigvals(rand(7,7))
k = .5
@time sys = rtf(k*fromroots(zer),fromroots(pol),var=:s)
@time sys1 = rtf(zer,pol,k,Ts = 0)
@test sys ≈ sys1

# operations
@time r1 = rtf(Polynomial(1), Polynomial([-1, 1])) # r1 =  1/(x-1)
@time r2 = rtf(Polynomial(-1.), Polynomial([1, 1])) # r2 = -1/(x+1)
@time r3 = rtf(Polynomial(1), Polynomial([-1, 1]), var = :s) # r3 =  1/(s-1)
@time r4 = rtf(Polynomial(1), Polynomial([-1, 1]), Ts=1, var = :z) # r4 =  1/(z-1)

@test_throws ArgumentError r3+r4
@test_throws ArgumentError r3*r4
@test_throws ArgumentError rtf(Polynomial(1), Polynomial(0))

@test r1+r2 ≈ rtf(Polynomial(2), Polynomial([-1, 0, 1]))
@test r1-r2 ≈ rtf(Polynomial([0, 2]), Polynomial([-1, 0, 1]))
@test r1*r2 ≈ rtf(Polynomial(-1), Polynomial([-1, 0, 1]))
@test r1/r2 ≈ inv(r2/r1) ≈ -rtf(Polynomial([1, 1]), Polynomial([-1, 1]))

# with polynomials
@time p1 = Polynomial([1, 1], :x)     # p1 = (x+1)
@time p2 = Polynomial([1, 1], :s)     # p2 = (s+1)
r1 = rtf(p1) # r1 = (x+1)/1

@test p1 == r1 == p1
@test p2 ≠ r1 ≠ p2

@test p1 ≈ r1 ≈ p1
@test !(p2 ≈ r1)
@test !(r1 ≈ p2)

@time r1 = rtf(Polynomial([0, 1, 1]), Polynomial([-1, 1]))

@test p1+r1 ≈ r1+p1 ≈ rtf(Polynomial([-1, 1, 2]), Polynomial([-1, 1]))
@test r1-p1 ≈ -(p1-r1) ≈ rtf(Polynomial([1, 1]),Polynomial([-1, 1]))
@test r1*p1 ≈ p1*r1 ≈ rtf(Polynomial([0, 1, 2, 1]), Polynomial([-1, 1]))
@test r1/p1 ≈ inv(p1/r1) ≈ rtf(Polynomial([0, 1]), Polynomial([-1, 1]))

# with numbers
n   = 3.
@time r1  = rtf(Polynomial(n*[1, 1]), Polynomial([1, 1]))

@test n == r1 == n
@test n+1 ≠ r1 ≠ n+2

@test n ≈ r1 ≈ n
@test n + 0*1im ≈ r1 ≈ n + 0*1im

@time r1  = rtf(Polynomial([-1, 1]), Polynomial([1, 1])) # r1 = (x-1)/(x+1)

@test r1+n ≈ n+r1 ≈ rtf(Polynomial([2, 4]), Polynomial([1, 1]))
@test r1-n ≈ -(n-r1) ≈ rtf(Polynomial([-4, -2]), Polynomial([1, 1]))
@test r1*n ≈ n*r1 ≈ rtf(Polynomial([-3, 3]), Polynomial([1, 1]))
@test r1/n ≈ inv(n/r1) ≈ rtf(Polynomial([-1, 1]), Polynomial([3, 3])) ≈ rtf(Polynomial([-1, 1]/3), Polynomial([1, 1]))

# bilinear transformation
a = 1; b = 2; c = 3; d = 4;
s = rtf(:s)
@time sys = rtf(Polynomial([b, a],:s), Polynomial([d, c],:s), Ts = 0)
@time sysi = rtf(Polynomial([-b, d],:s), Polynomial([a, -c],:s))
@test (a*sysi+b)/(c*sysi+d) == s
@test confmap(sys,sysi) == s

s = rtf('s'); z = rtf('z'); λ = rtf(:λ);
@time g, gi = rtfbilin("cayley");
@test confmap(g,gi) == s
@test confmap(gi,g) == z

@time g, gi = rtfbilin("tustin",Ts = 0.01);
@test confmap(g,gi) == s
@test confmap(gi,g) == rtf('z',Ts = 0.01)

@time g, gi = rtfbilin("tustin",Ts = 0.01);
@test confmap(g,gi) == s
@test confmap(gi,g) == rtf('z',Ts = 0.01)

g, gi = rtfbilin("euler",Ts = 0.01);
@test confmap(g,gi) == s
@test confmap(gi,g) == rtf('z',Ts = 0.01)

g, gi = rtfbilin("beuler",Ts = 0.01);
@test confmap(g,gi) == s
@test confmap(gi,g) == rtf('z',Ts = 0.01)

@time g, gi = rtfbilin("Moebius",Ts = 0.01,Tsi = 0.02, a = 1, b = 2, c = 3, d = 4);
@test confmap(g,gi) == rtf('z',Ts = 0.02)
@test confmap(gi,g) == rtf('z',Ts = 0.01)

@time g, gi = rtfbilin("Moebius", a = 1, b = 2, c = 3, d = 4);
@test confmap(g,gi) == s
@test confmap(gi,g) == s

@test_throws ErrorException rtfbilin("ddd")

# normalization
a = 1; b = 2; c = 3; d = 4;
sys = rtf(Polynomial([b, a],:s), Polynomial([d, c],:s), Ts = 0)
@time sysn = normalize(sys)
@test sysn ≈ sys && sysn.den.coeffs[end] == 1

# simplify
@time r = rtf(Polynomial(rand(3)),Polynomial(rand(4)))
@time t = simplify(r/r)
@test t ≈ 1

# manipulation of rational matrices

# s and z rational functions
s = rtf('s'); z = rtf('z');     # define the complex variables s and z
@time Gc = [s^2 s/(s+1); 0 1/s]
@time Gd = [z^2 z/(z-2); 0 1/z]
@test all(transpose(transpose(Gc)) .== Gc)
@test transpose(s) == s

# # s and z polynomials
# s = Polynomial([0, 1],'s'); z = Polynomial([0, 1],'z');
# @time Gc1 = [s^2 s/(s+1); 0 1/s]
# @time Gd1 = [z^2 z/(z-2); 0 1/z]

# @test all( Gc .== rtf.(Gc1,Ts = Gc[1].Ts,var=:s))
# @test all( Gd .== rtf.(Gd1,Ts = Gd[1].Ts,var=:z))

@time g, ginv = rtfbilin("c2d")
@test all(confmap.(confmap.(Gc,[g]),[ginv]) .== Gc)
@test all(confmap.(confmap.(Gd,[ginv]),[g]) .== Gd)

@test all(confmap(confmap(Gc,g),ginv) .== Gc)
@test all(confmap(confmap(Gd,ginv),g) .== Gd)

# concatenations
s = rtf('s'); z = rtf('z');     # define the complex variables s and z
@time Gc = [s^2 s/(s+1); 0 1/s]
@time Gd = [z^2 z/(z-2); 0 1/z]
@test all([s^2 s/(s+1)] .== Gc[1:1,:])
@test all([s/(s+1); 1/s] .== Gc[:,2:2])
# @time Gc1 = [s^2 s/(s+1) I]
# @time Gc2 = [s^2 s/(s+1) 1]
# @test all(Gc1 .== Gc2)
# @test all([s^2; s/(s+1); I] .==[s^2; s/(s+1); 1])
@test all([[s^2; s/(s+1)]; I] .==[s^2; s/(s+1); 1])

@test all([z^2 z/(z-2)] .== Gd[1:1,:])
@test all([[z^2 z/(z-2)] I; I] .== [Gd[1:1,:] I;I])
@test all([z/(z-2); 1/z] .== Gd[:,2:2])
@test all([[[z/(z-2); 1/z];I] I] .== [[Gd[:,2:2];I] I])

@test_throws ArgumentError [Gc Gd]
@test_throws ArgumentError [Gc; Gd]
@time [Gc Gc]; [Gd;Gd]

@time Rc=[Gc[:,2:2];I]
@time Rd=[Gd[:,2:2];I]
@test_throws ArgumentError [Rc Rd]
@test_throws ArgumentError [Rc; Rd]

@test_throws ArgumentError [Rc Gc]
@test_throws ArgumentError [Rc; Gc]

@time Rc=[Gc[:,2:2] I]
@time Rd=[Gd[2:2,:]; I]

try
   @time [Rc Gc I]
   @test true
catch
   @test false
end

try
   @time [Rd; Gd; I;]
   @test true
catch
   @test false
end

try
   s = sprint(print,rtf(1));
   s = sprint(print,rtf('s'));
   s = sprint(print,rtf('z'));
   s = sprint(print,rtf('z',Ts=2));
   λ = rtf('λ'); s = sprint(print,λ/(λ+1));
   λ = rtf('λ',Ts=1); s = sprint(print,(λ^2+1)/5);
   sprint(print,Gc);
   sprint(print,Gd);
   @test true
catch
   @test false
end

end #test

end #module


text:
f(x)=sqrt(x)
true

max_stars_repo_name: Alfheim/node-julia

max_issues_repo_name: Alfheim/node-julia

max_forks_repo_name: Alfheim/node-julia


text:
module VariantVisualization

using DataFrames
using PlotlyJS
using ORCA
using GeneticVariation
using ArgParse
using DelimitedFiles
using Test

export
    test_parse_main,
    format_reader,
    load_vcf,
    clean_column1!,
    load_siglist,
    #sig_list_vcf_filter,
    #genomic_range_vcf_filter,
    #load_sort_phenotype_matrix,
    #reorder_columns,
    select_columns,
    #genotype_heatmap2,
    genotype_heatmap2_new_legend,
    dp_heatmap2,
    avg_sample_dp_scatter,
    avg_variant_dp_line_chart,
    read_depth_threshhold,
    list_variant_positions_low_dp,
    list_sample_names_low_dp,
    avg_dp_variant,
    avg_dp_samples,
    jupyter_main,
    save_numerical_array,
    io_pass_filter,
    io_sig_list_vcf_filter,
    io_genomic_range_vcf_filter,
    generate_genotype_array,
    define_geno_dict,
    translate_genotype_to_num_array,
    translate_readdepth_strings_to_num_array,
    genotype_heatmap_with_groups,
    dp_heatmap2_with_groups,
    returnXY_column1!,
    pass_genomic_range_siglist_filter,
    pass_genomic_range_filter,
    pass_siglist_filter,
    genomic_range_siglist_filter,
    get_sample_names,
    sortcols_by_phenotype_matrix,
    find_group_label_indices,
    checkfor_outputdirectory,
    combined_all_genotype_array_functions,
    combined_all_read_depth_array_functions,
    combined_all_read_depth_array_functions_for_avg_dp,
    generate_chromosome_positions_for_hover_labels,
    clean_column1_chr,
    clean_column1_siglist!,
    process_plot_inputs,
    process_plot_inputs_for_grouped_data,
    returnXY_column1_siglist!,
    chromosome_label_generator,
    add_pheno_matrix_to_gt_data_for_plotting,
    add_pheno_matrix_to_dp_data_for_plotting,
    generate_hover_text_array,
    generate_hover_text_array_grouped,
    save_graphic,
    build_set_from_list,
    generate_legend_increments_ungrouped,
    generate_legend_increments_grouped

include("vcf_utils_complete.jl")
include("plot_utils.jl")
include("new_notebook_utils.jl")

end # module


text:
module TestReduce

using Test
using LazyGroupBy
using LazyGroupBy: foldxl, foldxt
using Transducers: Map

extrema_rf((min1, max1), (min2, max2)) = (min(min1, min2), max(max1, max2))

@testset "extrema" begin
    @test mapreduce.(
        x -> (x, x),
        extrema_rf,
        grouped(isodd, [0, 7, 3, 1, 5, 9, 4, 3, 0, 5]),
    ) == Dict(false => (0, 4), true => (1, 9))

    @testset for fold in [foldl, foldxl, foldxt, reduce]
        foldl = nothing

        @test fold.(
            extrema_rf,
            Map(x -> (x, x)),
            grouped(isodd, [0, 7, 3, 1, 5, 9, 4, 3, 0, 5]),
        ) == Dict(false => (0, 4), true => (1, 9))
    end
end

end  # module


text:
#! Here come all the functions for the model with trade in structured data and non-rival data

"""
    Solve the problem of an unconstrained firm that can choose D_S freely (should stay positive).
"""
function solve_firm_problem_unc(par::Pars_v3; type::Symbol)
    f(x) = residuals_unc(x; par, type)
    if type == :BS
        initial_x = [1.090961407034043,
            1.7305342659571796,
            0.8017602275570367,
            0.846302424724036]
    elseif type == :BU
        initial_x = [1.2797075968276403,
            1.7305374802420503,
            0.7399594318155204,
            0.7810681256207107]
    elseif type == :LS
        initial_x = [1.0909610262794553,
            -0.11836701171722647,
            0.8017601473223311,
            0.8463025640344691]
    elseif type == :LU
        initial_x = [1.2797075968276403,
            1.7305374802420503,
            0.7399594318155204,
            0.7810681256207107]
    end
    #res = optimize(f, initial_x, Optim.Options(g_tol = 1e-14))
    res = optimize(f, initial_x)
    # try
    #     @assert Optim.converged(res) == true
    # catch
    #     τ = par.τ
    #     println("Solution may be inaccurate for $τ and $type")
    # end
    return res
end

function solve_firm_problem_unc_nl(par::Pars_v3; type::Symbol, sec_period::Symbol=:no)
    res = zeros(4)
    f!(x) = residuals_unc_nl(x, res; par, type, sec_period)

    initial_x = [1.090961407034043,
        1.7305342659571796,
        0.8017602275570367,
        0.846302424724036]

    #res = optimize(f, initial_x, Optim.Options(g_tol = 1e-14))
    res = nlsolve(f!, initial_x)
    # try
    #     @assert Optim.converged(res) == true
    # catch
    #     τ = par.τ
    #     println("Solution may be inaccurate for $τ and $type")
    # end
    return res
end

function solve_firm_problem_unc_nl2(par::Pars_v3; type::Symbol, sec_period::Symbol=:no)
    if sec_period == :yes
        res = zeros(4)
        initial_x = [1.090961407034043,
            1.7305342659571796,
            0.8017602275570367,
            0.846302424724036]
    elseif sec_period == :no
        res = zeros(8)
        initial_x = [1.090961407034043,
            1.7305342659571796,
            0.8017602275570367,
            0.846302424724036,
            1.090961407034043,
            1.7305342659571796,
            0.8017602275570367,
            0.846302424724036]
    end
    f!(x) = residuals_unc_nl2(x, res; par, type, sec_period)
    res = nlsolve(f!, initial_x)
    return res
end

"""
    Solve the problem of a constrained firm with D_S = 0 (would like to buy internal data).
"""
function solve_firm_problem_con_nl(par::Pars_v3; type::Symbol)
    #f(x) = residuals_con(x; par, type)
    res = zeros(3)
    f!(x) = residuals_con_nl(x, res; par, type)
    initial_x = [0.5, 0.5, 0.5]
    #res = optimize(f, initial_x, Optim.Options(g_tol = 1e-14))
    res = nlsolve(f!, initial_x)
    #@assert Optim.converged(res) == true
    return res
end

function solve_firm_problem_con(par::Pars_v3; type::Symbol)
    f(x) = residuals_con(x; par, type)
    initial_x = [0.5, 0.5, 0.5]
    #res = optimize(f, initial_x, Optim.Options(g_tol = 1e-14))
    res = optimize(f, initial_x)
    #@assert Optim.converged(res) == true
    return res
end

"""
    Solve the problem of the firm.
"""
function solve_firm_problem(par::Pars_v3; type::Symbol)
    sol = solve_firm_problem_unc(par; type)
    l_P, l_G, l, D_E, D_I, D = get_solution(sol; par, type)

    D_S = sold_data(D_I, l_G; par, type)

    if D_S < 0
        sol_con = solve_firm_problem_con(par; type)
        return sol_con
    elseif D_S >= 0
        return sol
    else
        error("I never should be here")
    end
end

function solve_firm_problem_nl(par::Pars_v3; type::Symbol, sec_period::Symbol=:no)
    sol = solve_firm_problem_unc_nl(par; type, sec_period)
    l_P, l_G, l, D_E, D_I, D = get_solution_nl(sol; par, type)

    D_S = sold_data(D_I, l_G; par, type)
    #! Idea: Solve first unconstrained problem. If firm tries to buy internal data, constrain to not sharing data.
    if D_S < 0
        sol_con = solve_firm_problem_con_nl(par; type)
        return sol_con
    elseif D_S >= 0
        return sol
    else
        error("I never should be here")
    end
end

# function firm_profits(x; par::Pars_v3, type::Symbol)
#     l_P = x[1]^2
#     l_G = x[2]^2
#     l   = x[3]^2
#     D_E = x[4]^2
#     D_S = x[5]

#     # All other variables we can directly constrain to be

# end

"""
    Returns the labor market clearing condition residuals
"""
function labor_market_clearing(l_BS, l_LS, l_LU, l_BU; par)
    w_BS = weight_of_type(:BS; par)
    w_LS = weight_of_type(:LS; par)
    w_LU = weight_of_type(:LU; par)
    w_BU = weight_of_type(:BU; par)

    #! Labor of each type is fixed to one for the moment.
    return l_BS * w_BS + l_LS * w_LS + l_LU * w_LU + l_BU * w_BU - 1
end

"""
    Aggregate Production Function.
"""
function aggregate_output(Y_BS, Y_LS, Y_LU, Y_BU; par::Pars_v3)
    w_BS = weight_of_type(:BS; par)
    w_LS = weight_of_type(:LS; par)
    w_LU = weight_of_type(:LU; par)
    w_BU = weight_of_type(:BU; par)

    c = (par.σ - 1) / par.σ

    if par.compl_prod == :yes
        # Aggregat Production function with strategic complementarities in production
        # May lead to instability, could not solve the model so far with this specification.
        Y_agg = (w_BS * Y_BS^c + w_LS * Y_LS^c + w_LU * Y_LU^c + w_BU * Y_BU^c)^(1 / c)
    elseif par.compl_prod == :no
        # Aggregat Production function without strategic complementarities in production.
        # Demand for one good is independent of demand for other goods.
        Y_agg = (w_BS * Y_BS^c + w_LS * Y_LS^c + w_LU * Y_LU^c + w_BU * Y_BU^c) / c
    else
        error("I should never be here!")
    end
end

function aggregate_labor(l_BS, l_LS, l_LU, l_BU; par::Pars_v3)
    w_BS = weight_of_type(:BS; par)
    w_LS = weight_of_type(:LS; par)
    w_LU = weight_of_type(:LU; par)
    w_BU = weight_of_type(:BU; par)

    return w_BS * l_BS + w_LS * l_LS + w_LU * l_LU + w_BU * l_BU
end

sol_f_problem(par::Pars_v3; type::Symbol) = get_solution(solve_firm_problem(par; type); par, type)
sol_f_problem_nl(par::Pars_v3; type::Symbol) = get_solution_nl(solve_firm_problem_nl(par; type); par, type)

function get_ss_solution(sol)
    Y = sol.minimizer[1]^2
    w = sol.minimizer[2]^2
    w_P = sol.minimizer[3]^2
    w_G = sol.minimizer[4]^2
    p_D = sol.minimizer[5]^2

    return Y, w, w_P, w_G, p_D
end

function find_steady_state(par::Pars_v3; tol=1e-10, show_steps::Symbol=:yes)
    diff = 10
    counter = 0
    while diff > tol
        counter += 1
        p_D_sol = find_p_D(par)
        @set! par.p_D = p_D_sol

        l_P_BS, l_G_BS, l_BS, D_E_BS, D_I_BS, D_BS = sol_f_problem_nl(par; type=:BS)
        l_P_LS, l_G_LS, l_LS, D_E_LS, D_I_LS, D_LS = sol_f_problem_nl(par; type=:LS)
        l_P_LU, l_G_LU, l_LU, D_E_LU, D_I_LU, D_LU = sol_f_problem_nl(par; type=:LU)
        l_P_BU, l_G_BU, l_BU, D_E_BU, D_I_BU, D_BU = sol_f_problem_nl(par; type=:BU)

        l_P = [l_P_BS, l_P_LS, l_P_LU, l_P_BU]
        l_G = [l_G_BS, l_G_LS, l_G_LU, l_G_BU]
        l = [l_BS, l_LS, l_LU, l_BU]
        D_E = [D_E_BS, D_E_LS, D_E_LU, D_E_BU]
        D_I = [D_I_BS, D_I_LS, D_I_LU, D_I_BU]
        D = [D_BS, D_LS, D_LU, D_BU]

        #D_BS = bundle(D_I_BS, D_E_BS; par)
        #D_LS = bundle(D_I_LS, D_E_LS; par)
        #D_LU = bundle(D_I_LU, D_E_LU; par)
        #D_BU = bundle(D_I_BU, D_E_BU; par)

        Y_BS = firm_output(l_P_BS, D_BS, l_BS; par, type=:BS)
        Y_LS = firm_output(l_P_LS, D_LS, l_LS; par, type=:LS)
        Y_LU = firm_output(l_P_LU, D_LU, l_LU; par, type=:LU)
        Y_BU = firm_output(l_P_BU, D_BU, l_BU; par, type=:BU)

        Π_BS = firm_revenue(l_P_BS, D_BS, l_BS; par, type=:BS)
        Π_LS = firm_revenue(l_P_LS, D_LS, l_LS; par, type=:LS)
        Π_LU = firm_revenue(l_P_LU, D_LU, l_LU; par, type=:LU)
        Π_BU = firm_revenue(l_P_BU, D_BU, l_BU; par, type=:BU)

        K_BS = knowledge(l_P_BS, D_BS; par, type=:BS)
        K_LS = knowledge(l_P_LS, D_LS; par, type=:LS)
        K_LU = knowledge(l_P_LU, D_LU; par, type=:LU)
        K_BU = knowledge(l_P_BU, D_BU; par, type=:BU)

        Π = [Π_BS, Π_LS, Π_LU, Π_BU]

        Y = aggregate_output(Y_BS, Y_LS, Y_LU, Y_BU; par)
        w = wage_int_goods(Π; par)
        w_G = wage_gen_data(Π, D, l_G, D_I; par)
        w_P = wage_gen_proc(Π, D; par)

        if par.α_Y == 0
            res = zeros(3)

            res[1] = w - par.w
            res[2] = 10 * (w_G - par.w_G) # Make this constraint more important?
            res[3] = w_P - par.w_P

            diff = sum(res .^ 2)
        else
            par.α_Y != 0.0
            res = zeros(4)

            res[1] = Y - par.Y
            res[2] = w - par.w
            res[3] = 10 * (w_G - par.w_G) # Make this constraint more important?
            res[4] = w_P - par.w_P

            diff = sum(res .^ 2)
        end

        if show_steps == :yes
            #@show Y
            #@show w
            #@show w_G
            #@show w_P
            println("$counter. step, res =  $diff")
        end

        #! Try to make the convergence a bit more smooth, choose a small lambda
        @set! par.Y = par.λ * Y + (1 - par.λ) * par.Y
        @set! par.w = par.λ * w + (1 - par.λ) * par.w
        @set! par.w_G = par.λ * w_G + (1 - par.λ) * w_G
        @set! par.w_P = par.λ * w_P + (1 - par.λ) * w_P

        if diff < tol || counter > 100

            #! Check the solution with other equations! Did I write the equations for wages correctly?
            # try
            #     @assert abs(aggregate_labor(l_P_BS, l_P_LS, l_P_LU, l_P_BU; par) - par.P) < 1e-2
            #     @assert abs(aggregate_labor(l_G_BS, l_G_LS, l_G_LU, l_G_BU; par) - par.G) < 1e-2
            #     @assert abs(aggregate_labor(l_BS, l_LS, l_LU, l_BU; par) - par.L) < 1e-2
            # catch
            #     @show aggregate_labor(l_P_BS, l_P_LS, l_P_LU, l_P_BU; par)
            #     @show aggregate_labor(l_G_BS, l_G_LS, l_G_LU, l_G_BU; par)
            #     @show aggregate_labor(l_BS, l_LS, l_LU, l_BU; par)

            #     error("Labor Market Clearing did not work!")
            # end

            #! If convergence has been achieved, save the solution
            D_S_BS = sold_data(D_I_BS, l_G_BS; par, type=:BS)
            D_S_LS = sold_data(D_I_LS, l_G_LS; par, type=:LS)
            D_S_LU = sold_data(D_I_LU, l_G_LU; par, type=:LU)
            D_S_BU = sold_data(D_I_BU, l_G_BU; par, type=:BU)

            D_G_BS = data_gen(l_G_BS; par, type=:BS)
            D_G_LS = data_gen(l_G_LS; par, type=:LS)
            D_G_LU = data_gen(l_G_LU; par, type=:LU)
            D_G_BU = data_gen(l_G_BU; par, type=:BU)

            prof_BS = firm_profits(l_BS, l_P_BS, l_G_BS, D_S_BS, D_E_BS; par, type=:BS)
            prof_LS = firm_profits(l_LS, l_P_LS, l_G_LS, D_S_LS, D_E_LS; par, type=:LS)
            prof_LU = firm_profits(l_LU, l_P_LU, l_G_LU, D_S_LU, D_E_LU; par, type=:LU)
            prof_BU = firm_profits(l_BU, l_P_BU, l_G_BU, D_S_BU, D_E_BU; par, type=:BU)

            agg_D_S = D_S_BS + D_S_LS + D_S_LU + D_S_BU
            agg_D_G = D_G_BS + D_G_LS + D_G_LU + D_G_BU
            agg_D = D_BS + D_LS + D_LU + D_BU

            @set! par.D_G_SS = agg_D_G
            @set! par.D_S_SS = agg_D_S
            @set! par.Ω = data_multiplier(agg_D_G, agg_D_S; par)
            @set! par.Ω_bundle = data_multiplier_bundle(agg_D, agg_D_G)

            @set! par.prof_LU_ss = prof_LU # Labor good production small customer base, unsophisticated
            @set! par.prof_BU_ss = prof_BU # Labor good production big customer base, unsophisticated
            @set! par.prof_LS_ss = prof_LS # Labor good production small customer base, sophisticated
            @set! par.prof_BS_ss = prof_BS # Labor good production big customer base, sophisticated

            @set! par.l_LU_ss = l_LU # Labor good production small customer base, unsophisticated
            @set! par.l_BU_ss = l_BU # Labor good production big customer base, unsophisticated
            @set! par.l_LS_ss = l_LS # Labor good production small customer base, sophisticated
            @set! par.l_BS_ss = l_BS # Labor good production big customer base, sophisticated

            @set! par.l_P_LU_ss = l_P_LU # Labor good production small customer base, unsophisticated
            @set! par.l_P_BU_ss = l_P_BU # Labor good production big customer base, unsophisticated
            @set! par.l_P_LS_ss = l_P_LS # Labor good production small customer base, sophisticated
            @set! par.l_P_BS_ss = l_P_BS # Labor good production big customer base, sophisticated

            @set! par.l_G_LU_ss = l_G_LU # Labor good production small customer base, unsophisticated
            @set! par.l_G_BU_ss = l_G_BU # Labor good production big customer base, unsophisticated
            @set! par.l_G_LS_ss = l_G_LS # Labor good production small customer base, sophisticated
            @set! par.l_G_BS_ss = l_G_BS # Labor good production big customer base, sophisticated

            @set! par.D_S_LU_ss = D_S_LU # data sharing small customer base, unsophisticated
            @set! par.D_S_BU_ss = D_S_BU # data sharing big customer base, unsophisticated
            @set! par.D_S_LS_ss = D_S_LS # data sharing small customer base, sophisticated
            @set! par.D_S_BS_ss = D_S_BS # data sharing big customer base, sophisticated

            @set! par.D_G_LU_ss = D_G_LU # data generated small customer base, unsophisticated
            @set! par.D_G_BU_ss = D_G_BU # data generated big customer base, unsophisticated
            @set! par.D_G_LS_ss = D_G_LS # data generated small customer base, sophisticated
            @set! par.D_G_BS_ss = D_G_BS # data generated big customer base, sophisticated

            @set! par.D_LU_ss = D_LU # data generated small customer base, unsophisticated
            @set! par.D_BU_ss = D_BU # data generated big customer base, unsophisticated
            @set! par.D_LS_ss = D_LS # data generated small customer base, sophisticated
            @set! par.D_BS_ss = D_BS # data generated big customer base, sophisticated

            @set! par.D_E_LU_ss = D_E_LU # data generated small customer base, unsophisticated
            @set! par.D_E_BU_ss = D_E_BU # data generated big customer base, unsophisticated
            @set! par.D_E_LS_ss = D_E_LS # data generated small customer base, sophisticated
            @set! par.D_E_BS_ss = D_E_BS # data generated big customer base, sophisticated

            @set! par.D_I_LU_ss = D_I_LU # data generated small customer base, unsophisticated
            @set! par.D_I_BU_ss = D_I_BU # data generated big customer base, unsophisticated
            @set! par.D_I_LS_ss = D_I_LS # data generated small customer base, sophisticated
            @set! par.D_I_BS_ss = D_I_BS # data generated big customer base, sophisticated

            @set! par.K_LU_ss = K_LU # knowledge employed small customer base, unsophisticated
            @set! par.K_BU_ss = K_BU # knowledge employed big customer base, unsophisticated
            @set! par.K_LS_ss = K_LS # knowledge employed small customer base, sophisticated
            @set! par.K_BS_ss = K_BS # knowledge employed big customer base, sophisticated

            @set! par.Firm_sol = Save_Firm_Solution(l, l_P, l_G, D_I, D_E, D; par)
        end
    end
    return par
end

"""
    Solve all firm problems
"""
function find_steady_state_PE(par::Pars_v3; solve_price::Symbol = :yes)
    if solve_price == :yes
        p_D_sol = find_p_D(par)
        @set! par.p_D = p_D_sol
    end

    l_P_BS, l_G_BS, l_BS, D_E_BS, D_I_BS, D_BS = sol_f_problem_nl(par; type = :BS)
    l_P_LS, l_G_LS, l_LS, D_E_LS, D_I_LS, D_LS = sol_f_problem_nl(par; type = :LS)
    l_P_LU, l_G_LU, l_LU, D_E_LU, D_I_LU, D_LU = sol_f_problem_nl(par; type = :LU)
    l_P_BU, l_G_BU, l_BU, D_E_BU, D_I_BU, D_BU = sol_f_problem_nl(par; type = :BU)

    D_S_BS = sold_data(D_I_BS, l_G_BS; par, type = :BS)
    D_S_LS = sold_data(D_I_LS, l_G_LS; par, type = :LS)
    D_S_LU = sold_data(D_I_LU, l_G_LU; par, type = :LU)
    D_S_BU = sold_data(D_I_BU, l_G_BU; par, type = :BU)

    D_G_BS = data_gen(l_G_BS; par, type = :BS)
    D_G_LS = data_gen(l_G_LS; par, type = :LS)
    D_G_LU = data_gen(l_G_LU; par, type = :LU)
    D_G_BU = data_gen(l_G_BU; par, type = :BU)

    K_BS = knowledge(l_P_BS, D_BS; par, type = :BS)
    K_LS = knowledge(l_P_LS, D_LS; par, type = :LS)
    K_LU = knowledge(l_P_LU, D_LU; par, type = :LU)
    K_BU = knowledge(l_P_BU, D_BU; par, type = :BU)

    prof_BS = firm_profits(l_BS, l_P_BS, l_G_BS, D_S_BS, D_E_BS; par, type = :BS)
    prof_LS = firm_profits(l_LS, l_P_LS, l_G_LS, D_S_LS, D_E_LS; par, type = :LS)
    prof_LU = firm_profits(l_LU, l_P_LU, l_G_LU, D_S_LU, D_E_LU; par, type = :LU)
    prof_BU = firm_profits(l_BU, l_P_BU, l_G_BU, D_S_BU, D_E_BU; par, type = :BU)

    @set! par.l_LU_ss = l_LU # Labor good production small customer base, unsophisticated
    @set! par.l_BU_ss = l_BU # Labor good production big customer base, unsophisticated
    @set! par.l_LS_ss = l_LS # Labor good production small customer base, sophisticated
    @set! par.l_BS_ss = l_BS # Labor good production big customer base, sophisticated

    @set! par.prof_LU_ss = prof_LU # Labor good production small customer base, unsophisticated
    @set! par.prof_BU_ss = prof_BU # Labor good production big customer base, unsophisticated
    @set! par.prof_LS_ss = prof_LS # Labor good production small customer base, sophisticated
    @set! par.prof_BS_ss = prof_BS # Labor good production big customer base, sophisticated

    @set! par.l_P_LU_ss = l_P_LU # Labor good production small customer base, unsophisticated
    @set! par.l_P_BU_ss = l_P_BU # Labor good production big customer base, unsophisticated
    @set! par.l_P_LS_ss = l_P_LS # Labor good production small customer base, sophisticated
    @set! par.l_P_BS_ss = l_P_BS # Labor good production big customer base, sophisticated

    @set! par.l_G_LU_ss = l_G_LU # Labor good production small customer base, unsophisticated
    @set! par.l_G_BU_ss = l_G_BU # Labor good production big customer base, unsophisticated
    @set! par.l_G_LS_ss = l_G_LS # Labor good production small customer base, sophisticated
    @set! par.l_G_BS_ss = l_G_BS # Labor good production big customer base, sophisticated

    @set! par.D_S_LU_ss = D_S_LU # data sharing small customer base, unsophisticated
    @set! par.D_S_BU_ss = D_S_BU # data sharing big customer base, unsophisticated
    @set! par.D_S_LS_ss = D_S_LS # data sharing small customer base, sophisticated
    @set! par.D_S_BS_ss = D_S_BS # data sharing big customer base, sophisticated

    @set! par.D_G_LU_ss = D_G_LU # data generated small customer base, unsophisticated
    @set! par.D_G_BU_ss = D_G_BU # data generated big customer base, unsophisticated
    @set! par.D_G_LS_ss = D_G_LS # data generated small customer base, sophisticated
    @set! par.D_G_BS_ss = D_G_BS # data generated big customer base, sophisticated

    @set! par.D_LU_ss = D_LU # data generated small customer base, unsophisticated
    @set! par.D_BU_ss = D_BU # data generated big customer base, unsophisticated
    @set! par.D_LS_ss = D_LS # data generated small customer base, sophisticated
    @set! par.D_BS_ss = D_BS # data generated big customer base, sophisticated

    @set! par.D_E_LU_ss = D_E_LU # data generated small customer base, unsophisticated
    @set! par.D_E_BU_ss = D_E_BU # data generated big customer base, unsophisticated
    @set! par.D_E_LS_ss = D_E_LS # data generated small customer base, sophisticated
    @set! par.D_E_BS_ss = D_E_BS # data generated big customer base, sophisticated

    @set! par.D_I_LU_ss = D_I_LU # data generated small customer base, unsophisticated
    @set! par.D_I_BU_ss = D_I_BU # data generated big customer base, unsophisticated
    @set! par.D_I_LS_ss = D_I_LS # data generated small customer base, sophisticated
    @set! par.D_I_BS_ss = D_I_BS # data generated big customer base, sophisticated

    @set! par.K_LU_ss = K_LU # knowledge employed small customer base, unsophisticated
    @set! par.K_BU_ss = K_BU # knowledge employed big customer base, unsophisticated
    @set! par.K_LS_ss = K_LS # knowledge employed small customer base, sophisticated
    @set! par.K_BS_ss = K_BS # knowledge employed big customer base, sophisticated

    return par
end

"""
    Find first the price for data that clears the data market.
"""
function find_p_D(par::Pars_v3; verbose = false)
    f(p_D) = residuals_p_D(p_D^2; par)
    #f(p_D) = residuals_p_D(p_D; par)
    find_zero(f, 0.3, verbose = verbose)^2
    #find_zero(f, [0.001, 2], verbose = verbose, xatol = 1e-10, xrtol = 1e-10, atol = 1e-10, rtol = 1e-10)
end

function comp_statics(τ_set, ν_set, A_G_common_set, A_P_common_set; tol = 1e-10, ζ = 1.0, γ_S = 0.2, γ_U = 0.4)
    grid = gridmake(τ_set, ν_set, A_G_common_set, A_P_common_set)
    sol = Array{Any}(undef, 0)
    for i in 1:size(grid)[1]
        @show i
        @show grid[i, :]
        if i == 1
            push!(sol, find_steady_state(Pars_v3(τ = grid[i, 1], ν = grid[i, 2], A_G_common = grid[i, 3], A_P_common = grid[i, 4], ζ = ζ, γ_S = γ_S, γ_U = γ_U); tol))
        else
            #try
            push!(sol, find_steady_state(Pars_v3(τ = grid[i, 1], ν = grid[i, 2], A_G_common = grid[i, 3], A_P_common = grid[i, 4], Y = sol[i-1].Y, w = sol[i-1].w, w_P = sol[i-1].w_P, w_G = sol[i-1].w_G, ζ = ζ, γ_S = γ_S, γ_U = γ_S); tol))
            #catch
            #    push!(sol, find_steady_state(Pars_v3(τ = grid[i, 1], ν = grid[i, 2], A_G_common = grid[i, 3], A_P_common = grid[i, 4]); tol))
            #end
        end
    end
    return sol
end

function find_value_func(par::Pars_v3; τ)
    @set! par.τ = τ
    @set! par.second_period = :yes
    K_last_set = 0.0:0.1:2
    D_last_set = 0.0:0.1:2
    grid = gridmake(K_last_set, D_last_set)
    V = zeros(size(grid)[1], 4)
    par = find_steady_state(par, tol = 1e-12)
    i = 0
    f(K_last, D_last) = begin
        @set! par.D_last_LU = D_last # data generated small customer base, unsophisticated
        @set! par.D_last_BU = D_last # data generated big customer base, unsophisticated
        @set! par.D_last_LS = D_last # data generated small customer base, sophisticated
        @set! par.D_last_BS = D_last # data generated big customer base, sophisticated

        @set! par.K_last_LU = K_last # knowledge employed small customer base, unsophisticated
        @set! par.K_last_BU = K_last # knowledge employed big customer base, unsophisticated
        @set! par.K_last_LS = K_last # knowledge employed small customer base, sophisticated
        @set! par.K_last_BS = K_last # knowledge employed big customer base, sophisticated

        par = find_steady_state_PE(par, solve_price = :no)

        par.prof_BS_ss
    end

    A = [f(x, y) for x in K_last_set, y in D_last_set]

    # linear interpolation
    return LinearInterpolation((K_last_set, D_last_set), A)

    # cubic spline interpolation
    #interp_cubic = CubicSplineInterpolation((xs, ys), A)
    #DataFrame(K_last = grid[:, 1], D_last = grid[:, 2], V_BS = V[:, 1], V_LS = V[:, 2], V_LU = V[:, 3], V_BU = V[:, 4])
end


text:
using Base.Test

# see https://github.com/JuliaQuant/MarketData.jl
# and https://github.com/JuliaStats/TimeSeries.jl

using TALib
using TALib: basepath

using TimeSeries

@testset "level 2 timearrays" begin
    epsilon = 1e-10
    eps_price = 1e-6

    filename = joinpath(basepath(), "ford_2012.csv")
    ohlcv = readtimearray(filename)

    indic = MA(ohlcv)
    @test indic["Real"][end].values[1] ≈ 11.546 atol=eps_price

    indic = MA(ohlcv, price=:Open)
    @test indic["Real"][end].values[1] ≈ 11.4613333 atol=eps_price

    indic = BBANDS(ohlcv)
    @test indic["UpperBand"][end].values[1] ≈ 13.13191 atol=eps_price
    @test indic["MiddleBand"][end].values[1] ≈ 12.75400 atol=eps_price
    @test indic["LowerBand"][end].values[1] ≈ 12.37609 atol=eps_price

end


text:
# import ..PolaronMakie
# using ..PolaronMakie

"""
This is a playground for testing implementations of the X(u) function. Still need to find a detailed expansions of ReX(u) that is not horrific numerically, if such a thing exists!
"""

using BigCombinatorics
using Combinatorics
using SpecialFunctions: gamma, beta
using ArbNumerics, Readables
using Plots
using Struve
using MathLink
using QuadGK

function L(ν, z)
    ν = BigFloat(ν)
    z = BigFloat(z)
    coefficient = (z / 2)^(ν + 1)
    total_sum = 0.0
    k = 0
    prev_sum = 0.0
    while true
        next_sum = (z / 2)^(2 * k) / (gamma(3/2 + k) * gamma(3/2 + k + ν))
        if next_sum == prev_sum
            break
        end
        total_sum += next_sum
        k += 1
        prev_sum = next_sum
    end
    coefficient * total_sum
end

function I(ν, z)
    ν = BigInt(ν)
    z = BigFloat(z)
    total_sum = 0.0
    k = 0
    prev_sum = 0.0
    while true
        next_sum = (z / 2)^(2 * k + ν) / (gamma(k + 1) * gamma(k + ν + 1))
        if next_sum == prev_sum
            break
        end
        total_sum += next_sum
        k += 1
        prev_sum = next_sum
    end
    total_sum
end

function I_minus_L(n, x, a)
    setprecision(BigFloat, digits)
    n = BigFloat(n)
    x = BigFloat(x)
    a = BigFloat(a)

    d = weval(W"N"(W"Abs"(W"Subtract"(W"BesselI"(n+1, sign(x)^(n + 1) * x * a), W"StruveL"(-n-1, sign(x)^(n) * x * a))), digits)).value
    d_split = split(d, r"[`^]")
    if length(d_split) == 2
        return parse(BigFloat, d_split[1])
    else
        return parse(BigFloat, d_split[1] * "e" * d_split[3])
    end
end

function ℑχ_0(Ω, α, v, w, N = 10)

    R = (v^2 - w^2) / (w^2 * v)
    coefficient = 2 / 3 * α * (v / w)^3

    total_sum = 0.0
    for n in 0:Int(N)
        total_sum += - coefficient * sqrt(π) / (gamma(-n - 1/2) * gamma(n + 1)) * (-2 * R)^n / BigCombinatorics.doublefactorial(2 * n + 1) * abs(Ω - 1 - n * v)^(n + 1/2) * exp(-R * abs(Ω - 1 - n * v)) * (1 + sign(Ω - 1 - n * v))
    end
    return total_sum
end

function ℑχ(Ω, β, α, v, w, N = 10)

    Ω = BigFloat(Ω)
    β = BigFloat(β)
    α = BigFloat(α)
    v = BigFloat(v)
    w = BigFloat(w)

    R = (v^2 - w^2) / (w^2 * v)
    a = sqrt(β^2 / 4 + R * β * coth(β * v / 2))
    b = R * β / sinh(β * v / 2)

    coefficient = 2 * α / (3 * sqrt(π)) * (v / w)^3 * β^(3/2) * sinh(β * Ω / 2) / sinh(β / 2)

    @show(Ω)

    total_sum = BigFloat(0.0)

    for n in 0:Int(N)

        n_coefficient = -2 * beta(-1/2 - n, n + 1)^(-1) * (-1)^n * b^n * (4 * a)^(-n - 1) * sqrt(π) / gamma(n + 3/2)

        if isodd(n)
            for k in -1:Int((n - 1) / 2)
                k_coefficient_inverse = ((n + 1) * beta(n - k + 1, k + 1))

                y = [(Ω + 1 + v * (n - 2 * k)), (Ω + 1 - v * (n - 2 * k)), (Ω - 1 + v * (n - 2 * k)), (Ω - 1 - v * (n - 2 * k))]

                for x in y
                    total_sum +=  coefficient * n_coefficient * besselk(Int(n + 1), abs(Float64(x)) * Float64(a)) * abs(x)^(n + 1) / k_coefficient_inverse
                    # total_sum += coefficient * n_coefficient * k_coefficient * besselk(ArbReal(n + 1), ArbReal(abs(x) * a)) * abs(x)^(n + 1)
                end
            end
        end

        if iseven(n)
            for k in -1:Int(n / 2 - 1)

                k_coefficient_inverse = ((n + 1) * beta(n - k + 1, k + 1))
                y = [(Ω + 1 + v * (n - 2 * k)), (Ω + 1 - v * (n - 2 * k)), (Ω - 1 + v * (n - 2 * k)), (Ω - 1 - v * (n - 2 * k))]

                for x in y
                    total_sum +=  coefficient * n_coefficient * besselk(Int(n + 1), abs(Float64(x)) * Float64(a)) * abs(x)^(n + 1) / k_coefficient_inverse
                    # total_sum += coefficient * n_coefficient * k_coefficient * besselk(ArbReal(n + 1), ArbReal(abs(x) * a)) * abs(x)^(n + 1)
                end

                k_even_coefficient_inverse = ((n + 1) * beta(n/2 + 1, n/2 + 1))
                y_even = [(Ω + 1), (Ω - 1)]

                for x in y_even
                    total_sum +=  coefficient * n_coefficient * besselk(Int(n + 1), abs(Float64(x)) * Float64(a)) * abs(x)^(n + 1) / k_even_coefficient_inverse
                    # total_sum += coefficient * n_coefficient * k_even_coefficient * besselk(ArbReal(n + 1), ArbReal(abs(x) * a)) * abs(x)^(n + 1)
                end
            end
        end
    end
    println("/n")
    return total_sum
end

function ℜχ(Ω, β, α, v, w, N = 10)
    Ω = BigFloat(Ω)
    β = BigFloat(β)
    α = BigFloat(α)
    v = BigFloat(v)
    w = BigFloat(w)

    R = (v^2 - w^2) / (w^2 * v)
    a = sqrt(β^2 / 4 + R * β * coth(β * v / 2))
    b = R * β / sinh(β * v / 2)

    @show(Ω)

    coefficient = 2 * α / (3 * sqrt(π)) * (v / w)^3 * β^(3/2) / sinh(β / 2)

    # total_sum = BigFloat(0.0)
    #
    # for n in 0:Int(N)
    #
    #     n_coefficient = sinh(β * Ω / 2) * (1 / (-1/2 * beta(n + 1, -n - 1/2))) * π *  b^n * (1/4) * (1/a^(n + 1)) * Factorial(n + 1) / Factorial(2 * (n + 1))
    #
    #     if isodd(n)
    #         for k in -1:Int((n - 1) / 2)
    #
    #             k_coefficient = Binomial(n, k) * (1 - (-1)^n)
    #
    #             y = [(Ω + 1 + v * (n - 2 * k)), (Ω + 1 - v * (n - 2 * k)), (Ω - 1 + v * (n - 2 * k)), (Ω - 1 - v * (n - 2 * k))]
    #
    #             for x in y
    #                 total_sum += coefficient * n_coefficient * k_coefficient * I_minus_L(n, x, a) * x * abs(x)^n
    #             end
    #         end
    #     end
    #
    #     if iseven(n)
    #         for k in -1:Int(n / 2 - 1)
    #
    #             k_coefficient = Binomial(n, k) * (1 - (-1)^n)
    #
    #             y = [(Ω + 1 + v * (n - 2 * k)), (Ω + 1 - v * (n - 2 * k)), (Ω - 1 + v * (n - 2 * k)), (Ω - 1 - v * (n - 2 * k))]
    #
    #             for x in y
    #                 total_sum += coefficient * n_coefficient * k_coefficient * I_minus_L(n, x, a) * x * abs(x)^n
    #             end
    #
    #             k_even_coefficient = -Binomial(n, k) * (1 - (-1)^n) / ((n + 1) * beta(n/2 + 1, n/2 + 1))
    #
    #             y_even = [(Ω + 1), (Ω - 1)]
    #
    #             for x in y_even
    #                 total_sum += coefficient * n_coefficient * k_coefficient * I_minus_L(n, x, a) * x * abs(x)^n
    #             end
    #         end
    #     end
    # end

    first_integrand(t) = sinh(Ω * β / 2) * sin(Ω * t) * cos(t) / (a^2 + t^2 - b*cos(v * t))^(3 / 2)

    first_term = coefficient * quadgk(t -> first_integrand(t), 0, Inf)[1]

    second_integrand(t) = (1 - cosh(Ω * (t - β / 2))) * cosh(t) / (a^2 - t^2 - b * cosh(v * t))^(3 / 2)

    second_term = coefficient * quadgk(t -> second_integrand(t), 0, β / 2)[1]

    return first_term + second_term
end

function realX(Ω, β, α, v, w, N=10)

    R = (v^2 - w^2) / (w^2 * v)
    P = 1 / (exp(v * β) - 1)
    Q = 1 / (exp(β) - 1)

    ℜD(x) = w^2 / v^2 * (R * (1 - cos(v * x)) + 4 * R * P * (sin(v * x / 2))^2 + x^2 / β)

    ℑD(x) = -w^2 / v^2 * (R * sin(v * x) - x)

    RD(x) = sqrt(ℜD(x)^2 + ℑD(x)^2)
    θD(x) = atan(ℑD(x), ℜD(x))

    ℑS(x) = 2 * α / (3 * sqrt(π)) * ((1 + Q) * sin(x - 3 * θD(x) / 2) - Q * sin(x + 3 * θD(x) / 2)) / RD(x)^(3 / 2)

    integrand(x) = -(1 - cos(Ω * x)) * ℑS(x) / 150
    @show(Ω)
    return quadgk(x -> integrand(x), 0, Inf)[1]
end

function ℜχ_0(Ω, α, v, w, N = 10)
    R = (v^2 - w^2) / (w^2 * v)
    integrand(x, n) = ((n + 1/2) * x^(n - 1/2) * exp(-R * x) - R * x^(n + 1/2) * exp(-R * x)) * log(abs((1 + n * v + x)^2 / (Ω^2 - (1 + n * v + x)^2))^(1/2))
    total_sum = 0.0
    for n in 0:Int(N)
        total_sum += -(1 / gamma(n + 3/2)) * quadgk(x -> integrand(x, n), 0, Inf)[1]
    end
    return total_sum
end

# ImX(x, y, α, v, w) = ℑχ(x, y, α, v, w, 50)
#
# p = (:α, :v, :w)
# i = (5, 4.0, 2.1)
# pl = ((0.01, 10), (0.01, 10), (0.01, 10))
# al = ((0.01, 100), (0.0001, 100), (0.0, 200))
# a = Dict(
#     "ℜ_title" => "Real Χ(u)",
#     "ℜ_x" => "ℜ(u)",
#     "ℜ_y" => "ℑ(u)",
#     "ℜ_z" => "ℜ[D(u)]",
#     "ℑ_title" => "Imaginary Χ(u)",
#     "ℑ_x" => "ℜ(u)",
#     "ℑ_y" => "ℑ(u)",
#     "ℑ_z" => "ℑ[D(u)]",
# )
#
# vf = viewfunction(ImX, ImX, p, i, pl, a, axis_limits = al, len = 500)

β = 10
v = range(3.0, stop = 21.99, length = 100)
params = (5.0, 4.0, 2.1, 60)
# ImX = [ℑχ(t, β, params...) for t in v]
# ImX_0 = [ℑχ_0(t, params...) for t in v]
# plot(v, ImX_0, label = "β = ∞")
# plot!(v, ImX, label = "β = $(β)")
# xlabel!("Ω")
ReX_0 = [ℜχ_0(t, params...) for t in v]
# ReX = [ℜχ(t, β, params...) for t in v]
plot(v, ReX_0)
# plot!(v, ReX)
xlabel!("Ω")
ylabel!("Reχ(Ω)")


text:
## Esta función da la distancia de seguridad que debe de haber entre dos vehículos
function safe_distance(i, Highway1, N, delta, alfa::Array{Float64, 1} = Float64[0.8f0, 0.75f0])

    pvfo = Vehicle_Ahead(i, Highway1, N)
    d_ifo = (pvfo < N ? (pvfo-i-Highway1[pvfo].len) : delta*5 ) # 5 = vmax[2]
    tipo = ( Highway1[i].tipo == 1 ? 1 : 2)
    dfo = d_ifo + Int8( floor( (1.-alfa[tipo]) * Highway1[pvfo].speed + 0.5) )
    alfa = pvfo = d_ifo = 0
    return dfo
end

## Esta función hace que un vehículo cambie de carril
function Merge!(v0::Vehicle, Count0, C1::Highway1D, vmax::Array{Int64, 1} = [3, 5])

    i = v0.position
    if C1.highway[i].tipo == -1
        if (i < C1.N && C1.highway[i+1].tipo != 1) || i == C1.N

            pvbo = Vehicle_Behind(i, C1.highway, C1.N)
            d_bo = (pvbo < (C1.N) ? (i-pvbo-C1.highway[pvbo].len) : vmax[2] )
            vb = C1.highway[pvbo].speed
            if vb <= d_bo

                Change_Vehicle!( C1.highway[i], v0.speed, v0.position, v0.tipo, 1, v0.num, v0.len)
                Empty_Cell!( v0 )
                Count0 -= 1
                C1.count += 1
            end
        end
    end
end

## Esta función hace los cambios de carril desde un carril izquierdo hacia uno derecho
function Merge_Left_Right!(Ck, Ckminus1, k)

    for v in Ck.highway[end:-1:9]
        if (v.tipo == 1 || v.tipo == 2) && v.change == 0

            delta = (v.tipo == 2 ? Int8(3) : Int8(1))
            l = v.speed
            dfo = safe_distance(v.position, Ckminus1.highway, Ckminus1.N, delta)
            dfs = safe_distance(v.position, Ck.highway, Ck.N, delta)

            if v.tipo == 2 && dfo >= l*delta && (dfs >= delta*l || l > dfs)
                Merge!(v, Ck.count, Ckminus1)
            end

            if v.tipo == 1 #&& dfo >= l*delta && (dfs  >= delta*l || l > dfs)
                pos_behind = Vehicle_Behind(v.position, Ck.highway, Ck.N)
                sfd_b = safe_distance(pos_behind, Ck.highway, Ck.N, delta)

                if (dfs  >= delta*l || l > dfs || Ck.highway[pos_behind].speed > sfd_b) && dfo >= l*delta
                    Merge!(v, Ck.count, Ckminus1)
                end
            end
            l = dfo = dfs = delta = 0
        end
    end
end

## Esta función hace los cambios de carril desde un carril derecho hacia uno izquierdo
function Merge_Right_Left!(Ck, Ckplus1, k, kmax::Int64 = 3)

    for v in Ck.highway[end:-1:9]
        if (v.tipo == 2 || (v.tipo == 1 && k<kmax) ) && v.change == 0 #(v.tipo == 1 && k < 2)
        #if v.tipo == 2 && v.change == 0
            l = v.speed
            dfs = safe_distance(v.position, Ck.highway, Ck.N, Int8(1))

            if dfs < l
                if Ckplus1.highway[v.position].tipo == -1
                    dfo = safe_distance(v.position, Ckplus1.highway, Ckplus1.N, Int8(1))
                    if dfo >= l
                        Merge!(Ck.highway[v.position], Ck.count, Ckplus1)
                    end
                end
            end
            dfs = l = dfo = 0
        end
    end
end

## Esta función simula una rampa de entrada (tipo == 1) o de salida (tipo == 0)
function Rampa!( tipo, x0, lramp, pin, p, Ck, num, vmax::Array{Int8, 1} = Int8[3, 5])
    #vmax = [3, 5]

    # Rampa de salida, tipo == 0
    if tipo == 0
    ############################ cambio ###################################
        x = Pos_der(Ck.carretera[x0:x0+lramp])  # if longitud > 1 Ck.carretera[xi] && ... && Ck.carretera[xi+londitud] esten vacia
        if x <= x0+lramp && rand() < pin
            if Ck.cuenta > 0
                Ck.cuenta -= 1
            end
            Empty_Cell!(Ck.carretera[x])
        end

    # Rampa de entrada, tipo == 1
    elseif tipo == 1
    ############################ cambio ###################################
        x = Pos_siguiente(Ck.carretera[x0:x0+lramp])
        if x <= x0+lramp && rand() < pin
          Ck.cuenta += 1
          v_nueva = (rand() <= p ? vmax[1] : vmax[2])
          tipo_nuevo = (v_nueva == vmax[1] ? 1 : 2)
          Change_Vehicle!( Ck.carretera[x], v_nueva, x, tipo_nuevo, 1, num)
        end
        v_nueva = 0; tipo_nuevo = 0
    end
    x = 0
end

function Ramp!(x0, lramp, p2, p1, Ck, num, vmax::Array{Int8, 1} = Int8[2, 2], Len::Array{Int64, 1} = Int64[2, 1])

  ############# Ramp light vehicles #################

  ############# Entrance           #################
  if p2 > 0
    if rand() < p2
      intro = false
      x_end = x0+lramp
      x = Pos_next(Ck.highway[x0:x_end])
      while !intro && x > x0
        if Ck.highway[x+1].tipo != 1
          Ck.count += 1
          Change_Vehicle!( Ck.highway[x], vmax[2], x, 2, 1, num, Len[2])
          intro = true
        else
          x_end = x-1
          x = Pos_next(Ck.highway[x0:x_end])
        end
      end
    end
  end
  ############# Exit            #################
  if p2 < 0
     if rand() < -p2
       out = false
       x_end = x0+lramp
       x = Pos_right(Ck.highway[x0:x_end], Int8(2))
       while !out && x > x0
         if Ck.highway[x].tipo == 2
           Ck.count -= 1
           Empty_Cell!(Ck.highway[x])
           out = true
         else
           x_end = x-1
           x = Pos_right(Ck.highway[x0:x_end], Int8(2))
         end
       end
     end
   end

  ############# Ramp trucks    #################

  ############# Entrance         #################
  if p1 > 0
    if rand() < p1
      intro = false
      x_end = x0+lramp
      x = Pos_next(Ck.highway[x0:x_end])
      while !intro && x > x0
        if Ck.highway[x+1].tipo != 1 && Ck.highway[x-1].tipo == -1
          Ck.count += 1
          Change_Vehicle!(Ck.highway[x], vmax[1], x, 1, 1, num, Len[1])
          intro = true
        else
          x_end = x-1
          x = Pos_next(Ck.highway[x0:x_end])
        end
      end
    end
  end
    ############# Exit            #################
  if p1 < 0
    if rand() < -p1
      out = false
      x_end = x0+lramp
      x = Pos_right(Ck.highway[x0:x_end], Int8(1))
      while !out && x > x0
        if Ck.highway[x].tipo == 1
          Ck.count -= 1
          Empty_Cell!(Ck.highway[x])
          out = true
        else
          x_end = x-1
          x = Pos_right(Ck.highway[x0:x_end], Int8(2))
        end
      end
    end
  end
end


text:
@test Either{String, Int}(4) == Identity(4)
@test Either{String, Int}("left") == Const("left")
@test Either{Symbol}(true) == Identity(true)
@test Either{Symbol}(:fallback) == Const(:fallback)

@test @either(true ? 4 : "else") == Identity(4)
@test @either(false ? 4 : "else") == Const("else")
@test @either(if true
    4
else
    "else"
end) == Identity(4)
@test @either(if false
    4
else
    "else"
end) == Const("else")

@test iseither(Identity(4))
@test iseither(Const("hi"))
@test !iseither(Base.Some(nothing))

@test isright(either("hi", true, 3))
@test isright(Either{String}(3))
@test !isleft(Either{String}(3))

@test getright(Either{String}(3)) == 3
@test_throws MethodError getright(Either{String}("hi"))

@test getleft(Either{String}("hi")) == "hi"
@test_throws MethodError getleft(Either{String}(3))

@test getrightOption(Either{String}(3)) == Option(3)
@test getrightOption(Either{String}("fallback")) == Option()

@test getleftOption(Either{String}(3)) == Option()
@test getleftOption(Either{String}("fallback")) == Option("fallback")

@test getOption(Either{String}(3)) == Option(3)
@test getOption(Either{String}("fallback")) == Option()

@test eltype(Either{Int, String}) == String
@test eltype(Either{<:Any, String}) == String
@test eltype(Either) == Any

@test flip_left_right(Const(4)) == Identity(4)
@test flip_left_right(Identity("hi")) == Const("hi")
h = Identity(:hellp)
@test flip_left_right(flip_left_right(h)) == h