text:
export k_guessing_game, ambiguous_guessing_game

export maximum_likelihood_facet, anti_guessing_facet, ambiguous_guessing_facet, k_guessing_facet, non_negativity_facet
export coarse_grained_input_ambiguous_guessing_facet

"""
The ``k``-guessing game is a general family of Bell inequalities for signaling polytopes.
For any integer ``k\\in[0,Y]`` a ``k``-guessing game Bell inequality
``(\\mathbf{G}_{\\text{K}},\\gamma_{\\text{K}})``
bounds the signaling polytope ``\\mathcal{C}_d^{\\binom{Y}{k}\\to Y }``.
The columns of matrix ``\\mathbf{G}_{\\text{K}}\\in \\mathbf{R}^{Y \\times \\binom{Y}{k}}``
consist of all ``\\binom{Y}{k}`` ways to arrange ``k`` unit elements and ``(Y-k)``
null elements.
Hence each input ``x\\in[X]`` corresponds to a unique set of ``k`` correct
answers out of ``Y`` possible answers.
The upper bound is  ``\\gamma_{\\text{K}} = \\binom{Y}{k} - \\binom{Y-d}{k}``.

A ``k``-guessing game Bell inequality is constructed with the `k_guessing_game`
method,

```julia
k_guessing_game(Y :: Int64, d :: Int64, k :: Int64) :: BellGame
```

Parameters:
* `Y` - The number of outputs
* `d` - The signaling dimension
* `k` - The number of unit elements in each column. Must satisfy `Y ≥ k ≥ 0`.

For example,

```jldoctest k_guessing_game
using SignalingDimension

G_K = k_guessing_game(6,2,3)

# output

6×20 BellScenario.BellGame:
 1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  1  1  1  1  1  1  0  0  0  0
 1  0  0  0  1  1  1  0  0  0  1  1  1  0  0  0  1  1  1  0
 0  1  0  0  1  0  0  1  1  0  1  0  0  1  1  0  1  1  0  1
 0  0  1  0  0  1  0  1  0  1  0  1  0  1  0  1  1  0  1  1
 0  0  0  1  0  0  1  0  1  1  0  0  1  0  1  1  0  1  1  1
```

The upper bound is then

```jldoctest k_guessing_game
G_K.β

# output

16
```

The `k_guessing_game` can alternatively take a `BellScenario.LocalSignaling(X, Y, d)` scenario as input.

```julia
k_guessing_game( scenario :: LocalSignaling, k :: Int64 ) :: BellGame
```

A `DomainError` is thrown if `scenario.X !=  binomial(scenario.Y, k)`.
"""
function k_guessing_game(scenario :: LocalSignaling, k :: Int64) :: BellGame
    if scenario.X != binomial(scenario.Y, k)
        throw(DomainError(scenario.X, "Number of inputs `X` must be `binomial(Y, k)`."))
    end

    k_guessing_game(scenario.Y, scenario.d, k)
end
function k_guessing_game(Y :: Int64, d :: Int64, k :: Int64) :: BellGame
    if !(Y ≥ k ≥ 0)
        throw(DomainError(k, "Input k must satisfy `Y  ≥ k ≥ 0`"))
    end

    game = hcat(map( combo -> begin
        a = zeros(Int64, Y)
        a[combo] .= 1
        return a
    end, combinations(1:Y, k))...)

    max_score = binomial(Y,k) - binomial(Y-d,k) #sum(map(i -> binomial(Y-i, k-1), 1:d))

    BellGame(game, max_score)
end

"""
Ambiguous guessing games are a family of Bell inequalities for signaling polytopes.
An ambiguous guessing game, denoted by ``(\\mathbf{G}_{\\text{Q}},\\gamma_{\\text{Q}})``, bounds the
signaling polytope ``\\mathcal{C}_d^{X \\to Y}``.
The ``Y \\times X`` matrix ``\\mathbf{G}_{\\text{Q}}`` is described as having two types of rows:
* *Guessing Rows*: one column contains a non-zero element of value ``(X - d + 1)``;
* *Ambiguous Rows*: each column contains a ``1``.

For any integer ``k \\in [0,Y]`` An ambiguous guessing game is defined to have ``k`` guessing rows
and ``(Y-k)`` ammbiguous rows.
The upper bound on the  ambiguous guessing game Bell inequality is then ``\\gamma_{\\text{Q}} = d(X-d+1)``.

The `ambiguous_guessing_game` method constructs an ambiguous guessing game Bell inequality with `k`
guessing rows.

```julia
ambiguous_guessing_game(X::Int64, Y::Int64, d::Int64, k::Int64) :: BellGame
```

Parameters:
* `X` - The number of inputs
* `Y` - The number of outputs
* `d` - The signaling dimension
* `k` - Th number of guessing rows. Must satisfy `Y ≥ k ≥ 0`.

For example, the following constructed ambiguous guessing game has four guessing
rows and  3 ambiguous rows.

```jldoctest ambiguous_guessing_game
using SignalingDimension

(X, Y, d) = (6, 7, 3)

k = 4        # num guessing rows

G_Q = ambiguous_guessing_game(X, Y, d, k)

# output

7×6 BellScenario.BellGame:
 4  0  0  0  0  0
 0  4  0  0  0  0
 0  0  4  0  0  0
 0  0  0  4  0  0
 1  1  1  1  1  1
 1  1  1  1  1  1
 1  1  1  1  1  1
```

```jldoctest ambiguous_guessing_game
G_Q.β

# output

12
```

Alternatively, the `ambiguous_guessing_game` can accept a `BellScenario.LocalSignaling(X, Y, d)` scenario.

```julia
ambiguous_guessing_game(scenario :: LocalSignaling, k :: Int64) :: BellGame
```
"""
function ambiguous_guessing_game(scenario :: LocalSignaling, k :: Int64) :: BellGame
    ambiguous_guessing_game(scenario.X,scenario.Y,scenario.d,k)
end
function ambiguous_guessing_game(X::Int64, Y::Int64, d::Int64, k::Int64) :: BellGame
    if !(Y ≥ k ≥ 0)
        throw(DomainError(k, "num guessing rows `k` should satisfy `Y ≥ k ≥ 0`."))
    elseif !(min(X,Y) ≥ d ≥ 1)
        throw(DomainError(d, "`scenario.d` should satisfy `min(X, Y) ≥ d ≥ 1`."))
    end

    ambiguous_scalar = (X - d + 1)

    γ = d * ambiguous_scalar
    G = zeros(Int64, Y, X)

    for y in 1:k
        if y ≤ X
            G[y,y] = ambiguous_scalar
        else
            G[y,X] = ambiguous_scalar
        end
    end

    ambiguous_row = ones(Int64, X)
    for y in (k+1):Y
        G[y,:] =  ambiguous_row
    end

    BellGame(G, γ)
end

"""
The maximum likelihood facet ``(\\mathbf{G}_{\\text{ML}},d)`` tightly bounds the
``\\mathcal{C}_d^{N \\to N}``  signaling polytope for ``N > d > 1``.
The matrix ``\\mathbf{G}_{\\text{ML}}`` is a ``(k=1)``-guessing game (see [`k_guessing_game`](@ref)),
hence, ``\\mathbf{G}_{\\text{ML}}=\\mathbb{I}_{N}`` the ``N\\times N`` identity matrix.
The upper bound is  ``d \\geq \\langle \\mathbf{G}_{\\text{ML}}, \\mathbf{P}\\rangle``
for any channel ``\\mathbf{P}\\in\\mathcal{C}_d^{N \\to N}``.

```julia
maximum_likelihood_facet( N :: Int64, d :: Int64 ) :: BellGame
```

Construct the maximum likelihood facet for the ``\\mathcal{C}_d^{N \\to N}``
signaling polytope. Note that `N` specifies the number of inputs and outputs.
For example,

```jldoctest maximum_likelihood_facet
using SignalingDimension

N = 4    # number of inputs and outputs
d = 2    # bit signaling

G_ML = maximum_likelihood_facet(N, d)

# output

4×4 BellScenario.BellGame:
 1  0  0  0
 0  1  0  0
 0  0  1  0
 0  0  0  1
```

```jldoctest maximum_likelihood_facet
G_ML.β   # the upper bound

# output

2
```

A `DomainError` is thrown if the following requirements aren't satisfied:
* `N > 2
* `N > d > 1`
"""
function maximum_likelihood_facet(N :: Int64, d :: Int64) :: BellGame
    if N ≤ 2
        throw(DomainError(N, "Input N must satisfy `N > 2`"))
    elseif !(N > d > 1)
        throw(DomainError(d, "Input d must satisfy `N > d > 1`"))
    end

    BellGame(Matrix{Int64}(I,N,N), d)
end

"""
    anti_guessing_facet( N :: Int64, d :: Int64, ε :: Int64 ) :: BellGame

Constructs the anti-guessing facet of the ``\\mathcal{C}_d^{N \\to N}`` signaling polytope.
Input `ε` sets the size of the ``(k=N-1)``-guessing game block while the remaining
unit elements are arranged along the diagonal.
The upper bound of the anti-guessing game Bell inequality is ``\\gamma = \\varepsilon + d - 2``.
For example,

```jldoctest anti_guessing_game
using SignalingDimension

N = 6    # number of inputs and outputs
d = 2    # dit signaling
ε = 4    # anti-guessing block size

G_A = anti_guessing_facet(N,d,ε)

# output

6×6 BellScenario.BellGame:
 0  1  1  1  0  0
 1  0  1  1  0  0
 1  1  0  1  0  0
 1  1  1  0  0  0
 0  0  0  0  1  0
 0  0  0  0  0  1
```
```jldoctest anti_guessing_game
G_A.β    # upper bound

# output

4
```

Note in the above example that

A `DomainError` is thrown if the following requirements aren't satisfied:
* `N ≥ 4`
* `(N - 2) ≥ d ≥ 2`
* `(N - d + 1) ≥ ε ≥ 3`
"""
function anti_guessing_facet(N :: Int64, d :: Int64, ε :: Int64) :: BellGame
    if !(N ≥ 4)
        throw(DomainError(N, "Input N must satisfy `N ≥ 4`"))
    elseif !(N-2 ≥ d ≥ 2)
        throw(DomainError(d, "Input d must satisfy `(N - 2) ≥ d ≥ 2`"))
    elseif !(N-d+1 ≥ ε ≥ 3)
        throw(DomainError(ε, "Input ε must satisfy `(N - d + 1) ≥ ε ≥ 3`"))
    end

    err_game = ones(Int64, ε, ε) - I
    succ_game = Matrix{Int64}(I, (N-ε),(N-ε))

    BellGame(cat(
        cat(err_game, zeros(Int64, ε, (N-ε)), dims=2),
        cat(zeros(Int64, (N-ε), ε), succ_game, dims=2),
        dims = 1
    ), ε + d - 2)
end

"""
    ambiguous_guessing_facet( Y :: Int64, d :: Int64 ) :: BellGame

Constructs the ambiguous guessing facet for the ``\\mathcal{C}_d^{(Y-1)\\to Y} polytope.
This is a special case of the [`ambiguous_guessing_game`](@ref) where `Y = X - 1`.
For example,

```jldoctest ambiguous_guessing_facet
using SignalingDimension

Y = 5
d = 3

G_Q = ambiguous_guessing_facet(Y, d)

# output

5×4 BellScenario.BellGame:
 2  0  0  0
 0  2  0  0
 0  0  2  0
 0  0  0  2
 1  1  1  1
```
```jldoctest ambiguous_guessing_facet
G_Q.β

# output

6
```

A `DomainError` is thrown if the inputs don't satisfy the following requirements:
* `Y ≥ 4`
* `(Y - 2) ≥ d ≥ 2`
"""
function ambiguous_guessing_facet(Y :: Int64, d :: Int64) :: BellGame
    if !(Y ≥ 4)
        throw(DomainError(Y, "Input Y must satisfy `Y ≥ 4`"))
    elseif !(Y-2 ≥ d ≥ 2)
        throw(DomainError(d, "Input d must satisfy `(Y - 2) ≥ d ≥ 2`"))
    end

    succ_game = (Y-d)*Matrix{Int64}(I, Y-1, Y-1)

    BellGame( cat(succ_game, ones(Int64, 1, Y-1), dims = 1), d*(Y-d))
end

"""
    k_guessing_facet( Y :: Int64, d :: Int64, k :: Int64 ) :: BellGame

Constructs the ``k``-guessing facet for the ``\\mathcal{C}_d^{\\binom{Y}{k}\\to Y}``
signaling polytope. A ``k``-guessing game is tight when ``Y = d + k``.
Note that `k` is the number of unit elements in each column.

A `DomainError` is satisfied if the following requirements aren't satisfied:
* `Y == k + d`
* `Y ≥ 3`
"""
function k_guessing_facet(Y :: Int64, d :: Int64, k :: Int64) :: BellGame
    if Y != k + d
        throw(DomainError(d, "Input `d` must satisfy `Y == k + d`."))
    elseif !(Y ≥ 3)
        throw(DomainError(Y, "Input N must satisfy `Y ≥ 3`"))
    end

    k_guessing_game(Y, d, k)
end

"""
A non-negativity facet reflects the fact that ``P(y|x) \\geq 0``. These facets
bound all signaling polytopes.

    non_negativity_facet( X :: Int64, Y :: Int64 ) :: BellGame

Constructs the non-negativity game for a channel with `X` inputs and `Y` outputs.
For example

```jldoctest non_negativity_facet
using SignalingDimension

G = non_negativity_facet(3, 4)

# output

4×3 BellScenario.BellGame:
 1  0  0
 1  0  0
 1  0  0
 0  0  0
```
```jldoctest non_negativity_facet
G.β

# output

1
```

A `DomainError` is thrown if `X` or `Y` is not greater than 1.
"""
function non_negativity_facet(X :: Int64, Y :: Int64) :: BellGame
    if !(Y > 1)
        throw(DomainError(Y, "Inputs must satisfy `Y > 1`"))
    elseif !(X > 1)
        throw(DomainError(X, "Inputs must satisfy `X > 1`"))
    end

    m = zeros(Int64, Y, X)
    m[1:end-1,1] .= 1

    BellGame(m, 1)
end

"""
    coarse_grained_input_ambiguous_guessing_facet(
        Y :: Int64,
        d :: Int64
    ) :: BellGame

Constructs a canonical form for a input coarse-grained ambiguous game.
For example

```jldoctest
using SignalingDimension

Y = 4
d = 2

G = coarse_grained_input_ambiguous_guessing_facet(Y, d)

# output

4×4 BellScenario.BellGame:
 2  0  0  0
 0  2  0  0
 0  0  1  1
 1  1  1  0
```

A `DomainError` is thrown if the inputs don't satisfy the following requirements:
* `Y ≥ 4`
* `(Y - 2) ≥ d ≥ 2`
"""
function coarse_grained_input_ambiguous_guessing_facet(Y :: Int64, d :: Int64) :: BellGame
    G_Q = ambiguous_guessing_facet(Y, d)

    game = zeros(Int64,Y,Y)
    game[1:Y,1:end-1] = G_Q
    game[Y-1,Y] = game[Y-1,Y-1] - 1
    game[Y-1,Y-1] = 1

    BellGame(game, G_Q.β)
end


text:
using Test
using RandomV06

v6 = RandomV06

v6.srand(1)
@test v6.rand(1) == [0.23603334566204692]
r = v6.rand(5)
@test r[5] == 0.21096820215853596
@test v6.randperm(6) == [4; 1; 5; 3; 6; 2]

v6.srand(1)
@test v6.randcycle(5) == [4;1;5;3;2]

v6.srand(1)
@test v6.rand(1:100,5) == [99;47;60;70;74]

v6.randn(3,3)
v6.randn(3,3,3)

x = ones(5)
v6.rand!(x)
v6.randn!(x)
v6.randexp!(x)

v6.bitrand(3,4)
v6.randstring(4)
v6.randsubseq(collect(1:6),0.3)
x = []; v6.randsubseq!(x,collect(1:12),0.5)

x = collect(1:10)
v6.shuffle(x)
v6.shuffle!(x)

rng = v6.MersenneTwister(0);

V06 = RandomV06.V06
V07 = RandomV06.V07

seed_ver!(V06, 1)
rand_ver(V06, 3)
randn_ver(V06, 3)
rand_ver(V06, ["a"; "b"], 3, 4)
x = ones(4);
rand_ver!(V06, x)
randn_ver!(V06, x)
randexp_ver(V06, 3)
randexp_ver!(V06, x)
randstring_ver(V06, 3)
randperm_ver(V06, 7)
randcycle_ver(V06, 7)
shuffle_ver(V06, x)
shuffle_ver!(V06, x)

seed_ver!(V07, 1)
rand_ver(V07, 3)
randn_ver(V07, 3)
rand_ver(V07, ["a"; "b"], 3, 4)
x = ones(4);
rand_ver!(V07, x)
randn_ver!(V07, x)
randexp_ver(V07, 3)
randexp_ver!(V07, x)
randstring_ver(V07, 3)
randperm_ver(V07, 7)
randcycle_ver(V07, 7)
shuffle_ver(V07, x)
shuffle_ver!(V07, x)

ext: jl


text:
## Array views of other indexable containers
"""
    ViewArray(parent, indices)

Create an array `out` which is a lazy view of `parent`. The indices of `out` match the
indices of `indices`, and the values of `indices` are used to index `parent`. Unlike
`SubArray`, the parent type could be any indexable container (for instance, we can create an
`AbstractArray` view of the values stored in an `AbstractDict`).

See also `ViewDict`, `view`.

# Examples

```julia
julia> d = Dict(:a => 11, :b => 12, :c => 13)
Dict{Symbol,Int64} with 3 entries:
  :a => 11
  :b => 12
  :c => 13

julia> ViewArray(d, [:a, :c])
2-element ViewArray{Int64,1,Dict{Symbol,Int64},Array{Symbol,1}}:
 11
 13
```
"""
struct ViewArray{T, N, A, I <: AbstractArray{<:Any, N}} <: AbstractArray{T, N}
    parent::A
    indices::I
end

ViewArray(a, indices::AbstractArray) = ViewArray{_valtype(a)}(a, indices)
ViewArray{T}(a, indices::AbstractArray) where {T} = ViewArray{T, ndims(indices)}(a, indices)
ViewArray{T, N}(a, indices::AbstractArray) where {T, N} = ViewArray{T, N, typeof(a)}(a, indices)
ViewArray{T, N, A}(a, indices::AbstractArray) where {T, N, A} = ViewArray{T, N, A, typeof(indices)}(a, indices)

# vector & matrix versions
const ViewVector{T, A, I <: AbstractVector} = ViewArray{T, 1, A, I}
ViewVector(a, indices::AbstractVector) = ViewVector{_valtype(a)}(a, indices)

const ViewMatrix{T, A, I <: AbstractMatrix} = ViewArray{T, 2, A, I}
ViewMatrix(a, indices::AbstractMatrix) = ViewMatrix{_valtype(a)}(a, indices)

# `AbstractArray` implementation (no resizing)
IndexStyle(::Type{ViewArray{T, N, A, I}}) where {T, N, A, I} = IndexStyle(I)
axes(a::ViewArray) = axes(a.indices)
size(a::ViewArray) = size(a.indices)

@propagate_inbounds getindex(a::ViewArray, i::Int) = a.parent[a.indices[i]]
@propagate_inbounds getindex(a::ViewArray, is::Int...) = a.parent[a.indices[is...]]

@propagate_inbounds setindex!(a::ViewArray, v, i::Int) = (a.parent[a.indices[i]] = v)
@propagate_inbounds setindex!(a::ViewArray, v, is::Int...) = (a.parent[a.indices[is...]] = v)

# Connect with `view`
"""
    view(d::AbstractDict, inds::AbstractArray)

Create an array view over the dictionary `d` with the indices stored in `inds`.

See also `ViewArray`.

# Examples

```julia
julia> d = Dict(:a => 11, :b => 12, :c => 13)
Dict{Symbol,Int64} with 3 entries:
  :a => 11
  :b => 12
  :c => 13

julia> view(d, [:a, :c])
2-element ViewArray{Int64,1,Dict{Symbol,Int64},Array{Symbol,1}}:
 11
 13
```
"""
view(d::AbstractDict, inds::AbstractArray) = ViewArray{_valtype(d), ndims(inds)}(d, inds)
parent(a::ViewArray) = a.parent

## Dictionary views of indexable containers
"""
    ViewDict(parent, indices)

Create an array `out` which is a lazy view of `parent`. The indices of `out` match the
indices of `indices`, and the values of `indices` are used to index `parent`. The parent
must be an indexable type - for instance, an array or dictionary.

See also `ViewArray`, `view`.

# Examples

```julia
julia> d = Dict(:a => "Alice", :b => "Bob", :c => "Charlie")
Dict{Symbol,String} with 3 entries:
  :a => "Alice"
  :b => "Bob"
  :c => "Charlie"

julia> ViewDict(d, Dict(:aa => :a, :cc => :c))
ViewDict{Symbol,String,Dict{Symbol,String},Dict{Symbol,Symbol}} with 2 entries:
  :aa => "Alice"
  :cc => "Charlie"

julia> v = [11, 12, 13]
3-element Array{Int64,1}:
 11
 12
 13

julia> ViewDict(v, Dict(:a =>1 , :c => 3))
ViewDict{Symbol,Int64,Array{Int64,1},Dict{Symbol,Int64}} with 2 entries:
  :a => 11
  :c => 13
```
"""
struct ViewDict{K, V, A, I} <: AbstractDict{K, V}
    parent::A
    indices::I
end

ViewDict(a, indices) = ViewDict{keytype(indices)}(a, indices)
ViewDict{K}(a, indices) where {K} = ViewDict{K, _valtype(a)}(a, indices)
ViewDict{K, V}(a, indices) where {K, V} = ViewDict{K, V, typeof(a)}(a, indices)
ViewDict{K, V, A}(a, indices) where {K, V, A} = ViewDict{K, V, A, typeof(indices)}(a, indices)

# `AbstractDict` implementation (immutable keys)
keys(d::ViewDict) = keys(d.indices)
haskey(d::ViewDict, i) = haskey(d.indices, i)
length(d::ViewDict) = length(d.indices)

@propagate_inbounds getindex(d::ViewDict, i) = d.parent[d.indices[i]]
@propagate_inbounds setindex!(d::ViewDict, v, i) = (d.parent[d.indices[i]] = v)

function iterate(d::ViewDict, state...)
    tmp = iterate(d.indices, state...)
    if tmp === nothing
        return tmp
    end
    (kv, state2) = tmp
    return (kv.first => d.parent[kv.second], state2)
end

# Connect with `view`
"""
    view(parent, inds::AbstractDict)

Create a dictionary view over the indexable container `parent`. The indices of the output
are the indices of `inds`, and the values are the corresponding `parent[i]` for
`i in values(inds)`.

See also `ViewDict`.

# Examples

```julia
julia> d = Dict(:a => "Alice", :b => "Bob", :c => "Charlie")
Dict{Symbol,String} with 3 entries:
  :a => "Alice"
  :b => "Bob"
  :c => "Charlie"

julia> view(d, Dict(:aa => :a, :cc => :c))
ViewDict{Symbol,String,Dict{Symbol,String},Dict{Symbol,Symbol}} with 2 entries:
  :aa => "Alice"
  :cc => "Charlie"

julia> v = [11, 12, 13]
3-element Array{Int64,1}:
 11
 12
 13

julia> view(v, Dict(:a =>1 , :c => 3))
ViewDict{Symbol,Int64,Array{Int64,1},Dict{Symbol,Int64}} with 2 entries:
  :a => 11
  :c => 13
```
"""
view(d, inds::AbstractDict) = ViewDict(d, inds)
view(d::AbstractArray, inds::AbstractDict) = ViewDict(d, inds) # disambiguation
parent(d::ViewDict) = d.parent


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

using Markdown
using InteractiveUtils

# ╔═╡ dd4498f1-53c4-48e7-b651-70d48ae91379
using Pkg, DrWatson

# ╔═╡ ac129262-c0a8-47d2-af1a-1e1268096451
begin
	@quickactivate "NumericalMethodsforEngineers"
	using Symbolics, SymbolicUtils, Latexify
	using NumericalMethodsforEngineers, Test
end

# ╔═╡ 33c7c64a-7a2c-4beb-9150-74f6c8444f60
begin
	@variables t α σ(..) β[1:2]
end;

# ╔═╡ 39a8c494-5608-470a-89f6-9f4a68a9a527
@variables x y;

# ╔═╡ d38b1175-04c2-4f15-a100-ec7cda34c1de
begin
	@variables b[1:3]
	@variables w(..) x1(t) y z(t, α, x)
	b = [2, 5, 9]
end;

# ╔═╡ b70e72aa-32b9-4ae7-8b5f-83a770e52ac7
begin
	z1 = x^2 + y
	L(z1)
end

# ╔═╡ 2d01cb71-28b8-4c40-804b-68faf4a19890
begin
	expr = β[1]* x + y^α + σ(3) * (z1 - t) - β[2] * w(t - 1)
	L(expr)
end

# ╔═╡ fb931ec7-1704-4490-a82a-555a644e1bbc
@variables A[1:3, 1:3] B[1:3, 1:3];

# ╔═╡ f0c84794-6ecc-49ff-a769-e163ef8c44dc
x ~ y

# ╔═╡ 90cd94e1-602b-49af-9913-bfbcd134f520
x - y ~ 0

# ╔═╡ 28c312b7-b4a9-4d0d-8988-8ef8a16cb786
A .~ B

# ╔═╡ 076d6d44-a552-458c-b4ce-1351090f5353
A .~ 3x

# ╔═╡ 83feb017-91d9-4363-811b-9c2be8f650a2
A .~ [1 2 3;]' * x

# ╔═╡ 2202ac58-f7c6-49bf-beee-c32737a83b89
A .~ [1, 2, 3] * x

# ╔═╡ 7659d23b-58fc-4528-a42f-7a818b364958
L(A[3,2])

# ╔═╡ Cell order:
# ╠═dd4498f1-53c4-48e7-b651-70d48ae91379
# ╠═ac129262-c0a8-47d2-af1a-1e1268096451
# ╠═d38b1175-04c2-4f15-a100-ec7cda34c1de
# ╠═b70e72aa-32b9-4ae7-8b5f-83a770e52ac7
# ╠═2d01cb71-28b8-4c40-804b-68faf4a19890
# ╠═33c7c64a-7a2c-4beb-9150-74f6c8444f60
# ╠═39a8c494-5608-470a-89f6-9f4a68a9a527
# ╠═fb931ec7-1704-4490-a82a-555a644e1bbc
# ╠═f0c84794-6ecc-49ff-a769-e163ef8c44dc
# ╠═90cd94e1-602b-49af-9913-bfbcd134f520
# ╠═28c312b7-b4a9-4d0d-8988-8ef8a16cb786
# ╠═076d6d44-a552-458c-b4ce-1351090f5353
# ╠═83feb017-91d9-4363-811b-9c2be8f650a2
# ╠═2202ac58-f7c6-49bf-beee-c32737a83b89
# ╠═7659d23b-58fc-4528-a42f-7a818b364958


text: isdebug(group) = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, group, oneL0) !== nothing


text:
println("--------- statring LinearExpressions.jl tests ---------")
include("test/runtests.jl")
println("✓✓✓✓✓✓✓✓ all LinearExpressions.jl tests passed ✓✓✓✓✓✓✓✓")


text:
##### Array{T} #####
using Mmap

struct MMSer{T} end # sentinal type used in deserialize to switch to our
                    # custom memory-mapping deserializers

can_mmap(io::IOStream) = true
can_mmap(io::IO) = false

function padalign(io::IOStream, align=8)
    p = position(io)
    if p & (align-1) !== 0
        write(io, zeros(UInt8, (align-p % align)))
    end
end
padalign(io, align=8) = nothing

function seekpadded(io::IOStream, align=8)
    p = position(io)
    if p & (align-1) !== 0
        seek(io, p+(align-p%align))
    end
end
seekpadded(io, align=8) = nothing

function deserialize(s::AbstractSerializer, ::Type{MMSer{T}}) where T
    mmread(T, s, can_mmap(s.io))
end

mmwrite(io::IO, xs) = mmwrite(Serializer(io), xs)
mmwrite(io::AbstractSerializer, xs) = serialize(io, xs) # fallback

function mmwrite(io::AbstractSerializer, arr::A) where A<:Union{Array,BitArray}
    T = eltype(A)
    Serialization.serialize_type(io, MMSer{typeof(arr)})
    if T<:Union{} || T==Missing
        serialize(io, size(arr))
        return
    elseif isbitstype(T)
        serialize(io, size(arr))
        padalign(io.io, sizeof(eltype(arr)))
        write(io.io, arr)
        return
    end

    fl = fixedlength(T)
    if fl > 0
        serialize(io, size(arr))
        for x in arr
            fast_write(io.io, x)
        end
    else
        serialize(io, arr)
    end
end

function mmread(::Type{A}, io, mmap) where A<:Union{Array,BitArray}
    T = eltype(A)
    if T<:Union{} || T==Missing
        sz = deserialize(io)
        return Array{T}(undef, sz)
    elseif isbitstype(T)
        sz = deserialize(io)
        seekpadded(io.io, sizeof(T))

        if prod(sz) == 0
            return A(undef, sz...)
        end
        if mmap
            data = Mmap.mmap(io.io, A, sz, position(io.io))
            seek(io.io, position(io.io)+sizeof(data)) # move
            return data
        else
            data = Base.read!(io.io, A(undef, sz...))
            return data
        end
    end

    fl = fixedlength(T)
    if fl > 0
        sz = deserialize(io)
        arr = A(undef, sz...)
        @inbounds for i in eachindex(arr)
            arr[i] = fast_read(io.io, T)::T
        end
        return arr
    else
        return deserialize(io) # slow!!
    end
end

##### Array{String} #####

const UNDEF_LENGTH = typemax(UInt32) # if your string is exaclty 4GB you're out of luck

function mmwrite(io::AbstractSerializer, xs::Array{String})
    Serialization.serialize_type(io, MMSer{typeof(xs)})

    lengths = UInt32[]
    buffer = UInt8[]
    serialize(io, size(xs))
    # todo: write directly to buffer, but also mmap
    ptr = pointer(buffer)
    for i in 1:length(xs)
        if isassigned(xs, i)
            x = xs[i]
            l = sizeof(x)
            lb = length(buffer)
            push!(lengths, l)
            resize!(buffer, lb+l)
            unsafe_copyto!(pointer(buffer)+lb, pointer(x), l)
            ptr += l
        else
            push!(lengths, UNDEF_LENGTH)
        end
    end

    mmwrite(io, buffer)
    mmwrite(io, lengths)
end

function mmread(::Type{Array{String,N}}, io, mmap) where N
    sz = deserialize(io)
    buf = deserialize(io)
    lengths = deserialize(io)

   #@assert length(buf) == sum(filter(x->x>0, lengths))
   #@assert prod(sz) == length(lengths)

    ys = Array{String,N}(undef, (sz...,)) # output
    ptr = pointer(buf)
    @inbounds for i = 1:length(ys)
        l = lengths[i]
        l == UNDEF_LENGTH && continue
        ys[i] = unsafe_string(ptr, l)
        ptr += l
    end
    ys
end

## Optimized fixed length IO
## E.g. this is very good for `StaticArrays.MVector`s

function fixedlength(t::Type, cycles=IdDict())
    if isbitstype(t)
        return sizeof(t)
    elseif isa(t, UnionAll) || isabstracttype(t) || Base.isbitsunion(t)
        return -1
    end

    if haskey(cycles, t)
        return -1
    end
    cycles[t] = nothing
    lens = ntuple(i->fixedlength(fieldtype(t, i), copy(cycles)), fieldcount(t))
    if isempty(lens)
        # e.g. abstract type / array type
        return -1
    elseif any(x -> x < 0, lens)
        return -1
    else
        return sum(lens)
    end
end

fixedlength(::Type{>:Missing}, cycles=nothing) = -1
fixedlength(::Type{<:String}, cycles=nothing) = -1
fixedlength(::Type{Union{}}, cycles=nothing) = -1
fixedlength(::Type{<:Ptr}, cycles=nothing) = -1

function gen_writer(::Type{T}, expr) where T
    @assert fixedlength(T) >= 0 "gen_writer must be called for fixed length eltypes"
    if T<:Tuple && isbitstype(T)
        :(write(io, Ref{$T}($expr)))
    elseif length(T.types) > 0
        :(begin
              $([gen_writer(fieldtype(T, i), :(getfield($expr, $i))) for i=1:fieldcount(T)]...)
          end)
    elseif isbitstype(T) && sizeof(T) == 0
        return :(begin end)
    elseif isbitstype(T)
        return :(write(io, $expr))
    else
        error("Don't know how to serialize $T")
    end
end

function gen_reader(::Type{T}) where T
    @assert fixedlength(T) >= 0 "gen_reader must be called for fixed length eltypes"
    ex = if T<:Tuple
        if isbitstype(T)
            return :(read!(io, Ref{$T}())[])
        else
            exprs = [gen_reader(fieldtype(T, i)) for i=1:fieldcount(T)]
            return :(tuple($(exprs...)))
        end
    elseif length(T.types) > 0
        return :(ccall(:jl_new_struct, Any, (Any,Any...), $T, $([gen_reader(fieldtype(T, i)) for i=1:fieldcount(T)]...)))
    elseif isbitstype(T) && sizeof(T) == 0
        return :(ccall(:jl_new_struct, Any, (Any,Any...), $T))
    elseif isbitstype(T)
        return :($T; read(io, $T))
    else
        error("Don't know how to deserialize $T")
    end
    return :($T; $ex)
end

@generated function fast_write(io, x)
    gen_writer(x, :x)
end

@generated function fast_read(io, ::Type{T}) where T
    gen_reader(T)
end


text:
int_rules_1_2_2_5 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.2.5*P(x)*(a+b*x^2+c*x^4)^p =#
    @apply_utils Antiderivative(~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand(~Pq * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (PolyQ(~Pq, ~x) && IGtQ(~p, 0))
    @apply_utils Antiderivative(~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => Antiderivative(~x * PolynomialQuotient(~Pq, ~x, ~x) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (PolyQ(~Pq, ~x) && (EqQ(Coeff(~Pq, ~x, 0), 0) && Not(MatchQ(~Pq, (~x) ^ ~(m') * ~(u') <-- IntegerQ(m)))))
    @apply_utils Antiderivative(~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => Module([q = Expon(~Pq, ~x), k], Antiderivative(Sum(Coeff(~Pq, ~x, 2k) * (~x) ^ (2k), [k, 0, q / 2]) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x) + Antiderivative(~x * Sum(Coeff(~Pq, ~x, 2k + 1) * (~x) ^ (2k), [k, 0, (q - 1) / 2]) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (PolyQ(~Pq, ~x) && Not(PolyQ(~Pq, (~x) ^ 2)))
    @apply_utils Antiderivative(~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => With([d = Coeff(~Pq, ~x, 0), e = Coeff(~Pq, ~x, 2), f = Coeff(~Pq, ~x, 4)], (d * ~x * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1)) / ~a <-- EqQ(~a * e - ~b * d * (2 * ~p + 3), 0) && EqQ(~a * f - ~c * d * (4 * ~p + 5), 0)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (PolyQ(~Pq, (~x) ^ 2) && EqQ(Expon(~Pq, ~x), 4))
    @apply_utils Antiderivative(~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => With([d = Coeff(~Pq, ~x, 0), e = Coeff(~Pq, ~x, 2), f = Coeff(~Pq, ~x, 4), g = Coeff(~Pq, ~x, 6)], (~x * (3 * ~a * d + (~a * e - ~b * d * (2 * ~p + 3)) * (~x) ^ 2) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1)) / (3 * (~a) ^ 2) <-- EqQ(3 * (~a) ^ 2 * g - ~c * (4 * ~p + 7) * (~a * e - ~b * d * (2 * ~p + 3)), 0) && EqQ((3 * (~a) ^ 2 * f - 3 * ~a * ~c * d * (4 * ~p + 5)) - ~b * (2 * ~p + 5) * (~a * e - ~b * d * (2 * ~p + 3)), 0)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (PolyQ(~Pq, (~x) ^ 2) && EqQ(Expon(~Pq, ~x), 6))
    @apply_utils Antiderivative(~Pq / (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => Antiderivative(ExpandIntegrand(~Pq / (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (PolyQ(~Pq, (~x) ^ 2) && Expon(~Pq, (~x) ^ 2) > 1)
    @apply_utils Antiderivative(~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => ((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ FracPart(~p) / ((4 * ~c) ^ IntPart(~p) * (~b + 2 * ~c * (~x) ^ 2) ^ (2 * FracPart(~p)))) * Antiderivative(~Pq * (~b + 2 * ~c * (~x) ^ 2) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (PolyQ(~Pq, (~x) ^ 2) && (Expon(~Pq, (~x) ^ 2) > 1 && EqQ((~b) ^ 2 - 4 * ~a * ~c, 0)))
    @apply_utils Antiderivative(~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => With([d = Coeff(PolynomialRemainder(~Pq, ~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4, ~x), ~x, 0), e = Coeff(PolynomialRemainder(~Pq, ~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4, ~x), ~x, 2)], (~x * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1) * ((~a * ~b * e - d * ((~b) ^ 2 - 2 * ~a * ~c)) - ~c * (~b * d - 2 * ~a * e) * (~x) ^ 2)) / (2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / (2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1) * ExpandToSum((((2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) * PolynomialQuotient(~Pq, ~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4, ~x) + (~b) ^ 2 * d * (2 * ~p + 3)) - 2 * ~a * ~c * d * (4 * ~p + 5)) - ~a * ~b * e) + ~c * (4 * ~p + 7) * (~b * d - 2 * ~a * e) * (~x) ^ 2, ~x), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (PolyQ(~Pq, (~x) ^ 2) && (Expon(~Pq, (~x) ^ 2) > 1 && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && LtQ(~p, -1))))
    @apply_utils Antiderivative(~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => With([q = Expon(~Pq, (~x) ^ 2), e = Coeff(~Pq, (~x) ^ 2, Expon(~Pq, (~x) ^ 2))], (e * (~x) ^ (2q - 3) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1)) / (~c * (2q + 4 * ~p + 1)) + (1 / (~c * (2q + 4 * ~p + 1))) * Antiderivative((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p * ExpandToSum(((~c * (2q + 4 * ~p + 1) * ~Pq - ~a * e * (2q - 3) * (~x) ^ (2q - 4)) - ~b * e * ((2q + 2 * ~p) - 1) * (~x) ^ (2q - 2)) - ~c * e * (2q + 4 * ~p + 1) * (~x) ^ (2q), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (PolyQ(~Pq, (~x) ^ 2) && (Expon(~Pq, (~x) ^ 2) > 1 && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(LtQ(~p, -1)))))
    @apply_utils Antiderivative(~Pq * (~Q4) ^ ~p, ~x) => With([a = Coeff(~Q4, ~x, 0), b = Coeff(~Q4, ~x, 1), c = Coeff(~Q4, ~x, 2), d = Coeff(~Q4, ~x, 3), e = Coeff(~Q4, ~x, 4)], Subst(Antiderivative(SimplifyIntegrand(ReplaceAll(~Pq, (~x->begin
                                                #= none:1 =#
                                                #= none:3 =#
                                                -d / (4 * e) + ~x
                                            end)) * (((a + d ^ 4 / (256 * e ^ 3)) - (b * d) / (8 * e)) + (c - (3 * d ^ 2) / (8 * e)) * (~x) ^ 2 + e * (~x) ^ 4) ^ ~p, ~x), ~x), ~x, d / (4 * e) + ~x) <-- EqQ((d ^ 3 - 4 * c * d * e) + 8 * b * e ^ 2, 0) && NeQ(d, 0)) <-- FreeQ(~p, ~x) && (PolyQ(~Pq, ~x) && (PolyQ(~Q4, ~x, 4) && Not(IGtQ(~p, 0))))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.2 Quartic/1.2.2.5 P(x) (a+b x^2+c x^4)^p.jl


text:
# some tests are mysteriously broken with certain hardware/software.
# use a horrible macro to mark those tests as "potentially broken"
@eval Test begin
    export @test_maybe_broken

    macro test_maybe_broken(ex, kws...)
        test_expr!("@test_maybe_broken", ex, kws...)
        orig_ex = Expr(:inert, ex)
        result = get_test_result(ex, __source__)
        quote
            x = $result
            if x.value
                do_test(x, $orig_ex)
            else
                do_broken_test(x, $orig_ex)
            end
        end
    end
end


text:
## Create a list of flights based on the Flight id, this output should include the passenger Id,
# relevant IATA/FAA codes, the departure time, the arrival time
# (times to be converted toHH:MM:SS format), and the flight times.
NMR.clean_split_dir()
jobid = 2
NMR.runjob(NMR.Nmr(jobid, ["AComp_Passenger_data.csv"], NMR.mapper_parserecordacomp, "acompdata.json"))
output = NMR.runcombiner("acompdata.json", NMR.combiner_parsejson)

# store output
NMR.save_job_output(jobid, "", output)


text:
@recipe function f(mds::MDS{<:Real}; mds_axes=(1,2))
    length(mds_axes) in [2,3] || throw(ArgumentError("Can only accept 2 or 3 mds axes"))
    xax = mds_axes[1]
    yax = mds_axes[2]
    ev = eigvals(mds)
    var_explained = [v / sum(ev) for v in ev]
    tfm = collect(transform(mds)')

    xlabel --> "MDS$xax"
    ylabel --> "MDS$yax"
    seriestype := :scatter
    aspect_ratio --> 1

    if length(mds_axes) == 3
        zax = mds_axes[3]
        zlabel --> "MDS$zax"
        tfm[:,xax], tfm[:,yax], tfm[:,zax]
    else
        tfm[:,xax], tfm[:,yax]
    end
end

#= This needs to wait on a different PCA API in MultivariateStats.jl
@recipe function f(pca::PCA{<:Real}; pca_axes=(1,2))
end
=#


text:
using LightGraphs
using MetaGraphs

#=
definition.jl
this script contains the basic definitions and constructors for the LPBasic problem
=#
export Line, ProblemInstance, LPBasic, VerifyInstance

#structs declarations:
#defines what a line is
"""
    Line
A basic data structure for model construction. Represents a directed walk through a transportation network

# Constructor:
    Line(id::Union{Int64, String}, walk::Array{Tuple{Int64, Int64}, 1}, frequency::Int64, cost::Float64)

# Fields:
- id::Union{Int64, String} - an Int64 or String to help distinguis this lines from other lines. Uniqueness is not enforced
- walk::Array{Tuple{Int64, Int64}, 1} - an array of edge indices, representing a directed walk through the transportation network
- frequency::Int64 - the base frequency for this line
- cost::Float64 - the cost per frequency unit of this line
"""
struct Line
    id::Union{Int64, String}
    #a walk is defined as a sequence of *edges*
    walk::Array{Tuple{Int64, Int64}, 1}
    frequency::Int64
    #the cost *per frequency unit*
    cost::Float64

    Line(id, walk, frequency, cost) =

        if length(walk) < 1 error("Invalid line: walk must contain at least one edge. Id:" * string(id))
        elseif frequency <= 0 error("Invalid Line: frequency must be strictly positive: " * string(frequency))
        elseif cost < 0 error("Invalid Line: cost must be positive: " * string(cost))
        else new(id, walk, frequency, cost) end

end

"""
    VerifyInstance(network::MetaDiGraph, lines::Array{Line,1})
    VerifyInstance(instance::ProblemInstance)

Verify if a given problem instance is valid per the LPBasic requirements. Raise error if there are inconsistensies, otherwise return true.

This includes checking if:
- all lines are unique(unique id, unique walk, unique frequency)
- each edge has defined its expected properties
- there are no invalid values for frequencies and weights

# to-do:
- is the network/graph is connected?(is there is a path between every pair of vertices?)
- redundant information matches one another? (:line_coverage)

# return value:
- true - indicating that the instance is valid, or
- *raise error* - indicating that the instance is not valid and the reason why
"""
function VerifyInstance(network::MetaDiGraph, lines::Array{Line,1})

    #verifies that any redundant information on the LPBasic instance actually matches. This should help detect any errors in the reading of input data
    #lines listed for a given edge as 'covering' actually cover the given edge (to-do, line number is not unique yet)

    #verifies that at least one line is defined
    if length(lines) <= 0 error("Invalid instance: 0 lines are defined") end

    #verifies that all lines are unique
    setLines = Set{Line}()
    for line in lines
        if in(line, setLines)
            error("Invalid instance: repeated lines detected")
        else
            push!(setLines, line)
        end
    end
    #ps: validation of each line is done elsewhere(on struct constructor)

    #network validation
    #verifies that at least two nodes have been defined
    if nv(network) < 2
        error("Invalid instance: network must contain at least two nodes")
    end
    #verifies that at least one edge exists
    if ne(network) < 1
        error("Invalid instance: network must contain at least one edge")
    end

    #verifies that each edge has defined its expected properties
    #weight, f_min, f_max
    for edge in edges(network)
        base_freq, f_min, f_max, weight,line_coverage = nothing, nothing, nothing, nothing, nothing
        try
            base_freq = get_prop(network, edge, :base_freq)
            f_min = get_prop(network, edge, :f_min)
            f_max = get_prop(network, edge, :f_max)
            weight = get_prop(network, edge, :weight)
            line_coverage = get_prop(network, edge, :line_coverage)
        catch e
            #at least one of the properties was not defined
            error("Expected properties (base_freq, f_min, f_max, weight) were not defined as expected for edge " + edge)
        end

        #if they have been defined, type check them
        #to-do: should f_min and f_max really be restricted to integer values?
        if !(base_freq isa Integer) error("base_freq was not defined as an integer in " * repr(edge)) end
        if !(f_min isa Integer) error("f_min was not defined as an integer in " * repr(edge)) end
        if !(f_max isa Integer) error("f_max was not defined as an integer in " * repr(edge)) end
        if !(weight isa Number) error("weight was not defined as a number in " * repr(edge)) end
        if !(line_coverage isa Dict{Int64,Int64}) error("line_coverage was not defined as a Dict{Int64,Int64} in " * repr(edge)) end

    end

    #checks redundancy on :line_coverage (to-do)

    #if no problems were detected, instance is considered valid
    return true
end

"""
    ProblemInstance
Struct containing all data of a single LPBasic problem's instance.

# Constructor:
    ProblemInstance(network::MetaDiGraph, lines::Array{Line,1})

# Fields:
- network::MetaDiGraph - representation the transportation network
- lines::Array{Line,1} - array containing all lines to be considered in the optimization
"""
struct ProblemInstance
	#a metadigraph with edges having the properties:
	#weight, f_min, f_max, (coverage_list)
    network::MetaDiGraph
    #an array of lines
    lines::Array{Line,1}

    ProblemInstance(network, lines) = VerifyInstance(network, lines) && new(network,lines) #short circuit: returns the second value if the first is true

end

#second method for VerifyInstance function
"""
    VerifyInstance(instance::ProblemInstance)

Do the same as VerifyInstance(instance.network, instance.lines)
"""
function VerifyInstance(instance::ProblemInstance)
    VerifyInstance(instance.network, instance.lines)
end


text:
module SpringCollab2020DiagonalWingedStrawberry

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/diagonalWingedStrawberry" DiagonalWingedStrawberry(x::Integer, y::Integer, order::Integer=-1, checkpointID::Integer=-1)

const placements = Ahorn.PlacementDict(
	"Diagonal Winged Strawberry (Spring Collab 2020)" => Ahorn.EntityPlacement(
		DiagonalWingedStrawberry
	)
)

function Ahorn.selection(entity::DiagonalWingedStrawberry)
	x, y = Ahorn.position(entity)

	return Ahorn.getSpriteRectangle("collectables/strawberry/wings01", x, y)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::DiagonalWingedStrawberry, room::Maple.Room)
	x, y = Ahorn.position(entity)

	Ahorn.drawSprite(ctx, "collectables/strawberry/wings01", x, y)
end

end


text:
# testing PIT calculation
using myJuliaUtils
# using StatsBase
using Test
# using BenchmarkTools

a = randn(10, 10000)
b = zeros(10)
Q = pspit(a, b)
@test length(Q.z) == length(b)
@test sum(Q.zhist.weights) == 10

nbins = rand(1:20)
a = randn(100, 10000)
b = zeros(100)
Q = pspit(a, b, nbins = nbins)
@test length(Q.z) == length(b)
@test length(Q.zhist.weights) == nbins
@test sum(Q.zhist.weights) == 100
# @btime pspit(a,b)

nbins = rand(1:20)
d = rand(100:200)
ndraws = rand(10000:20000)
a = randn(d, ndraws)
b = randn(d)
Q = pspit(a, b, nbins = nbins)
@test length(Q.z) == length(b)
@test length(Q.zhist.weights) == nbins
@test sum(Q.zhist.weights) == d

# Q = pspit(a, b, nbins = nbins, returnZ=true)

# using DelimitedFiles
# AAA=readdlm("/Users/zymantas/Downloads/testPIT.txt", ',');
# pspit(AAA[:,1:end-1], AAA[:,end])

a = randn(100, 10000)
b = zeros(100)
Q = pscrps(a, b)
@test length(Q.crps_individual) == length(b)

max_stars_repo_name: zymbuzz/myJuliaUtils

max_issues_repo_name: zymbuzz/myJuliaUtils

max_forks_repo_name: zymbuzz/myJuliaUtils


text:
module AWSClusterManagers

using AWSBatch: JobQueue, max_vcpus, run_batch
using Dates: Dates, Period, Minute, Second, Millisecond
using Distributed: Distributed, ClusterManager, WorkerConfig, cluster_cookie, start_worker
using JSON: JSON
using Memento: Memento, getlogger, warn, notice, info, debug
using Mocking: Mocking, @mock
using Sockets: IPAddr, IPv4, @ip_str, accept, connect, getipaddr, listen, listenany

export AWSBatchManager, AWSBatchNodeManager, DockerManager, start_batch_node_worker

const LOGGER = getlogger(@__MODULE__)

function __init__()
    # https://invenia.github.io/Memento.jl/latest/faq/pkg-usage/
    return Memento.register(LOGGER)
end

include("socket.jl")
include("container.jl")
include("batch.jl")
include("batch_node.jl")
include("docker.jl")

end  # module


text:
include("rob_lib.jl")
# and so this talse begins we mark all thr mark_corners
function mark_corners!(r::Robot)
    num_verticale = moves!(r, Sud)
    num_horizont = moves!(r, West)

    for i in 3:6
        putmarker!(r)
        movements!(r, HorizonSide(i % 4)) # that is out pioritat cuuse  1 0 3 2
    end

    Back!(r, Nord, num_verticale)
    Back!(r, Ost, num_horizont)
end

function moves!(r::Robot, side::HorizonSide)
    num_steps = 0
    while move_if_possible!(r, side)
        num_steps += 1
    end
    return num_steps
end

max_stars_repo_name: NikitaDestroyer/Basics

max_issues_repo_name: NikitaDestroyer/Basics

max_forks_repo_name: NikitaDestroyer/Basics


text:
export FunctionObservable, make

import Base: show

"""
    FunctionObservable{D}(f::F)

    FunctionObservable wraps a callback `f`, which is called when the Observable is initially subscribed to.
    This function is given an Actor, to which new values can be nexted (with `next!(actor, data)`),
    or an `error!` method can be called to raise an error, or `complete!` can be called to notify of a successful completion.

    # Arguments
    - `f::F`: function to be invoked on subscription

    See also: [`Subscribable`](@ref), [`make`](@ref)
"""
struct FunctionObservable{D, F} <: Subscribable{D}
    f :: F
end

function on_subscribe!(observable::FunctionObservable{D}, actor::A) where { D, A }
    wrapper = FunctionObservableActorWrapper{D, A}(false, actor)
    subscription = FunctionObservableSubscription{D, A}(wrapper)
    observable.f(wrapper)
    return subscription
end

mutable struct FunctionObservableActorWrapper{D, A} <: Actor{D}
    is_unsubscribed :: Bool
    actor           :: A
end

function on_next!(wrapper::FunctionObservableActorWrapper{D}, data::D) where D
    if !wrapper.is_unsubscribed
        next!(wrapper.actor, data)
    end
end

function on_error!(wrapper::FunctionObservableActorWrapper, err)
    if !wrapper.is_unsubscribed
        error!(wrapper.actor, err)
    end
end

function on_complete!(wrapper::FunctionObservableActorWrapper)
    if !wrapper.is_unsubscribed
        complete!(wrapper.actor)
    end
end

struct FunctionObservableSubscription{D, A} <: Teardown
    wrapper :: FunctionObservableActorWrapper{D, A}
end

as_teardown(::Type{<:FunctionObservableSubscription}) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::FunctionObservableSubscription)
    subscription.wrapper.is_unsubscribed = true
    return nothing
end

"""
    make(f::Function, type::Type{D})

Creation operator for the `FunctionObservable`.

# Arguments
- `f`: function to be invoked on subscription
- `type`: type of data in observable

# Examples
```jldoctest
using Rocket

source = make(Int) do actor
    next!(actor, 0)
    complete!(actor)
end

subscription = subscribe!(source, logger());
unsubscribe!(subscription)
;

# output

[LogActor] Data: 0
[LogActor] Completed

```

```jldoctest
using Rocket

source = make(Int) do actor
    next!(actor, 0)
    setTimeout(100) do
        next!(actor, 1)
        complete!(actor)
    end
end

subscription = subscribe!(source, logger())
unsubscribe!(subscription)
;

# output

[LogActor] Data: 0
```

See also: [`FunctionObservable`](@ref), [`subscribe!`](@ref), [`logger`](@ref)
"""
make(f::F, type::Type{D}) where { D, F } = FunctionObservable{D, F}(f)

Base.show(io::IO, ::FunctionObservable{D}) where D  = print(io, "FunctionObservable($D)")
Base.show(io::IO, ::FunctionObservableSubscription) = print(io, "FunctionObservableSubscription()")


text:
using LinearAlgebra

function gradient_descent(P, q, x₀; α=0.1,
                maxiter=1000, ϵ=1e-5)
    x = copy(x₀)
    ∇f = x -> P * x + q
    Δx = -∇f(x)
    iter = 0
    while norm(Δx) > ϵ || iter ≤ maxiter
        iter += 1
        x .+= α * Δx
        Δx .= - ∇f(x)
    end
    return x
end

P = [10.0 -1.0;
    -1.0 1.0];
q = [0; -10.0];
x₀ = zeros(2);


text:
module PD

using ..Cos

include("PDDoc.jl")
include("PDPageElement.jl")
include("PDFontMetrics.jl")
include("PDFonts.jl")
include("PDXObject.jl")
include("PDDocImpl.jl")
include("PDPage.jl")
include("PDFontTables.jl")
include("PDOutline.jl")
include("PDSignature.jl")

end


text:
using GslibIO
using FileIO
using Meshes
using GeoStatsBase
using Test

# environment settings
datadir = joinpath(@__DIR__,"data")

@testset "GslibIO.jl" begin
  @testset "Basic" begin
    fname = tempname()*".gslib"

    props3D = rand(10,10,10), rand(10,10,10)
    props2D = rand(10,10), rand(10,10)

    for (prop1, prop2) in [props2D, props3D]
      FileIO.save(fname, [prop1,prop2])
      grid = FileIO.load(fname)
      @test grid[:prop1] == vec(prop1)
      @test grid[:prop2] == vec(prop2)

      FileIO.save(fname, grid)
      grid = FileIO.load(fname)
      @test grid[:prop1] == vec(prop1)
      @test grid[:prop2] == vec(prop2)

      FileIO.save(fname, prop1)
      grid = FileIO.load(fname)
      @test grid[:prop1] == vec(prop1)
    end

    rm(fname)
  end

  @testset "LegacyParser" begin
    fname = joinpath(datadir, "legacy_grid.gslib")

    spec = GslibIO.parse_legacy(fname)

    @test isequal(spec.data[:, 1], [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8])
    @test isequal(spec.data[:, 2], [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0])
    @test isequal(spec.data[:, 3], [0.8, 0.7, 0.8, 0.7, 0.8, 0.7, 0.8, -999.0])

    @test spec.varnames == [:Porosity, :Lithology, Symbol("Water Saturation")]

    fname = joinpath(datadir,"legacy_pset.gslib")

    spec = GslibIO.parse_legacy(fname)

    @test spec.varnames == [:East, :North, :Elevation, :Porosity]

    @test isequal(spec.data[:, 1], [10.0, 20.0, 30.0, 40.0])
    @test isequal(spec.data[:, 2], [11.0, 21.0, 31.0, 41.0])
    @test isequal(spec.data[:, 3], [12.0, 22.0, 32.0, 42.0])
    @test isequal(spec.data[:, 4], [0.1, 0.2, 0.3, 0.4])
  end

  @testset "LegacyGrid" begin
    fname = joinpath(datadir,"legacy_grid.gslib")

    sdata = GslibIO.load_legacy(fname, (2,2,2))

    @test size(domain(sdata)) == (2,2,2)
    @test minimum(domain(sdata)) == Point(0.,0.,0.)
    @test spacing(domain(sdata)) == [1.,1.,1.]

    por = sdata[:Porosity]
    lit = sdata[:Lithology]
    sat = sdata[Symbol("Water Saturation")]
    @test isequal(por, [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8])
    @test isequal(lit, [1,2,3,4,5,6,7,8])
    @test isequal(sat, [0.8,0.7,0.8,0.7,0.8,0.7,0.8,NaN])

    # test if storing/loading recovers data
    fname = tempname()*".gslib"
    GslibIO.save_legacy(fname, sdata, na=-999)
    ndata = GslibIO.load_legacy(fname, (2,2,2), na=-999)
    @test isequal(domain(sdata), domain(ndata))
    @test isequal(values(sdata), values(ndata))

    rm(fname)
  end

  @testset "LegacyPointSet" begin
    fname = joinpath(datadir,"legacy_pset.gslib")

    sdata = GslibIO.load_legacy(fname, (:East, :North, :Elevation))

    sdomain = domain(sdata)
    @test embeddim(sdomain) == 3
    @test nelements(sdomain) == 4

    por = sdata[:Porosity]
    @test centroid(sdomain, 1) == Point(10.0, 11.0, 12.0)
    @test centroid(sdomain, 2) == Point(20.0, 21.0, 22.0)
    @test centroid(sdomain, 3) == Point(30.0, 31.0, 32.0)
    @test centroid(sdomain, 4) == Point(40.0, 41.0, 42.0)
    @test isequal(por, [0.1, 0.2, 0.3, 0.4])

    # test when coordnames are not in varnames
    @test_throws AssertionError GslibIO.load_legacy(fname, (:x, :y, :Elevation))

    # test if storing/loading recovers data
    fname = tempname()*".gslib"
    GslibIO.save_legacy(fname, sdata, coordnames=(:x, :y, :Elevation))
    ndata = GslibIO.load_legacy(fname, (:x, :y, :Elevation))
    @test sdata == ndata

    rm(fname)
  end

  @testset "LegacyInvalid" begin
    fname = joinpath(datadir,"legacy_invalid.gslib")

    @test_throws MethodError GslibIO.load_legacy(fname, (:x, :y, :z))
  end
end


text:
# Linear covariance function

@inline dotijp(X::MatF64, i::Int, j::Int, p::Int) = X[p,i]*X[p,j]
@inline function dotij(X::MatF64, i::Int, j::Int, dim::Int)
	s=zero(eltype(X))
	@inbounds @simd for p in 1:dim
		s+=dotijp(X,i,j,p)
	end
	return s
end
include("lin_iso.jl")
include("lin_ard.jl")

"""
    Lin(ll::Union{Float64,Vector{Float64}})

Create linear kernel with length scale `exp.(ll)`.

See also [`LinIso`](@ref) and [`LinArd`](@ref).
"""
Lin(ll::Float64) = LinIso(ll)
Lin(ll::Vector{Float64}) = LinArd(ll)


text:
##
# Authentication mechanisms
# Only SASL-Plain supported for now
type AccumuloAuth
    mechanism::String
    callback::Function

    AccumuloAuth(mechanism::String=SASL_MECH_PLAIN, callback::Function=Thrift.sasl_callback_default) = new(mechanism, callback)
end

function AccumuloAuthSASLPlain(uid::String, passwd::String; zid::String="")
    function callback(part::Symbol)
        (part == :authcid) && (return uid)
        (part == :passwd) && (return passwd)
        (part == :show) && (return uid)
        (part == :mechanism) && (return "SASL-Plain")
        return zid
    end
    AccumuloAuth(SASL_MECH_PLAIN, callback)
end

function show(io::IO, auth::AccumuloAuth)
    uid = auth.callback(:show)
    mech = auth.callback(:mechanism)
    println(io, "AccumuloAuth ($mech): $uid")
    nothing
end

##
# AccumuloConn holds the thrift connection and protocol objects.
# It also holds the accumulo session handle for this connection.
type AccumuloConn
    transport::TTransport
    protocol::TProtocol
    client::AccumuloProxyClient
    handle::Vector{UInt8}
    connstr::String

    function AccumuloConn(host::String, port::Integer, auth::AccumuloAuth; tprotocol::Symbol=:compact)
        transport = TFramedTransport(TSocket(host, port))
        if tprotocol === :binary
            protocol = TBinaryProtocol(transport, true)
        elseif tprotocol === :compact
            protocol = TCompactProtocol(transport)
        else
            error("Unsupported protocol: $tprotocol")
        end
        client = AccumuloProxyClient(protocol)
        uid = auth.callback(:show)
        connstr = "accumulo://$(uid)@$(host):$port"
        new(transport, protocol, client, connect(transport, client, auth), connstr)
    end

    function connect(transport::TTransport, client::AccumuloProxyClient, auth::AccumuloAuth)
        open(transport)
        uid = auth.callback(:authcid)
        creds = Dict{String,String}("password" => auth.callback(:passwd))
        login(client, uid, creds)
    end
end

function show(io::IO, conn::AccumuloConn)
    println(io, conn.connstr)
    nothing
end

function close(conn::AccumuloConn)
    close(conn.transport)
    nothing
end

#
# AccumuloSession holds a connection and session status
type AccumuloSession
    conn::AccumuloConn

    function AccumuloSession(host::String="localhost", port::Integer=42424, auth::AccumuloAuth=AccumuloAuth(); tprotocol::Symbol=:compact)
        new(AccumuloConn(host, port, auth; tprotocol=tprotocol))
    end
end

function show(io::IO, sess::AccumuloSession)
    print(io, "AccumuloSession: ")
    show(io, sess.conn)
end
close(session::AccumuloSession) = close(session.conn)
client(session::AccumuloSession) = session.conn.client
handle(session::AccumuloSession) = session.conn.handle


text:
#===============================================================================

    Routines for computing the log-likelihood functions and solve!'ing
    the path from the Wiener path and computing the log-likelihood at
    the same time.

===============================================================================#

@doc raw"""
    loglikhd([::IntegrationRule=::LeftRule], P::GuidProp, X::Trajectory; skip=0)

Compute path contribution to the log-likelihood function, i.e.:
```math
\int_0^T G(t, X_t) dt,
```
where
```math
G(t,x):=\left[
    (b-\tilde{b})^T\tilde{r}
    + 0.5 tr\left\{
        (a-\tilde{a})(\tilde{r}\tilde{r}^T-H)
    \right\}
\right](t,x),
```
and
$\tilde{r}(t,x):=\nabla\log\tilde{\rho}(t,x)$.

    loglikhd(
        [::IntegrationRule=::LeftRule],
        PP::AbstractArray{<:GuidProp}, XX::AbstractArray{<:Trajectory};
        skip=0
    )
Compute path contribution to the log-likelihood function for a sequence of
segments.
"""
function loglikhd end

function loglikhd(P::GuidProp, X::Trajectory; skip=0)
    loglikhd(LeftRule(), P, X, P.guiding_term_solver; skip=skip)
end

function loglikhd(ir::IntegrationRule, P::GuidProp, X::Trajectory; skip=0)
    loglikhd(ir, P, X, P.guiding_term_solver; skip=skip)
end

function loglikhd(
        ::LeftRule,
        P::GuidProp,
        X::Trajectory,
        ::AbstractGuidingTermSolver{:outofplace};
        skip=0
    )
    tt, xx = X.t, X.x
    ll = 0.0
    N = length(tt)-1-skip

    for i in 1:N
        x, s, dt = xx[i], tt[i], tt[i+1]-tt[i]
        r_i = ∇logρ(i, x, P)
        b_i = DD.b(s, x, P.P_target)
        btil_i = DD.b(s, x, P.P_aux)

        ll += dot(b_i-btil_i, r_i) * dt

        if !DD.constdiff(P)
            H_i = H(i, x, P)
            a_i = DD.a(s, x, P.P_target)
            atil_i = DD.a(s, x, P.P_aux)
            ll += 0.5 * tr( (a_i - atil_i)*(r_i*r_' - H_i) ) * dt
        end
    end
    ll
end

# NOTE worry about this later...
function loglikhd(
        ::LeftRule,
        P::GuidProp,
        X::Trajectory,
        ::AbstractGuidingTermSolver{:inplace};
        skip=0
    )
    tt, xx = X.t, x.x
    ll = 0.0
    N = length(tt)-1-skip

    for i in 1:N
        x, s, dt = xx[i], tt[i], tt[i+1]-tt[i]
        ∇logρ!(P.buffer.r_i, i, x, P)
        DD.b(P.buffer.b_i, s, x, P.P_target)
        DD.b(s.buffer.btil_i, x, P.P_aux)
        for j in 1:length(P.buffer.b_i)
            P.buffer.b_i[j] -= P.buffer.btil_i[j]
        end

        ll += dot(P.buffer.b_i, P.buffer.r_i) * dt

        if !DD.constdiff(P)
            H_i = H(i, x, P)
            DD.a(P.buffer.a_i, s, x, P.P_target)
            DD.a(P.buffer.atil_i, s, x, P.P_aux)
            ll -= 0.5*sum( (a_i - atil_i).*H_i ) * dt
            ll += 0.5*( r_i'*(a_i - atil_i)*r_i ) * dt
        end
    end
    ll
end

function loglikhd(
        PP::AbstractArray{<:GuidProp},
        XX::AbstractArray{<:Trajectory};
        skip=0
    )
    loglikhd(LeftRule(), PP, XX; skip=skip)
end

function loglikhd(
        ir::IntegrationRule,
        PP::AbstractArray{<:GuidProp},
        XX::AbstractArray{<:Trajectory};
        skip=0
    )
    ll_tot = loglikhd_obs(PP[1], XX[1].x[1])
    for i in eachindex(PP)
        ll_tot += loglikhd(ir, PP[i], XX[i], PP[i].guiding_term_solver; skip=skip)
    end
    ll_tot
end

@doc raw"""
    loglikhd_obs(P::GuidProp, x0)

Compute the contribution of end-points to the log-likelihood function, i.e.:
$\log\rho(t,x)$.
"""
loglikhd_obs(P::GuidProp, x0) = loglikhd_obs(P, x0, P.guiding_term_solver)

function loglikhd_obs(P::GuidProp, x0, ::AbstractGuidingTermSolver{:outofplace})
    - 0.5 * ( x0'*H(P, 1)*x0 - 2.0*dot(F(P, 1), x0) ) - c(P, 1)
end

# worry about it later
function loglikhd_obs(P::GuidProp, x0, ::AbstractGuidingTermSolver{:inplace})
    - 0.5 * ( x0'*H(P, 1)*x0 - 2.0*dot(F(P, 1), x₀) ) - c(P, 1)
end

#===============================================================================
            solve! with simultaneous computation of log-likelihood
===============================================================================#
"""
    solve_and_ll!(X, W, P, y1)

Compute the trajectory under the law `P` for a given Wiener noise `W` and
a starting point `y1`. Store the trajectory in `XX`. Compute the log-likelihood
(path contribution) along the way. Return `success_flag` and log-likelihood.
`success_flag` is set to false only if sampling was prematurely halted due to
`XX` violating assumptions about state space.
"""
function solve_and_ll!(X, W, P::GuidProp, y1; skip=0)
    solve_and_ll!(X, W, P, P.guiding_term_solver, y1; skip=skip)
end

function solve_and_ll!(XX, WW, PP::AbstractArray{<:GuidProp}, y1; skip=0)
    ll_tot = loglikhd_obs(PP[1], y1)
    for i in eachindex(PP)
        success, ll = solve_and_ll!(XX[i], WW[i], PP[i], PP[i].guiding_term_solver, y1; skip=skip)
        success || return false, ll
        ll_tot += ll
        y1 = XX[i].x[end]
    end
    true, ll_tot
end

function solve_and_ll!(
        XX,
        WW,
        P::GuidProp,
        ::AbstractGuidingTermSolver{:outofplace},
        y1;
        skip=0
    )
    yy, ww, tt = XX.x, WW.x, XX.t
    N = length(XX)

    yy[1] = DD.value(y1)
    x = y1
    ll = 0.0
    for i in 1:(N-1)
        add_to_ll = (i < N-skip)
        s = tt[i]
        dt = tt[i+1] - tt[i]
        dW = ww[i+1] - ww[i]

        r_i = ∇logρ(i, x, P)
        b_i = DD.b(s, x, P.P_target)
        btil_i = DD.b(s, x, P.P_aux)

        σ_i = DD.σ(s, x, P.P_target)
        a_i = σ_i*σ_i'

        add_to_ll && (ll += dot(b_i-btil_i, r_i) * dt)

        if !DD.constdiff(P)
            H_i = H(i, x, P)
            atil_i = DD.a(s, x, P.P_aux)
            add_to_ll && (ll += 0.5*tr( (a_i - atil_i)*(r_i*r_i'-H_i') ) * dt)
        end

        x = x + (a_i*r_i + b_i)*dt + σ_i*dW

        yy[i+1] = DD.value(x)

        DD.bound_satisfied(P, yy[i+1]) || return false, -Inf
    end
    true, ll
end

#NOTE worry about it later
function solve_and_ll!(
        XX,
        WW,
        P::GuidProp,
        ::AbstractGuidingTermSolver{:inplace},
        y1;
        skip=0
    )
end


text:
#!/usr/bin/env julia
# Copyright (c) 2016, Jost von Hardenberg - ISAC-CNR, Italy
using RainFARM
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--radius", "-r"
            help = "Smoothing radius (in grid units or pixels (if negative))"
            arg_type = Float64
            default = 0.
        "--gaussian", "-g"
            help = "Use a Gaussian kernel instead of a top hat. The radius parameter mean sigma in that case."
            action = :store_true
        "--varname", "-v"
            help = "Input variable name (in orofile)"
            arg_type = AbstractString
            default = ""
        "infile"
            help = "The file to downscale"
            arg_type = AbstractString
            required = true
        "outfile"
            help = "The output file name"
            arg_type = AbstractString
            required = true
    end

    s.description = "Smooth input file using convolution with a circle of fixed radius"
    s.version = "0.1"
    s.add_version = true

    return parse_args(s)
end

args = parse_commandline()
radius = args["radius"]
fgauss = args["gaussian"]
filein = args["infile"]
fileout = args["outfile"]
varname = args["varname"]

(tin0, lonl0, latl0) = read_netcdf2d(filein, varname)
if length(size(lonl0))>1
    dxl = max(lonl0[2,1]-lonl0[1,1], lonl0[1,2]-lonl0[1,1])
else
    dxl = lonl0[2]-lonl0[1]
end

(tin, lonl, latl, varname) = read_netcdf2d(filein, varname)
nx = size(tin, 1)
ny = size(tin, 2)
nt = size(tin, 3)

if fgauss && (nx!=ny)
    println("Domain is ", nx, " by ", ny, " pixel")
    println("If Gaussian smoothing is used then the domain should be square!")
    exit(1)
end

if radius>0
    nf2 = div(radius, dxl)
else
    nf2 = -radius
end

println("Original resolution: ", dxl)

if fgauss
    nas = div(nx, nf2*2)
    println("Gaussian smoothing, nas=", nas)
    println("Smoothing with sigma ", nf2, " pixel")
    for i = 1:nt
        tin[:,:,i] = smoothconv(tin[:,:,i], nas)
    end
else
   println("Smoothing with radius ",nf2 ," pixel")
   for i=1:nt
       tin[:,:,i]=smooth(tin[:,:,i],nf2)
    end
end

run(`cp $filein $fileout`)
overwrite_netcdf2d(fileout, tin ,varname)


text:
# Primitive cells used for testing

# a=1 for fcc, cubic and bcc
cellfcc = [-0.5 0 0.5; 0 0.5 0.5; -0.5 0.5 0];
cellcubic = [1.0 0 0; 0 1.0 0; 0 0 1.0];
#cellbcc = [-0.5 0.5 0.5; 0.5 -0.5 0.5; 0.5 0.5 -0.5];
# Through some trial and error testing, I determed that the "bcc" cell used in Table 3 of Militzer must have contained a typo. The real bcc cell was not working (table values did not match up). After seeing other cells work and scouring the code, I could not find the answer. Finally, I checked an invalid primitive cell and I think I found the issue. All the numbers (rmax, dmin, cell angles...) match up when the basis vectors are created with one basis vector missing a negative sign for one. This creates an primitive cell that has atoms overlapping when the basis vectors are translated in space, but it's technically still a valid primitive cell and can be used to check the algorithm's correctness.
cellbcc = [0.5 0.5 0.5; 0.5 -0.5 0.5; 0.5 0.5 -0.5];
hcp_ca = √(8.0/3.0)  # ideal c/a ratio
# In principle the lattice parameter should be one, but the test values don't
# work out otherwise. Perhaps it's possible that the paper makes a mistake of
# not dividing by the lattice parameter. The dmin reported is 0.561 and by
# guess and check with the data 0.5612 gives reproduces the dmin data the best
# with rounding.
hcp_a = 0.5612
cellhcp = hcp_a.*[1 0 0; -1/2 √(3.0)/2 0; 0 0 hcp_ca];


text:
import PeriodicTable

"""Return the data directory with pseudopotential files"""
datadir_psp() = joinpath(get(ENV, "DFTK_DATADIR", DFTK_DATADIR), "psp")

"""
    load_psp(key; datadir_psp)

Load a pseudopotential file from the library of pseudopotentials.
The file is searched in the directory `datadir_psp` and by the `key`.
If the `key` is a path to a valid file, the extension is used to determine
the type of the pseudopotential file format and a respective class is returned.
"""
function load_psp(key::AbstractString; datadir_psp=datadir_psp(), identifier="")
    if startswith(key, "hgh/") || endswith(key, ".hgh")
        parser = parse_hgh_file
        extension = ".hgh"
    else
        error("Could not determine pseudopotential family of '$key'")
    end
    isfile(key) && return parser(key, identifier=identifier)

    # Not a file: Threat as identifier, add extension if needed
    @assert identifier == ""
    identifier = lowercase(key)
    fullpath = joinpath(datadir_psp, identifier)
    isfile(fullpath) || (fullpath = fullpath * extension)

    if isfile(fullpath)
        parser(fullpath, identifier=identifier)
    else
        error("Could not find pseudopotential for identifier " *
              "'$identifier' in directory '$datadir_psp'")
    end
end

function load_psp(element::Union{Symbol,Integer}; family="hgh", core=:fullcore, kwargs...)
    list = list_psp(element; family=family, core=core, kwargs...)
    if length(list) != 1
        error("Parameters passed to load_psp do not uniquely identify a PSP file.")
    end
    load_psp(list[1].path, identifier=list[1].identifier)
end


text:
module Transpiler

using BetweenFlags

include(joinpath("intrinsics","FortranIntrinsics.jl"))
include(joinpath("intrinsics","PythonIntrinsics.jl"))
include("types.jl")
include("remove_comments.jl")
include("file_managemnet.jl")

extensions_from(::Union{Python,Cython}) = (".pxd", ".py")
extension_from_append(::Union{Python,Cython}) = ".pyx"
extension_from(::Julia) = ".jl"
extension_to(::Python) = ".py"
folder_filters(::Union{Python,Cython}) = "__pycache__"
extension_to(::Julia) = ".jl"

include(joinpath("python_to_julia","python_to_julia.jl"))

end # module


text:
@export Pb(m::AbstractLedger) = m[Pb][1].pb::Float64

@export H_magnitudes(ls::Vector{<:AbstractLedger}) = map(x -> singleton(x, H).magnitude, ls)
@export H_vecs(ls::Vector{<:AbstractLedger}) = map(x -> singleton(x, H).v, ls)

@export Gd_angles(ls::Vector{<:AbstractLedger}, i) = map(x -> x[Spin][i].ϕ, ls)

@export Gd_entities(l::Ledger) = collect(@entities_in(l[Spin]))
@export L_entities(l::Ledger)  = collect(@entities_in(l[L]))

@export L_angles(ls::Vector{<:AbstractLedger}, i) = map(x -> x[L][i].ϕ, ls)
@export L_vecs(ls::Vector{<:AbstractLedger}, i) = map(x -> x[L][i].v, ls)

@export spins(ls::Vector{<:AbstractLedger}, i) = map(x->x[Spin][i].v, ls)

@export Pb(ls::Vector{<:AbstractLedger}) = [m[Pb][1].pb for m in ls]

@export Etot(ls::Vector{<:AbstractLedger}, esym=:E_all) = map(x->calc_etot(x, esym).e, ls)

@export function calc_etot(m::AbstractLedger, sym=:E_all)
    m[Etot].data[1] = Etot(0.0)
    update(stage(m, sym), m)
    return m[Etot].data[1]
end

@export function ϕ_parallel_perpendicular(l::Ledger)
    h = singleton(l, H)
    H_direction = Gd.polar2vec(h.ϕ, 0.0)
    perp_direction = Gd.polar2vec(h.ϕ + π/2, 0.0)
    ϕ_parallel = zeros(8)
    ϕ_perpendicular = zeros(8)
    for i=1:8
        ϕ_parallel[i] = l[Spin][i].v ⋅ H_direction
        ϕ_perpendicular[i] = l[Spin][i].v ⋅ perp_direction
    end
end

@export function M_spins(l::AbstractLedger, is...)
    h = singleton(l, H)
    p = singleton(l, Model)
    H_direction = polar2vec(h.ϕ, 0.0)
    return p.χ_Gd * sum([H_direction ⋅ l[Spin][i].v for i in is])
end

M_spins(l::Ledger) = M_spins(l, 1:length(l[Spin])...)

@export function M_L(l::AbstractLedger, is...)
    h = singleton(l, H)
    p = singleton(l, Model)
    H_direction = polar2vec(h.ϕ, 0.0)
    return mapreduce(x -> - 2*p.χ_L * H_direction ⋅ (l[L][x].v * (l[L][x].v ⋅ h.v)), +, is)
end
M_L(l::Ledger) = M_L(l, 1,2)
@export M(l::AbstractLedger) = M_spins(l) + M_L(l)

@export J_parallel(l::Ledger)      = J_parallel(singleton(l, Model))
@export J_perpendicular(l::Ledger) = J_perpendicular(singleton(l, Model))

@export function all_angles!(angle_vector::AbstractVector{Float64}, m::AbstractLedger)
    for (i, e) in enumerate(@entities_in(m[Optimize]))
        @inbounds if e ∈ m[Spin]
            angle_vector[i] = m[Spin][e].ϕ
        elseif e ∈ m[L]
            angle_vector[i] = m[L][e].ϕ
        else
            angle_vector[i] = m[Pb][e].pb
        end
    end
    return angle_vector
end

@export all_angles(l::AbstractLedger) = all_angles!(zeros(length(l[Optimize])), l)

@export function set_all_angles!(l::AbstractLedger, angles)
    for i=1:10:length(angles)
        for j = 0:7
            l[Spin].data[i+j] = Spin(ϕ=angles[i+j])
        end
        l[L].data[div(i,10)+1] = L(ϕ=angles[i+8])
        l[L].data[div(i,10)+2] = L(ϕ=angles[i+9])
    end
    return l
end

@export function set_angle!(m::AbstractLedger, e::Entity, ϕ)
    if e in m[L]
        m[L][e] = L(ϕ = ϕ)
    else
        m[Spin][e] = Spin(ϕ=ϕ)
    end
end

@export function field_on_L(l)
    h = singleton(l, H)
    p = singleton(l, Model)
    Ls = l[L]
    easy = l[EasyAxis]
    Gd_spins = l[Spin]
    J1 = p.J1
    J2 = p.J2
    field = (H = [-p.χ_L * h.v * (Ls[1].v ⋅ h.v), -p.χ_L * h.v * (Ls[2].v ⋅ h.v)],
             L = [-p.γ * Ls[2].v * (Ls[1].v ⋅ Ls[2].v),-p.γ * Ls[1].v * (Ls[2].v ⋅ Ls[1].v)],
             easy   = [p.K_L * easy[9].v * (easy[9].v⋅Ls[1].v), p.K_L * easy[10].v * (easy[10].v⋅Ls[2].v)],
             Gd     = -1 .* [J2 * Gd_spins[1].v +
                             J1 * Gd_spins[2].v -
                             J1 * Gd_spins[3].v -
                             J2 * Gd_spins[4].v -
                             J2 * Gd_spins[5].v -
                             J1 * Gd_spins[6].v +
                             J1 * Gd_spins[7].v +
                             J2 * Gd_spins[8].v,
                             J1 * Gd_spins[1].v +
                             J2 * Gd_spins[2].v +
                             J2 * Gd_spins[3].v +
                             J1 * Gd_spins[4].v -
                             J1 * Gd_spins[5].v -
                             J2 * Gd_spins[6].v -
                             J2 * Gd_spins[7].v -
                             J1 * Gd_spins[8].v])
    return field
end

@export Gd_field_on_L1L2(Gd_spins, J1, J2) = -1 .* (J2 * Gd_spins[1].v +
                                   J1 * Gd_spins[2].v -
                                   J1 * Gd_spins[3].v -
                                   J2 * Gd_spins[4].v -
                                   J2 * Gd_spins[5].v -
                                   J1 * Gd_spins[6].v +
                                   J1 * Gd_spins[7].v +
                                   J2 * Gd_spins[8].v,
                                   J1 * Gd_spins[1].v +
                                   J2 * Gd_spins[2].v +
                                   J2 * Gd_spins[3].v +
                                   J1 * Gd_spins[4].v -
                                   J1 * Gd_spins[5].v -
                                   J2 * Gd_spins[6].v -
                                   J2 * Gd_spins[7].v -
                                   J1 * Gd_spins[8].v)

@export Gd_field_on_L1L2(l) = Gd_field_on_L1L2(l[Spin].data, singleton(l, Model).J1, singleton(l, Model).J2)

@export function field_on_Gds(l)
    h = singleton(l, H)
    p = singleton(l, Model)
    field = (H = fill(p.χ_Gd * h.v, 8),
             dipole = zeros(Vec3, 8),
             L      = zeros(Vec3, 8),
             easy   = zeros(Vec3, 8),
             sum    = zeros(Vec3, 8))

    spin = l[Spin]
    easy = l[EasyAxis]
    for (i, s) in enumerate(@entities_in(spin))
        d = spin[s] ⋅ easy[s]
        field.easy[i] += 2 * p.K_Gd * d * easy[s].v
    end
    J1 = p.J1
    J2 = p.J2
    L1 = l[L][1].v
    L2 = l[L][2].v
    @inbounds for i=1:8:length(spin)
        field.L[i]   += (-J1 * L2 - J2 * L1)
        field.L[i+1] += (-L1 * J1 - J2 * L2)
        field.L[i+2] += ( J1 * L1 - J2 * L2)
        field.L[i+3] += (-J1 * L2 + J2 * L1)
        field.L[i+4] -= (-J1 * L2 - J2 * L1)
        field.L[i+5] -= (-L1 * J1 - J2 * L2)
        field.L[i+6] -= ( J1 * L1 - J2 * L2)
        field.L[i+7] -= (-J1 * L2 + J2 * L1)
    end
    bonds = l[Bond]

    for b in bonds
        pr = b.r
        s1 = spin[b.i1].v
        s2 = spin[b.i1].v
        field.dipole[b.i1] += p.dipole_strength * (3 * (s2 ⋅ pr) * pr - b.nr2*s2)*b.inv_nr5
        field.dipole[b.i2] += p.dipole_strength * (3 * (s1 ⋅ pr) * pr - b.nr2*s1)*b.inv_nr5
    end
    for i = 1:8
        field.sum[i] = field.H[i] + field.dipole[i] + field.L[i] + field.easy[i]
    end
    return field
end

function fill_bonds!(l::AbstractLedger)
    Bond ∈ l && empty!(l[Bond])
    JBond ∈ l && empty!(l[JBond])

    p = singleton(l, Model)
    dip_bonds = Bond[]
    Jbonds = JBond[]
    a = Vec3(7.353099822998047, 0.0, 0)
    b = Vec3(0.0, 8.537099794497102, 0)
    c = Vec3(0.0, 0.0, 5.6807)
    for e1 in @entities_in(l[Spin] && l[Position])
        for e2 in @entities_in(l[Spin] && l[Position])
            p1 = l[Position][e1].p
            tp2 = l[Position][e2].p
            for ib=-1:1:1, ic = -1:1:1, ia=-2:2:2
                p2 = tp2 + ia*a + ib*b + ic * c
                n = norm(p2-p1)
                if e1 ∈ l[Gd] && e2 ∈ l[Gd] && 0.1 < n < p.dipole_cutoff
                    push!(dip_bonds, Bond(e1, e2, p2-p1))
                else
                    J = 0.0
                    if n - 3.91 < 0.01
                        J = 0.4
                    elseif n - 4.56 < 0.01
                        J = 0.25
                    elseif n - 3.57 < 0.01
                        J = 1.26
                    elseif n - 3.373 < 0.01
                        J = -12.8
                    elseif n - 3.544 < 0.01
                        J = -9.2
                    elseif n - 2.843 < 0.01
                        J = 13.67
                    elseif n - 3.800 < 0.01
                        J = -0.6
                    elseif n - 3.30315 < 0.01
                        J = 1.69
                    elseif n - 3.9577 < 0.01
                        J = -0.8
                    elseif n - 3.35751 < 0.01
                        J = -0.1
                    end
                    if J != 0
                        push!(Jbonds, JBond(e1, e2, J))
                    end
                end
            end
        end
    end
    for b in unique(Jbonds)
        Entity(l, b)
    end
    for b in unique(dip_bonds)
        Entity(l, b)
    end
end

function Overseer.Ledger(p::Model)
    na = 2
    nb = 1
    nc = 1
    m = Ledger(Stage(:E_Gd, [E_Gd_H(), E_Gd_easy()]),
               Stage(:E_L, [E_LL(), E_L_H(), E_L_easy()]),
               Stage(:E_Gd_L, [E_Gd_L()]),
               Stage(:E_all, [E_Gd_H(), E_L_H(), E_LL(), E_L_easy(), E_Gd_easy(), E_Gd_L(), E_Dipole(), E_L2H2(), E_Pb()]),
               Stage(:E_no_Pb, [E_Gd_H(), E_L_H(), E_LL(), E_L_easy(), E_Gd_easy(), E_Gd_L(), E_Dipole(), E_L2H2()]),
               Stage(:E_dipole, [E_Dipole()]))

    a = Vec3(7.353099822998047, 0.0, 0)
    b = Vec3(0.0, 8.537099794497102, 0)
    c = Vec3(0.0, 0.0, 5.6807)

    t = [Vec3(2.6507924861907957, 5.733516221984254, 0.0),
         Vec3(1.0257574253082276, 1.4649663247357028, 0.0),
         Vec3(6.3273423976898195, 7.0721334697614, 0.0),
         Vec3(4.702307336807251,  2.8035835725128484, 0.0)]

    Gd_ϕs = p.Gd_easy

    Gd_entities = Entity[]
    for ib=1:nb, ic=1:nc, ia = 1:na
        for (ϕ, easy, pos) in zip(Gd_ϕs, p.Gd_easy, t)
            push!(Gd_entities, Entity(m, Spin(ϕ=ϕ), EasyAxis(ϕ=easy), Position(pos + (ia-1)*a + (ib-1)*b + (ic-1) * c), Gd()))
        end
    end
    L1e = Entity(m, L(ϕ=p.L1_easy_ϕ), EasyAxis(ϕ=p.L1_easy_ϕ))
    L2e = Entity(m, L(ϕ=p.L2_easy_ϕ), EasyAxis(ϕ=p.L2_easy_ϕ))

    pe  = Entity(m, p, Etot(0.0), H(ϕ=0.0, magnitude=0.0), E_field())
    fill_bonds!(m)
    Entity(m, Pb())
    m.components[Optimize] = Component{Optimize}()
    optimize!(m, Spin, L, Pb)
    m.components[Optimize] = Component{Optimize}()

    return m
end

@export function JLedger(p::Model, na=2, nb=1, nc=1)
    m = Ledger(Stage(:E_Gd, [E_Gd_H(), E_Gd_easy()]),
               Stage(:E_L, [E_LL(), E_L_H(), E_L_easy()]),
               Stage(:E_Gd_L, [E_Gd_L()]),
               Stage(:E_all, [E_Gd_H(), E_Gd_easy(), E_Dipole(), E_Gd_easy4th(), E_Pb_J(),E_J()]),
               Stage(:E_dipole, [E_Dipole()]))

    a = Vec3(7.353099822998047, 0.0, 0)
    b = Vec3(0.0, 8.537099794497102, 0)
    c = Vec3(0.0, 0.0, 5.6807)

    t = [Vec3(2.6507924861907957, 5.733516221984254, 0.0),
         Vec3(1.0257574253082276, 1.4649663247357028, 0.0),
         Vec3(6.3273423976898195, 7.0721334697614, 0.0),
         Vec3(4.702307336807251,  2.8035835725128484, 0.0)]

    t_Mn = [Vec3(3.6765499114990234, 0.0,0.0),
            Vec3(0.0, 4.268549897248551, 0.0),
            Vec3(3.0280065071105957, 3.0033517077040806, 0.0),
            Vec3(6.704556418609619, 1.2651981895444706, 0.0),
            Vec3(4.325093315887451, 5.533748086793022, 0.0),
            Vec3(0.6485434043884277, 7.271901604952632, 0.0)]

    Gd_ϕs = p.Gd_easy

    Gd_entities = Entity[]
    for ib=1:nb, ic=1:nc, ia = 1:na
        for (phi, easy, pos) in zip(Gd_ϕs, p.Gd_easy, t)
            offset = iseven(ia) ? deg2rad(0) : deg2rad(180)
            push!(Gd_entities, Entity(m, Spin(ϕ=phi+offset), EasyAxis(ϕ=easy), Position(pos + (ia-1)*a + (ib-1)*b + (ic-1) * c), Gd()))
        end
    end

    Mn_easy = [EasyAxis(ϕ=deg2rad(23.4)), EasyAxis(ϕ=deg2rad(-23.4)), EasyAxis(ϕ=deg2rad(-23.4+180)), EasyAxis(ϕ=deg2rad(23.4+180)), EasyAxis(ϕ=deg2rad(-23.4)), EasyAxis(ϕ= deg2rad(23.4+180))]
    Mn_entities = Entity[]
    for ib=1:nb, ic=1:nc, ia = 1:na
        offset = iseven(ia) ? deg2rad(180) : 0
        for (easy, pos) in zip(Mn_easy, t_Mn)
            push!(Mn_entities, Entity(m, Spin(ϕ=easy.ϕ+offset), easy, Position(pos + (ia-1)*a + (ib-1)*b + (ic-1) * c), Mn()))
        end
    end
    pe  = Entity(m, p, Etot(0.0), H(ϕ=0.0, magnitude=0.0), E_field())
    fill_bonds!(m)
    Entity(m, Pb())
    m.components[Optimize] = Component{Optimize}()
    optimize!(m, Spin, Pb)
    m.components[Optimize] = Component{Optimize}()

    return m
end

@export function find_states(l::Ledger, max_tries_without_new=100)
    best_minimum = optimize(l, Spin, L).minimum
    angle_sets = Vector{Float64}[]
    tries = 0
    while length(angle_sets) < 4
        res = bboptimize(l, Spin, L)
        if best_fitness(res) <= 0.98*best_minimum
            new_angles = mod.(optim_values(l, (Spin, L)), 2π)
            if !any(x -> abs(sum(new_angles .- x)) <= 5e-4, angle_sets)
                tries = 0
                push!(angle_sets, new_angles)
            end
            if best_fitness(res) < best_minimum
                best_minimum = best_fitness(res)
            end
        end
        tries += 1
    end
    return angle_sets
end

@export function generate_Hsweep_states(l::Ledger, Hϕ, Hrange)
    H_ledger = copy(l)
    ledgers = Ledger[]
    for i in 1:length(Hrange)
        hmag = Hrange[i]
        if isempty(H_ledger[H])
            Entity(H_ledger, H(ϕ=deg2rad(Hϕ), magnitude=hmag))
        else
            H_ledger[H].data[1] = H(ϕ=deg2rad(Hϕ), magnitude=hmag)
        end
        res = optimize(H_ledger)
        set_all_angles!(H_ledger, res.minimizer[1:end-1])
        H_ledger[Pb].data[1] = Pb(res.minimizer[end])
        push!(ledgers, copy(H_ledger))
    end
    return ledgers
end

@export function generate_Esweep_states(l::Ledger, Erange)
    H_ledger = copy(l)
    ledgers = Ledger[]
    for i in 1:length(Erange)
        hmag = Erange[i]
        if isempty(H_ledger[E_field])
            Entity(H_ledger, E_field(hmag))
        else
            H_ledger[E_field].data[1] = E_field(hmag)
        end
        res = optimize(H_ledger)
        set_all_angles!(H_ledger, res.minimizer[1:end-1])
        H_ledger[Pb].data[1] = Pb(res.minimizer[end])
        push!(ledgers, copy(H_ledger))
    end
    return ledgers
end

function find_path_ranges(sweep_ledgers, start_H, stop_H)
    min, startid = findmin(map(x -> abs(start_H - singleton(x, H).magnitude), sweep_ledgers))
    min, stopid  = findmin(map(x -> abs(stop_H - singleton(x, H).magnitude), sweep_ledgers))

    real_start_H = singleton(sweep_ledgers[startid], H).magnitude
    real_stop_H  = singleton(sweep_ledgers[stopid], H).magnitude

    findnext_start_H_id = i -> findnext(x -> singleton(x, H).magnitude == real_start_H, sweep_ledgers, i)
    findnext_stop_H_id  = i -> findnext(x -> singleton(x, H).magnitude == real_stop_H, sweep_ledgers, i)
    path_ranges = StepRange[startid:stopid]
    for i=2:5
        lastid = last(path_ranges[i-1])
        r = i%2 == 0 ? reverse(findnext_stop_H_id(last(path_ranges[i-1])+1):findnext_start_H_id(last(path_ranges[i-1])+1)) : (findnext_start_H_id(first(path_ranges[i-1])+1):findnext_stop_H_id(first(path_ranges[i-1])+1))
        push!(path_ranges, r)
    end
    return path_ranges
end

function neb_barrier_ledgers(sweep_ledgers, path_ranges, c, n_interpolations)
    path_ledgers = [Ledger[] for i=1:length(path_ranges[1])]
    Threads.@threads for i = 1:length(path_ranges[1])
        r1 = path_ranges[1][i]
        r2 = path_ranges[2][i]
        r3 = path_ranges[3][i]
        r4 = path_ranges[4][i]
        r5 = path_ranges[5][i]
        ledgers1 = smart_interpolate_states(sweep_ledgers[r1:r2], n_interpolations)
        ledgers2 = smart_interpolate_states(sweep_ledgers[r2:r3], n_interpolations)
        ledgers3 = smart_interpolate_states(sweep_ledgers[r3:r4], n_interpolations)
        ledgers4 = smart_interpolate_states(sweep_ledgers[r4:r5], n_interpolations)
        for opt_ledgers in (ledgers1, ledgers2, ledgers3, ledgers4)
            for l in opt_ledgers
                l[H].data[1] = opt_ledgers[1][H].data[1]
            end
        end
        nebbed_ledgers1 = neb(ledgers1, c)
        nebbed_ledgers2 = neb(ledgers2, c)
        nebbed_ledgers3 = neb(ledgers3, c)
        nebbed_ledgers4 = neb(ledgers4, c)
        nebbed_ledgers1 = neb(nebbed_ledgers1, c/2)
        nebbed_ledgers2 = neb(nebbed_ledgers2, c/2)
        nebbed_ledgers3 = neb(nebbed_ledgers3, c/2)
        nebbed_ledgers4 = neb(nebbed_ledgers4, c/2)
        path_ledgers[i] = [nebbed_ledgers1..., nebbed_ledgers2..., nebbed_ledgers3..., nebbed_ledgers4...]
    end
    return path_ledgers
end

"Calculates the energy barriers for a set of sweep states (full cycle) between start_H and stop_H, performing NEB with stiffness c"
@export function energy_barriers(sweep_ledgers, start_H, stop_H, c, n_interpolations=10)

    path_ranges = find_path_ranges(sweep_ledgers, start_H, stop_H)
    for r in path_ranges
        sweep_ledgers[first(r)][H].data[1] = H(magnitude=start_H, sweep_ledgers[first(r)][H].data[1])
        optimize(sweep_ledgers[first(r)], Spin, L, Pb)
        sweep_ledgers[last(r)][H].data[1] = H(magnitude=stop_H, sweep_ledgers[first(r)][H].data[1])
        optimize(sweep_ledgers[last(r)], Spin, L, Pb)
    end
    barrier_ledgers = neb_barrier_ledgers(sweep_ledgers, path_ranges, c, n_interpolations)
    return [getfield.(calc_etot.(ls, :E_all), :e) for ls in barrier_ledgers], H_magnitudes(sweep_ledgers[path_ranges[1]])
end

# "Calculates the energy barriers for a set of sweep states (full cycle) between start_H and stop_H, performing NEB with stiffness c"
# @export function energy_barriers_new(sweep_ledgers, c, n_interpolations=10)
#     Hstart, Hstop = find_hysteresis(sweep_ledgers)
#     path_ranges = find_path_ranges(sweep_ledgers, Hstart, Hstop)
#     sweep_segments = map(x -> sweep_ledgers[x], path_ranges)

#     @assert all(x -> length(x) == length(sweep_segments[1]), sweep_segments) "Sweep segments don't have the same length"
#     # Hmags = H_magnitudes(sweep_ledgers)
#     energies = [Float64[] for i = 1:length(sweep_segments[1])-1]
#     for is = 1:length(sweep_segments)-1
#         for i = 1:length(sweep_segments[is])-1
#             ledgers_to_interpolate = [sort(sweep_segments[is], by=x->singleton(x, H).magnitude)[i:end]..., reverse(sort(sweep_segments[is+1], by=x->singleton(x, H).magnitude))[1:end-(i-1)]...]
#             set_optim_values!(ledgers_to_interpolate[1], optimize(ledgers_to_interpolate[1], Spin, L, Pb).minimizer)
#             set_optim_values!(ledgers_to_interpolate[end], optimize(ledgers_to_interpolate[end], Spin, L, Pb).minimizer)
#             ledgers_to_neb = smart_interpolate_states(ledgers_to_interpolate, n_interpolations)
#             for l in ledgers_to_neb
#                 l[H].data[1] = singleton(ledgers_to_neb[1], H)
#             end
#             nebbed_ledgers = neb(ledgers_to_neb, c)
#             nebbed_ledgers = neb(nebbed_ledgers, c/2)
#             append!(energies[i], getfield.(calc_etot.(nebbed_ledgers), :e))
#         end
#     end
#     return energies
# end

"""
    smart_interpolate_states(ledgers, n)

Takes a set of states, spaces out the angle interpolation such that it's uniform,
and adds or removes states until the set number n is reached.
"""
function smart_interpolate_states(ledgers, n)
    out_ledgers = deepcopy.(ledgers)
    space_out_states!(out_ledgers)
    n_to_add = n - length(out_ledgers)
    interval = 1
    while n_to_add > 0
        if interval + 1 >= length(out_ledgers)
            interval = 1
        end
        l1 = out_ledgers[interval]
        l2 = out_ledgers[interval+1]
        l_to_add = deepcopy(l1)
        set_optim_values!(l_to_add, (all_angles(l1) .+ all_angles(l2))./2)
        insert!(out_ledgers, interval + 1, l_to_add)
        interval += 2
        n_to_add -= 1
    end
    interval = length(out_ledgers) - 1
    while n_to_add < 0
        if interval <= 1
            interval = length(out_ledgers)-1
        end
        deleteat!(out_ledgers, interval)
        interval -= 2
        n_to_add += 1
    end

    space_out_states!(out_ledgers)
    return out_ledgers
end

@export function find_sweep_segments(hsweep)
    Hrange = H_magnitudes(hsweep)

    Hsegments = Vector{AbstractLedger}[]
    t_hsegment = AbstractLedger[]
    cursign = sign(Hrange[2]-Hrange[1])

    for i = 2:length(Hrange)
        tsign = sign(Hrange[i]-Hrange[i-1])
        if tsign != cursign
            cursign = tsign
            push!(t_hsegment, hsweep[i-1])
            push!(Hsegments, copy(t_hsegment))
            empty!(t_hsegment)
        else
            push!(t_hsegment, hsweep[i-1])
        end
    end
    push!(t_hsegment, hsweep[end])
    push!(Hsegments, t_hsegment)
    return Hsegments
end

@export function find_hysteresis(hsweep,threshold=0.001)
    Hrange = H_magnitudes(hsweep)
    Hsegments = find_sweep_segments(hsweep)
    pbs = [Pb.(h) for h in Hsegments]
    Hrs = [H_magnitudes(h) for h in Hsegments]

    diffs_min = [norm(pbs[1][i] - pbs[4][findfirst(x->x==Hrs[1][i], Hrs[4])]) for i=1:length(Hrs[1])-1]

    diffs_max = [norm(pbs[1][i] - pbs[2][findfirst(x->x==Hrs[1][i], Hrs[2])]) for i=1:length(Hrs[1])-1]

    imin = findfirst(x->x>threshold, diffs_min)-1
    imax = findlast(x -> x>threshold, diffs_max)+1
    return Hrange[imin], Hrange[imax]
end


text:
# generalTools.jl

export logsumexp, rDirichlet, rGenDirichlet, lmvbeta, shape_Dir2GenDir;

"""
    logsumexp(x[, usemax])

Computes `log(sum(exp(x)))` in a stable manner.

### Example
```julia
  x = rand(5)
  logsumexp(x)
  log(sum(exp.(x)))
```
"""
function logsumexp(x::Array{Float64}, usemax::Bool=true)
  if usemax
    m = maximum(x)
  else
    m = minimum(x)
  end

  return m + log(sum(exp.(x .- m)))
end

"""
    logsumexp(x, region[, usemax])

Computes `log(sum(exp(x)))` in a stable manner along dimensions specified.

### Example
```julia
  x = reshape(collect(1:24)*1.0, (2,3,4))
  logsumexp(x, 2)
```
"""
function logsumexp(x::Array{Float64}, region, usemax::Bool=true)
  if usemax
    ms = maximum(x, dims=region)
  else
    ms = minimum(x, dims=region)
  end
  bc_xminusms = broadcast(-, x, ms)

  expxx = exp.(bc_xminusms)
  sumexpxx = sum(expxx, dims=region)

  return log.(sumexpxx) .+ ms
end

"""
    rDirichlet(α[; logout=false, logcompute=true])

  Single draw from Dirichlet distribution, option for log scale, computation on log scale.

  ### Example
  ```julia
  rDirichlet(ones(5),true)
  ```
"""
function rDirichlet(α::Array{Float64, 1}; logout::Bool=false, logcompute::Bool=true)
  all(α .> 0.0) || error("All shape parameters must be positive.")

  k = length(α)
  xx = Vector{Float64}(undef, k) # allows changes to elements
  s = 0.0

  if logcompute
    for i in 1:k
      xx[i] = log(rand(Gamma(α[i], 1.0)))
    end
    s = logsumexp(xx)
    out = xx .- s
    if !logout
        out = exp.(out)
    end
  else
    for i in 1:k
      xx[i] = rand(Gamma(α[i], 1.0))
      s += xx[i]
    end
    out = xx / s
    if logout
        out = log.(out)
    end
  end

  return out
end

"""
    rGenDirichlet(a, b[, logout])

  Single draw from the generalized Dirichlet distribution (Connor and Mosimann '69), option for log scale.

  ### Example
  ```julia
  a = ones(5)
  b = ones(5)
  rGenDirichlet(a, b)
  ```
"""
function rGenDirichlet(a::Vector{Float64}, b::Vector{Float64}; logout::Bool=false)
    n = length(a)
    K = n + 1
    length(b) == n || error("Dimension mismatch between a and b.")
    all(a .> 0.0) || error("All elements of a must be positive.")
    all(b .> 0.0) || error("All elements of b must be positive.")

    lz = Vector{Float64}(undef, n)
    loneminusz = Vector{Float64}(undef, n)
    for i = 1:n
        lx1 = log( rand( Distributions.Gamma( a[i] ) ) )
        lx2 = log( rand( Distributions.Gamma( b[i] ) ) )
        lxm = max(lx1, lx2)
        lxsum = lxm + log( exp(lx1 - lxm) + exp(lx2 - lxm) ) # logsumexp
        lz[i] = lx1 - lxsum
        loneminusz[i] = lx2 - lxsum
    end

    ## break the Stick
    lw = Vector{Float64}(undef, K)
    lwhatsleft = 0.0

    for i in 1:n
        lw[i] = lz[i] + lwhatsleft
        # lwhatsleft += log( 1.0 - exp(lw[i] - lwhatsleft) ) # logsumexp (not numerically stable)
        lwhatsleft += loneminusz[i]
    end
    lw[K] = copy(lwhatsleft)

    if logout
        out = (lw, lz)
    else
        out = (exp.(lw), exp.(lz))
    end
    return out
end

"""
    lmvbeta(x)

Computes the natural log of ``∏(Γ(x)) / Γ(sum(x))``.

### Example
```julia
  x = rand(5)
  lmvbeta(x)
```
"""
function lmvbeta(x::Array{Float64})
    return sum(loggamma.(x)) - loggamma(sum(x))
end

"""
    shape_Dir2GenDir(α)

Converts Dirichlet shape parameter vector to set of shape parameters for latent beta variables in the ggeneralized Dirichlet distribution (Connor and Mosimann '69).

### Example
```julia
  α = ones(5)
  shape_Dir2GenDir(α)
```
"""
function shape_Dir2GenDir(α::Vector{Float64})
    all(α .> 0.0) || error("All shape parameters must be positive.")

    K = length(α)
    n = K - 1
    rcrα = cumsum( α[range(K, step=-1, length=n)] )[range(n, step=-1, length=n)]

    return (α[1:n], rcrα)
end